changeset 3:59c387dcb2d2 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mutect_wrapper.py	Thu Feb 13 09:27:12 2014 -0500
@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+
+"""
+Runs muTect on normal/tumor bam pair.
+
+usage: mutect_wrapper.py [options]
+
+See below for options
+"""
+
+import optparse, os, shutil, subprocess, sys, tempfile
+GALAXY_EXT_TO_GATK_EXT = { 'gatk_interval':'intervals', 'bam_index':'bam.bai', 'gatk_dbsnp':'dbSNP', 'picard_interval_list':'interval_list' } #items not listed here will use the galaxy extension as-is
+GALAXY_EXT_TO_GATK_FILE_TYPE = GALAXY_EXT_TO_GATK_EXT #for now, these are the same, but could be different if needed
+DEFAULT_GATK_PREFIX = "gatk_file"
+CHUNK_SIZE = 2**20 #1mb
+
+def stop_err( msg ):
+    sys.stderr.write( '%s\n' % msg )
+    sys.exit()
+
+
+def __main__():
+    tmp_dir = tempfile.mkdtemp( prefix='tmp-muTect-' )
+    print tmp_dir
+    #Parse Command Line
+    parser = optparse.OptionParser()
+    parser.add_option( '-R', '--ref', dest='ref', help='The reference genome to use ' )
+    parser.add_option( '-n', '--input_normal', dest='normal',  help='The Normal Tissue BAM file' )
+    ##parser.add_option( '','--index_normal',dest='index_normal',help='index of normal bam file')
+    parser.add_option( '-t', '--input_tumor', dest='tumor',  help='The Tumor Tissue BAM file' )
+    ##parser.add_option( '','--index_tumor',dest='index_tumor',help='index of tumor bam file')
+    parser.add_option( '', '--callstats', dest='callstats', help='The file to save the call statistics (txt format)' )
+    parser.add_option( '', '--coverage', dest='coverage', help='The file to save the coverage wig (wig format)' )
+    parser.add_option( '', '--vcf', dest='vcf',help='output VCF file.')
+    parser.add_option( '', '--cosmic', dest='cosmic',help='COSMIC VCF file')
+    parser.add_option( '', '--dbsnp',dest='dbsnp',help='dbSNP VCF file')
+    parser.add_option( '-L','--intervals',dest='intervals',help='Interval file (-L)')
+    parser.add_option( '-T', '--analysis_type', dest='analysis', help='Analysis Type(default = MuTect)' )
+    parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)',default='pre_set' )
+    parser.add_option( '-j', '--jar', dest='jar', help='path to the mutect jar file.')
+    parser.add_option( '', '--artifact', dest='artifact', help='Disable HC filter')
+    parser.add_option( '', '--pon', dest='pon', help='Panel of Normal (VCF file)')
+    
+    (options, args) = parser.parse_args()
+
+
+    # make temp directory for placement of indices
+    #tmp_index_dir = tempfile.mkdtemp()
+    #tmp_dir = tempfilte.mkdtemp()
+    # set runtime arguments
+    # input
+    normalpath = "%s/normal.bam" %(tmp_dir)
+    normalidx = "%s/normal.bai" %( tmp_dir)
+    tumorpath = "%s/tumor.bam" %(tmp_dir)
+    tumoridx = "%s/tumor.bai" %(tmp_dir)
+    os.symlink(options.normal, normalpath)
+    os.symlink(options.tumor, tumorpath)
+    #os.symlink(options.index_normal, normalidx)
+    #os.symlink(options.index_tumor, tumoridx)
+    try:
+	ic = "cd %s && samtools index normal.bam" % (tmp_dir)
+	os.system(ic)
+        ic = "cd %s && samtools index tumor.bam" % (tmp_dir)
+	os.system(ic)
+    except:
+	raise "indexing of bam files failed. "+str(e)
+
+    arguments = '-T %s --enable_extended_output --input_file:normal %s --input_file:tumor %s -R %s --dbsnp %s' % ( options.analysis, normalpath, tumorpath, options.ref, options.dbsnp)
+    if options.cosmic:
+	arguments += ' --cosmic %s' %(options.cosmic)
+    if options.intervals:
+	if os.path.isfile(options.intervals):
+		intervalpath = "%s/intervals.bed" %(tmp_dir)
+		os.symlink(options.intervals, intervalpath)
+		arguments += ' -L %s' %(intervalpath)
+	elif options.intervals != 'None' :
+		arguments += ' -L %s' %(options.intervals)
+    ## disable HC filters
+    if options.artifact:
+	arguments =+ ' --artifact_detection_mode'
+    ## enable PON
+    if options.pon:
+	arguments =+ ' --normal_panel %s ' %(options.pon)
+    # output
+    arguments += ' --out %s --coverage_file %s --vcf %s' % (options.callstats, options.coverage,options.vcf)
+
+    # all parameters set?
+    #if options.params != 'pre_set':
+	
+    # final command
+    command = 'java -Xmx2G -jar %s %s ' % (options.jar, arguments)
+    #print command 
+    try:
+	os.system(command)
+    except:
+	raise "muTect Failed" + str(e)
+    #if tmp_dir and os.path.exists( tmp_dir ):
+    #    shutil.rmtree( tmp_dir )
+
+	 
+if __name__=="__main__": __main__()