Mercurial > repos > artbio > probecoverage

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