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