Mercurial > repos > artbio > probecoverage

annotate probecoverage.r @ 0:3c0451ca266e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
author artbio
date Tue, 07 Jan 2020 11:08:31 +0000
parents
children 9eb4a7000c1e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
1 ## Setup R error handling to go to stderr
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
2 options( show.error.messages=F,
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
3 error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
4 warnings()
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
5 library(optparse)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
6 library(ggplot2)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
7 library(reshape2)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
8
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
9 option_list <- list(
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
10 make_option(c("-i", "--input"), type="character", help="Path to dataframe"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
11 make_option(c("-t", "--title"), type="character", help="Main Title"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
12 make_option("--xlab", type = "character", help="X-axis legend"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
13 make_option("--ylab", type = "character", help="Y-axis legend"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
14 make_option("--sample", type = "character", help="a space separated of sample labels"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
15 make_option("--method", type = "character", help="bedtools or pysam"),
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
16 make_option(c("-o", "--output"), type = "character", help="path to the pdf plot")
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
17 )
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
18
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
19 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
20 args = parse_args(parser)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
21 samples = substr(args$sample, 2, nchar(args$sample)-2)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
22 samples = strsplit(samples, ", ")
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
23
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
24 # data frames implementation
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
25
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
26 Table <- read.delim(args$input, header=F)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
27 headers = c("chromosome", "start", "end", "id")
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
28 for (i in seq(1, length(Table)-4)) {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
29 headers <- c(headers, samples[[1]][i])
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
30 colnames(Table) <- headers
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
31 }
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
32
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
33 ## function
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
34 if (args$method == 'bedtools') {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
35 cumul <- function(x,y) sum(Table[,y]/(Table$end-Table$start) > x)/length(Table$chromosome)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
36 } else {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
37 cumul <- function(x,y) sum(Table[,y] > x)/length(Table$chromosome)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
38 }
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
39 scaleFUN <- function(x) sprintf("%.3f", x)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
40
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
41 ## end of function
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
42 ## let's do a dataframe before plotting
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
43 if (args$method == 'bedtools') {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
44 maxdepth <- trunc(max(Table[,5:length(Table)]/(Table$end-Table$start))) + 20
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
45 } else {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
46 maxdepth <- trunc(max(Table[,5:length(Table)])) + 20
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
47 }
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
48
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
49 graphpoints <- data.frame(1:maxdepth)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
50 i <- 5
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
51 for (colonne in colnames(Table)[5:length(colnames(Table))]) {
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
52 graphpoints <- cbind(graphpoints, mapply(cumul, 1:maxdepth, rep(i, maxdepth)))
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
53 i <- i + 1
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
54 }
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
55 colnames(graphpoints) <- c("Depth", colnames(Table)[5:length(Table)])
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
56 maxfrac = max(graphpoints[,2:length(graphpoints)])
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
57
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
58 graphpoints <- melt(graphpoints, id.vars="Depth", variable.name="Samples", value.name="sample_value")
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
59
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
60 ## GRAPHS
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
61
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
62 pdf(file=args$output)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
63 ggplot(data=graphpoints, aes(x=Depth, y=sample_value, colour=Samples)) +
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
64 geom_line(size=1) +
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
65 scale_x_continuous(trans='log2', breaks = 2^(seq(0,log(maxdepth, 2)))) +
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
66 scale_y_continuous(breaks = seq(0, maxfrac, by=maxfrac/10), labels=scaleFUN) +
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
67 labs(x=args$xlab, y=args$ylab, title=args$title) +
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
68 theme(legend.position="top", legend.title=element_blank(), legend.text=element_text(colour="blue", size=7))
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
69
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
70
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
71 ## facet_wrap(~Samples, ncol=2)
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
72
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
73 devname=dev.off()
3c0451ca266e "planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/probecoverage commit 62a8b073b9ac98b0231641e5266768e7f8b80b89"
artbio
parents:
diff changeset
74