details: http://www.bx.psu.edu/hg/galaxy/rev/0cb12b9d9d4d changeset: 2658:0cb12b9d9d4d user: Kelly Vincent <kpvincent@bx.psu.edu> date: Tue Sep 01 17:01:44 2009 -0400 description: Added several options to the Samtools pileup tool 2 file(s) affected in this change: tools/samtools/sam_pileup.py tools/samtools/sam_pileup.xml diffs (107 lines): diff -r a2d69372d699 -r 0cb12b9d9d4d tools/samtools/sam_pileup.py --- a/tools/samtools/sam_pileup.py Tue Sep 01 14:20:49 2009 -0400 +++ b/tools/samtools/sam_pileup.py Tue Sep 01 17:01:44 2009 -0400 @@ -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 @@ 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 -r a2d69372d699 -r 0cb12b9d9d4d tools/samtools/sam_pileup.xml --- a/tools/samtools/sam_pileup.xml Tue Sep 01 14:20:49 2009 -0400 +++ b/tools/samtools/sam_pileup.xml Tue Sep 01 17:01:44 2009 -0400 @@ -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" />