[2] | 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 | |
---|