comparison mutect_wrapper.py @ 3:59c387dcb2d2 draft

Uploaded

3:59c387dcb2d2

#!/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__()