/tools/solid_tools/maq_cs_wrapper.py

  1#!/usr/bin/env python
  2#Guruprasad Ananda
  3#MAQ mapper for SOLiD colourspace-reads
  4
  5import sys, os, zipfile, tempfile, subprocess
  6
  7def stop_err( msg ):
  8    sys.stderr.write( "%s\n" % msg )
  9    sys.exit()
 10 
 11def __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
267if __name__=="__main__":
268    __main__()
269        
270