Mercurial > repos > jpayne > seqsero_v2

comparison libs/mapping_and_assembly_hybrid.py @ 0:4ff2aee11e5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author jpayne
date Tue, 06 Nov 2018 09:45:57 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:4ff2aee11e5b
1 import os,sys,glob,time,itertools,subprocess
2 from collections import OrderedDict
3 from Initial_Conditions import phase1
4 from Initial_Conditions import phase2
5 from Initial_Conditions import phaseO
6 from Initial_Conditions import sero
7 from distutils.version import LooseVersion
8
9 import csv
10 from subprocess import check_output
11
12
13 def xml_parse_score_comparision_seqsero(xmlfile):
14 #used to do seqsero xml analysis
15 from Bio.Blast import NCBIXML
16 handle=open(xmlfile)
17 handle=NCBIXML.parse(handle)
18 handle=list(handle)
19 List=[]
20 List_score=[]
21 List_ids=[]
22 for i in range(len(handle)):
23 if len(handle[i].alignments)>0:
24 for j in range(len(handle[i].alignments)):
25 score=0
26 ids=0
27 List.append(handle[i].query.strip()+"___"+handle[i].alignments[j].hit_def)
28 for z in range(len(handle[i].alignments[j].hsps)):
29 if "last" in handle[i].query or "first" in handle[i].query:
30 score+=handle[i].alignments[j].hsps[z].bits
31 ids+=float(handle[i].alignments[j].hsps[z].identities)/handle[i].query_length
32 else:
33 if handle[i].alignments[j].hsps[z].align_length>=30:
34 #for the long alleles, filter noise parts
35 score+=handle[i].alignments[j].hsps[z].bits
36 ids+=float(handle[i].alignments[j].hsps[z].identities)/handle[i].query_length
37 List_score.append(score)
38 List_ids.append(ids)
39 temp=zip(List,List_score,List_ids)
40 Final_list=sorted(temp, key=lambda d:d[1], reverse = True)
41 return Final_list
42
43
44 def Uniq(L,sort_on_fre="none"): #return the uniq list and the count number
45 Old=L
46 L.sort()
47 L = [L[i] for i in range(len(L)) if L[i] not in L[:i]]
48 count=[]
49 for j in range(len(L)):
50 y=0
51 for x in Old:
52 if L[j]==x:
53 y+=1
54 count.append(y)
55 if sort_on_fre!="none":
56 d=zip(*sorted(zip(count, L)))
57 L=d[1]
58 count=d[0]
59 return (L,count)
60
61
62 def judge_fliC_or_fljB_from_head_tail_for_one_contig(nodes_vs_score_list):
63 #used to predict it's fliC or fljB for one contig, based on tail and head score, but output the score difference,if it is very small, then not reliable, use blast score for whole contig to test
64 #this is mainly used for
65 a=nodes_vs_score_list
66 fliC_score=0
67 fljB_score=0
68 for z in a:
69 if "fliC" in z[0]:
70 fliC_score+=z[1]
71 elif "fljB" in z[0]:
72 fljB_score+=z[1]
73 if fliC_score>=fljB_score:
74 role="fliC"
75 else:
76 role="fljB"
77 return (role,abs(fliC_score-fljB_score))
78
79 def judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(node_name,Final_list_passed):
80 #used to predict contig is fliC or fljB, if the differnce score value on above head_and_tail is less than 10 (quite small)
81 #also used when no head or tail got blasted score for the contig
82 role=""
83 for z in Final_list_passed:
84 if node_name in z[0]:
85 role=z[0].split("_")[0]
86 break
87 return role
88
89
90 def fliC_or_fljB_judge_from_head_tail_sequence(nodes_list,tail_head_list,Final_list_passed):
91 #nodes_list is the c created by c,d=Uniq(nodes) in below function
92 first_target=""
93 role_list=[]
94 for x in nodes_list:
95 a=[]
96 role=""
97 for y in tail_head_list:
98 if x in y[0]:
99 a.append(y)
100 if len(a)==4:
101 #compare two heads (37 > 30)
102 #four contigs, most perfect assembly, high quality
103 """
104 for z in a:
105 if "fliC_first_37" in z[0]:
106 t1=z[1]
107 elif "fljB_first_37" in z[0]:
108 t2=z[1]
109 if t1>=t2:
110 role="fliC"
111 else:
112 role="fljB"
113 """
114 role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
115 if diff<20:
116 role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
117 elif len(a)==3:
118 """
119 #compare the number, because hybrid problem
120 temp=[]
121 for z in a:
122 temp.append(z[0].split("_")[0])
123 m,n=Uniq(temp)#only two choices in m or n
124 if n[0]>n[1]:
125 role=m[0]
126 else:
127 role=m[1]
128 """
129 ###however, if the one with highest score is the fewer one, compare their accumulation score
130 role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
131 if diff<20:
132 role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
133 ###end of above score comparison
134 elif len(a)==2:
135 #must on same node, if not, then decide with unit blast score, blast-score/length_of_special_sequence(30 or 37)
136 temp=[]
137 for z in a:
138 temp.append(z[0].split("_")[0])
139 m,n=Uniq(temp)#should only have one choice, but weird situation might occur too
140 if len(m)==1:
141 pass
142 else:
143 pass
144 #print "head and tail not belong to same role, now let's guess based on maximum likelihood"
145 role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
146 if diff<20:
147 role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
148 """
149 max_unit_score=0
150 for z in a:
151 unit_score=z[-1]/int(z[0].split("__")[1])
152 if unit_score>=max_unit_score:
153 role=z[0].split("_")[0]
154 max_unit_score=unit_score
155 """
156 ###need to desgin a algorithm to guess most possible situation for nodes_list, See the situations of test evaluation
157 elif len(a)==1:
158 #that one
159 role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
160 if diff<20:
161 role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
162 #role=a[0][0].split("_")[0]
163 #need to evaluate, in future, may set up a cut-off, if not met, then just find Final_list_passed best match,like when "a==0"
164 else:#a==0
165 #use Final_list_passed best match
166 for z in Final_list_passed:
167 if x in z[0]:
168 role=z[0].split("_")[0]
169 break
170 #print x,role,len(a)
171 role_list.append((role,x))
172 if len(role_list)==2:
173 if role_list[0][0]==role_list[1][0]:#this is the most cocmmon error, two antigen were assigned to same phase
174 #just use score to do a final test
175 role_list=[]
176 for x in nodes_list:
177 role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
178 role_list.append((role,x))
179 return role_list
180
181 def decide_contig_roles_for_H_antigen(Final_list):
182 #used to decide which contig is FliC and which one is fljB
183 contigs=[]
184 Final_list_passed=[x for x in Final_list if float(x[0].split("_cov_")[1])>=3.5 and (x[1]>=int(x[0].split("__")[1]) or x[1]>=int(x[0].split("___")[1].split("_")[3]))]
185 nodes=[]
186 for x in Final_list_passed:
187 if x[0].startswith("fl") and "last" not in x[0] and "first" not in x[0]:
188 nodes.append(x[0].split("___")[1].strip())
189 c,d=Uniq(nodes)#c is node_list
190 #print c
191 tail_head_list=[x for x in Final_list if ("last" in x[0] or "first" in x[0])]
192 roles=fliC_or_fljB_judge_from_head_tail_sequence(c,tail_head_list,Final_list_passed)
193 return roles
194
195 def Combine(b,c):
196 fliC_combinations=[]
197 fliC_combinations.append(",".join(c))
198 temp_combinations=[]
199 for i in range(len(b)):
200 for x in itertools.combinations(b,i+1):
201 temp_combinations.append(",".join(x))
202 for x in temp_combinations:
203 temp=[]
204 for y in c:
205 temp.append(y)
206 temp.append(x)
207 temp=",".join(temp)
208 temp=temp.split(",")
209 temp.sort()
210 temp=",".join(temp)
211 fliC_combinations.append(temp)
212 return fliC_combinations
213
214 def decide_O_type_and_get_special_genes(Final_list):
215 #decide O based on Final_list
216 O_choice="?"
217 O_list=[]
218 special_genes=[]
219 Final_list_passed=[x for x in Final_list if float(x[0].split("_cov_")[1])>=3.5 and (x[1]>=int(x[0].split("__")[1]) or x[1]>=int(x[0].split("___")[1].split("_")[3]))]
220 nodes=[]
221 for x in Final_list_passed:
222 if x[0].startswith("O-"):
223 nodes.append(x[0].split("___")[1].strip())
224 elif not x[0].startswith("fl"):
225 special_genes.append(x)
226 #print "special_genes:",special_genes
227 c,d=Uniq(nodes)
228 #print "potential O antigen contig",c
229 final_O=[]
230 O_nodes_list=[]
231 for x in c:#c is the list for contigs
232 temp=0
233 for y in Final_list_passed:
234 if x in y[0] and y[0].startswith("O-"):
235 final_O.append(y)
236 break
237 ### O contig has the problem of two genes on same contig, so do additional test
238 potenial_new_gene=""
239 for x in final_O:
240 pointer=0 #for genes merged or not
241 #not consider O-1,3,19_not_in_3,10, too short compared with others
242 if "O-1,3,19_not_in_3,10" not in x[0] and int(x[0].split("__")[1].split("___")[0])+800 <= int(x[0].split("length_")[1].split("_")[0]):#gene length << contig length; for now give 300*2 (for secureity can use 400*2) as flank region
243 pointer=x[0].split("___")[1].strip()#store the contig name
244 print pointer
245 if pointer!=0:#it has potential merge event
246 for y in Final_list:
247 if pointer in y[0] and y not in final_O and (y[1]>=int(y[0].split("__")[1].split("___")[0])*1.5 or (y[1]>=int(y[0].split("__")[1].split("___")[0])*y[2] and y[1]>=400)):#that's a realtively strict filter now; if passed, it has merge event and add one more to final_O
248 potenial_new_gene=y
249 print potenial_new_gene
250 break
251 if potenial_new_gene!="":
252 print "two differnt genes in same contig, fix it for O antigen"
253 final_O.append(potenial_new_gene)
254 ### end of the two genes on same contig test
255 if len(final_O)==0:
256 #print "$$$No Otype, due to no hit"#may need to be changed
257 O_choice="-"
258 else:
259 O_list=[]
260 for x in final_O:
261 O_list.append(x[0].split("__")[0])
262 if not "O-1,3,19_not_in_3,10__130" in x[0]:#O-1,3,19_not_in_3,10 is too small, which may affect further analysis
263 O_nodes_list.append(x[0].split("___")[1])
264 ### special test for O9,46 and O3,10 family
265 if "O-9,46_wbaV" in O_list:#not sure should use and float(O9_wbaV)/float(num_1) > 0.1
266 if "O-9,46_wzy" in O_list:#and float(O946_wzy)/float(num_1) > 0.1
267 O_choice="O-9,46"
268 #print "$$$Most possilble Otype: O-9,46"
269 elif "O-9,46,27_partial_wzy" in O_list:#and float(O94627)/float(num_1) > 0.1
270 O_choice="O-9,46,27"
271 #print "$$$Most possilble Otype: O-9,46,27"
272 else:
273 O_choice="O-9"#next, detect O9 vs O2?
274 O2=0
275 O9=0
276 for z in special_genes:
277 if "tyr-O-9" in z[0]:
278 O9=z[1]
279 elif "tyr-O-2" in z[0]:
280 O2=z[1]
281 if O2>O9:
282 O_choice="O-2"
283 elif O2<O9:
284 pass
285 else:
286 pass
287 #print "$$$No suitable one, because can't distinct it's O-9 or O-2, but O-9 has a more possibility."
288 elif ("O-3,10_wzx" in O_list) and ("O-9,46_wzy" in O_list):#and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1
289 if "O-3,10_not_in_1,3,19" in O_list:#and float(O310_no_1319)/float(num_1) > 0.1
290 O_choice="O-3,10"
291 #print "$$$Most possilble Otype: O-3,10 (contain O-3,10_not_in_1,3,19)"
292 else:
293 O_choice="O-1,3,19"
294 #print "$$$Most possilble Otype: O-1,3,19 (not contain O-3,10_not_in_1,3,19)"
295 ### end of special test for O9,46 and O3,10 family
296 else:
297 try:
298 max_score=0
299 for x in final_O:
300 if x[1]>=max_score:
301 max_score=x[1]
302 O_choice=x[0].split("_")[0]
303 if O_choice=="O-1,3,19":
304 O_choice=final_O[1][0].split("_")[0]
305 #print "$$$Most possilble Otype: ",O_choice
306 except:
307 pass
308 #print "$$$No suitable Otype, or failure of mapping (please check the quality of raw reads)"
309 #print "O:",O_choice,O_nodes_list
310 return O_choice,O_nodes_list,special_genes,final_O
311
312 def seqsero_from_formula_to_serotypes(Otype,fliC,fljB,special_gene_list):
313 #like test_output_06012017.txt
314 #can add more varialbles like sdf-type, sub-species-type in future (we can conclude it into a special-gene-list)
315 from Initial_Conditions import phase1
316 from Initial_Conditions import phase2
317 from Initial_Conditions import phaseO
318 from Initial_Conditions import sero
319 seronames=[]
320 for i in range(len(phase1)):
321 fliC_combine=[]
322 fljB_combine=[]
323 if phaseO[i]==Otype:
324 ### for fliC, detect every possible combinations to avoid the effect of "["
325 if phase1[i].count("[")==0:
326 fliC_combine.append(phase1[i])
327 elif phase1[i].count("[")>=1:
328 c=[]
329 b=[]
330 if phase1[i][0]=="[" and phase1[i][-1]=="]" and phase1[i].count("[")==1:
331 content=phase1[i].replace("[","").replace("]","")
332 fliC_combine.append(content)
333 fliC_combine.append("-")
334 else:
335 for x in phase1[i].split(","):
336 if "[" in x:
337 b.append(x.replace("[","").replace("]",""))
338 else:
339 c.append(x)
340 fliC_combine=Combine(b,c) #Combine will offer every possible combinations of the formula, like f,[g],t: f,t f,g,t
341 ### end of fliC "[" detect
342 ### for fljB, detect every possible combinations to avoid the effect of "["
343 if phase2[i].count("[")==0:
344 fljB_combine.append(phase2[i])
345 elif phase2[i].count("[")>=1:
346 d=[]
347 e=[]
348 if phase2[i][0]=="[" and phase2[i][-1]=="]" and phase2[i].count("[")==1:
349 content=phase2[i].replace("[","").replace("]","")
350 fljB_combine.append(content)
351 fljB_combine.append("-")
352 else:
353 for x in phase2[i].split(","):
354 if "[" in x:
355 d.append(x.replace("[","").replace("]",""))
356 else:
357 e.append(x)
358 fljB_combine=Combine(d,e)
359 ### end of fljB "[" detect
360 new_fliC=fliC.split(",") #because some antigen like r,[i] not follow alphabetical order, so use this one to judge and can avoid missings
361 new_fliC.sort()
362 new_fliC=",".join(new_fliC)
363 new_fljB=fljB.split(",")
364 new_fljB.sort()
365 new_fljB=",".join(new_fljB)
366 if (new_fliC in fliC_combine or fliC in fliC_combine) and (new_fljB in fljB_combine or fljB in fljB_combine):
367 seronames.append(sero[i])
368 #analyze seronames
369 if len(seronames)==0:
370 seronames=["N/A (The predicted antigenic profile does not exist in the White-Kauffmann-Le Minor scheme)"]
371 star=""
372 star_line=""
373 if len(seronames)>1:#there are two possible predictions for serotypes
374 star="*"
375 star_line="The predicted serotypes share the same general formula:\t"+Otype+":"+fliC+":"+fljB+"\n"##
376 print "\n"
377 predict_form=Otype+":"+fliC+":"+fljB#
378 predict_sero=(" or ").join(seronames)
379 ###special test for Enteritidis
380 if predict_form=="9:g,m:-":
381 sdf="-"
382 for x in special_gene_list:
383 if x[0].startswith("sdf"):
384 sdf="+"
385 predict_form=predict_form+"\nSdf prediction:"+sdf
386 if sdf=="-":
387 star="*"
388 star_line="Additional characterization is necessary to assign a serotype to this strain. Commonly circulating strains of serotype Enteritidis are sdf+, although sdf- strains of serotype Enteritidis are known to exist. Serotype Gallinarum is typically sdf- but should be quite rare. Sdf- strains of serotype Enteritidis and serotype Gallinarum can be differentiated by phenotypic profile or genetic criteria.\n"#+##
389 predict_sero="See comments below"
390 ###end of special test for Enteritidis
391 elif predict_form=="4:i:-":
392 predict_sero="potential monophasic variant of Typhimurium"
393 elif predict_form=="4:r:-":
394 predict_sero="potential monophasic variant of Heidelberg"
395 elif predict_form=="4:b:-":
396 predict_sero="potential monophasic variant of Paratyphi B"
397 elif predict_form=="8:e,h:1,2":
398 predict_sero="Newport"
399 star="*"
400 star_line="Serotype Bardo shares the same antigenic profile with Newport, but Bardo is exceedingly rare."
401 claim="The serotype(s) is/are the only serotype(s) with the indicated antigenic profile currently recognized in the Kauffmann White Scheme. New serotypes can emerge and the possibility exists that this antigenic profile may emerge in a different subspecies. Identification of strains to the subspecies level should accompany serotype determination; the same antigenic profile in different subspecies is considered different serotypes."##
402 if "N/A" in predict_sero:
403 claim=""
404 if "Typhimurium" in predict_sero or predict_form=="4:i:-":
405 normal=0
406 mutation=0
407 for x in special_gene_list:
408 if "oafA-O-4_full" in x[0]:
409 normal=x[1]
410 elif "oafA-O-4_5-" in x[0]:
411 mutation=x[1]
412 if normal>mutation:
413 #print "$$$Typhimurium"
414 pass
415 elif normal<mutation:
416 predict_sero=predict_sero.strip()+"(O5-)"
417 star="*"#
418 star_line="Detected the deletion of O5-."
419 #print "$$$Typhimurium_O5-"
420 else:
421 #print "$$$Typhimurium, even no 7 bases difference"
422 pass
423 return predict_form,predict_sero,star,star_line,claim
424
425 def main():
426 database=sys.argv[1]#used to extract reads
427 mapping_mode=sys.argv[2]#mem or sampe
428 threads=sys.argv[3]
429 for_fq=sys.argv[4]
430 rev_fq=sys.argv[5]
431 current_time=time.strftime("%Y_%m_%d_%H_%M_%S", time.localtime())
432 sam=for_fq+".sam"
433 bam=for_fq+".bam"
434 sorted_bam=for_fq+"_sorted.bam"
435 mapped_fq1=for_fq+"_mapped.fq"
436 mapped_fq2=rev_fq+"_mapped.fq"
437 combined_fq=for_fq+"_combined.fq"
438 for_sai=for_fq+".sai"
439 rev_sai=rev_fq+".sai"
440 print "building database..."
441 #os.system("bwa index "+database+ " 2> /dev/null")
442 os.system("bwa index "+database+ " 2>> data_log.txt ")
443 print "mapping..."
444 if mapping_mode=="mem":
445 os.system("bwa mem -t "+threads+" "+database+" "+for_fq+" "+rev_fq+" > "+sam+ " 2>> data_log.txt")
446 elif mapping_mode=="sam":
447 os.system("bwa aln -t "+threads+" "+database+" "+for_fq+" > "+for_sai+ " 2>> data_log.txt")
448 os.system("bwa aln -t "+threads+" "+database+" "+rev_fq+" > "+rev_sai+ " 2>> data_log.txt")
449 os.system("bwa sampe "+database+" "+for_sai+" "+ rev_sai+" "+for_fq+" "+rev_fq+" > "+sam+ " 2>> data_log.txt")
450 os.system("samtools view -@ "+threads+" -F 4 -Sbh "+sam+" > "+bam)
451 os.system("samtools view -@ "+threads+" -h -o "+sam+" "+bam)
452 ### check the version of samtools then use differnt commands
453 samtools_version=subprocess.Popen(["samtools"],stdout=subprocess.PIPE,stderr=subprocess.PIPE)
454 out, err = samtools_version.communicate()
455 version = err.split("ersion:")[1].strip().split(" ")[0].strip()
456 print "check samtools version:",version
457 if LooseVersion(version)<=LooseVersion("1.2"):
458 os.system("samtools sort -@ "+threads+" -n "+bam+" "+for_fq+"_sorted")
459 else:
460 os.system("samtools sort -@ "+threads+" -n "+bam+" >"+sorted_bam)
461 ### end of samtools version check and its analysis
462 os.system("bamToFastq -i "+sorted_bam+" -fq "+combined_fq)
463 os.system("bamToFastq -i "+sorted_bam+" -fq "+mapped_fq1+" -fq2 "+mapped_fq2 + " 2>> data_log.txt")#2> /dev/null if want no output
464 outdir=current_time+"_temp"
465 print "assembling..."
466 if int(threads)>4:
467 t="4"
468 else:
469 t=threads
470 os.system("spades.py --careful --pe1-s "+combined_fq+" --pe1-1 "+mapped_fq1+" --pe1-2 "+mapped_fq2+" -t "+t+" -o "+outdir+ " >> data_log.txt 2>&1")
471 new_fasta=for_fq+"_"+database+"_"+mapping_mode+".fasta"
472 os.system("mv "+outdir+"/contigs.fasta "+new_fasta+ " 2> /dev/null")
473 #os.system("mv "+outdir+"/scaffolds.fasta "+new_fasta+ " 2> /dev/null") contigs.fasta
474 os.system("rm -rf "+outdir+ " 2> /dev/null")
475 ### begin blast
476 print "blasting..."
477 print "\n"
478 xmlfile=for_fq+"-extracted_vs_"+database+"_"+mapping_mode+".xml"
479 os.system('makeblastdb -in '+new_fasta+' -out '+new_fasta+'_db '+'-dbtype nucl >> data_log.txt 2>&1') #temp.txt is to forbid the blast result interrupt the output of our program###1/27/2015
480 print check_output("blastn -word_size 10 -query "+database+" -db "+new_fasta+"_db -out "+xmlfile+" -outfmt 5 >> data_log.txt 2>&1", shell=True)###1/27/2015
481 Final_list=xml_parse_score_comparision_seqsero(xmlfile)
482 Final_list_passed=[x for x in Final_list if float(x[0].split("_cov_")[1])>=3.5 and (x[1]>=int(x[0].split("__")[1]) or x[1]>=int(x[0].split("___")[1].split("_")[3]))]
483 fliC_choice="-"
484 fljB_choice="-"
485 fliC_contig="NA"
486 fljB_contig="NA"
487 fliC_length=0 #can be changed to coverage in future
488 fljB_length=0 #can be changed to coverage in future
489 O_choice=""#no need to decide O contig for now, should be only one
490 O_choice,O_nodes,special_gene_list,O_nodes_roles=decide_O_type_and_get_special_genes(Final_list)#decide the O antigen type and also return special-gene-list for further identification
491 O_choice=O_choice.split("-")[-1].strip()
492 H_contig_roles=decide_contig_roles_for_H_antigen(Final_list)#decide the H antigen contig is fliC or fljB
493 log_file=open("SeqSero_hybrid_assembly_log.txt","a")
494 print "O_contigs:"
495 log_file.write("O_contigs:\n")
496 for x in O_nodes_roles:
497 if "O-1,3,19_not_in_3,10" not in x[0]:#O-1,3,19_not_in_3,10 is just a small size marker
498 print x[0].split("___")[-1],x[0].split("__")[0],"blast score:",x[1],"identity%:",str(round(x[2]*100,2))+"%"
499 log_file.write(x[0].split("___")[-1]+" "+x[0].split("__")[0]+" "+"blast score: "+str(x[1])+"identity%:"+str(round(x[2]*100,2))+"%"+"\n")
500 print "H_contigs:"
501 log_file.write("H_contigs:\n")
502 H_contig_stat=[]
503 for i in range(len(H_contig_roles)):
504 x=H_contig_roles[i]
505 a=0
506 for y in Final_list_passed:
507 if x[1] in y[0] and y[0].startswith(x[0]):
508 if "first" in y[0] or "last" in y[0]: #this is the final filter to decide it's fliC or fljB, if can't pass, then can't decide
509 for y in Final_list_passed: #it's impossible to has the "first" and "last" allele as prediction, so re-do it
510 if x[1] in y[0]:#it's very possible to be third phase allele, so no need to make it must be fliC or fljB
511 print x[1],"can't_decide_fliC_or_fljB",y[0].split("_")[1],"blast_score:",y[1],"identity%:",str(round(y[2]*100,2))+"%"
512 log_file.write(x[1]+" "+x[0]+" "+y[0].split("_")[1]+" "+"blast_score: "+str(y[1])+" identity%:"+str(round(y[2]*100,2))+"%"+"\n")
513 H_contig_roles[i]="can't decide fliC or fljB, may be third phase"
514 break
515 else:
516 print x[1],x[0],y[0].split("_")[1],"blast_score:",y[1],"identity%:",str(round(y[2]*100,2))+"%"
517 log_file.write(x[1]+" "+x[0]+" "+y[0].split("_")[1]+" "+"blast_score: "+str(y[1])+" identity%:"+str(round(y[2]*100,2))+"%"+"\n")
518 break
519 for x in H_contig_roles:
520 #if multiple choices, temporately select the one with longest length for now, will revise in further change
521 if "fliC" == x[0] and int(x[1].split("_")[3])>=fliC_length and x[1] not in O_nodes:#remember to avoid the effect of O-type contig, so should not in O_node list
522 fliC_contig=x[1]
523 fliC_length=int(x[1].split("_")[3])
524 elif "fljB" == x[0] and int(x[1].split("_")[3])>=fljB_length and x[1] not in O_nodes:
525 fljB_contig=x[1]
526 fljB_length=int(x[1].split("_")[3])
527 for x in Final_list_passed:
528 if fliC_choice=="-" and "fliC_" in x[0] and fliC_contig in x[0] :
529 fliC_choice=x[0].split("_")[1]
530 elif fljB_choice=="-" and "fljB_" in x[0] and fljB_contig in x[0]:
531 fljB_choice=x[0].split("_")[1]
532 elif fliC_choice!="-" and fljB_choice!="-":
533 break
534 # print "\n"
535 # print "SeqSero Input files:",for_fq,rev_fq
536 # print "Most possible O antigen:",O_choice
537 # print "Most possible H1 antigen:",fliC_choice
538 # print "Most possible H2 antigen:",fljB_choice
539
540
541
542 #print Final_list
543 ###output
544 predict_form,predict_sero,star,star_line,claim=seqsero_from_formula_to_serotypes(O_choice,fliC_choice,fljB_choice,special_gene_list)
545
546 result = OrderedDict(sample_name=for_fq.split('_')[0],
547 O_antigen_prediction=O_choice,
548 H1_antigen_prediction=fliC_choice,
549 H2_antigen_prediction=fljB_choice,
550 predicted_antigenic_profile=predict_form,
551 predicted_serotypes=predict_sero)
552
553 print result
554
555 with open("Seqsero_result.tsv", 'w') as results_file:
556 wtr = csv.DictWriter(delimiter='\t', fieldnames=result.keys())
557
558 # new_file=open("Seqsero_result.txt","w")
559 # new_file.write("Input files:\t"+for_fq+" "+rev_fq+"\n"+"O antigen prediction:\t"+"O-"+O_choice+"\n"+"H1 antigen prediction(fliC):\t"+fliC_choice+"\n"+"H2 antigen prediction(fljB):\t"+fljB_choice+"\n"+"Predicted antigenic profile:\t"+predict_form+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+star+star_line+claim+"\n")#+##
560 # new_file.close()
561 # os.system("cat Seqsero_result.txt")
562
563 if __name__ == '__main__':
564 main()