Mercurial > repos > xuebing > sharplabtool

comparison tools/solid_tools/maq_cs_wrapper.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 #!/usr/bin/env python
2 #Guruprasad Ananda
3 #MAQ mapper for SOLiD colourspace-reads
4
5 import sys, os, zipfile, tempfile, subprocess
6
7 def stop_err( msg ):
8 sys.stderr.write( "%s\n" % msg )
9 sys.exit()
10
11 def __main__():
12
13 out_fname = sys.argv[1].strip()
14 out_f2 = open(sys.argv[2].strip(),'r+')
15 ref_fname = sys.argv[3].strip()
16 f3_read_fname = sys.argv[4].strip()
17 f3_qual_fname = sys.argv[5].strip()
18 paired = sys.argv[6]
19 if paired == 'yes':
20 r3_read_fname = sys.argv[7].strip()
21 r3_qual_fname = sys.argv[8].strip()
22 min_mapqual = int(sys.argv[9].strip())
23 max_mismatch = int(sys.argv[10].strip())
24 out_f3name = sys.argv[11].strip()
25 subprocess_dict = {}
26
27 ref_csfa = tempfile.NamedTemporaryFile()
28 ref_bfa = tempfile.NamedTemporaryFile()
29 ref_csbfa = tempfile.NamedTemporaryFile()
30 cmd2_1 = 'maq fasta2csfa %s > %s 2>&1' %(ref_fname,ref_csfa.name)
31 cmd2_2 = 'maq fasta2bfa %s %s 2>&1' %(ref_csfa.name,ref_csbfa.name)
32 cmd2_3 = 'maq fasta2bfa %s %s 2>&1' %(ref_fname,ref_bfa.name)
33 try:
34 os.system(cmd2_1)
35 os.system(cmd2_2)
36 os.system(cmd2_3)
37 except Exception, erf:
38 stop_err(str(erf)+"Error processing reference sequence")
39
40 if paired == 'yes': #paired end reads
41 tmpf = tempfile.NamedTemporaryFile() #forward reads
42 tmpr = tempfile.NamedTemporaryFile() #reverse reads
43 tmps = tempfile.NamedTemporaryFile() #single reads
44 tmpffastq = tempfile.NamedTemporaryFile()
45 tmprfastq = tempfile.NamedTemporaryFile()
46 tmpsfastq = tempfile.NamedTemporaryFile()
47
48 cmd1 = "solid2fastq_modified.pl 'yes' %s %s %s %s %s %s %s 2>&1" %(tmpf.name,tmpr.name,tmps.name,f3_read_fname,f3_qual_fname,r3_read_fname,r3_qual_fname)
49 try:
50 os.system(cmd1)
51 os.system('gunzip -c %s >> %s' %(tmpf.name,tmpffastq.name))
52 os.system('gunzip -c %s >> %s' %(tmpr.name,tmprfastq.name))
53 os.system('gunzip -c %s >> %s' %(tmps.name,tmpsfastq.name))
54
55 except Exception, eq:
56 stop_err("Error converting data to fastq format." + str(eq))
57
58 #make a temp directory where the split fastq files will be stored
59 try:
60 split_dir = tempfile.mkdtemp()
61 split_file_prefix_f = tempfile.mktemp(dir=split_dir)
62 split_file_prefix_r = tempfile.mktemp(dir=split_dir)
63 splitcmd_f = 'split -a 2 -l %d %s %s' %(32000000,tmpffastq.name,split_file_prefix_f) #32M lines correspond to 8M reads
64 splitcmd_r = 'split -a 2 -l %d %s %s' %(32000000,tmprfastq.name,split_file_prefix_r) #32M lines correspond to 8M reads
65
66 os.system(splitcmd_f)
67 os.system(splitcmd_r)
68 os.chdir(split_dir)
69 ii = 0
70 for fastq in os.listdir(split_dir):
71 if not fastq.startswith(split_file_prefix_f.split("/")[-1]):
72 continue
73 fastq_r = split_file_prefix_r + fastq.split(split_file_prefix_f.split("/")[-1])[1] #find the reverse strand fastq corresponding to formward strand fastq
74 tmpbfq_f = tempfile.NamedTemporaryFile()
75 tmpbfq_r = tempfile.NamedTemporaryFile()
76 cmd3 = 'maq fastq2bfq %s %s 2>&1; maq fastq2bfq %s %s 2>&1; maq map -c %s.csmap %s %s %s 1>/dev/null 2>&1; maq mapview %s.csmap > %s.txt' %(fastq,tmpbfq_f.name,fastq_r,tmpbfq_r.name,fastq,ref_csbfa.name,tmpbfq_f.name,tmpbfq_r.name,fastq,fastq)
77 subprocess_dict['sp'+str(ii+1)] = subprocess.Popen([cmd3],shell=True,stdout=subprocess.PIPE)
78 ii += 1
79 while True:
80 all_done = True
81 for j,k in enumerate(subprocess_dict.keys()):
82 if subprocess_dict['sp'+str(j+1)].wait() != 0:
83 err = subprocess_dict['sp'+str(j+1)].communicate()[1]
84 if err != None:
85 stop_err("Mapping error: %s" %err)
86 all_done = False
87 if all_done:
88 break
89 cmdout = "for map in *.txt; do cat $map >> %s; done" %(out_fname)
90 os.system(cmdout)
91
92 tmpcsmap = tempfile.NamedTemporaryFile()
93 cmd_cat_csmap = "for csmap in *.csmap; do cat $csmap >> %s; done" %(tmpcsmap.name)
94 os.system(cmd_cat_csmap)
95
96 tmppileup = tempfile.NamedTemporaryFile()
97 cmdpileup = "maq pileup -m %s -q %s %s %s > %s" %(max_mismatch,min_mapqual,ref_bfa.name,tmpcsmap.name,tmppileup.name)
98 os.system(cmdpileup)
99 tmppileup.seek(0)
100 print >> out_f2, "#chr\tposition\tref_nt\tcoverage\tSNP_count\tA_count\tT_count\tG_count\tC_count"
101 for line in file(tmppileup.name):
102 elems = line.strip().split()
103 ref_nt = elems[2].capitalize()
104 read_nt = elems[4]
105 coverage = int(elems[3])
106 a,t,g,c = 0,0,0,0
107 ref_nt_count = 0
108 for ch in read_nt:
109 ch = ch.capitalize()
110 if ch not in ['A','T','G','C',',','.']:
111 continue
112 if ch in [',','.']:
113 ch = ref_nt
114 ref_nt_count += 1
115 try:
116 nt_ind = ['A','T','G','C'].index(ch)
117 if nt_ind == 0:
118 a+=1
119 elif nt_ind == 1:
120 t+=1
121 elif nt_ind == 2:
122 g+=1
123 else:
124 c+=1
125 except ValueError, we:
126 print >>sys.stderr, we
127 print >> out_f2, "%s\t%s\t%s\t%s\t%s\t%s" %("\t".join(elems[:4]),coverage-ref_nt_count,a,t,g,c)
128 except Exception, er2:
129 stop_err("Encountered error while mapping: %s" %(str(er2)))
130
131
132 else: #single end reads
133 tmpf = tempfile.NamedTemporaryFile()
134 tmpfastq = tempfile.NamedTemporaryFile()
135 cmd1 = "solid2fastq_modified.pl 'no' %s %s %s %s %s %s %s 2>&1" %(tmpf.name,None,None,f3_read_fname,f3_qual_fname,None,None)
136 try:
137 os.system(cmd1)
138 os.system('gunzip -c %s >> %s' %(tmpf.name,tmpfastq.name))
139 tmpf.close()
140 except:
141 stop_err("Error converting data to fastq format.")
142
143 #make a temp directory where the split fastq files will be stored
144 try:
145 split_dir = tempfile.mkdtemp()
146 split_file_prefix = tempfile.mktemp(dir=split_dir)
147 splitcmd = 'split -a 2 -l %d %s %s' %(32000000,tmpfastq.name,split_file_prefix) #32M lines correspond to 8M reads
148 os.system(splitcmd)
149 os.chdir(split_dir)
150 for i,fastq in enumerate(os.listdir(split_dir)):
151 tmpbfq = tempfile.NamedTemporaryFile()
152 cmd3 = 'maq fastq2bfq %s %s 2>&1; maq map -c %s.csmap %s %s 1>/dev/null 2>&1; maq mapview %s.csmap > %s.txt' %(fastq,tmpbfq.name,fastq,ref_csbfa.name,tmpbfq.name,fastq,fastq)
153 subprocess_dict['sp'+str(i+1)] = subprocess.Popen([cmd3],shell=True,stdout=subprocess.PIPE)
154
155 while True:
156 all_done = True
157 for j,k in enumerate(subprocess_dict.keys()):
158 if subprocess_dict['sp'+str(j+1)].wait() != 0:
159 err = subprocess_dict['sp'+str(j+1)].communicate()[1]
160 if err != None:
161 stop_err("Mapping error: %s" %err)
162 all_done = False
163 if all_done:
164 break
165
166 cmdout = "for map in *.txt; do cat $map >> %s; done" %(out_fname)
167 os.system(cmdout)
168
169 tmpcsmap = tempfile.NamedTemporaryFile()
170 cmd_cat_csmap = "for csmap in *.csmap; do cat $csmap >> %s; done" %(tmpcsmap.name)
171 os.system(cmd_cat_csmap)
172
173 tmppileup = tempfile.NamedTemporaryFile()
174 cmdpileup = "maq pileup -m %s -q %s %s %s > %s" %(max_mismatch,min_mapqual,ref_bfa.name,tmpcsmap.name,tmppileup.name)
175 os.system(cmdpileup)
176 tmppileup.seek(0)
177 print >> out_f2, "#chr\tposition\tref_nt\tcoverage\tSNP_count\tA_count\tT_count\tG_count\tC_count"
178 for line in file(tmppileup.name):
179 elems = line.strip().split()
180 ref_nt = elems[2].capitalize()
181 read_nt = elems[4]
182 coverage = int(elems[3])
183 a,t,g,c = 0,0,0,0
184 ref_nt_count = 0
185 for ch in read_nt:
186 ch = ch.capitalize()
187 if ch not in ['A','T','G','C',',','.']:
188 continue
189 if ch in [',','.']:
190 ch = ref_nt
191 ref_nt_count += 1
192 try:
193 nt_ind = ['A','T','G','C'].index(ch)
194 if nt_ind == 0:
195 a+=1
196 elif nt_ind == 1:
197 t+=1
198 elif nt_ind == 2:
199 g+=1
200 else:
201 c+=1
202 except:
203 pass
204 print >> out_f2, "%s\t%s\t%s\t%s\t%s\t%s" %("\t".join(elems[:4]),coverage-ref_nt_count,a,t,g,c)
205 except Exception, er2:
206 stop_err("Encountered error while mapping: %s" %(str(er2)))
207
208 #Build custom track from pileup
209 chr_list=[]
210 out_f2.seek(0)
211 fcov = tempfile.NamedTemporaryFile()
212 fout_a = tempfile.NamedTemporaryFile()
213 fout_t = tempfile.NamedTemporaryFile()
214 fout_g = tempfile.NamedTemporaryFile()
215 fout_c = tempfile.NamedTemporaryFile()
216 fcov.write('''track type=wiggle_0 name="Coverage track" description="Coverage track (from Galaxy)" color=0,0,0 visibility=2\n''')
217 fout_a.write('''track type=wiggle_0 name="Track A" description="Track A (from Galaxy)" color=255,0,0 visibility=2\n''')
218 fout_t.write('''track type=wiggle_0 name="Track T" description="Track T (from Galaxy)" color=0,255,0 visibility=2\n''')
219 fout_g.write('''track type=wiggle_0 name="Track G" description="Track G (from Galaxy)" color=0,0,255 visibility=2\n''')
220 fout_c.write('''track type=wiggle_0 name="Track C" description="Track C (from Galaxy)" color=255,0,255 visibility=2\n''')
221
222 for line in out_f2:
223 if line.startswith("#"):
224 continue
225 elems = line.split()
226 chr = elems[0]
227
228 if chr not in chr_list:
229 chr_list.append(chr)
230 if not (chr.startswith('chr') or chr.startswith('scaffold')):
231 chr = 'chr'
232 header = "variableStep chrom=%s" %(chr)
233 fcov.write("%s\n" %(header))
234 fout_a.write("%s\n" %(header))
235 fout_t.write("%s\n" %(header))
236 fout_g.write("%s\n" %(header))
237 fout_c.write("%s\n" %(header))
238 try:
239 pos = int(elems[1])
240 cov = int(elems[3])
241 a = int(elems[5])
242 t = int(elems[6])
243 g = int(elems[7])
244 c = int(elems[8])
245 except:
246 continue
247 fcov.write("%s\t%s\n" %(pos,cov))
248 try:
249 a_freq = a*100./cov
250 t_freq = t*100./cov
251 g_freq = g*100./cov
252 c_freq = c*100./cov
253 except ZeroDivisionError:
254 a_freq=t_freq=g_freq=c_freq=0
255 fout_a.write("%s\t%s\n" %(pos,a_freq))
256 fout_t.write("%s\t%s\n" %(pos,t_freq))
257 fout_g.write("%s\t%s\n" %(pos,g_freq))
258 fout_c.write("%s\t%s\n" %(pos,c_freq))
259
260 fcov.seek(0)
261 fout_a.seek(0)
262 fout_g.seek(0)
263 fout_t.seek(0)
264 fout_c.seek(0)
265 os.system("cat %s %s %s %s %s | cat > %s" %(fcov.name,fout_a.name,fout_t.name,fout_g.name,fout_c.name,out_f3name))
266
267 if __name__=="__main__":
268 __main__()
269
270