|
Commit message:
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/coverage_report commit 1b647bb088f62c1369d963f030caecaa9ee003c7 |
|
added:
CoverageReport.pl CoverageReport.xml README test-data/answer.pdf test-data/ecoli.bam test-data/ecoli_6.bed |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 CoverageReport.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CoverageReport.pl Wed Oct 25 12:37:35 2017 -0400 |
| [ |
| b'@@ -0,0 +1,851 @@\n+#!/usr/bin/env perl\n+# load modules\n+use warnings;\n+use Getopt::Std;\n+use File::Basename;\n+use Number::Format;\n+use Cwd;\n+\n+my %opts;\n+\n+# number format\n+my $de = new Number::Format(-thousands_sep =>\',\',-decimal_point => \'.\');\n+\n+##########\n+## opts ##\n+##########\n+## input files\n+# b : path to input (b)am file\n+# t : path to input (t)arget regions in BED format\n+## output files\n+# o : report pdf (o)utput file\n+## entries in the report\n+# r : Coverage per (r)egion (boolean)\n+# s : (s)ubregion coverage if average < specified (plots for positions along target region) (boolean)\n+# S : (S)ubregion coverage for ALL failed exons => use either s OR S or you will have double plots. \n+# A : (A)ll exons will be plotted. \n+# L : (L)ist failed exons instead of plotting\n+# m : (m)inimal Coverage threshold\n+# f : fraction of average as threshold\n+# n : sample (n)ame.\n+ \n+\n+getopts(\'b:t:o:rsSALm:n:f:\', \\%opts) ;\n+\n+my $tmp = getcwd();\n+# make output directory in (tmp) working dir\n+our $wd = "$tmp/Coverage.".int(rand(1000));\n+while (-d $wd) {\n+\t$wd = "$tmp/Coverage.".int(rand(1000));\n+}\n+system("mkdir $wd");\n+\n+## variables\n+our %commandsrun = ();\n+my ($thresh,$frac,$pdffile,$samplename,$totalmapped);\n+\n+\n+if (!exists($opts{\'b\'}) || !-e $opts{\'b\'}) {\n+\tdie(\'Bam File not found\');\n+}\n+if (!exists($opts{\'t\'}) || !-e $opts{\'t\'}) {\n+\tdie(\'Target File (BED) not found\');\n+}\n+\n+if (exists($opts{\'m\'})) {\n+\t$thresh = $opts{\'m\'};\n+}\n+else {\n+\t$thresh = 40;\n+}\n+\n+if (exists($opts{\'f\'})) {\n+\t$frac = $opts{\'f\'};\n+}\n+else {\n+\t$frac = 0.2;\n+}\n+\n+if (exists($opts{\'o\'})) {\n+\t$pdffile = $opts{\'o\'};\n+}\n+else {\n+\t$pdffile = "$wd/CoverageReport.pdf";\n+}\n+\n+\n+# 1. Global Summary => default\n+&GlobalSummary($opts{\'b\'}, $opts{\'t\'});\n+\n+# 2. Coverage per position\n+&SubRegionCoverage($opts{\'b\'}, $opts{\'t\'});\n+our %filehash;\n+if (exists($opts{\'s\'}) || exists($opts{\'S\'}) || exists($opts{\'A\'}) || exists($opts{\'L\'})) {\n+\tsystem("mkdir $wd/SplitFiles");\n+\t## get position coverages\n+\t## split input files\n+\topen IN, "$wd/Targets.Position.Coverage";\n+\tmy $fileidx = 0;\n+\tmy $currreg = \'\';\n+\twhile (<IN>) {\n+\t\tmy $line = $_;\n+\t\tchomp($line);\n+\t\tmy @p = split(/\\t/,$line);\n+\t\tmy $reg = $p[0].\'-\'.$p[1].\'-\'.$p[2]; #.$p[3];\n+\t\tmy $ex = $p[3];\n+\t\tif ($reg ne $currreg) {\n+\t\t\t## new exon open new outfile\n+\t\t\tif ($currreg ne \'\') {\n+\t\t\t\t## filehandle is open. close it\n+\t\t\t\tclose OUT; \n+\t\t\t}\n+\t\t\tif (!exists($filehash{$reg})) {\n+\t\t\t\t$fileidx++;\n+\t\t\t\t$filehash{$reg}{\'idx\'} = $fileidx;\n+\t\t\t\t$filehash{$reg}{\'exon\'} = $ex;\n+\t\t\t\topen OUT, ">> $wd/SplitFiles/File_$fileidx.txt";\n+\t\t\t\t$currreg = $reg;\n+\t\t\t}\n+\t\t\telse {\n+\t\t\t\topen OUT, ">> $wd/SplitFiles/File_".$filehash{$reg}{\'idx\'}.".txt";\n+\t\t\t\t$currreg = $reg;\n+\t\t\t}\n+\t\t}\n+\t\t## print the line to the open filehandle. \t\n+\t\tprint OUT "$line\\n";\n+\t}\n+\tclose OUT;\n+\tclose IN;\n+\n+}\n+\n+## sort output files according to targets file\n+if (exists($opts{\'r\'}) ) {\n+\tmy %hash = ();\n+\topen IN, "$wd/Targets.Global.Coverage";\n+\twhile (<IN>) {\n+ chomp;\n+ my @p = split(/\\t/,$_) ;\n+ $hash{$p[3]} = $_;\n+\t}\n+\tclose IN;\n+\topen OUT, ">$wd/Targets.Global.Coverage";\n+\topen IN, $opts{\'t\'};\n+\twhile (<IN>) {\n+ chomp;\n+ my @p = split(/\\t/,$_) ;\n+ print OUT $hash{$p[3]} . "\\n";\n+\t}\n+\tclose IN;\n+\tclose OUT;\n+}\n+\n+\n+####################################\n+## PROCESS RESULTS & CREATE PLOTS ##\n+####################################\n+system("mkdir $wd/Report");\n+\n+system("mkdir $wd/Rout");\n+system("mkdir $wd/Plots");\n+\n+$samplename = $opts{\'n\'};\n+$samplename =~ s/_/\\\\_/g;\n+$samplename =~ s/\\..+$//g;\n+\n+# 0. Preamble\n+## compose preamble\n+open OUT, ">$wd/Report/Report.tex";\n+print OUT \'\\documentclass[a4paper,10pt]{article}\'."\\n";\n+print OUT \'\\usepackage[left=2cm,top=1.5cm,right=1.5cm,bottom=2.5cm,nohead]{geometry}\'."\\n";\n+print OUT \'\\usepackage{longtable}\'."\\n";\n+print OUT \'\\usepackage{fancyhdr}\'."\\n";\n+print OUT \'\\usepackage{fontspec}\' . "\\n";\n+print OUT \'\\usepackage{color}\'."\\n";\n+print OUT \'\\definecolo'..b'esize\\begin{longtable}[l]{@{\\extracolsep{\\fill}}llll}\'."\\n".\'\\hline\'."\\n";\n+\tprint TEX \'\\textbf{Target Name} & \\textbf{Genomic Position} & \\textbf{Avg.Coverage} & \\textbf{Min.Coverage} \\\\\\\\\'."\\n".\'\\hline\'."\\n";\n+\tprint TEX \'\\endhead\'."\\n";\n+\tprint TEX \'\\hline \'."\\n".\'\\multicolumn{4}{r}{{\\textsl{\\footnotesize Continued on next page}}} \\\\\\\\ \'."\\n".\'\\hline\' ."\\n". \'\\endfoot\' . "\\n". \'\\endlastfoot\' . "\\n";\n+ \n+\t$col = 1;\n+\topen IN, "$wd/Targets.Global.Coverage";\n+\twhile (<IN>) {\n+\t\tchomp($_);\n+\t\tmy @p = split(/\\t/,$_);\n+\t\tmy $region = $p[0].\'-\'.$p[1].\'-\'.$p[2];\n+\t\tmy $exon = $filehash{$region}{\'exon\'};\n+\t\t# grep exon to tmp file\n+\t\tmy $exonfile = "$wd/SplitFiles/File_".$filehash{$region}{\'idx\'}.".txt";\n+\t\t## determine transcript orientation.\n+ my $firstline = `head -n 1 $exonfile`;\n+ my @firstcols = split(/\\t/,$firstline);\n+ my $orient = $firstcols[5];\n+\t\tmy $genomicchr = $firstcols[0];\n+\t\tmy $genomicstart = $firstcols[1];\n+\t\tmy $genomicstop = $firstcols[2];\n+\n+\t\tif ($orient eq \'+\') {\n+\t\t\tmy $bps = $genomicstop - $genomicstart + 1;\n+\t\t\t$subtitle = "$genomicchr:".$de->format_number($genomicstart)."+".$de->format_number($genomicstop);\n+\n+\t\t}\n+\t\telse {\n+\t\t\tmy $bps = $genomicstop - $genomicstart + 1;\n+\t\t\t$subtitle = "$genomicchr:".$de->format_number($genomicstart)."-".$de->format_number($genomicstop);\n+\t\t}\n+\n+\t\t# check if failed\n+\t\tmy $cs = `cut -f $covcol \'$exonfile\' `;\n+\t\tmy @c = split(/\\n/,$cs);\n+\t\tmy ($avg,$med,$min,$max,$first,$third,$ontarget) = arraystats(@c);\n+\n+\t\tif ($min >= $thresh) {\n+\t\t\t# lowest coverage > threshold => skip\n+\t\t\tnext;\n+\t\t}\n+\n+\t\t# print to .tex table \n+\t\tif (length($exon) > 30) {\n+\t\t\t$exon = substr($exon,0,27) . \'...\';\n+\t\t}\n+\t\t$exon =~ s/_/ /g;\n+\t\t$exon =~ s/\\|/ /g;\n+\n+\t\tprint TEX "$exon & $subtitle & ".int($avg)." & $min ".\'\\\\\\\\\'."\\n"; \n+\t}\n+\tclose IN;\n+\tprint TEX \'\\hline\'."\\n";\n+\tprint TEX \'\\end{longtable}}\'."\\n";\n+\tclose TEX;\n+}\n+\n+\n+## Close document\n+open OUT, ">>$wd/Report/Report.tex";\n+print OUT \'\\label{endofdoc}\'."\\n";\n+print OUT \'\\end{document}\'."\\n";\n+close OUT;\n+system("cd $wd/Report && tectonic Report.tex") == 0 || die "Could not run tectonic with following error \'$!\'\\n";\n+\n+## mv report to output file\n+system("cp -f $wd/Report/Report.pdf \'$pdffile\'");\n+##create tar.gz file\n+system("mkdir $wd/Results");\n+system("cp -Rf $wd/Plots $wd/Results/");\n+system("cp -Rf $wd/Report/ $wd/Results/");\n+if (-e "$wd/Targets.Global.Coverage") {\n+\tsystem("cp -Rf $wd/Targets.Global.Coverage $wd/Results/");\n+}\n+if (-e "$wd/Targets.Position.Coverage") {\n+\tsystem("cp -Rf $wd/Targets.Position.Coverage $wd/Results/");\n+}\n+\n+\n+exit;\n+\n+###############\n+## FUNCTIONS ##\n+###############\n+sub arraystats{\n+\tmy @array = @_;\n+\tmy $count = scalar(@array);\n+\t@array = sort { $a <=> $b } @array;\n+\t# median\n+\tmy $median = 0;\n+\tif ($count % 2) { \n+\t\t$median = $array[int($count/2)]; \n+\t} else { \n+\t\t$median = ($array[$count/2] + $array[$count/2 - 1]) / 2; \n+\t} \n+\t# average \n+\tmy $sum = 0;\n+\tforeach (@array) { $sum += $_; } \n+\tmy $average = $sum / $count;\n+\t# quantiles (rounded)\n+\tmy $quart = int($count/4) ;\t\n+\tmy $first = $array[$quart];\n+\tmy $third = $array[($quart*3)];\n+\tmy $min = $array[0];\n+\tmy $max = $array[($count-1)];\n+\treturn ($average,$median,$min,$max,$first,$third,$sum);\n+}\n+\n+sub GlobalSummary {\n+\tmy ($bam,$targets) = @_;\n+\n+\tmy $command = "cd $wd && coverageBed -abam $bam -b $targets > $wd/Targets.Global.Coverage";\n+\tif (exists($commandsrun{$command})) {\n+\t\treturn;\n+\t}\n+\tsystem($command);\n+\t$commandsrun{$command} = 1;\n+}\n+\n+sub CoveragePerRegion {\n+\tmy ($bam,$targets) = @_;\n+\tmy $command = "cd $wd && coverageBed -abam $bam -b $targets > $wd/Targets.Global.Coverage";\n+\tif (exists($commandsrun{$command})) {\n+\t\treturn;\n+\t}\n+\tsystem($command);\n+\t$commandsrun{$command} = 1;\n+}\n+\n+sub SubRegionCoverage {\n+\tmy ($bam,$targets) = @_;\n+\tmy $command = "cd $wd && coverageBed -abam $bam -b $targets -d > $wd/Targets.Position.Coverage";\n+\tsystem($command);\n+\t$commandsrun{$command} = 1;\n+}\n+\n' |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 CoverageReport.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CoverageReport.xml Wed Oct 25 12:37:35 2017 -0400 |
| [ |
| @@ -0,0 +1,105 @@ +<tool id="CoverageReport2" name="Panel Coverage Report" version="0.0.4"> + <description>Provide mapping statistics</description> + <requirements> + <requirement type="package" version="1.75">perl-number-format</requirement> + <requirement type="package" version="3.4.1">R</requirement> + <requirement type="package" version="2.17.0">bedtools</requirement> + <requirement type="package" version="0.1.18">samtools</requirement> + <requirement type="package" version="0.1.6">tectonic</requirement> + </requirements> + <command detect_errors="exit_code"><![CDATA[ + perl '$__tool_directory__/CoverageReport.pl' + ## input files + -b '$input1' + -t '$input2' + + ## output files + -o '$output1' + + ## run parameters + $perGene + $PositionLevel + -m $threshold + -f $frac + ## sample name + -n '$input1.name' + ]]></command> + <inputs> + <param name="input1" type="data" format="bam" label="BAM file" help="BAM file of mapped reads" /> + <param name="input2" type="data" format="bed" label="Target Regions BED" help="BED file containing regions of interest. See below for format" /> + <param name="threshold" type="integer" value="40" label="Minimal Coverage Threshold" help="Default: 40" /> + <param name="frac" type="float" value="0.2" label="Fraction of Average Coverage for usage in plot" help="Default: 0.2" /> + <param name="perGene" type="boolean" truevalue="-r" falsevalue="" checked="True" label="Plot exon coverages for all genes in targets"/> + <param name="PositionLevel" type="select" label="Perform Per Exon Analysis" help="Only Failed: Only those exons not reaching global coverage above threshold, or 100%. All Exons: This can take a very long time for large panels! Select all failed to check all exons for local failures." > + <option value='' selected="true">None</option> + <option value='-s'>Plot Only Globally Failed</option> + <option value='-S'>Plot All Failed Exons</option> + <option value='-A'>Plot All Exons</option> + <option value='-L'>List All Failed Exons</option> + </param> + </inputs> + <outputs> + <data format="pdf" name="output1" label="${tool.name} on ${on_string}: PDF Report"/> + </outputs> + <tests> + <test> + <param name="input1" value="ecoli.bam" /> + <param name="input2" value="ecoli_6.bed" /> + <output name="output1" value="answer.pdf" compare="sim_size"/> + </test> + </tests> + <help><![CDATA[ +**What it does** + +This tool creates a coverage report for QC purposes. By default, average coverage statistics are provided, taken from samtools flagstats. If specified, it can also create overviews per gene in the BED file, and sub-exon plots for failed exons. + +------ + +**BED format** + +The BED file containing targets of interest has very specific format requirements. You **must** use the following format:: + + Column 1: Chromosome : Use the same syntax as the references used by Galaxy. Check your sam-headers for the correct format. ('chr1' vs '1') + Column 2: Start Position + Column 3: End Position + Column 4: Target Name. Use : "GENE-NAME<space>Exon_number" : This is split on the space after 'GeneName' for correct grouping. + Column 5: Score : ignored, use '0' + Column 6: Strand: ignored,'+' or '-' + +.. class:: infomark + +Note: The exons for the plots will be ordered in the same way as the exons in the BED file. + +------ + +**Input formats** + +BAM file for reads, BED file for targets. + +------ + +**Outputs** + +The output files are a PDF report. +The output tables are (tab seperated txt files): + +**Targets.Global.Coverage** : Original BED file + following columns:: + - Total coverage in target + - Bases in target with coverage + - Length of target + - Percent of target covered + +**Targets.Position.Coverage** : Original BED file + following columns:: + - Position in target region + - Coverage at position + + ]]></help> + <citations> + <citation type="bibtex">@ARTICLE{a1, + title = {This tool creates a coverage report for QC purposes.}, + author = {Philip Mabon}, + url = {https://github.com/galaxyproject/tools-iuc/} + } + }</citation> + </citations> +</tool> |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 README --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README Wed Oct 25 12:37:35 2017 -0400 |
| b |
| @@ -0,0 +1,14 @@ +Panel Resequencing Report +========================= + +Original wrapper was created by Geert Vandeweyer. It is free to use, adapt and redistribute, but please mention the original source. + +Current author is Philip Mabon (Takadonet) + + +History +======= +0.0.1 - Initial release, no dependency handling +0.0.2 - Updated release, following toolshed recommendations. +0.0.4 - Changing over to bioconda for samtools,bedtools,tectonic and perl modules. + |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 test-data/answer.pdf |
| b |
| Binary file test-data/answer.pdf has changed |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 test-data/ecoli.bam |
| b |
| Binary file test-data/ecoli.bam has changed |
| b |
| diff -r 000000000000 -r 30f8f85e3f98 test-data/ecoli_6.bed --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/ecoli_6.bed Wed Oct 25 12:37:35 2017 -0400 |
| b |
| @@ -0,0 +1,1 @@ +gi|49175990|ref|NC_000913.2|_Escherichia_coli_str._K-12_substr._MG1655,_complete_genome,_cropped_to_first_1000_nucleotides 1 1000 gi|49175990|ref|NC_000913.2|_Escherichia_coli_str._K-12_substr._MG1655,_complete_genome,_cropped_to_first_1000_nucleotides 5 6 |