Mercurial > repos > drosofff > mississippi_gcc

changeset 4:2025f1b016be draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mvdbeek
date Tue, 24 Jun 2014 09:49:51 -0400
parents 5dd9f42f72db
children 457746c01695
files mississippi_gcc/readmap.xml
diffstat 1 files changed, 225 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mississippi_gcc/readmap.xml	Tue Jun 24 09:49:51 2014 -0400
@@ -0,0 +1,225 @@
+<tool id="Readmap" name="Generate readmap and histograms from alignment files" version="0.9.2">
+  <description>from sRbowtie aligment</description>
+  <requirements><requirement type='package'>bowtie-inspect</requirement></requirements>
+  <parallelism method="basic"></parallelism>
+<command interpreter="python">
+        readmap.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_readmap
+		  $readmap_dataframe
+		  --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
+
+</command>
+  <inputs>
+       <conditional name="refGenomeSource">
+           <param name="genomeSource" type="select" label="Will you select a reference genome from your history or use a built-in index?" help="Built-ins were indexed using default options">
+               <option value="indexed">Use a built-in index</option>
+               <option value="history">Use one from the history</option>
+           </param>
+           <when value="indexed">
+	     <repeat name="series" title="Add alignment files">
+	       <param name="input" type="data" label="Select multiple alignments to parse">
+                  <validator type="dataset_metadata_in_data_table" table_name="bowtie_indexes" metadata_name="dbkey" metadata_column="0" message="database not set for this bowtie output. Select the database(=genome used for matching) manually, or select a reference fasta from your history."/>
+               </param>
+	       <param name="norm" type="float" value="1" label="Indicate a normalization factor to compare multiple aligments"/>
+	     </repeat>
+           </when>
+           <when value="history">
+	     <repeat name="series" title="Add alignment files">
+	       <param name="input" type="data" label="Select multiple alignments to parse"/>
+	       <param name="norm" type="integer" value="1" label="Indicate a normalization factor to compare multiple aligments"/>
+	     </repeat>
+	   </when>
+       </conditional>
+                <param name="gff" type="data" optional="true" label="Optional: select a GFF to investigate regions of interest" help="GFF must match genome build"/>
+                 <!-- <validator type="dataset_metadata_in_data_table" table_name="bowtie_indexes" metadata_name="dbkey" metadata_column="0" message="GFF database and alignment file databse do not match!"/> -->
+                <param name="minquery" type="integer" size="3" value="18" label="Min size of query small RNAs" help="'18' = 18 nucleotides"/>
+                <param name="maxquery" type="integer" size="3" value="28" label="Max size of query small RNAs" help="'28' = 28 nucleotides"/>
+                <param name="title" type="text" size="15" value= "Readmaps and size distributions" label="Main Titles"/>
+                <param name="xlabel" type="text" size="15" value="Coordinates/read size" label="x axis label"/>
+                <param name="ylabel" type="text" size="15" value="Number of reads" label="y axis label"/>
+                <param name="rows_per_page" type="text" size="9" value="8" label="How many items to display per page?">
+		  <validator type="in_range" min="6" max="20" message="Select between 6 and 20 rows, as the readability will suffer otherwise."/>
+                </param>
+  </inputs>
+   <configfiles>
+     <configfile name="plotCode">
+      ## 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
+
+      rm=read.delim("${readmap_dataframe}", header=T, row.names=NULL)
+      pdf(file="${readmap_PDF}", paper="special", height=11.69, width=8.2677)
+      n_samples=length(unique(rm\$sample))
+
+      genes=unique(levels(rm\$gene))
+      per_gene_readmap=lapply(genes, function(x) subset(rm, gene==x))
+      n_genes=length(per_gene_readmap)
+      
+
+      par.settings.readmap=list(layout.heights=list(top.padding=0, bottom.padding=-3), fontsize = list(text=96/${rows_per_page}, points=8))
+      par.settings.size=list(layout.heights=list(top.padding=-1, bottom.padding=-3), fontsize = list(text=96/${rows_per_page}, points=8))
+      par.settings.combination.readmap=list(layout.heights=list(top.padding=0, bottom.padding=-3), fontsize = list(text=96/${rows_per_page}, points=8))
+      par.settings.combination.size=list(layout.heights=list(top.padding=-2, bottom.padding=-0.5), fontsize = list(text=96/${rows_per_page}, points=8))
+      
+
+      plot_readmap=function(df, ...) {
+      combineLimits(xyplot(count~coord|factor(sample, levels=unique(sample))+reorder(gene, count, function(x) -sum(abs(x))), 
+      data=df, 
+      type='h', 
+      scales= list(relation="free", x=list(rot=0, cex=0.75, axs="i", tck=0.5), y=list(tick.number=4, rot=90, cex=0.75)),
+      xlab=NULL, main=NULL, ylab=NULL, 
+      as.table=T, 
+      origin = 0, 
+      horizontal=FALSE, 
+      group=polarity,
+      col=c("red","blue"),
+      ...))
+      }
+
+      plot_size_distribution= function(df, ...) {
+          smR.prepanel=function(x,y,...){; yscale=c(-max(abs(y)), max(abs(y)));list(ylim=yscale);}
+         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'),
+          cex=0.75,
+          scales=list(y=list(tick.number=4, rot=90, relation="free"), cex=0.75),
+          prepanel=smR.prepanel,
+          xlab = NULL,
+          ylab = NULL,
+#          par.settings=list(layout.heights=list(top.padding=-2, bottom.padding=-3), fontsize = list(text=8, points=8)),
+          main = NULL , as.table=TRUE, newpage = T, ...)
+          combineLimits(bc)
+          }
+
+      for (i in seq(1,n_genes,${rows_per_page})) {
+        start=i
+        end=i+${rows_per_page}-1
+        if (end>n_genes) {end=n_genes}
+	readmap_plot.list=lapply(per_gene_readmap[start:end], function(x) plot_readmap(x, par.settings=par.settings.readmap))
+	args.list=c(readmap_plot.list, list(nrow=${rows_per_page}, ncol=1, main="${title}", left="${ylabel}", sub="readmap coordinate"))
+        do.call(grid.arrange, args.list)
+      }
+
+      devname=dev.off()
+
+      size=read.delim("${size_distribution_dataframe}", header=T, row.names=NULL)
+      per_gene_size=lapply(genes, function(x) subset(size, gene==x))
+
+      pdf(file="${size_PDF}", paper="special", height=11.69, width=8.2677)
+
+      for (i in seq(1,n_genes,${rows_per_page})) {
+        start=i
+        end=i+${rows_per_page}-1
+        if (end>n_genes) {end=n_genes}
+        plot.list=lapply(per_gene_size[start:end], function(x) plot_size_distribution(x, par.settings=par.settings.size))
+        args.list=c(plot.list, list(nrow=${rows_per_page}, ncol=1, main="${title}", left="${ylabel}", sub="readsize in nucleotides"))
+        do.call(grid.arrange, args.list)
+      }
+
+      devname=dev.off()
+
+      pdf(file="${combi_PDF}", paper="special", height=11.69, width=8.2677)
+
+      for (i in seq(1,n_genes,${rows_per_page}/2)) {
+        start=i
+        end=i+${rows_per_page}/2-1
+        if (end>n_genes) {end=n_genes}
+	read_plot.list=lapply(per_gene_readmap[start:end], function(x) plot_readmap(x, par.settings=par.settings.combination.readmap))
+        size_plot.list=lapply(per_gene_size[start:end], function(x) plot_size_distribution(x, strip=FALSE, par.settings=par.settings.combination.size))
+	plot.list=rbind(read_plot.list, size_plot.list )
+        args.list=c(plot.list, list(nrow=${rows_per_page}, ncol=1, main="${title}", left="${ylabel}", sub="${xlabel}"))
+        do.call(grid.arrange, args.list)
+      }
+
+      devname=dev.off()
+
+
+     </configfile>
+   </configfiles>
+
+   <outputs>
+   <data format="tabular" name="readmap_dataframe" label="Readmap dataframe"/>
+   <data format="tabular" name="size_distribution_dataframe" label="Size distributionn dataframe"/>
+   <data format="pdf" name="readmap_PDF" label="Readmaps"/>
+   <data format="pdf" name="size_PDF" label="Size distribution"/>
+   <data format="pdf" name="combi_PDF" label="Size distributiona and Readmaps"/>
+   </outputs>
+<help>
+
+**What it does**
+
+Takes one or more alignment files (BAM, SAM or tabular bowtie output) as input and produces a "Readmap", 
+where by default for each "chromosome" the position of the read is recorded on the x-axis, and the y-axis indicates 
+the number of reads per position. Reads that map in sense are on the top, reads that map antisense are on the bottom.
+
+
+.. 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/readmap.png 
+    :height: 800 
+    :width: 500
+
+</help>
+</tool>