Mercurial > repos > artbio > probecoverage
changeset 8:0adb846ca056 draft default tip
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 | |
| files | multicov.py probecoverage.r probecoverage.xml test-data/graph.pdf test-data/graph_pysam.pdf |
| diffstat | 5 files changed, 46 insertions(+), 35 deletions(-) [+] |
line wrap: on
line diff
--- a/multicov.py Thu Jun 03 23:59:31 2021 +0000 +++ b/multicov.py Tue Feb 20 08:34:16 2024 +0000 @@ -53,6 +53,7 @@ line_counter += 1 suffix = '\t'.join(crossline) print('%s\t%s' % (prefix, suffix)) + F.close() if __name__ == "__main__":
--- a/probecoverage.r Thu Jun 03 23:59:31 2021 +0000 +++ b/probecoverage.r Tue Feb 20 08:34:16 2024 +0000 @@ -1,9 +1,11 @@ ## Setup R error handling to go to stderr -options(show.error.messages = F, - error = function() { - cat(geterrmessage(), file = stderr()); q("no", 1, F) - } - ) +options( + show.error.messages = FALSE, + error = function() { + cat(geterrmessage(), file = stderr()) + q("no", 1, FALSE) + } +) warnings() library(optparse) library(ggplot2) @@ -17,44 +19,44 @@ make_option("--sample", type = "character", help = "a space separated of sample labels"), make_option("--method", type = "character", help = "bedtools or pysam"), make_option(c("-o", "--output"), type = "character", help = "path to the pdf plot") - ) - +) + parser <- OptionParser(usage = "%prog [options] file", option_list = option_list) args <- parse_args(parser) samples <- substr(args$sample, 2, nchar(args$sample) - 2) samples <- strsplit(samples, ", ") - + # data frames implementation -table <- read.delim(args$input, header = F) +table <- read.delim(args$input, header = FALSE) headers <- c("chromosome", "start", "end", "id") for (i in seq(1, length(table) - 4)) { headers <- c(headers, samples[[1]][i]) -colnames(table) <- headers + colnames(table) <- headers } ## function if (args$method == "bedtools") { cumul <- function(x, y) sum(table[, y] / (table$end - table$start) > x) / length(table$chromosome) - } else { +} else { cumul <- function(x, y) sum(table[, y] > x) / length(table$chromosome) - } +} scale_fun <- function(x) sprintf("%.3f", x) ## end of function ## let's do a dataframe before plotting if (args$method == "bedtools") { maxdepth <- trunc(max(table[, 5:length(table)] / (table$end - table$start))) + 20 - } else { +} else { maxdepth <- trunc(max(table[, 5:length(table)])) + 20 - } +} graphpoints <- data.frame(1:maxdepth) i <- 5 for (colonne in colnames(table)[5:length(colnames(table))]) { - graphpoints <- cbind(graphpoints, mapply(cumul, 1:maxdepth, rep(i, maxdepth))) + graphpoints <- cbind(graphpoints, mapply(cumul, 1:maxdepth, rep(i, maxdepth))) i <- i + 1 - } +} colnames(graphpoints) <- c("Depth", colnames(table)[5:length(table)]) maxfrac <- max(graphpoints[, 2:length(graphpoints)]) @@ -64,11 +66,13 @@ pdf(file = args$output) ggplot(data = graphpoints, aes(x = Depth, y = sample_value, colour = Samples)) + - geom_line(size = 1) + - scale_x_continuous(trans = "log2", breaks = 2 ^ (seq(0, log(maxdepth, 2)))) + - scale_y_continuous(breaks = seq(0, maxfrac, by = maxfrac / 10), labels = scale_fun) + - labs(x = args$xlab, y = args$ylab, title = args$title) + - theme(legend.position = "top", legend.title = element_blank(), - legend.text = element_text(colour = "blue", size = 7)) + geom_line(size = 1) + + scale_x_continuous(trans = "log2", breaks = 2^(seq(0, log(maxdepth, 2)))) + + scale_y_continuous(breaks = seq(0, maxfrac, by = maxfrac / 10), labels = scale_fun) + + labs(x = args$xlab, y = args$ylab, title = args$title) + + theme( + legend.position = "top", legend.title = element_blank(), + legend.text = element_text(colour = "blue", size = 7) + ) devname <- dev.off()
--- a/probecoverage.xml Thu Jun 03 23:59:31 2021 +0000 +++ b/probecoverage.xml Tue Feb 20 08:34:16 2024 +0000 @@ -1,19 +1,25 @@ -<tool id="probecoverage" name="Probe Coverage" version="0.7.0"> +<tool id="probecoverage" name="Probe Coverage" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="@PROFILE@"> <description></description> - <requirements> - <requirement type="package" version="1.10=h9402c20_2">samtools</requirement> - <requirement type="package" version="2.29.2=hc088bd4_0">bedtools</requirement> - <requirement type="package" version="1.6.4=r36h6115d3f_0">r-optparse</requirement> - <requirement type="package" version="3.2.1=r36h6115d3f_0">r-ggplot2</requirement> - <requirement type="package" version="1.17.3=py37h95a1406_0">numpy</requirement> - <requirement type="package" version="0.16.0.1">pysam</requirement> - </requirements> + <macros> + <token name="@TOOL_VERSION@">0.22.0</token> + <token name="@VERSION_SUFFIX@">0</token> + <token name="@PROFILE@">23.0</token> + </macros> + <requirements> + <requirement type="package" version="@TOOL_VERSION@">pysam</requirement> + <requirement type="package" version="2.31.1">bedtools</requirement> + <requirement type="package" version="1.7.4">r-optparse</requirement> + <requirement type="package" version="3.4.4">r-ggplot2</requirement> + <requirement type="package" version="1.4.4">r-reshape2</requirement> + <requirement type="package" version="1.26.4">numpy</requirement> + </requirements> <stdio> <exit_code range="1:" level="fatal" description="Tool exception" /> </stdio> <command detect_errors="exit_code"><![CDATA[ - #for $file in $inputs - samtools index '$file' && + #for $i, $file in enumerate($inputs): + ln -f -s '$file' ${i}.bam && + ln -f -s '$file.metadata.bam_index' ${i}.bam.bai && #end for #if $method == 'pysam': python $__tool_directory__/multicov.py @@ -21,8 +27,8 @@ bedtools multicov #end if -bams - #for $file in $inputs - '$file' + #for $i, $file in enumerate($inputs): + '${i}.bam' #end for -bed '$bed' > $coverage_dataframe && Rscript '$__tool_directory__'/probecoverage.r