jpayne@0: import os,sys,glob,time,itertools,subprocess
jpayne@0: from collections import OrderedDict
jpayne@0: from Initial_Conditions import phase1
jpayne@0: from Initial_Conditions import phase2
jpayne@0: from Initial_Conditions import phaseO
jpayne@0: from Initial_Conditions import sero
jpayne@0: from distutils.version import LooseVersion
jpayne@0: 
jpayne@0: import csv
jpayne@0: from subprocess import check_output
jpayne@0: 
jpayne@0: 
jpayne@0: def xml_parse_score_comparision_seqsero(xmlfile):
jpayne@0:   #used to do seqsero xml analysis
jpayne@0:   from Bio.Blast import NCBIXML
jpayne@0:   handle=open(xmlfile)
jpayne@0:   handle=NCBIXML.parse(handle)
jpayne@0:   handle=list(handle)
jpayne@0:   List=[]
jpayne@0:   List_score=[]
jpayne@0:   List_ids=[]
jpayne@0:   for i in range(len(handle)):
jpayne@0:     if len(handle[i].alignments)>0:
jpayne@0:       for j in range(len(handle[i].alignments)):
jpayne@0:         score=0
jpayne@0:         ids=0
jpayne@0:         List.append(handle[i].query.strip()+"___"+handle[i].alignments[j].hit_def)
jpayne@0:         for z in range(len(handle[i].alignments[j].hsps)):
jpayne@0:           if "last" in handle[i].query or "first" in handle[i].query:
jpayne@0:             score+=handle[i].alignments[j].hsps[z].bits
jpayne@0:             ids+=float(handle[i].alignments[j].hsps[z].identities)/handle[i].query_length
jpayne@0:           else:
jpayne@0:             if handle[i].alignments[j].hsps[z].align_length>=30:
jpayne@0:               #for the long alleles, filter noise parts
jpayne@0:               score+=handle[i].alignments[j].hsps[z].bits
jpayne@0:               ids+=float(handle[i].alignments[j].hsps[z].identities)/handle[i].query_length
jpayne@0:         List_score.append(score)
jpayne@0:         List_ids.append(ids)
jpayne@0:   temp=zip(List,List_score,List_ids)
jpayne@0:   Final_list=sorted(temp, key=lambda d:d[1], reverse = True)
jpayne@0:   return Final_list
jpayne@0: 
jpayne@0: 
jpayne@0: def Uniq(L,sort_on_fre="none"): #return the uniq list and the count number
jpayne@0:   Old=L
jpayne@0:   L.sort()
jpayne@0:   L = [L[i] for i in range(len(L)) if L[i] not in L[:i]]
jpayne@0:   count=[]
jpayne@0:   for j in range(len(L)):
jpayne@0:     y=0
jpayne@0:     for x in Old:
jpayne@0:       if L[j]==x:
jpayne@0:         y+=1
jpayne@0:     count.append(y)
jpayne@0:   if sort_on_fre!="none":
jpayne@0:     d=zip(*sorted(zip(count, L)))
jpayne@0:     L=d[1]
jpayne@0:     count=d[0]
jpayne@0:   return (L,count)
jpayne@0: 
jpayne@0: 
jpayne@0: def judge_fliC_or_fljB_from_head_tail_for_one_contig(nodes_vs_score_list):
jpayne@0:   #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:   #this is mainly used for 
jpayne@0:   a=nodes_vs_score_list
jpayne@0:   fliC_score=0
jpayne@0:   fljB_score=0
jpayne@0:   for z in a:
jpayne@0:     if "fliC" in z[0]:
jpayne@0:       fliC_score+=z[1]
jpayne@0:     elif "fljB" in z[0]:
jpayne@0:       fljB_score+=z[1]
jpayne@0:   if fliC_score>=fljB_score:
jpayne@0:     role="fliC"
jpayne@0:   else:
jpayne@0:     role="fljB"
jpayne@0:   return (role,abs(fliC_score-fljB_score))
jpayne@0: 
jpayne@0: def judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(node_name,Final_list_passed):
jpayne@0:   #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:   #also used when no head or tail got blasted score for the contig
jpayne@0:   role=""
jpayne@0:   for z in Final_list_passed:
jpayne@0:     if node_name in z[0]:
jpayne@0:       role=z[0].split("_")[0]
jpayne@0:       break
jpayne@0:   return role
jpayne@0: 
jpayne@0: 
jpayne@0: def fliC_or_fljB_judge_from_head_tail_sequence(nodes_list,tail_head_list,Final_list_passed):
jpayne@0:   #nodes_list is the c created by c,d=Uniq(nodes) in below function
jpayne@0:   first_target=""
jpayne@0:   role_list=[]
jpayne@0:   for x in nodes_list:
jpayne@0:     a=[]
jpayne@0:     role=""
jpayne@0:     for y in tail_head_list:
jpayne@0:       if x in y[0]:
jpayne@0:         a.append(y)
jpayne@0:     if len(a)==4:
jpayne@0:       #compare two heads (37 > 30)
jpayne@0:       #four contigs, most perfect assembly, high quality
jpayne@0:       """
jpayne@0:       for z in a:
jpayne@0:         if "fliC_first_37" in z[0]:
jpayne@0:           t1=z[1]
jpayne@0:         elif "fljB_first_37" in z[0]:
jpayne@0:           t2=z[1]
jpayne@0:       if t1>=t2:
jpayne@0:         role="fliC"
jpayne@0:       else:
jpayne@0:         role="fljB"
jpayne@0:       """
jpayne@0:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@0:       if diff<20:
jpayne@0:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
jpayne@0:     elif len(a)==3:
jpayne@0:       """
jpayne@0:       #compare the number, because hybrid problem
jpayne@0:       temp=[]
jpayne@0:       for z in a:
jpayne@0:         temp.append(z[0].split("_")[0])
jpayne@0:       m,n=Uniq(temp)#only two choices in m or n
jpayne@0:       if n[0]>n[1]:
jpayne@0:         role=m[0]
jpayne@0:       else:
jpayne@0:         role=m[1]
jpayne@0:       """
jpayne@0:       ###however, if the one with highest score is the fewer one, compare their accumulation score
jpayne@0:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@0:       if diff<20:
jpayne@0:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
jpayne@0:       ###end of above score comparison
jpayne@0:     elif len(a)==2:
jpayne@0:       #must on same node, if not, then decide with unit blast score, blast-score/length_of_special_sequence(30 or 37)
jpayne@0:       temp=[]
jpayne@0:       for z in a:
jpayne@0:         temp.append(z[0].split("_")[0])
jpayne@0:       m,n=Uniq(temp)#should only have one choice, but weird situation might occur too
jpayne@0:       if len(m)==1:
jpayne@0:         pass
jpayne@0:       else:
jpayne@0:         pass
jpayne@0:         #print "head and tail not belong to same role, now let's guess based on maximum likelihood"
jpayne@0:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@0:       if diff<20:
jpayne@0:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
jpayne@0:         """
jpayne@0:         max_unit_score=0
jpayne@0:         for z in a:
jpayne@0:           unit_score=z[-1]/int(z[0].split("__")[1])
jpayne@0:           if unit_score>=max_unit_score:
jpayne@0:             role=z[0].split("_")[0]
jpayne@0:             max_unit_score=unit_score
jpayne@0:         """
jpayne@0:         ###need to desgin a algorithm to guess most possible situation for nodes_list, See the situations of test evaluation
jpayne@0:     elif len(a)==1:
jpayne@0:       #that one
jpayne@0:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@0:       if diff<20:
jpayne@0:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
jpayne@0:       #role=a[0][0].split("_")[0]
jpayne@0:       #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:     else:#a==0
jpayne@0:       #use Final_list_passed best match
jpayne@0:       for z in Final_list_passed:
jpayne@0:         if x in z[0]:
jpayne@0:           role=z[0].split("_")[0]
jpayne@0:           break
jpayne@0:     #print x,role,len(a)
jpayne@0:     role_list.append((role,x))
jpayne@0:   if len(role_list)==2:
jpayne@0:     if role_list[0][0]==role_list[1][0]:#this is the most cocmmon error, two antigen were assigned to same phase
jpayne@0:       #just use score to do a final test
jpayne@0:       role_list=[]
jpayne@0:       for x in nodes_list: 
jpayne@0:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list_passed)
jpayne@0:         role_list.append((role,x))
jpayne@0:   return role_list
jpayne@0: 
jpayne@0: def decide_contig_roles_for_H_antigen(Final_list):
jpayne@0:   #used to decide which contig is FliC and which one is fljB
jpayne@0:   contigs=[]
jpayne@0:   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:   nodes=[]
jpayne@0:   for x in Final_list_passed:
jpayne@0:     if x[0].startswith("fl") and "last" not in x[0] and "first" not in x[0]:
jpayne@0:       nodes.append(x[0].split("___")[1].strip())
jpayne@0:   c,d=Uniq(nodes)#c is node_list
jpayne@0:   #print c
jpayne@0:   tail_head_list=[x for x in Final_list if ("last" in x[0] or "first" in x[0])]
jpayne@0:   roles=fliC_or_fljB_judge_from_head_tail_sequence(c,tail_head_list,Final_list_passed)
jpayne@0:   return roles
jpayne@0: 
jpayne@0: def Combine(b,c):
jpayne@0:   fliC_combinations=[]
jpayne@0:   fliC_combinations.append(",".join(c))
jpayne@0:   temp_combinations=[]
jpayne@0:   for i in range(len(b)):
jpayne@0:     for x in itertools.combinations(b,i+1):
jpayne@0:       temp_combinations.append(",".join(x))
jpayne@0:   for x in temp_combinations:
jpayne@0:     temp=[]
jpayne@0:     for y in c:
jpayne@0:       temp.append(y)
jpayne@0:     temp.append(x)
jpayne@0:     temp=",".join(temp)
jpayne@0:     temp=temp.split(",")
jpayne@0:     temp.sort()
jpayne@0:     temp=",".join(temp)
jpayne@0:     fliC_combinations.append(temp)
jpayne@0:   return fliC_combinations
jpayne@0: 
jpayne@0: def decide_O_type_and_get_special_genes(Final_list):
jpayne@0:   #decide O based on Final_list
jpayne@0:   O_choice="?"
jpayne@0:   O_list=[]
jpayne@0:   special_genes=[]
jpayne@0:   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:   nodes=[]
jpayne@0:   for x in Final_list_passed:
jpayne@0:     if x[0].startswith("O-"):
jpayne@0:       nodes.append(x[0].split("___")[1].strip())
jpayne@0:     elif not x[0].startswith("fl"):
jpayne@0:       special_genes.append(x)
jpayne@0:   #print "special_genes:",special_genes
jpayne@0:   c,d=Uniq(nodes)
jpayne@0:   #print "potential O antigen contig",c
jpayne@0:   final_O=[]
jpayne@0:   O_nodes_list=[]
jpayne@0:   for x in c:#c is the list for contigs
jpayne@0:     temp=0
jpayne@0:     for y in Final_list_passed:
jpayne@0:       if x in y[0] and y[0].startswith("O-"):
jpayne@0:         final_O.append(y)
jpayne@0:         break
jpayne@0:   ### O contig has the problem of two genes on same contig, so do additional test
jpayne@0:   potenial_new_gene=""
jpayne@0:   for x in final_O:
jpayne@0:     pointer=0 #for genes merged or not
jpayne@0:     #not consider O-1,3,19_not_in_3,10, too short compared with others
jpayne@0:     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:       pointer=x[0].split("___")[1].strip()#store the contig name
jpayne@0:       print pointer
jpayne@0:     if pointer!=0:#it has potential merge event
jpayne@0:       for y in Final_list:
jpayne@0:         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:           potenial_new_gene=y
jpayne@0:           print potenial_new_gene
jpayne@0:           break
jpayne@0:   if potenial_new_gene!="":
jpayne@0:     print "two differnt genes in same contig, fix it for O antigen"
jpayne@0:     final_O.append(potenial_new_gene)
jpayne@0:   ### end of the two genes on same contig test
jpayne@0:   if len(final_O)==0:
jpayne@0:     #print "$$$No Otype, due to no hit"#may need to be changed
jpayne@0:     O_choice="-"
jpayne@0:   else:
jpayne@0:     O_list=[]
jpayne@0:     for x in final_O:
jpayne@0:       O_list.append(x[0].split("__")[0])
jpayne@0:       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:         O_nodes_list.append(x[0].split("___")[1])
jpayne@0:     ### special test for O9,46 and O3,10 family
jpayne@0:     if "O-9,46_wbaV" in O_list:#not sure should use and float(O9_wbaV)/float(num_1) > 0.1
jpayne@0:       if "O-9,46_wzy" in O_list:#and float(O946_wzy)/float(num_1) > 0.1
jpayne@0:         O_choice="O-9,46"
jpayne@0:         #print "$$$Most possilble Otype:  O-9,46"
jpayne@0:       elif "O-9,46,27_partial_wzy" in O_list:#and float(O94627)/float(num_1) > 0.1
jpayne@0:         O_choice="O-9,46,27"
jpayne@0:         #print "$$$Most possilble Otype:  O-9,46,27"
jpayne@0:       else:
jpayne@0:         O_choice="O-9"#next, detect O9 vs O2?
jpayne@0:         O2=0
jpayne@0:         O9=0
jpayne@0:         for z in special_genes:
jpayne@0:           if "tyr-O-9" in z[0]:
jpayne@0:             O9=z[1]
jpayne@0:           elif "tyr-O-2" in z[0]:
jpayne@0:             O2=z[1]
jpayne@0:         if O2>O9:
jpayne@0:           O_choice="O-2"
jpayne@0:         elif O2<O9:
jpayne@0:           pass
jpayne@0:         else:
jpayne@0:           pass
jpayne@0:           #print "$$$No suitable one, because can't distinct it's O-9 or O-2, but O-9 has a more possibility."
jpayne@0:     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:       if "O-3,10_not_in_1,3,19" in O_list:#and float(O310_no_1319)/float(num_1) > 0.1
jpayne@0:         O_choice="O-3,10"
jpayne@0:         #print "$$$Most possilble Otype:  O-3,10 (contain O-3,10_not_in_1,3,19)"
jpayne@0:       else:
jpayne@0:         O_choice="O-1,3,19"
jpayne@0:         #print "$$$Most possilble Otype:  O-1,3,19 (not contain O-3,10_not_in_1,3,19)"
jpayne@0:     ### end of special test for O9,46 and O3,10 family
jpayne@0:     else:
jpayne@0:       try: 
jpayne@0:         max_score=0
jpayne@0:         for x in final_O:
jpayne@0:           if x[1]>=max_score:
jpayne@0:             max_score=x[1]
jpayne@0:             O_choice=x[0].split("_")[0]
jpayne@0:         if O_choice=="O-1,3,19":
jpayne@0:           O_choice=final_O[1][0].split("_")[0]
jpayne@0:         #print "$$$Most possilble Otype: ",O_choice
jpayne@0:       except:
jpayne@0:         pass
jpayne@0:         #print "$$$No suitable Otype, or failure of mapping (please check the quality of raw reads)"
jpayne@0:   #print "O:",O_choice,O_nodes_list
jpayne@0:   return O_choice,O_nodes_list,special_genes,final_O
jpayne@0: 
jpayne@0: def seqsero_from_formula_to_serotypes(Otype,fliC,fljB,special_gene_list):
jpayne@0:   #like test_output_06012017.txt
jpayne@0:   #can add more varialbles like sdf-type, sub-species-type in future (we can conclude it into a special-gene-list)
jpayne@0:   from Initial_Conditions import phase1
jpayne@0:   from Initial_Conditions import phase2
jpayne@0:   from Initial_Conditions import phaseO
jpayne@0:   from Initial_Conditions import sero
jpayne@0:   seronames=[]
jpayne@0:   for i in range(len(phase1)):
jpayne@0:     fliC_combine=[]
jpayne@0:     fljB_combine=[]
jpayne@0:     if phaseO[i]==Otype:
jpayne@0:       ### for fliC, detect every possible combinations to avoid the effect of "["
jpayne@0:       if phase1[i].count("[")==0:
jpayne@0:         fliC_combine.append(phase1[i])
jpayne@0:       elif phase1[i].count("[")>=1:
jpayne@0:         c=[]
jpayne@0:         b=[]
jpayne@0:         if phase1[i][0]=="[" and phase1[i][-1]=="]" and phase1[i].count("[")==1:
jpayne@0:           content=phase1[i].replace("[","").replace("]","")
jpayne@0:           fliC_combine.append(content)
jpayne@0:           fliC_combine.append("-")
jpayne@0:         else:
jpayne@0:           for x in phase1[i].split(","):
jpayne@0:             if "[" in x:
jpayne@0:               b.append(x.replace("[","").replace("]",""))
jpayne@0:             else:
jpayne@0:               c.append(x)
jpayne@0:           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:       ### end of fliC "[" detect
jpayne@0:       ### for fljB, detect every possible combinations to avoid the effect of "["
jpayne@0:       if phase2[i].count("[")==0:
jpayne@0:         fljB_combine.append(phase2[i])
jpayne@0:       elif phase2[i].count("[")>=1:
jpayne@0:         d=[]
jpayne@0:         e=[]
jpayne@0:         if phase2[i][0]=="[" and phase2[i][-1]=="]" and phase2[i].count("[")==1:
jpayne@0:           content=phase2[i].replace("[","").replace("]","")
jpayne@0:           fljB_combine.append(content)
jpayne@0:           fljB_combine.append("-")
jpayne@0:         else:
jpayne@0:           for x in phase2[i].split(","):
jpayne@0:             if "[" in x:
jpayne@0:               d.append(x.replace("[","").replace("]",""))
jpayne@0:             else:
jpayne@0:               e.append(x)
jpayne@0:           fljB_combine=Combine(d,e)
jpayne@0:       ### end of fljB "[" detect
jpayne@0:       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:       new_fliC.sort()
jpayne@0:       new_fliC=",".join(new_fliC)
jpayne@0:       new_fljB=fljB.split(",")
jpayne@0:       new_fljB.sort()
jpayne@0:       new_fljB=",".join(new_fljB)
jpayne@0:       if (new_fliC in fliC_combine or fliC in fliC_combine) and (new_fljB in fljB_combine or fljB in fljB_combine):
jpayne@0:         seronames.append(sero[i])
jpayne@0:   #analyze seronames
jpayne@0:   if len(seronames)==0:
jpayne@0:     seronames=["N/A (The predicted antigenic profile does not exist in the White-Kauffmann-Le Minor scheme)"]
jpayne@0:   star=""
jpayne@0:   star_line=""
jpayne@0:   if len(seronames)>1:#there are two possible predictions for serotypes
jpayne@0:     star="*"
jpayne@0:     star_line="The predicted serotypes share the same general formula:\t"+Otype+":"+fliC+":"+fljB+"\n"##
jpayne@0:   print "\n"
jpayne@0:   predict_form=Otype+":"+fliC+":"+fljB#
jpayne@0:   predict_sero=(" or ").join(seronames)
jpayne@0:   ###special test for Enteritidis
jpayne@0:   if predict_form=="9:g,m:-":
jpayne@0:     sdf="-"
jpayne@0:     for x in special_gene_list:
jpayne@0:       if x[0].startswith("sdf"):
jpayne@0:         sdf="+"
jpayne@0:     predict_form=predict_form+"\nSdf prediction:"+sdf
jpayne@0:     if sdf=="-":
jpayne@0:       star="*"
jpayne@0:       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:       predict_sero="See comments below"
jpayne@0:   ###end of special test for Enteritidis
jpayne@0:   elif predict_form=="4:i:-":
jpayne@0:     predict_sero="potential monophasic variant of Typhimurium"
jpayne@0:   elif predict_form=="4:r:-":
jpayne@0:     predict_sero="potential monophasic variant of Heidelberg"
jpayne@0:   elif predict_form=="4:b:-":
jpayne@0:     predict_sero="potential monophasic variant of Paratyphi B"
jpayne@0:   elif predict_form=="8:e,h:1,2":
jpayne@0:     predict_sero="Newport"
jpayne@0:     star="*"
jpayne@0:     star_line="Serotype Bardo shares the same antigenic profile with Newport, but Bardo is exceedingly rare."
jpayne@0:   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:   if "N/A" in predict_sero:
jpayne@0:     claim=""
jpayne@0:   if "Typhimurium" in predict_sero or predict_form=="4:i:-":
jpayne@0:     normal=0
jpayne@0:     mutation=0
jpayne@0:     for x in special_gene_list:
jpayne@0:       if "oafA-O-4_full" in x[0]:
jpayne@0:         normal=x[1]
jpayne@0:       elif "oafA-O-4_5-" in x[0]:
jpayne@0:         mutation=x[1]
jpayne@0:     if normal>mutation:
jpayne@0:       #print "$$$Typhimurium"
jpayne@0:       pass
jpayne@0:     elif normal<mutation:
jpayne@0:       predict_sero=predict_sero.strip()+"(O5-)"
jpayne@0:       star="*"#
jpayne@0:       star_line="Detected the deletion of O5-."
jpayne@0:       #print "$$$Typhimurium_O5-"
jpayne@0:     else:
jpayne@0:       #print "$$$Typhimurium, even no 7 bases difference"
jpayne@0:       pass
jpayne@0:   return predict_form,predict_sero,star,star_line,claim
jpayne@0: 
jpayne@0: def main():
jpayne@0:   database=sys.argv[1]#used to extract reads
jpayne@0:   mapping_mode=sys.argv[2]#mem or sampe
jpayne@0:   threads=sys.argv[3]
jpayne@0:   for_fq=sys.argv[4]
jpayne@0:   rev_fq=sys.argv[5]
jpayne@0:   current_time=time.strftime("%Y_%m_%d_%H_%M_%S", time.localtime())
jpayne@0:   sam=for_fq+".sam"
jpayne@0:   bam=for_fq+".bam"
jpayne@0:   sorted_bam=for_fq+"_sorted.bam"
jpayne@0:   mapped_fq1=for_fq+"_mapped.fq"
jpayne@0:   mapped_fq2=rev_fq+"_mapped.fq"
jpayne@0:   combined_fq=for_fq+"_combined.fq"
jpayne@0:   for_sai=for_fq+".sai"
jpayne@0:   rev_sai=rev_fq+".sai"
jpayne@0:   print "building database..."
jpayne@0:   #os.system("bwa index "+database+ " 2> /dev/null")
jpayne@0:   os.system("bwa index "+database+ " 2>> data_log.txt ")
jpayne@0:   print "mapping..."
jpayne@0:   if mapping_mode=="mem":
jpayne@0:     os.system("bwa mem -t "+threads+" "+database+" "+for_fq+" "+rev_fq+" > "+sam+ " 2>> data_log.txt")
jpayne@0:   elif mapping_mode=="sam":
jpayne@0:     os.system("bwa aln -t "+threads+" "+database+" "+for_fq+" > "+for_sai+ " 2>> data_log.txt")
jpayne@0:     os.system("bwa aln -t "+threads+" "+database+" "+rev_fq+" > "+rev_sai+ " 2>> data_log.txt")
jpayne@0:     os.system("bwa sampe "+database+" "+for_sai+" "+ rev_sai+" "+for_fq+" "+rev_fq+" > "+sam+ " 2>> data_log.txt")
jpayne@0:   os.system("samtools view -@ "+threads+" -F 4 -Sbh "+sam+" > "+bam)
jpayne@0:   os.system("samtools view -@ "+threads+" -h -o "+sam+" "+bam)
jpayne@0:   ### check the version of samtools then use differnt commands
jpayne@0:   samtools_version=subprocess.Popen(["samtools"],stdout=subprocess.PIPE,stderr=subprocess.PIPE)
jpayne@0:   out, err = samtools_version.communicate()
jpayne@0:   version = err.split("ersion:")[1].strip().split(" ")[0].strip()
jpayne@0:   print "check samtools version:",version
jpayne@0:   if LooseVersion(version)<=LooseVersion("1.2"):
jpayne@0:     os.system("samtools sort -@ "+threads+" -n "+bam+" "+for_fq+"_sorted")
jpayne@0:   else:
jpayne@0:     os.system("samtools sort -@ "+threads+" -n "+bam+" >"+sorted_bam)
jpayne@0:   ### end of samtools version check and its analysis
jpayne@0:   os.system("bamToFastq -i "+sorted_bam+" -fq "+combined_fq)
jpayne@0:   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:   outdir=current_time+"_temp"
jpayne@0:   print "assembling..."
jpayne@0:   if int(threads)>4:
jpayne@0:     t="4"
jpayne@0:   else:
jpayne@0:     t=threads
jpayne@0:   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:   new_fasta=for_fq+"_"+database+"_"+mapping_mode+".fasta"
jpayne@0:   os.system("mv "+outdir+"/contigs.fasta "+new_fasta+ " 2> /dev/null")
jpayne@0:   #os.system("mv "+outdir+"/scaffolds.fasta "+new_fasta+ " 2> /dev/null") contigs.fasta
jpayne@0:   os.system("rm -rf "+outdir+ " 2> /dev/null")
jpayne@0:   ### begin blast
jpayne@0:   print "blasting..."
jpayne@0:   print "\n"
jpayne@0:   xmlfile=for_fq+"-extracted_vs_"+database+"_"+mapping_mode+".xml"
jpayne@0:   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:   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:   Final_list=xml_parse_score_comparision_seqsero(xmlfile)
jpayne@0:   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:   fliC_choice="-"
jpayne@0:   fljB_choice="-"
jpayne@0:   fliC_contig="NA"
jpayne@0:   fljB_contig="NA"
jpayne@0:   fliC_length=0 #can be changed to coverage in future
jpayne@0:   fljB_length=0 #can be changed to coverage in future
jpayne@0:   O_choice=""#no need to decide O contig for now, should be only one
jpayne@0:   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:   O_choice=O_choice.split("-")[-1].strip()
jpayne@0:   H_contig_roles=decide_contig_roles_for_H_antigen(Final_list)#decide the H antigen contig is fliC or fljB
jpayne@0:   log_file=open("SeqSero_hybrid_assembly_log.txt","a")
jpayne@0:   print "O_contigs:"
jpayne@0:   log_file.write("O_contigs:\n")
jpayne@0:   for x in O_nodes_roles:
jpayne@0:     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:       print x[0].split("___")[-1],x[0].split("__")[0],"blast score:",x[1],"identity%:",str(round(x[2]*100,2))+"%"
jpayne@0:       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:   print "H_contigs:"
jpayne@0:   log_file.write("H_contigs:\n")
jpayne@0:   H_contig_stat=[]
jpayne@0:   for i in range(len(H_contig_roles)):
jpayne@0:     x=H_contig_roles[i]
jpayne@0:     a=0
jpayne@0:     for y in Final_list_passed:
jpayne@0:       if x[1] in y[0] and y[0].startswith(x[0]):
jpayne@0:         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:           for y in Final_list_passed: #it's impossible to has the "first" and "last" allele as prediction, so re-do it
jpayne@0:             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:               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:               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:               H_contig_roles[i]="can't decide fliC or fljB, may be third phase"
jpayne@0:               break
jpayne@0:         else:
jpayne@0:           print x[1],x[0],y[0].split("_")[1],"blast_score:",y[1],"identity%:",str(round(y[2]*100,2))+"%"
jpayne@0:           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:         break
jpayne@0:   for x in H_contig_roles:
jpayne@0:     #if multiple choices, temporately select the one with longest length for now, will revise in further change
jpayne@0:     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:       fliC_contig=x[1]
jpayne@0:       fliC_length=int(x[1].split("_")[3])
jpayne@0:     elif "fljB" == x[0] and int(x[1].split("_")[3])>=fljB_length and x[1] not in O_nodes:
jpayne@0:       fljB_contig=x[1]
jpayne@0:       fljB_length=int(x[1].split("_")[3])
jpayne@0:   for x in Final_list_passed:
jpayne@0:     if fliC_choice=="-" and "fliC_" in x[0] and fliC_contig in x[0] :
jpayne@0:       fliC_choice=x[0].split("_")[1]
jpayne@0:     elif fljB_choice=="-" and "fljB_" in x[0] and fljB_contig in x[0]:
jpayne@0:       fljB_choice=x[0].split("_")[1]
jpayne@0:     elif fliC_choice!="-" and fljB_choice!="-":
jpayne@0:       break
jpayne@0:   # print "\n"
jpayne@0:   # print "SeqSero Input files:",for_fq,rev_fq
jpayne@0:   # print "Most possible O antigen:",O_choice
jpayne@0:   # print "Most possible H1 antigen:",fliC_choice
jpayne@0:   # print "Most possible H2 antigen:",fljB_choice
jpayne@0: 
jpayne@0: 
jpayne@0: 
jpayne@0:   #print Final_list
jpayne@0:   ###output
jpayne@0:   predict_form,predict_sero,star,star_line,claim=seqsero_from_formula_to_serotypes(O_choice,fliC_choice,fljB_choice,special_gene_list)
jpayne@0:   
jpayne@0:   result = OrderedDict(sample_name=for_fq.split('_')[0],
jpayne@0:                 O_antigen_prediction=O_choice,
jpayne@0:                 H1_antigen_prediction=fliC_choice,
jpayne@0:                 H2_antigen_prediction=fljB_choice,
jpayne@0:                 predicted_antigenic_profile=predict_form,
jpayne@0:                 predicted_serotypes=predict_sero)
jpayne@0: 
jpayne@0:   print result
jpayne@0: 
jpayne@0:   with open("Seqsero_result.tsv", 'w') as results_file:
jpayne@0:     wtr = csv.DictWriter(delimiter='\t', fieldnames=result.keys())
jpayne@0:   
jpayne@0:   # new_file=open("Seqsero_result.txt","w")
jpayne@0:   # 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:   # new_file.close()
jpayne@0:   # os.system("cat Seqsero_result.txt")
jpayne@0: 
jpayne@0: if __name__ == '__main__':
jpayne@0:   main()