/tools/solid_tools/maq_cs_wrapper.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 270 lines · 238 code · 26 blank · 6 comment · 53 complexity · 9593f2f461ff43d0e32e268dbbb88dc0 MD5 · raw file

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