-
Notifications
You must be signed in to change notification settings - Fork 0
/
Texomer.py
executable file
·365 lines (356 loc) · 14.4 KB
/
Texomer.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
'''
Created on May 23, 2017
@author: FWang9
'''
import sys,os,time
import subprocess
from optparse import OptionParser
from tempfile import gettempdir
def VarscanTrans(varscan,case):
data=open(varscan)
germlineout=open(case+".germline","w")
somaticout=open(case+".somatic","w")
out="chr"+"\t"+"position"+"\t"+"ref"+"\t"+"alt"+"\t"+"refNumN"+'\t'+"altNumN"+"\t"+"refNumT"+"\t"+"altNumT"
print >> germlineout,out
print >> somaticout,out
next(data)
for line in data:
line=line[0:-1].split("\t")
if "chr" not in line[0]:
chr="chr"+line[0]
else:
chr=line[0]
position=line[1]
ref=line[2]
var=line[3]
refNumN=line[4]
altNumN=line[5]
refNumT=line[8]
altNumT=line[9]
gt=line[7]
type=line[12]
pgerm=float(line[13])
psoma=float(line[14])
if type == "Germline" and gt not in ["A","C","G","T"] and pgerm < 0.01:
out=chr+"\t"+position+"\t"+ref+"\t"+var+"\t"+refNumN+"\t"+altNumN+"\t"+refNumT+"\t"+altNumT
print >> germlineout,out
if type == "LOH" and psoma < 0.01:
out=chr+"\t"+position+"\t"+ref+"\t"+var+"\t"+refNumN+"\t"+altNumN+"\t"+refNumT+"\t"+altNumT
print >> germlineout,out
if type == "Somatic" and psoma < 0.01:
out=chr+"\t"+position+"\t"+ref+"\t"+var+"\t"+refNumN+"\t"+altNumN+"\t"+refNumT+"\t"+altNumT
print >> somaticout,out
def SNVexp(germline,somatic,samfile,outputname,case,bedtools):
chro=range(23)
chro.remove(0)
chrom=[]
for ele in chro:
chrom.append("chr"+str(ele))
chrom.append("chrX")
chrom.append("chrY")
SNV={}
chrindex=0
SNPdata=open(germline)
next(SNPdata)
for line in SNPdata:
line=line[0:-1].split("\t")
if "chr" not in line[0]:
chr="chr"+line[0]
else:
chr=line[0]
pos=int(line[1])
ref=line[2]
alt=line[3]
SNV.setdefault(chr,{})[pos]="germline:"+ref+":"+alt
SNPdata.close()
somaticdata=open(somatic)
next(somaticdata)
for line in somaticdata:
line=line[0:-1].split("\t")
if "chr" not in line[0]:
chr="chr"+line[0]
else:
chr=line[0]
pos=int(line[1])
ref=line[2]
alt=line[3]
SNV.setdefault(chr,{})[pos]="somatic:"+ref+":"+alt
somaticdata.close()
sam=open(samfile)
tempout=open(case+".bed","w")
for line in sam:
line=line[0:-1].split("\t")
readid=line[0]
if "@" not in readid:
if "chr" not in line[2]:
chr="chr"+line[2]
else:
chr=line[2]
start=line[3]
seq=line[9]
out=chr+"\t"+start+"\t"+str(int(start)+len(seq)-1)+"\t"+readid+"\t"+seq
print >> tempout,out
sam.close()
tempout.close()
output=open(outputname,'w')
out="chr\tpos\tref\talt\trefNum\taltNum\ttype"
print >> output,out
snvpos=open(case+".snptemp.bed","w")
chrom1=list(set(chrom).intersection(SNV.keys()))
for chromosome in chrom1:
temp=SNV[chromosome]
for key in sorted(temp.keys()):
type=temp[key].split(":")[0]
ref=temp[key].split(":")[1]
alt=temp[key].split(":")[2]
out=chromosome+"\t"+str(key)+"\t"+str(key)+"\t"+ref+"\t"+alt+"\t"+type
print >> snvpos,out
snvpos.close()
os.system(bedtools+" intersect -a "+case+".bed -b "+case+".snptemp.bed -wa -wb > "+case+".snp.read.bed")
snpdata=open(case+".snp.read.bed")
read={}
SNPpos={}
for line in snpdata:
line=line[0:-1].split("\t")
chr=line[0]
start=int(line[1])
end=int(line[2])
readid=line[3]
sequen=line[4]
pos=int(line[6])
ref=line[8]
alt=line[9]
if sequen[pos-start] == ref:
SNPpos.setdefault(chr,[]).append(pos)
read.setdefault(chr+":"+str(pos),[]).append("REF\t"+readid)
elif sequen[pos-start] == alt:
SNPpos.setdefault(chr,[]).append(pos)
read.setdefault(chr+":"+str(pos),[]).append("ALT\t"+readid)
snpdata.close()
for chromosome in chrom:
temppos=sorted(list(set(SNPpos[str(chromosome)])))
temp=SNV[str(chromosome)]
for key in temppos:
type=temp[key].split(":")[0]
ref=temp[key].split(":")[1]
alt=temp[key].split(":")[2]
nref=0
nalt=0
for id in list(set(read[str(chromosome)+":"+str(key)])):
if "REF\t" in id:
nref=nref+1
elif "ALT\t" in id:
nalt=nalt+1
out=str(chromosome)+"\t"+str(key)+"\t"+ref+"\t"+alt+"\t"+str(nref)+"\t"+str(nalt)+"\t"+type
print >> output,out
output.close()
def datatrans(input,output):
data=open(input)
outputdata=open(output,"w")
line=next(data)[0:-1]
print >> outputdata,line
for line in data:
print >> outputdata,"chr"+line[0:-1]
data.close()
outputdata.close()
def SNPcalling(Tumor,Normal,samtools,picard,ref,Varscan):
os.system(samtools+" index "+Tumor)
os.system(samtools+" index "+Normal)
tumorpile="tumor.pileup"
normalpile="normal.pileup"
commond=samtools+" mpileup -f "+ref+" "+Tumor+" > "+tumorpile
os.system(commond)
commond=samtools+" mpileup -f "+ref+" "+Normal+" > "+normalpile
os.system(commond)
commond="java -jar "+Varscan+" somatic "+normalpile+" "+tumorpile+" "+"case.varscan"
os.system(commond)
VarscanTrans(varscan="case.varscan.snp",case="case")
def getPath(path):
path1=path.split("/")
if path1[0] == ".":
if (len(path1)==1):
newpath=os.getcwd()
else:
newpath=os.getcwd()
for ele in path1:
if ele !=".":
newpath=newpath+"/"+ele
elif path1[0]=="..":
i = 0
for ele in path1:
if ele == "..":
i=i+1
path2=os.getcwd()
path2=path2.split("/")
newpath="/"+path2[0]
if len(path2)-i > 1:
for j in range(1,len(path2)-i):
newpath=newpath+"/"+path2[j]
else:
newpath=path
return newpath
def main():
usage = "usage: python %prog [-t <tumor bam file>] [-n <normal bam file>] [-r <RNA bam file>] [-v <Varscan mutation calling output>] [-g <Defiend germline mutation input file>] [-s <Defined somatic mutation input file>] [-o <output path>] [-u <Optimization>] [-e <Defined expression file of mutation>] -p <Texomer path> -I <input form> "
description = "[options] is optional. The path of Texomer and the form of input file are required in Texomer. Three forms of input file are allowed: BAM, Varscan and Defined. If -I BAM is selected, -t and -n are required. Texomer requires bam file is alligned based on GRCH38. If -I Varscan is selected, -v is required. If -I Defined is selected, -g and -s are required. RNA-seq bam file is optional for Texomer. Texomer would do estimation only at DNA level if without RNA data input."
op = OptionParser(version="%prog 0.1",description=description,usage=usage,add_help_option=False)
op.add_option("-h","--help",action="help",
help="Show this help message and exit.")
op.add_option("-p","--Texomer",dest="Texomer",type="str",
help="the path of Texomer")
op.add_option("-I","--Input",dest="Input",type="str",
help="the form of input file: BAM, Varscan, and Defined")
op.add_option("-t","--Tumor",dest="Tumor",type="str",
help="Tumor WES bam file.")
op.add_option("-n","--Normal",dest="Normal",type="str",
help="Normal WES bam file.")
op.add_option("-o","--outpath",dest="outpath",type="str",
help="the output path.")
op.add_option("-r","--RNA",dest="RNA",type="str",
help="RNA-seq bam file")
op.add_option("-g","--germline",dest="germline",type="str",
help="You can input your own germline mutation file together with somatic mutation file (-s). The file include 8 columns: chromosome, position, RefAllele, Altallele, read counts of RefAllele in normal, read counts of Altallele in normal, read counts of RefAllele in tumor and read counts of Altallele in tumor with header seperated by tab. \t")
op.add_option("-s","--somatic",dest="somatic",type="str",
help="You can input your own somatic mutation file together with germline mutation file (-g). The file include 8 columns: chromosome, position, RefAllele, Altallele, read counts of RefAllele in normal, read counts of Altallele in normal, read counts of RefAllele in tumor and read counts of Altallele in tumor with header seperated by tab. \t")
op.add_option("-v","--varscan",dest="varscan",type="str",
help="the output of Varscan2 based on somatic calling")
op.add_option("-e","--snvexpress",dest="snvexpress",type="str",
help="The allelic read count of mutation from RNA-seq including 7 columns: chromosome, position, ref, alt, refnum, altnum and type (germline or somatic)")
op.add_option("-u","--iter",dest="iter",type="int",
help="Optimization using somatic mutation, 0 corresponding to no optimization and 1 is optimization. The default = 1")
(options,args) = op.parse_args()
if not options.Texomer or not options.Input:
op.print_help()
sys.exit(1)
if options.Input not in ["BAM","Varscan","Defined"]:
print "Please select input form from BAM, varscan and Defined!"
op.print_help()
sys.exit(1)
currentPath=os.getcwd()
Texomerpath=options.Texomer
Texomerpath=getPath(Texomerpath)
samtools=Texomerpath+"/samtools"
bedtools=Texomerpath+"/bedtools"
picard=Texomerpath+"/picard.jar"
ref=Texomerpath+"/reference/GRCH38.fa"
Varscan=Texomerpath+"/VarScan.v2.3.9.jar"
if not options.outpath:
outpath=os.getcwd()
else:
outpath=options.outpath
outpath=getPath(outpath)
outpath1=outpath.split('/')
path1="/"
for i in range(1,(len(outpath1)-1)):
path1=path1+outpath1[i]+"/"
if outpath1[len(outpath1)-1] not in os.listdir(path1):
os.chdir(path1)
os.system("mkdir "+outpath1[len(outpath1)-1])
os.chdir(outpath)
os.system("mkdir temp")
os.chdir(currentPath)
#Rscript = options.Rscript
if options.Input == "Varscan":
if not options.varscan:
print "Please input mutation file matched with Varscan output through -v."
sys.exit(1)
else:
varscan=options.varscan
print "Transfer varscan2 output"
varscan=getPath(varscan)
case=varscan.split("/")[len(varscan.split("/"))-1]
os.chdir(outpath+"/temp")
VarscanTrans(varscan,case)
germlinename = os.getcwd()+"/"+case+".germline"
somaticname = os.getcwd()+"/"+case+".somatic"
elif options.Input == "BAM":
if options.Tumor and options.Normal:
Tumor=options.Tumor
Tumor=getPath(Tumor)
Normal=options.Normal
Normal=getPath(Normal)
case="case"
os.chdir(outpath+"/temp")
SNPcalling(Tumor=Tumor,Normal=Normal,samtools=samtools,picard=picard,ref=ref,Varscan=Varscan)
germlinename = os.getcwd()+"/"+case+".germline"
somaticname = os.getcwd()+"/"+case+".somatic"
os.system("rm "+Tumor+".bai")
os.system("rm "+Normal+".bai")
else:
print "Please input bam file through -t and -n."
sys.exit(1)
elif options.Input == "Defined":
if options.germline and options.somatic:
germlinename = options.germline
germlinename=getPath(germlinename)
somaticname = options.somatic
somaticname=getPath(somaticname)
else:
print "Please input defined germline and somatic mutation files through -g and -s."
sys.exit(1)
os.chdir(outpath+"/temp")
data=open(somaticname)
next(data)
line=next(data)
data.close()
if "chr" not in line:
case=somaticname.split("/")[len(somaticname.split("/"))-1]
somaticoutname = os.getcwd()+"/"+case+".somatic"
datatrans(input=somaticname,output=somaticoutname)
somaticname=somaticoutname
data=open(germlinename)
next(data)
line=next(data)
data.close()
if "chr" not in line:
case=germlinename.split("/")[len(germlinename.split("/"))-1]
germlineoutname = os.getcwd()+"/"+case+".germline"
datatrans(input=germlinename,output=germlineoutname)
germlinename=germlineoutname
os.chdir(currentPath)
if not options.RNA:
if not options.snvexpress:
print "No RNA data! Only estimate at DNA level!"
RNAname="NA"
else:
RNAname=options.snvexpress
RNAname=getPath(RNAname)
data=open(RNAname)
next(data)
line=next(data)[0:-1].split("\t")
data.close()
if "chr" not in line[0]:
case=outputname.split("/")[len(outputname.split("/"))-1]
os.chdir(outpath+"/temp")
datatrans(input=outputname,output=case)
#outputname=case
else:
RNAname = options.RNA
RNAname=getPath(RNAname)
case=RNAname.split("/")[len(RNAname.split("/"))-1]
print "Transfer bam to sam"
os.chdir(outpath+"/temp")
os.system(samtools+ " index "+RNAname)
os.system(samtools+" view "+RNAname+" > "+RNAname+".sam")
samfilename=RNAname+".sam"
outputname=case+".SNV"
if os.path.exists(outputname):
os.remove(outputname)
print "Align mutations to RNA data"
SNVexp(germline=germlinename,somatic=somaticname,samfile=samfilename,outputname=outputname,case=case,bedtools=bedtools)
os.system("rm "+RNAname+".sam")
os.system("rm "+RNAname+".bai")
RNAname=outputname
print "Infering at DNA level"
if not options.iter:
optindex=1
else:
optindex=options.iter
os.system("Rscript "+Texomerpath+"/Texomer.R "+germlinename+" "+somaticname+" "+RNAname+" "+Texomerpath+" "+outpath+" "+str(optindex))
os.chdir(outpath)
os.system("rm -r temp")
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt me, see you!\n")
sys.exit(0)