Mercurial > repos > artbio > probecoverage

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