Mercurial > repos > artbio > probecoverage

annotate probecoverage.r @ 1:9eb4a7000c1e draft default tip

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