# HG changeset patch # User mvdbeek # Date 1403617842 14400 # Node ID 241eb82a2e7b0dff37008ccf7fd8f2ff8a11ca72 # Parent 457746c0169556a03907a5ba8053b49a7fa6a651 Uploaded diff -r 457746c01695 -r 241eb82a2e7b mississippi_gcc/size_histogram.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mississippi_gcc/size_histogram.xml Tue Jun 24 09:50:42 2014 -0400 @@ -0,0 +1,200 @@ + + from sRbowtie aligment + bowtie-inspect + + + size_histogram.py + #if $refGenomeSource.genomeSource == "history": + --reference_fasta ## sys.argv[2] + $refGenomeSource.ownFile ## index source + #else: + #silent reference= filter( lambda x: str( x[0] ) == str( $refGenomeSource.series[0].input.dbkey ), $__app__.tool_data_tables[ 'bowtie_indexes' ].get_fields() )[0][-1] + --reference_bowtie_index + $reference + #end if + --rcode + $plotCode + --output_size_distribution + $size_distribution_dataframe + --minquery + $minquery + --maxquery + $maxquery + --input + #for $i in $refGenomeSource.series + $i.input + #end for + --ext + #for $i in $refGenomeSource.series + $i.input.ext + #end for + --label + #for $i in $refGenomeSource.series + "$i.input.name" + #end for + --normalization_factor + #for $i in $refGenomeSource.series + $i.norm + #end for + #if $gff: + --gff + $gff + #end if + #if $global.value == 'yes': + --global_size + #end if + #if $collapsestrands.value == 'yes': + --collapse + #end if + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ## Setup R error handling to go to stderr + options( show.error.messages=F, + error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) + library(RColorBrewer) + library(lattice) + library(latticeExtra) + library(grid) + library(gridExtra) + ##cheetahtemplate data frame implementation + + size=read.delim("${size_distribution_dataframe}", header=T, row.names=NULL) + + n_samples=length(unique(size\$sample)) + genes=unique(levels(size\$gene)) + n_genes=length(genes) + + par.settings.size=list(layout.heights=list(top.padding=-1, bottom.padding=-3, strip = .75), fontsize = list(text=96/${rows_per_page}, points=8)) + smR.prepanel=function(x,y,...){; yscale=c(-max(abs(y)), max(abs(y)));list(ylim=yscale);} + + plot_size_distribution= function(df, ...) { + bc= barchart(count~as.factor(size)|factor(sample, levels=unique(sample))+gene, data = df, origin = 0, + horizontal=FALSE, + group=polarity, + stack=TRUE, + col=c('red', 'blue'), + strip = strip.custom(par.strip.text = list(cex = 0.5)), + cex=0.75, + scales=list(y=list(tick.number=4, rot=90, relation="free"), cex=0.75), + xlab = "readsize in nucleotides", + ylab = "${ylabel}", + main="${title}" , + as.table=TRUE, newpage = T, ...) + combineLimits(update(useOuterStrips(bc), layout=c(n_samples,${rows_per_page})), margin.x=F, margin.y=1) + } + + per_gene_size=lapply(genes, function(x) subset(size, gene==x)) + + global = "no" + #if $global.value == 'yes': + global = "yes" + #end if + + if (global=="no") { + options(warn=-1) + pdf(file="${size_PDF}", paper="special", height=11.69, width=8.2677*n_samples/4) + plot_size_distribution(size, par.settings=par.settings.size, prepanel=smR.prepanel) + } else { + pdf(file="${size_PDF}", paper="special", height=11.69, width=8.2677) + bc= barchart(count~as.factor(size)|factor(sample, levels=unique(sample)), data = size, origin = 0, + horizontal=FALSE, + group=polarity, + stack=TRUE, + col=c('red', 'blue'), + cex=0.75, + par.settings=list(fontsize = list(text=8, points=8)), + scales=list(y=list(tick.number=4, rot=90, relation="same"), cex=0.75), + xlab = "readsize in nucleotides", + ylab = "${ylabel}", + main="${title}" , as.table=TRUE, newpage = T, + aspect=0.5) + #layout=c(n_samples, ${rows_per_page})) + bc + } + devname=dev.off() + + + + + + + + + + +**What it does** + +Takes one or more alignment files (BAM, SAM or tabular bowtie output) as input and produces a histogram of read sizes, +where by default for each "chromosome" a histogram of read sizes is drawn. +Reads that map in sense are on the top (red), reads that map antisense are on the bottom (blue). + + +.. class:: warningmark + +'''TIP''' The input data can be produced using the sRbowtie tool. + +---- + +'''Example''' + +Query sequence:: +For a SAM file as the following: + + 5 16 2L_79 24393 255 17M * 0 0 CCTTCATCTTTTTTTTT IIIIIIIIIIIIIIIII XA:i:0 MD:Z:17 NM:i:0 + + 11 0 2R_1 12675 255 21M * 0 0 AAAAAAAACGCGTCCTTGTGC IIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:21 NM:i:0 + + 2 16 2L_5 669 255 23M * 0 0 TGTTGCTGCATTTCTTTTTTTTT IIIIIIIIIIIIIIIIIIIIIII XA:i:0 MD:Z:23 NM:i:0 + +produce a plot like this: + +---- + +.. image:: static/images/size_histogram.png + :height: 800 + :width: 500 + + +