comparison probecoverage.r @ 0:dbeb4a0abfc6 draft

planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 102edc0c210d94e9d72f913e2d18c19220c4167c
author artbio
date Sat, 23 Sep 2017 05:11:49 -0400
parents
children ebe5ec2e244d
comparison
equal deleted inserted replaced
-1:000000000000 0:dbeb4a0abfc6
1 ## Setup R error handling to go to stderr
2 options( show.error.messages=F,
3 error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
4 warnings()
5 library(optparse)
6 library(ggplot2)
7 library(reshape2)
8
9 option_list <- list(
10 make_option(c("-i", "--input"), type="character", help="Path to dataframe"),
11 make_option(c("-t", "--title"), type="character", help="Main Title"),
12 make_option("--xlab", type = "character", help="X-axis legend"),
13 make_option("--ylab", type = "character", help="Y-axis legend"),
14 make_option("--sample", type = "character", help="a space separated of sample labels"),
15 make_option(c("-o", "--output"), type = "character", help="path to the pdf plot")
16 )
17
18 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
19 args = parse_args(parser)
20 samples = substr(args$sample, 2, nchar(args$sample)-2)
21 samples = strsplit(samples, ", ")
22
23 # data frames implementation
24
25 Table <- read.delim(args$input, header=F)
26 headers = c("chromosome", "start", "end", "id")
27 for (i in seq(1, length(Table)-4)) {
28 headers <- c(headers, samples[[1]][i])
29 colnames(Table) <- headers
30 }
31
32 ## function
33
34 cumul <- function(x,y) sum(Table[,y]/(Table$end-Table$start) > x)/length(Table$chromosome)
35
36 ## end of function
37 ## let's do a dataframe before plotting
38 maxdepth <- trunc(max(Table[,5:length(Table)]/(Table$end-Table$start))) + 20
39 graphpoints <- data.frame(1:maxdepth)
40 i <- 5
41 for (colonne in colnames(Table)[5:length(colnames(Table))]) {
42 graphpoints <- cbind(graphpoints, mapply(cumul, 1:maxdepth, rep(i, maxdepth)))
43 i <- i + 1
44 }
45 colnames(graphpoints) <- c("Depth", colnames(Table)[5:length(Table)])
46 maxfrac = max(graphpoints[,2:length(graphpoints)])
47
48 graphpoints <- melt(graphpoints, id.vars="Depth", variable.name="Samples", value.name="sample_value")
49
50 ## GRAPHS
51
52 pdf(file=args$output)
53 ggplot(data=graphpoints, aes(x=Depth, y=sample_value, group=Samples, colour=Samples)) +
54 geom_line(size=1) +
55 scale_x_continuous(trans='log2', breaks = 2^(seq(0,log(maxdepth, 2)))) +
56 scale_y_continuous(breaks = seq(0, maxfrac, by=maxfrac/10)) +
57 labs(x=args$xlab, y=args$ylab, title=args$title) +
58 theme(legend.position="top", legend.title=element_blank(), legend.text=element_text(colour="blue", size=7))
59
60 devname=dev.off()
61