From 1008d8692f5c6a6f7d13d11f8eacb3babbbf7dc8 Mon Sep 17 00:00:00 2001
From: Kelly Vincent <kpvincent@bx.psu.edu>
Date: Tue, 1 Sep 2009 17:01:44 -0400
Subject: [PATCH] Added several options to the Samtools pileup tool

---
 tools/samtools/sam_pileup.py  | 22 +++++++++++++++------
 tools/samtools/sam_pileup.xml | 37 +++++++++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+), 6 deletions(-)

diff --git a/tools/samtools/sam_pileup.py b/tools/samtools/sam_pileup.py
index f84bd90d6e5..4e11f48d61d 100644
--- a/tools/samtools/sam_pileup.py
+++ b/tools/samtools/sam_pileup.py
@@ -4,15 +4,22 @@
 Creates a pileup file from a bam file and a reference.
 
 usage: %prog [options]
-   -i, --input1=i: bam file
+   -p, --input1=p: bam file
    -o, --output1=o: Output pileup
-   -r, --ref=r: Reference file type
+   -R, --ref=R: Reference file type
    -n, --ownFile=n: User-supplied fasta reference file
    -d, --dbkey=d: dbkey of user-supplied file
    -x, --indexDir=x: Index directory
    -b, --bamIndex=b: BAM index file
+   -s, --lastCol=s: Print the mapping quality as the last column
+   -i, --indels=i: Only output lines containing indels
+   -M, --mapCap=M: Cap mapping quality
+   -c, --consensus=c: Call the consensus sequence using MAQ consensu model
+   -T, --theta=T: Theta paramter (error dependency coefficient)
+   -N, --hapNum=N: Number of haplotypes in sample
+   -r, --fraction=r: Expected fraction of differences between a pair of haplotypes
+   -I, --phredProb=I: Phred probability of an indel in sequencing/prep
 
-usage: %prog input1 output1 ref_type refFile ownFile dbkey index_dir bam_index
 """
 
 import os, sys, tempfile
@@ -49,17 +56,20 @@ def __main__():
     tmpf0bambai = '%s.bam.bai' % tmpf0.name
     tmpf1 = tempfile.NamedTemporaryFile(dir=tmp_dir)
     tmpf1fai = '%s.fai' % tmpf1.name
+    opts = '%s %s -M %s' % (('','-s')[options.lastCol=='yes'], ('','-i')[options.indels=='yes'], options.mapCap)
+    if options.consensus == 'yes':
+        opts += ' -c -T %s -N %s -r %s -I %s' % (options.theta, options.hapNum, options.fraction, options.phredProb)
     cmd1 = None
     cmd2 = 'cp %s %s; cp %s %s' % (options.input1, tmpf0bam, options.bamIndex, tmpf0bambai)
-    cmd3 = 'samtools pileup -f %s %s > %s 2> /dev/null'
+    cmd3 = 'samtools pileup %s -f %s %s > %s 2> /dev/null'
     if options.ref =='indexed':
         full_path = "%s.fai" % seq_path 
         if not os.path.exists( full_path ):
             stop_err( "No sequences are available for '%s', request them by reporting this error." % options.dbkey )
-        cmd3 = cmd3 % (seq_path, tmpf0bam, options.output1)
+        cmd3 = cmd3 % (opts, seq_path, tmpf0bam, options.output1)
     elif options.ref == 'history':
         cmd1 = 'cp %s %s; cp %s.fai %s' % (options.ownFile, tmpf1.name, options.ownFile, tmpf1fai)
-        cmd3 = cmd3 % (tmpf1.name, tmpf0bam, options.output1)
+        cmd3 = cmd3 % (opts, tmpf1.name, tmpf0bam, options.output1)
     # index reference if necessary
     if cmd1:
         try:
diff --git a/tools/samtools/sam_pileup.xml b/tools/samtools/sam_pileup.xml
index 4054d1f1d56..90e0f5844ec 100644
--- a/tools/samtools/sam_pileup.xml
+++ b/tools/samtools/sam_pileup.xml
@@ -13,6 +13,21 @@
 	    --dbkey=${input1.metadata.dbkey}
         --indexDir=${GALAXY_DATA_INDEX_DIR}
         --bamIndex=${input1.metadata.bam_index}
+        --lastCol=$lastCol
+        --indels=$indels
+        --mapCap=$mapCap
+        --consensus=$c.consensus
+        #if $c.consensus == "yes":
+         --theta=$c.theta
+         --hapNum=$c.hapNum
+         --fraction=$c.fraction
+         --phredProb=$c.phredProb
+        #else:
+         --theta="None"
+         --hapNum="None"
+         --fraction="None"
+         --phredProb="None"
+        #end if
   </command>
   <inputs>
     <conditional name="refOrHistory">
@@ -31,6 +46,28 @@
         <param name="ownFile" type="data" format="fasta" metadata_name="dbkey" label="Select a reference genome" />
       </when>
     </conditional>	
+    <param name="lastCol" type="select" label="Whether or not to print the mapping quality as the last column" help="Makes the output easier to parse, but is space inefficient">
+      <option value="no">Do not print the mapping quality as the last column</option>
+      <option value="yes">Print the mapping quality as the last column</option>
+    </param>
+    <param name="indels" type="select" label="Whether or not to print only output pileup lines containing indels">
+      <option value="no">Print all lines</option>
+      <option value="yes">Print only lines containing indels</option>
+    </param>
+    <param name="mapCap" type="integer" value="60" label="Where to cap mapping quality" />
+    <conditional name="c">
+      <param name="consensus" type="select" label="Whether or not to call the consensus sequence using the MAQ consensus model">
+        <option value="no">Don't use MAQ consensus model</option>
+        <option value="yes">Use the MAQ consensus model</option>
+      </param> 
+      <when value="no" />
+      <when value="yes">
+        <param name="theta" type="float" value="0.85" label="Theta paramter (error dependency coefficient) in the MAQ consensus calling model" />
+        <param name="hapNum" type="integer" value="2" label="Number of haplotypes in the sample" help="Greater than or equal to 2" />
+        <param name="fraction" type="float" value="0.001" label="Expected fraction of differences between a pair of haplotypes" />
+        <param name="phredProb" type="integer" value="40" label="Phred probability of an indel in sequencing/prep" />
+      </when>
+    </conditional>
   </inputs>
   <outputs>
     <data format="tabular" name="output1" />