|
0
|
1
|
|
14
|
2 import argparse, os, shutil, subprocess, sys, tempfile, time, shlex
|
|
0
|
3
|
|
|
4 parser = argparse.ArgumentParser(description='')
|
|
8
|
5 parser.add_argument('-r', dest='inputFastaFile', required=True, help='the reference file')
|
|
|
6 parser.add_argument('-b', dest='inputBamFile', required=True, help='the bam file')
|
|
14
|
7 parser.add_argument('-s', dest='insert_size', type=int, default=None, required=False, help='the insert size')
|
|
|
8 parser.add_argument('-t', dest='sampleTag', default='default', required=False, help='the sample tag')
|
|
8
|
9 # parser.add_argument('-o1', dest='outputRaw', help='the output raw', required=True)
|
|
|
10 parser.add_argument('-o2', dest='outputVcfFile', help='the output vcf', required=True)
|
|
0
|
11 parser.add_argument('--number_of_threads', dest='number_of_threads', type=int, default='1')
|
|
|
12 parser.add_argument('--window_size', dest='window_size', type=int, default='5')
|
|
|
13 parser.add_argument('--sequencing_error_rate', dest='sequencing_error_rate', type=float, default='0.01')
|
|
|
14 parser.add_argument('--sensitivity', dest='sensitivity', default='0.95', type=float)
|
|
|
15 parser.add_argument('--report_long_insertions', dest='report_long_insertions', action='store_true', default=False)
|
|
|
16 parser.add_argument('--report_duplications', dest='report_duplications', action='store_true', default=False)
|
|
|
17 parser.add_argument('--report_inversions', dest='report_inversions', action='store_true', default=False)
|
|
|
18 parser.add_argument('--report_breakpoints', dest='report_breakpoints', action='store_true', default=False)
|
|
|
19
|
|
9
|
20 parser.add_argument('-u', '--maximum_allowed_mismatch_rate', dest='maximum_allowed_mismatch_rate', type=float, default='0.02')
|
|
8
|
21 parser.add_argument('--report_close_mapped_reads', dest='report_close_mapped_reads', action='store_true', default=False)
|
|
|
22 parser.add_argument('--report_only_close_mapped_reads', dest='report_only_close_mapped_reads', action='store_true', default=False)
|
|
|
23 parser.add_argument('--report_interchromosomal_events', dest='report_interchromosomal_events', action='store_true', default=False)
|
|
|
24 parser.add_argument('--IndelCorrection', dest='IndelCorrection', action='store_true', default=False)
|
|
|
25 parser.add_argument('--NormalSamples', dest='NormalSamples', action='store_true', default=False)
|
|
9
|
26 parser.add_argument('-a', '--additional_mismatch', dest='additional_mismatch', type=int, default='1')
|
|
|
27 parser.add_argument('-m', '--min_perfect_match_around_BP', dest='min_perfect_match_around_BP', type=int, default='3')
|
|
|
28 parser.add_argument('-v', '--min_inversion_size', dest='min_inversion_size', type=int, default='50')
|
|
|
29 parser.add_argument('-d', '--min_num_matched_bases', dest='min_num_matched_bases', type=int, default='30')
|
|
|
30 parser.add_argument('-B', '--balance_cutoff', dest='balance_cutoff', type=int, default='0')
|
|
|
31 parser.add_argument('-A', '--anchor_quality', dest='anchor_quality', type=int, default='0')
|
|
|
32 parser.add_argument('-M', '--minimum_support_for_event', dest='minimum_support_for_event', type=int, default='3')
|
|
14
|
33 parser.add_argument('-n', '--NM', dest='NM', type=int, default='2')
|
|
8
|
34 parser.add_argument('--detect_DD', dest='detect_DD', action='store_true', default=False)
|
|
|
35 parser.add_argument('--MAX_DD_BREAKPOINT_DISTANCE', dest='MAX_DD_BREAKPOINT_DISTANCE', type=int, default='350')
|
|
|
36 parser.add_argument('--MAX_DISTANCE_CLUSTER_READS', dest='MAX_DISTANCE_CLUSTER_READS', type=int, default='100')
|
|
|
37 parser.add_argument('--MIN_DD_CLUSTER_SIZE', dest='MIN_DD_CLUSTER_SIZE', type=int, default='3')
|
|
|
38 parser.add_argument('--MIN_DD_BREAKPOINT_SUPPORT', dest='MIN_DD_BREAKPOINT_SUPPORT', type=int, default='3')
|
|
|
39 parser.add_argument('--MIN_DD_MAP_DISTANCE', dest='MIN_DD_MAP_DISTANCE', type=int, default='8000')
|
|
|
40 parser.add_argument('--DD_REPORT_DUPLICATION_READS', dest='DD_REPORT_DUPLICATION_READS', action='store_true', default=False)
|
|
0
|
41
|
|
14
|
42 parser.add_argument('-z', '--input_SV_Calls_for_assembly', dest='input_SV_Calls_for_assembly', action='store_true', default=False)
|
|
0
|
43
|
|
8
|
44
|
|
|
45 def execute(cmd, output=None):
|
|
|
46 import subprocess, sys, shlex
|
|
|
47 # function to execute a cmd and report if an error occur
|
|
|
48 print(cmd)
|
|
0
|
49 try:
|
|
8
|
50 process = subprocess.Popen(args=shlex.split(cmd), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
|
|
0
|
51 process.wait()
|
|
8
|
52 stdout,stderr = process.communicate()
|
|
|
53 except Exception, e: # une erreur de ma commande : stderr
|
|
|
54 sys.stderr.write("problem doing : %s\n%s\n" %(cmd, e))
|
|
|
55 return
|
|
|
56 if output:
|
|
|
57 output = open(output, 'w')
|
|
|
58 output.write(stdout)
|
|
|
59 output.close()
|
|
|
60 if stderr != '': # une erreur interne au programme : stdout (sinon, souvent des warning arrete les programmes)
|
|
|
61 sys.stdout.write("warning or error while doing : %s\n-----\n%s-----\n\n" %(cmd, stderr))
|
|
0
|
62
|
|
|
63
|
|
8
|
64 def indexBam(inputFastaFile, inputBamFile):
|
|
10
|
65 cmd = "samtools index %s" %(inputBamFile)
|
|
0
|
66 execute(cmd)
|
|
8
|
67 cmd = "samtools faidx %s" %(inputFastaFile)
|
|
0
|
68 execute(cmd)
|
|
|
69
|
|
|
70
|
|
14
|
71 def config(inputBamFile, meanInsertSize, tag, tempDir):
|
|
8
|
72 configFile = tempDir+"/pindel_configFile"
|
|
|
73 fil = open(configFile, 'w')
|
|
14
|
74 fil.write("%s\t%s\t%s" %(inputBamFile, meanInsertSize, tag))
|
|
0
|
75 fil.close()
|
|
8
|
76 return configFile
|
|
0
|
77
|
|
|
78
|
|
8
|
79 def pindel(reference, configFile, args, tempDir):
|
|
|
80 pindelTempDir=tempDir+"/pindel"
|
|
|
81 cmd = "pindel -f %s -i %s -o %s " %(reference, configFile, pindelTempDir)
|
|
|
82 cmd += " --number_of_threads %d --window_size %d --sequencing_error_rate %f --sensitivity %f" %(args.number_of_threads, args.window_size, args.sequencing_error_rate, args.sensitivity)
|
|
|
83 cmd += " -u %f -n %d -a %d -m %d -v %d -d %d -B %d -A %d -M %d " %(args.maximum_allowed_mismatch_rate, args.NM, args.additional_mismatch, args.min_perfect_match_around_BP, args.min_inversion_size, args.min_num_matched_bases, args.balance_cutoff, args.anchor_quality, args.minimum_support_for_event)
|
|
|
84
|
|
0
|
85 if args.report_long_insertions:
|
|
|
86 cmd += ' --report_long_insertions '
|
|
|
87 if args.report_duplications:
|
|
|
88 cmd += ' --report_duplications '
|
|
|
89 if args.report_inversions:
|
|
|
90 cmd += ' --report_inversions '
|
|
|
91 if args.report_breakpoints:
|
|
|
92 cmd += ' --report_breakpoints '
|
|
8
|
93
|
|
|
94 if args.report_close_mapped_reads:
|
|
|
95 cmd += ' --report_close_mapped_reads '
|
|
|
96 if args.report_only_close_mapped_reads:
|
|
|
97 cmd += ' --report_only_close_mapped_reads '
|
|
|
98 if args.report_interchromosomal_events:
|
|
|
99 cmd += ' --report_interchromosomal_events '
|
|
|
100 if args.IndelCorrection:
|
|
|
101 cmd += ' --IndelCorrection '
|
|
|
102 if args.NormalSamples:
|
|
|
103 cmd += ' --NormalSamples '
|
|
|
104 if args.input_SV_Calls_for_assembly:
|
|
|
105 cmd += ' --input_SV_Calls_for_assembly %s ' %(args.input_SV_Calls_for_assembly)
|
|
|
106
|
|
|
107 if args.detect_DD:
|
|
|
108 cmd += ' -q '
|
|
|
109 cmd += ' --MAX_DD_BREAKPOINT_DISTANCE '+str(args.MAX_DD_BREAKPOINT_DISTANCE)
|
|
|
110 cmd += ' --MAX_DISTANCE_CLUSTER_READS '+str(args.MAX_DISTANCE_CLUSTER_READS)
|
|
|
111 cmd += ' --MIN_DD_CLUSTER_SIZE '+str(args.MIN_DD_CLUSTER_SIZE)
|
|
|
112 cmd += ' --MIN_DD_BREAKPOINT_SUPPORT '+str(args.MIN_DD_BREAKPOINT_SUPPORT)
|
|
|
113 cmd += ' --MIN_DD_MAP_DISTANCE '+str(args.MIN_DD_MAP_DISTANCE)
|
|
|
114 if args.DD_REPORT_DUPLICATION_READS:
|
|
|
115 cmd += ' --DD_REPORT_DUPLICATION_READS '
|
|
|
116
|
|
|
117 execute(cmd)
|
|
|
118 return pindelTempDir
|
|
0
|
119
|
|
|
120
|
|
8
|
121 def move(avant, apres):
|
|
0
|
122 if os.path.exists(avant):
|
|
|
123 execute("mv %s %s" %(avant, apres))
|
|
|
124
|
|
|
125
|
|
8
|
126 def pindel2vcf(outputVcf, inputFastaFile, name, pindelTempDir):
|
|
0
|
127 date = str(time.strftime('%d/%m/%y',time.localtime()))
|
|
8
|
128 cmd = "pindel2vcf -P %s -r %s -R %s -d %s" %(pindelTempDir, inputFastaFile, name, date)
|
|
|
129 print("commande = " + cmd)
|
|
0
|
130 execute(cmd)
|
|
8
|
131 if os.path.exists(pindelTempDir+".vcf"):
|
|
|
132 execute("cp %s %s" %(pindelTempDir+".vcf", outputVcf))
|
|
|
133
|
|
|
134
|
|
14
|
135 # def getMeanInsertSize(bamFile, tempDir):
|
|
|
136 # samFile = tempDir+"/sam"
|
|
|
137 # histogramFile = tempDir+"/histogram"
|
|
|
138 # outputFile = tempDir+"/output"
|
|
|
139 # cmd = "samtools view -h %s -o %s" % (bamFile, samFile)
|
|
|
140 # execute(cmd, )
|
|
|
141 # cmd = "picard-tools CollectInsertSizeMetrics H=%s O=%s I=%s" % (histogramFile, outputFile, samFile)
|
|
|
142 # execute(cmd)
|
|
|
143 # try:
|
|
|
144 # f = open(outputFile, 'r')
|
|
|
145 # lignes = f.readlines()
|
|
|
146 # f.close()
|
|
|
147 # meanInsertSize = lignes[7].split('\t')[0]
|
|
|
148 # except:
|
|
|
149 # sys.stderr.write("problem while opening file %s " %(outputFile))
|
|
|
150 # return 0
|
|
|
151 # return meanInsertSize
|
|
0
|
152
|
|
|
153
|
|
|
154 def __main__():
|
|
|
155
|
|
|
156 args = parser.parse_args()
|
|
8
|
157 tempDir = tempfile.mkdtemp()
|
|
0
|
158
|
|
|
159 try:
|
|
8
|
160 indexBam(args.inputFastaFile, args.inputBamFile)
|
|
14
|
161 # if args.insert_size==None:
|
|
|
162 # meanInsertSize = int(getMeanInsertSize(args.inputBamFile, tempDir))
|
|
|
163 # else:
|
|
|
164 meanInsertSize=args.insert_size
|
|
|
165 configFile = config(args.inputBamFile, meanInsertSize, args.sampleTag, tempDir)
|
|
8
|
166 pindelTempDir = pindel(args.inputFastaFile, configFile, args, tempDir)
|
|
14
|
167 pindel2vcf(args.outputVcfFile, args.inputFastaFile, args.sampleTag, pindelTempDir)
|
|
|
168
|
|
0
|
169 finally:
|
|
8
|
170 if os.path.exists(tempDir):
|
|
|
171 shutil.rmtree(tempDir)
|
|
0
|
172
|
|
|
173 if __name__=="__main__":
|
|
|
174 __main__()
|
