|
3
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 """
|
|
|
4 Runs muTect on normal/tumor bam pair.
|
|
|
5
|
|
|
6 usage: mutect_wrapper.py [options]
|
|
|
7
|
|
|
8 See below for options
|
|
|
9 """
|
|
|
10
|
|
|
11 import optparse, os, shutil, subprocess, sys, tempfile
|
|
|
12 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
|
|
|
13 GALAXY_EXT_TO_GATK_FILE_TYPE = GALAXY_EXT_TO_GATK_EXT #for now, these are the same, but could be different if needed
|
|
|
14 DEFAULT_GATK_PREFIX = "gatk_file"
|
|
|
15 CHUNK_SIZE = 2**20 #1mb
|
|
|
16
|
|
|
17 def stop_err( msg ):
|
|
|
18 sys.stderr.write( '%s\n' % msg )
|
|
|
19 sys.exit()
|
|
|
20
|
|
|
21
|
|
|
22 def __main__():
|
|
|
23 tmp_dir = tempfile.mkdtemp( prefix='tmp-muTect-' )
|
|
|
24 print tmp_dir
|
|
|
25 #Parse Command Line
|
|
|
26 parser = optparse.OptionParser()
|
|
|
27 parser.add_option( '-R', '--ref', dest='ref', help='The reference genome to use ' )
|
|
|
28 parser.add_option( '-n', '--input_normal', dest='normal', help='The Normal Tissue BAM file' )
|
|
|
29 ##parser.add_option( '','--index_normal',dest='index_normal',help='index of normal bam file')
|
|
|
30 parser.add_option( '-t', '--input_tumor', dest='tumor', help='The Tumor Tissue BAM file' )
|
|
|
31 ##parser.add_option( '','--index_tumor',dest='index_tumor',help='index of tumor bam file')
|
|
|
32 parser.add_option( '', '--callstats', dest='callstats', help='The file to save the call statistics (txt format)' )
|
|
|
33 parser.add_option( '', '--coverage', dest='coverage', help='The file to save the coverage wig (wig format)' )
|
|
|
34 parser.add_option( '', '--vcf', dest='vcf',help='output VCF file.')
|
|
|
35 parser.add_option( '', '--cosmic', dest='cosmic',help='COSMIC VCF file')
|
|
|
36 parser.add_option( '', '--dbsnp',dest='dbsnp',help='dbSNP VCF file')
|
|
|
37 parser.add_option( '-L','--intervals',dest='intervals',help='Interval file (-L)')
|
|
|
38 parser.add_option( '-T', '--analysis_type', dest='analysis', help='Analysis Type(default = MuTect)' )
|
|
|
39 parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)',default='pre_set' )
|
|
|
40 parser.add_option( '-j', '--jar', dest='jar', help='path to the mutect jar file.')
|
|
|
41 parser.add_option( '', '--artifact', dest='artifact', help='Disable HC filter')
|
|
|
42 parser.add_option( '', '--pon', dest='pon', help='Panel of Normal (VCF file)')
|
|
|
43
|
|
|
44 (options, args) = parser.parse_args()
|
|
|
45
|
|
|
46
|
|
|
47 # make temp directory for placement of indices
|
|
|
48 #tmp_index_dir = tempfile.mkdtemp()
|
|
|
49 #tmp_dir = tempfilte.mkdtemp()
|
|
|
50 # set runtime arguments
|
|
|
51 # input
|
|
|
52 normalpath = "%s/normal.bam" %(tmp_dir)
|
|
|
53 normalidx = "%s/normal.bai" %( tmp_dir)
|
|
|
54 tumorpath = "%s/tumor.bam" %(tmp_dir)
|
|
|
55 tumoridx = "%s/tumor.bai" %(tmp_dir)
|
|
|
56 os.symlink(options.normal, normalpath)
|
|
|
57 os.symlink(options.tumor, tumorpath)
|
|
|
58 #os.symlink(options.index_normal, normalidx)
|
|
|
59 #os.symlink(options.index_tumor, tumoridx)
|
|
|
60 try:
|
|
|
61 ic = "cd %s && samtools index normal.bam" % (tmp_dir)
|
|
|
62 os.system(ic)
|
|
|
63 ic = "cd %s && samtools index tumor.bam" % (tmp_dir)
|
|
|
64 os.system(ic)
|
|
|
65 except:
|
|
|
66 raise "indexing of bam files failed. "+str(e)
|
|
|
67
|
|
|
68 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)
|
|
|
69 if options.cosmic:
|
|
|
70 arguments += ' --cosmic %s' %(options.cosmic)
|
|
|
71 if options.intervals:
|
|
|
72 if os.path.isfile(options.intervals):
|
|
|
73 intervalpath = "%s/intervals.bed" %(tmp_dir)
|
|
|
74 os.symlink(options.intervals, intervalpath)
|
|
|
75 arguments += ' -L %s' %(intervalpath)
|
|
|
76 elif options.intervals != 'None' :
|
|
|
77 arguments += ' -L %s' %(options.intervals)
|
|
|
78 ## disable HC filters
|
|
|
79 if options.artifact:
|
|
|
80 arguments =+ ' --artifact_detection_mode'
|
|
|
81 ## enable PON
|
|
|
82 if options.pon:
|
|
|
83 arguments =+ ' --normal_panel %s ' %(options.pon)
|
|
|
84 # output
|
|
|
85 arguments += ' --out %s --coverage_file %s --vcf %s' % (options.callstats, options.coverage,options.vcf)
|
|
|
86
|
|
|
87 # all parameters set?
|
|
|
88 #if options.params != 'pre_set':
|
|
|
89
|
|
|
90 # final command
|
|
|
91 command = 'java -Xmx2G -jar %s %s ' % (options.jar, arguments)
|
|
|
92 #print command
|
|
|
93 try:
|
|
|
94 os.system(command)
|
|
|
95 except:
|
|
|
96 raise "muTect Failed" + str(e)
|
|
|
97 #if tmp_dir and os.path.exists( tmp_dir ):
|
|
|
98 # shutil.rmtree( tmp_dir )
|
|
|
99
|
|
|
100
|
|
|
101 if __name__=="__main__": __main__()
|
