jpayne@1: #!/usr/bin/env python3
jpayne@1: 
jpayne@1: import sys
jpayne@1: import time
jpayne@1: import random
jpayne@1: import os
jpayne@1: import subprocess
jpayne@1: import gzip
jpayne@1: import io
jpayne@1: import pickle
jpayne@1: import argparse
jpayne@1: import itertools
jpayne@1: from distutils.version import LooseVersion
jpayne@1: 
jpayne@1: from collections import OrderedDict
jpayne@1: from csv import DictWriter
jpayne@1: 
jpayne@1: dirpath = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
jpayne@1: 
jpayne@1: ### SeqSero Kmer
jpayne@1: def parse_args():
jpayne@1:     "Parse the input arguments, use '-h' for help."
jpayne@1:     parser = argparse.ArgumentParser(usage='SeqSero2_package.py -t <data_type> -m <mode> -i <input_data> [-p <number of threads>] [-b <BWA_algorithm>]\n\nDevelopper: Shaokang Zhang (zskzsk@uga.edu), Hendrik C Den-Bakker (Hendrik.DenBakker@uga.edu) and Xiangyu Deng (xdeng@uga.edu)\n\nContact email:seqsero@gmail.com')#add "-m <data_type>" in future
jpayne@1:     parser.add_argument("-i",nargs="+",help="<string>: path/to/input_data")
jpayne@1:     parser.add_argument("-t",choices=['1','2','3','4','5','6'],help="<int>: '1'(pair-end reads, interleaved),'2'(pair-end reads, seperated),'3'(single-end reads), '4'(assembly),'5'(nanopore fasta),'6'(nanopore fastq)")
jpayne@7:     parser.add_argument("-b",choices=['sam','mem'],default="mem",help="<string>: mode used for mapping in allele workflow, 'sam'(bwa samse/sampe), 'mem'(bwa mem), default=mem") 
jpayne@7:     parser.add_argument("-p",default="1",help="<int>: threads used for allele workflow, if p>4, only 4 threads will be used for assembly step, default=1") 
jpayne@7:     parser.add_argument("-m",choices=['k','a'],default="a",help="<string>: 'k'(kmer workflow), 'a'(allele workflow), default=a") 
jpayne@1:     parser.add_argument("-d",help="<string>: output directory name, if not set, the output directory would be 'SeqSero_result_'+time stamp+one random number")
jpayne@1:     parser.add_argument("-c",action="store_true",help="<flag>: if '-c' was flagged, SeqSero2 will use clean mode and only output serotyping prediction, the directory containing log files will be deleted")
jpayne@1:     return parser.parse_args()
jpayne@1: 
jpayne@1: def reverse_complement(sequence):
jpayne@1:     complement = {
jpayne@1:         'A': 'T',
jpayne@1:         'C': 'G',
jpayne@1:         'G': 'C',
jpayne@1:         'T': 'A',
jpayne@1:         'N': 'N',
jpayne@1:         'M': 'K',
jpayne@1:         'R': 'Y',
jpayne@1:         'W': 'W',
jpayne@1:         'S': 'S',
jpayne@1:         'Y': 'R',
jpayne@1:         'K': 'M',
jpayne@1:         'V': 'B',
jpayne@1:         'H': 'D',
jpayne@1:         'D': 'H',
jpayne@1:         'B': 'V'
jpayne@1:     }
jpayne@1:     return "".join(complement[base] for base in reversed(sequence))
jpayne@1: 
jpayne@1: 
jpayne@1: def createKmerDict_reads(list_of_strings, kmer):
jpayne@1:     kmer_table = {}
jpayne@1:     for string in list_of_strings:
jpayne@1:         sequence = string.strip('\n')
jpayne@1:         for i in range(len(sequence) - kmer + 1):
jpayne@1:             new_mer = sequence[i:i + kmer].upper()
jpayne@1:             new_mer_rc = reverse_complement(new_mer)
jpayne@1:             if new_mer in kmer_table:
jpayne@1:                 kmer_table[new_mer.upper()] += 1
jpayne@1:             else:
jpayne@1:                 kmer_table[new_mer.upper()] = 1
jpayne@1:             if new_mer_rc in kmer_table:
jpayne@1:                 kmer_table[new_mer_rc.upper()] += 1
jpayne@1:             else:
jpayne@1:                 kmer_table[new_mer_rc.upper()] = 1
jpayne@1:     return kmer_table
jpayne@1: 
jpayne@1: 
jpayne@1: def multifasta_dict(multifasta):
jpayne@1:     multifasta_list = [
jpayne@1:         line.strip() for line in open(multifasta, 'r') if len(line.strip()) > 0
jpayne@1:     ]
jpayne@1:     headers = [i for i in multifasta_list if i[0] == '>']
jpayne@1:     multifasta_dict = {}
jpayne@1:     for h in headers:
jpayne@1:         start = multifasta_list.index(h)
jpayne@1:         for element in multifasta_list[start + 1:]:
jpayne@1:             if element[0] == '>':
jpayne@1:                 break
jpayne@1:             else:
jpayne@1:                 if h[1:] in multifasta_dict:
jpayne@1:                     multifasta_dict[h[1:]] += element
jpayne@1:                 else:
jpayne@1:                     multifasta_dict[h[1:]] = element
jpayne@1:     return multifasta_dict
jpayne@1: 
jpayne@1: 
jpayne@1: def multifasta_single_string(multifasta):
jpayne@1:     multifasta_list = [
jpayne@1:         line.strip() for line in open(multifasta, 'r')
jpayne@1:         if (len(line.strip()) > 0) and (line.strip()[0] != '>')
jpayne@1:     ]
jpayne@1:     return ''.join(multifasta_list)
jpayne@1: 
jpayne@1: 
jpayne@1: def chunk_a_long_sequence(long_sequence, chunk_size=60):
jpayne@1:     chunk_list = []
jpayne@1:     steps = len(long_sequence) // 60  #how many chunks
jpayne@1:     for i in range(steps):
jpayne@1:         chunk_list.append(long_sequence[i * chunk_size:(i + 1) * chunk_size])
jpayne@1:     chunk_list.append(long_sequence[steps * chunk_size:len(long_sequence)])
jpayne@1:     return chunk_list
jpayne@1: 
jpayne@1: 
jpayne@1: def target_multifasta_kmerizer(multifasta, k, kmerDict):
jpayne@1:     forward_length = 300  #if find the target, put forward 300 bases
jpayne@1:     reverse_length = 2200  #if find the target, put backward 2200 bases
jpayne@1:     chunk_size = 60  #it will firstly chunk the single long sequence to multiple smaller sequences, it controls the size of those smaller sequences
jpayne@1:     target_mers = []
jpayne@1:     long_single_string = multifasta_single_string(multifasta)
jpayne@1:     multifasta_list = chunk_a_long_sequence(long_single_string, chunk_size)
jpayne@1:     unit_length = len(multifasta_list[0])
jpayne@1:     forward_lines = int(forward_length / unit_length) + 1
jpayne@1:     reverse_lines = int(forward_length / unit_length) + 1
jpayne@1:     start_num = 0
jpayne@1:     end_num = 0
jpayne@1:     for i in range(len(multifasta_list)):
jpayne@1:         if i not in range(start_num, end_num):  #avoid computational repetition
jpayne@1:             line = multifasta_list[i]
jpayne@1:             start = int((len(line) - k) // 2)
jpayne@1:             s1 = line[start:k + start]
jpayne@1:             if s1 in kmerDict:  #detect it is a potential read or not (use the middle part)
jpayne@1:                 if i - forward_lines >= 0:
jpayne@1:                     start_num = i - forward_lines
jpayne@1:                 else:
jpayne@1:                     start_num = 0
jpayne@1:                 if i + reverse_lines <= len(multifasta_list) - 1:
jpayne@1:                     end_num = i + reverse_lines
jpayne@1:                 else:
jpayne@1:                     end_num = len(multifasta_list) - 1
jpayne@1:                 target_list = [
jpayne@1:                     x.strip() for x in multifasta_list[start_num:end_num]
jpayne@1:                 ]
jpayne@1:                 target_line = "".join(target_list)
jpayne@1:                 target_mers += [
jpayne@1:                     k1 for k1 in createKmerDict_reads([str(target_line)], k)
jpayne@1:                 ]  ##changed k to k1, just want to avoid the mixes of this "k" (kmer) to the "k" above (kmer length)
jpayne@1:         else:
jpayne@1:             pass
jpayne@1:     return set(target_mers)
jpayne@1: 
jpayne@1: 
jpayne@1: def target_read_kmerizer(file, k, kmerDict):
jpayne@1:     i = 1
jpayne@1:     n_reads = 0
jpayne@1:     total_coverage = 0
jpayne@1:     target_mers = []
jpayne@1:     if file.endswith(".gz"):
jpayne@1:         file_content = io.BufferedReader(gzip.open(file))
jpayne@1:     else:
jpayne@1:         file_content = open(file, "r").readlines()
jpayne@1:     for line in file_content:
jpayne@1:         start = int((len(line) - k) // 2)
jpayne@1:         if i % 4 == 2:
jpayne@1:             if file.endswith(".gz"):
jpayne@1:                 s1 = line[start:k + start].decode()
jpayne@1:                 line = line.decode()
jpayne@1:             else:
jpayne@1:                 s1 = line[start:k + start]
jpayne@1:             if s1 in kmerDict:  #detect it is a potential read or not (use the middle part)
jpayne@1:                 n_reads += 1
jpayne@1:                 total_coverage += len(line)
jpayne@1:                 target_mers += [
jpayne@1:                     k1 for k1 in createKmerDict_reads([str(line)], k)
jpayne@1:                 ]  #changed k to k1, just want to avoid the mixes of this "k" (kmer) to the "k" above (kmer length)
jpayne@1:         i += 1
jpayne@1:         if total_coverage >= 4000000:
jpayne@1:             break
jpayne@1:     return set(target_mers)
jpayne@1: 
jpayne@1: 
jpayne@1: def minion_fasta_kmerizer(file, k, kmerDict):
jpayne@1:     i = 1
jpayne@1:     n_reads = 0
jpayne@1:     total_coverage = 0
jpayne@1:     target_mers = {}
jpayne@1:     for line in open(file):
jpayne@1:         if i % 2 == 0:
jpayne@1:             for kmer, rc_kmer in kmers(line.strip().upper(), k):
jpayne@1:                 if (kmer in kmerDict) or (rc_kmer in kmerDict):
jpayne@1:                     if kmer in target_mers:
jpayne@1:                         target_mers[kmer] += 1
jpayne@1:                     else:
jpayne@1:                         target_mers[kmer] = 1
jpayne@1:                     if rc_kmer in target_mers:
jpayne@1:                         target_mers[rc_kmer] += 1
jpayne@1:                     else:
jpayne@1:                         target_mers[rc_kmer] = 1
jpayne@1:         i += 1
jpayne@1:     return set([h for h in target_mers])
jpayne@1: 
jpayne@1: 
jpayne@1: def minion_fastq_kmerizer(file, k, kmerDict):
jpayne@1:     i = 1
jpayne@1:     n_reads = 0
jpayne@1:     total_coverage = 0
jpayne@1:     target_mers = {}
jpayne@1:     for line in open(file):
jpayne@1:         if i % 4 == 2:
jpayne@1:             for kmer, rc_kmer in kmers(line.strip().upper(), k):
jpayne@1:                 if (kmer in kmerDict) or (rc_kmer in kmerDict):
jpayne@1:                     if kmer in target_mers:
jpayne@1:                         target_mers[kmer] += 1
jpayne@1:                     else:
jpayne@1:                         target_mers[kmer] = 1
jpayne@1:                     if rc_kmer in target_mers:
jpayne@1:                         target_mers[rc_kmer] += 1
jpayne@1:                     else:
jpayne@1:                         target_mers[rc_kmer] = 1
jpayne@1:         i += 1
jpayne@1:     return set([h for h in target_mers])
jpayne@1: 
jpayne@1: 
jpayne@1: def multifasta_single_string2(multifasta):
jpayne@1:     single_string = ''
jpayne@1:     with open(multifasta, 'r') as f:
jpayne@1:         for line in f:
jpayne@1:             if line.strip()[0] == '>':
jpayne@1:                 pass
jpayne@1:             else:
jpayne@1:                 single_string += line.strip()
jpayne@1:     return single_string
jpayne@1: 
jpayne@1: 
jpayne@1: def kmers(seq, k):
jpayne@1:     rev_comp = reverse_complement(seq)
jpayne@1:     for start in range(1, len(seq) - k + 1):
jpayne@1:         yield seq[start:start + k], rev_comp[-(start + k):-start]
jpayne@1: 
jpayne@1: 
jpayne@1: def multifasta_to_kmers_dict(multifasta,k_size):#used to create database kmer set
jpayne@1:     multi_seq_dict = multifasta_dict(multifasta)
jpayne@1:     lib_dict = {}
jpayne@1:     for h in multi_seq_dict:
jpayne@1:         lib_dict[h] = set(
jpayne@1:             [k for k in createKmerDict_reads([multi_seq_dict[h]], k_size)])
jpayne@1:     return lib_dict
jpayne@1: 
jpayne@1: 
jpayne@1: def Combine(b, c):
jpayne@1:     fliC_combinations = []
jpayne@1:     fliC_combinations.append(",".join(c))
jpayne@1:     temp_combinations = []
jpayne@1:     for i in range(len(b)):
jpayne@1:         for x in itertools.combinations(b, i + 1):
jpayne@1:             temp_combinations.append(",".join(x))
jpayne@1:     for x in temp_combinations:
jpayne@1:         temp = []
jpayne@1:         for y in c:
jpayne@1:             temp.append(y)
jpayne@1:         temp.append(x)
jpayne@1:         temp = ",".join(temp)
jpayne@1:         temp = temp.split(",")
jpayne@1:         temp.sort()
jpayne@1:         temp = ",".join(temp)
jpayne@1:         fliC_combinations.append(temp)
jpayne@1:     return fliC_combinations
jpayne@1: 
jpayne@1: 
jpayne@1: def seqsero_from_formula_to_serotypes(Otype, fliC, fljB, special_gene_list,subspecies):
jpayne@1:     #like test_output_06012017.txt
jpayne@1:     #can add more varialbles like sdf-type, sub-species-type in future (we can conclude it into a special-gene-list)
jpayne@4:     from Initial_Conditions import phase1,phase2,phaseO,sero,subs,remove_list,rename_dict
jpayne@4:     rename_dict_not_anymore=[rename_dict[x] for x in rename_dict]
jpayne@4:     rename_dict_all=rename_dict_not_anymore+list(rename_dict) #used for decide whether to 
jpayne@1:     seronames = []
jpayne@1:     seronames_none_subspecies=[]
jpayne@1:     for i in range(len(phase1)):
jpayne@1:         fliC_combine = []
jpayne@1:         fljB_combine = []
jpayne@1:         if phaseO[i] == Otype: # no VII in KW, but it's there
jpayne@1:             ### for fliC, detect every possible combinations to avoid the effect of "["
jpayne@1:             if phase1[i].count("[") == 0:
jpayne@1:                 fliC_combine.append(phase1[i])
jpayne@1:             elif phase1[i].count("[") >= 1:
jpayne@1:                 c = []
jpayne@1:                 b = []
jpayne@1:                 if phase1[i][0] == "[" and phase1[i][-1] == "]" and phase1[i].count(
jpayne@1:                         "[") == 1:
jpayne@1:                     content = phase1[i].replace("[", "").replace("]", "")
jpayne@1:                     fliC_combine.append(content)
jpayne@1:                     fliC_combine.append("-")
jpayne@1:                 else:
jpayne@1:                     for x in phase1[i].split(","):
jpayne@1:                         if "[" in x:
jpayne@1:                             b.append(x.replace("[", "").replace("]", ""))
jpayne@1:                         else:
jpayne@1:                             c.append(x)
jpayne@1:                     fliC_combine = Combine(
jpayne@1:                         b, c
jpayne@1:                     )  #Combine will offer every possible combinations of the formula, like f,[g],t: f,t  f,g,t
jpayne@1:             ### end of fliC "[" detect
jpayne@1:             ### for fljB, detect every possible combinations to avoid the effect of "["
jpayne@1:             if phase2[i].count("[") == 0:
jpayne@1:                 fljB_combine.append(phase2[i])
jpayne@1:             elif phase2[i].count("[") >= 1:
jpayne@1:                 d = []
jpayne@1:                 e = []
jpayne@1:                 if phase2[i][0] == "[" and phase2[i][-1] == "]" and phase2[i].count(
jpayne@1:                         "[") == 1:
jpayne@1:                     content = phase2[i].replace("[", "").replace("]", "")
jpayne@1:                     fljB_combine.append(content)
jpayne@1:                     fljB_combine.append("-")
jpayne@1:                 else:
jpayne@1:                     for x in phase2[i].split(","):
jpayne@1:                         if "[" in x:
jpayne@1:                             d.append(x.replace("[", "").replace("]", ""))
jpayne@1:                         else:
jpayne@1:                             e.append(x)
jpayne@1:                     fljB_combine = Combine(d, e)
jpayne@1:             ### end of fljB "[" detect
jpayne@1:             new_fliC = fliC.split(
jpayne@1:                 ","
jpayne@1:             )  #because some antigen like r,[i] not follow alphabetical order, so use this one to judge and can avoid missings
jpayne@1:             new_fliC.sort()
jpayne@1:             new_fliC = ",".join(new_fliC)
jpayne@1:             new_fljB = fljB.split(",")
jpayne@1:             new_fljB.sort()
jpayne@1:             new_fljB = ",".join(new_fljB)
jpayne@1:             if (new_fliC in fliC_combine
jpayne@1:                     or fliC in fliC_combine) and (new_fljB in fljB_combine
jpayne@1:                                                   or fljB in fljB_combine):
jpayne@4:               ######start, remove_list,rename_dict, added on 11/11/2018
jpayne@4:               if sero[i] not in remove_list:
jpayne@4:                 temp_sero=sero[i]
jpayne@4:                 if temp_sero in rename_dict:
jpayne@4:                   temp_sero=rename_dict[temp_sero] #rename if in the rename list
jpayne@4:                 if temp_sero not in seronames:#the new sero may already included, if yes, then not consider
jpayne@4:                   if subs[i] == subspecies:
jpayne@4:                     seronames.append(temp_sero)
jpayne@4:                   seronames_none_subspecies.append(temp_sero)
jpayne@4:                 else:
jpayne@4:                   pass
jpayne@4:               else:
jpayne@4:                 pass
jpayne@4:               ######end, added on 11/11/2018
jpayne@1:     #analyze seronames
jpayne@1:     subspecies_pointer=""
jpayne@1:     if len(seronames) == 0 and len(seronames_none_subspecies)!=0:
jpayne@1:       seronames=seronames_none_subspecies
jpayne@1:       subspecies_pointer="1"
jpayne@1:     if len(seronames) == 0:
jpayne@1:         seronames = [
jpayne@1:             "N/A (The predicted antigenic profile does not exist in the White-Kauffmann-Le Minor scheme)"
jpayne@1:         ]
jpayne@1:     star = ""
jpayne@1:     star_line = ""
jpayne@1:     if len(seronames) > 1:  #there are two possible predictions for serotypes
jpayne@1:         star = "*"
jpayne@7:         #changed 04072019
jpayne@7:         #star_line = "The predicted serotypes share the same general formula:\t" + Otype + ":" + fliC + ":" + fljB + "\n"
jpayne@1:     if subspecies_pointer=="1" and len(seronames_none_subspecies)!=0:
jpayne@1:       star="*"
jpayne@7:       star_line=" The predicted O and H antigens correspond to serotype '"+(" or ").join(seronames)+"' in the Kauffmann-White scheme. The predicted subspecies by SalmID (github.com/hcdenbakker/SalmID) may not be consistent with subspecies designation in the Kauffmann-White scheme." + star_line
jpayne@7:       #star_line="The formula with this subspieces prediction can't get a serotype in KW manual, and the serotyping prediction was made without considering it."+star_line
jpayne@1:     if  Otype=="":
jpayne@1:       Otype="-"
jpayne@1:     predict_form = Otype + ":" + fliC + ":" + fljB
jpayne@1:     predict_sero = (" or ").join(seronames)
jpayne@1:     ###special test for Enteritidis
jpayne@1:     if predict_form == "9:g,m:-":
jpayne@1:         sdf = "-"
jpayne@1:         for x in special_gene_list:
jpayne@1:             if x.startswith("sdf"):
jpayne@1:                 sdf = "+"
jpayne@7:                 #star_line="Detected sdf gene, a marker to differentiate Gallinarum and Enteritidis"
jpayne@7:                 star_line=" sdf gene detected." # ed_SL_04152019: new output format
jpayne@7:         #predict_form = predict_form + " Sdf prediction:" + sdf
jpayne@7:         predict_form = predict_form #changed 04072019
jpayne@1:         if sdf == "-":
jpayne@1:             star = "*"
jpayne@7:             #star_line="Didn't detected sdf gene, a marker to differentiate Gallinarum and Enteritidis"
jpayne@7:             star_line=" sdf gene not detected." # ed_SL_04152019: new output format
jpayne@7:             #changed in 04072019, for new output
jpayne@7:             #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@7:             #predict_sero = "Gallinarum/Enteritidis" #04132019, for new output requirement
jpayne@7:             predict_sero = "Gallinarum or Enteritidis" # ed_SL_04152019: new output format
jpayne@1:     ###end of special test for Enteritidis
jpayne@1:     elif predict_form == "4:i:-":
jpayne@7:         predict_sero = "4:i:-"
jpayne@1:     elif predict_form == "4:r:-":
jpayne@7:         predict_sero = "4:r:-"
jpayne@1:     elif predict_form == "4:b:-":
jpayne@7:         predict_sero = "4:b:-"
jpayne@4:     #elif predict_form == "8:e,h:1,2": #removed after official merge of newport and bardo
jpayne@4:         #predict_sero = "Newport"
jpayne@4:         #star = "*"
jpayne@4:         #star_line = "Serotype Bardo shares the same antigenic profile with Newport, but Bardo is exceedingly rare."
jpayne@7:     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.\n"
jpayne@1:     if "N/A" in predict_sero:
jpayne@1:         claim = ""
jpayne@1:     #special test for Typhimurium
jpayne@1:     if "Typhimurium" in predict_sero or predict_form == "4:i:-":
jpayne@1:         normal = 0
jpayne@1:         mutation = 0
jpayne@1:         for x in special_gene_list:
jpayne@1:             if "oafA-O-4_full" in x:
jpayne@1:                 normal = float(special_gene_list[x])
jpayne@1:             elif "oafA-O-4_5-" in x:
jpayne@1:                 mutation = float(special_gene_list[x])
jpayne@1:         if normal > mutation:
jpayne@1:             pass
jpayne@1:         elif normal < mutation:
jpayne@7:             #predict_sero = predict_sero.strip() + "(O5-)"
jpayne@7:             predict_sero = predict_sero.strip() #diable special sero for new output requirement, 04132019
jpayne@1:             star = "*"
jpayne@7:             #star_line = "Detected the deletion of O5-."
jpayne@7:             star_line = " Detected a deletion that causes O5- variant of Typhimurium." # ed_SL_04152019: new output format
jpayne@1:         else:
jpayne@1:             pass
jpayne@1:     #special test for Paratyphi B
jpayne@1:     if "Paratyphi B" in predict_sero or predict_form == "4:b:-":
jpayne@1:         normal = 0
jpayne@1:         mutation = 0
jpayne@1:         for x in special_gene_list:
jpayne@1:             if "gntR-family-regulatory-protein_dt-positive" in x:
jpayne@1:                 normal = float(special_gene_list[x])
jpayne@1:             elif "gntR-family-regulatory-protein_dt-negative" in x:
jpayne@1:                 mutation = float(special_gene_list[x])
jpayne@1:         #print(normal,mutation)
jpayne@1:         if normal > mutation:
jpayne@7:             #predict_sero = predict_sero.strip() + "(dt+)" #diable special sero for new output requirement, 04132019
jpayne@7:             predict_sero = predict_sero.strip()+' var. L(+) tartrate+' if "Paratyphi B" in predict_sero else predict_sero.strip() # ed_SL_04152019: new output format
jpayne@1:             star = "*"
jpayne@7:             #star_line = "Didn't detect the SNP for dt- which means this isolate is a Paratyphi B variant L(+) tartrate(+)."
jpayne@7:             star_line = " The SNP that causes d-Tartrate nonfermentating phenotype was not detected. " # ed_SL_04152019: new output format
jpayne@1:         elif normal < mutation:
jpayne@7:             #predict_sero = predict_sero.strip() + "(dt-)" #diable special sero for new output requirement, 04132019
jpayne@7:             predict_sero = predict_sero.strip()
jpayne@1:             star = "*"
jpayne@7:             #star_line = "Detected the SNP for dt- which means this isolate is a systemic pathovar of Paratyphi B."
jpayne@7:             star_line = " Detected the SNP for d-Tartrate nonfermenting phenotype." # ed_SL_04152019: new output format
jpayne@1:         else:
jpayne@1:             star = "*"
jpayne@7:             star_line = " Failed to detect the SNP for dt-, can't decide it's a Paratyphi B variant L(+) tartrate(+) or not."
jpayne@1:     #special test for O13,22 and O13,23
jpayne@1:     if Otype=="13":
jpayne@1:       #ex_dir = os.path.dirname(os.path.realpath(__file__))
jpayne@1:       f = open(dirpath + '/special.pickle', 'rb')
jpayne@1:       special = pickle.load(f)
jpayne@1:       O22_O23=special['O22_O23']
jpayne@4:       if predict_sero.split(" or ")[0] in O22_O23[-1] and predict_sero.split(" or ")[0] not in rename_dict_all:#if in rename_dict_all, then it means already merged, no need to analyze 
jpayne@1:         O22_score=0
jpayne@1:         O23_score=0
jpayne@1:         for x in special_gene_list:
jpayne@1:             if "O:22" in x:
jpayne@1:                 O22_score = O22_score+float(special_gene_list[x])
jpayne@1:             elif "O:23" in x:
jpayne@1:                 O23_score = O23_score+float(special_gene_list[x])
jpayne@1:         #print(O22_score,O23_score)
jpayne@1:         for z in O22_O23[0]:
jpayne@1:           if predict_sero.split(" or ")[0] in z:
jpayne@1:             if O22_score > O23_score:
jpayne@1:               star = "*"
jpayne@7:               #star_line = "Detected O22 specific genes to further differenciate '"+predict_sero+"'." #diabled for new output requirement, 04132019
jpayne@1:               predict_sero = z[0]
jpayne@1:             elif O22_score < O23_score:
jpayne@1:               star = "*"
jpayne@7:               #star_line = "Detected O23 specific genes to further differenciate '"+predict_sero+"'." #diabled for new output requirement, 04132019
jpayne@1:               predict_sero = z[1]
jpayne@1:             else:
jpayne@1:               star = "*"
jpayne@7:               #star_line = "Fail to detect O22 and O23 differences." #diabled for new output requirement, 04132019
jpayne@7:     if " or " in predict_sero:
jpayne@7:       star_line = star_line + " The predicted serotypes share the same general formula:\t" + Otype + ":" + fliC + ":" + fljB + "\n"
jpayne@1:     #special test for O6,8 
jpayne@4:     #merge_O68_list=["Blockley","Bovismorbificans","Hadar","Litchfield","Manhattan","Muenchen"] #remove 11/11/2018, because already in merge list
jpayne@4:     #for x in merge_O68_list:
jpayne@4:     #  if x in predict_sero:
jpayne@4:     #    predict_sero=x
jpayne@4:     #    star=""
jpayne@4:     #    star_line=""
jpayne@1:     #special test for Montevideo; most of them are monophasic
jpayne@4:     #if "Montevideo" in predict_sero and "1,2,7" in predict_form: #remove 11/11/2018, because already in merge list
jpayne@4:       #star="*"
jpayne@4:       #star_line="Montevideo is almost always monophasic, having an antigen called for the fljB position may be a result of Salmonella-Salmonella contamination."
jpayne@1:     return predict_form, predict_sero, star, star_line, claim
jpayne@1: ### End of SeqSero Kmer part
jpayne@1: 
jpayne@1: ### Begin of SeqSero2 allele prediction and output
jpayne@1: def xml_parse_score_comparision_seqsero(xmlfile):
jpayne@1:   #used to do seqsero xml analysis
jpayne@1:   from Bio.Blast import NCBIXML
jpayne@1:   handle=open(xmlfile)
jpayne@1:   handle=NCBIXML.parse(handle)
jpayne@1:   handle=list(handle)
jpayne@1:   List=[]
jpayne@1:   List_score=[]
jpayne@1:   List_ids=[]
jpayne@1:   List_query_region=[]
jpayne@1:   for i in range(len(handle)):
jpayne@1:     if len(handle[i].alignments)>0:
jpayne@1:       for j in range(len(handle[i].alignments)):
jpayne@1:         score=0
jpayne@1:         ids=0
jpayne@1:         cover_region=set() #fixed problem that repeated calculation leading percentage > 1
jpayne@1:         List.append(handle[i].query.strip()+"___"+handle[i].alignments[j].hit_def)
jpayne@1:         for z in range(len(handle[i].alignments[j].hsps)):
jpayne@1:           hsp=handle[i].alignments[j].hsps[z]
jpayne@1:           temp=set(range(hsp.query_start,hsp.query_end))
jpayne@1:           if len(cover_region)==0:
jpayne@1:             cover_region=cover_region|temp
jpayne@1:             fraction=1
jpayne@1:           else:
jpayne@1:             fraction=1-len(cover_region&temp)/float(len(temp))
jpayne@1:             cover_region=cover_region|temp
jpayne@1:           if "last" in handle[i].query or "first" in handle[i].query:
jpayne@1:             score+=hsp.bits*fraction
jpayne@1:             ids+=float(hsp.identities)/handle[i].query_length*fraction
jpayne@1:           else:
jpayne@1:             score+=hsp.bits*fraction
jpayne@1:             ids+=float(hsp.identities)/handle[i].query_length*fraction
jpayne@1:         List_score.append(score)
jpayne@1:         List_ids.append(ids)
jpayne@1:         List_query_region.append(cover_region)
jpayne@1:   temp=zip(List,List_score,List_ids,List_query_region)
jpayne@1:   Final_list=sorted(temp, key=lambda d:d[1], reverse = True)
jpayne@1:   return Final_list
jpayne@1: 
jpayne@1: 
jpayne@1: def Uniq(L,sort_on_fre="none"): #return the uniq list and the count number
jpayne@1:   Old=L
jpayne@1:   L.sort()
jpayne@1:   L = [L[i] for i in range(len(L)) if L[i] not in L[:i]]
jpayne@1:   count=[]
jpayne@1:   for j in range(len(L)):
jpayne@1:     y=0
jpayne@1:     for x in Old:
jpayne@1:       if L[j]==x:
jpayne@1:         y+=1
jpayne@1:     count.append(y)
jpayne@1:   if sort_on_fre!="none":
jpayne@1:     d=zip(*sorted(zip(count, L)))
jpayne@1:     L=d[1]
jpayne@1:     count=d[0]
jpayne@1:   return (L,count)
jpayne@1: 
jpayne@1: def judge_fliC_or_fljB_from_head_tail_for_one_contig(nodes_vs_score_list):
jpayne@1:   #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@1:   #this is mainly used for 
jpayne@1:   a=nodes_vs_score_list
jpayne@1:   fliC_score=0
jpayne@1:   fljB_score=0
jpayne@1:   for z in a:
jpayne@1:     if "fliC" in z[0]:
jpayne@1:       fliC_score+=z[1]
jpayne@1:     elif "fljB" in z[0]:
jpayne@1:       fljB_score+=z[1]
jpayne@1:   if fliC_score>=fljB_score:
jpayne@1:     role="fliC"
jpayne@1:   else:
jpayne@1:     role="fljB"
jpayne@1:   return (role,abs(fliC_score-fljB_score))
jpayne@1: 
jpayne@1: def judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(node_name,Final_list,Final_list_passed):
jpayne@1:   #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@1:   #also used when no head or tail got blasted score for the contig
jpayne@1:   role=""
jpayne@1:   for z in Final_list_passed:
jpayne@1:     if node_name in z[0]:
jpayne@1:       role=z[0].split("_")[0]
jpayne@1:       break
jpayne@1:   return role
jpayne@1: 
jpayne@1: def fliC_or_fljB_judge_from_head_tail_sequence(nodes_list,tail_head_list,Final_list,Final_list_passed):
jpayne@1:   #nodes_list is the c created by c,d=Uniq(nodes) in below function
jpayne@1:   first_target=""
jpayne@1:   role_list=[]
jpayne@1:   for x in nodes_list:
jpayne@1:     a=[]
jpayne@1:     role=""
jpayne@1:     for y in tail_head_list:
jpayne@1:       if x in y[0]:
jpayne@1:         a.append(y)
jpayne@1:     if len(a)==4:
jpayne@1:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@1:       if diff<20:
jpayne@1:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list,Final_list_passed)
jpayne@1:     elif len(a)==3:
jpayne@1:       ###however, if the one with highest score is the fewer one, compare their accumulation score
jpayne@1:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@1:       if diff<20:
jpayne@1:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list,Final_list_passed)
jpayne@1:       ###end of above score comparison
jpayne@1:     elif len(a)==2:
jpayne@1:       #must on same node, if not, then decide with unit blast score, blast-score/length_of_special_sequence(30 or 37)
jpayne@1:       temp=[]
jpayne@1:       for z in a:
jpayne@1:         temp.append(z[0].split("_")[0])
jpayne@1:       m,n=Uniq(temp)#should only have one choice, but weird situation might occur too
jpayne@1:       if len(m)==1:
jpayne@1:         pass
jpayne@1:       else:
jpayne@1:         pass
jpayne@1:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@1:       if diff<20:
jpayne@1:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list,Final_list_passed)
jpayne@1:         ###need to desgin a algorithm to guess most possible situation for nodes_list, See the situations of test evaluation
jpayne@1:     elif len(a)==1:
jpayne@1:       #that one
jpayne@1:       role,diff=judge_fliC_or_fljB_from_head_tail_for_one_contig(a)
jpayne@1:       if diff<20:
jpayne@1:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list,Final_list_passed)
jpayne@1:       #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@1:     else:#a==0
jpayne@1:       #use Final_list_passed best match
jpayne@1:       for z in Final_list_passed:
jpayne@1:         if x in z[0]:
jpayne@1:           role=z[0].split("_")[0]
jpayne@1:           break
jpayne@1:     #print x,role,len(a)
jpayne@1:     role_list.append((role,x))
jpayne@1:   if len(role_list)==2:
jpayne@1:     if role_list[0][0]==role_list[1][0]:#this is the most cocmmon error, two antigen were assigned to same phase
jpayne@1:       #just use score to do a final test
jpayne@1:       role_list=[]
jpayne@1:       for x in nodes_list: 
jpayne@1:         role=judge_fliC_or_fljB_from_whole_contig_blast_score_ranking(x,Final_list,Final_list_passed)
jpayne@1:         role_list.append((role,x))
jpayne@1:   return role_list
jpayne@1: 
jpayne@1: def decide_contig_roles_for_H_antigen(Final_list,Final_list_passed):
jpayne@1:   #used to decide which contig is FliC and which one is fljB
jpayne@1:   contigs=[]
jpayne@1:   nodes=[]
jpayne@1:   for x in Final_list_passed:
jpayne@1:     if x[0].startswith("fl") and "last" not in x[0] and "first" not in x[0]:
jpayne@1:       nodes.append(x[0].split("___")[1].strip())
jpayne@1:   c,d=Uniq(nodes)#c is node_list
jpayne@1:   #print c
jpayne@1:   tail_head_list=[x for x in Final_list if ("last" in x[0] or "first" in x[0])]
jpayne@1:   roles=fliC_or_fljB_judge_from_head_tail_sequence(c,tail_head_list,Final_list,Final_list_passed)
jpayne@1:   return roles
jpayne@1: 
jpayne@1: def decide_O_type_and_get_special_genes(Final_list,Final_list_passed):
jpayne@1:   #decide O based on Final_list
jpayne@1:   O_choice="?"
jpayne@1:   O_list=[]
jpayne@1:   special_genes={}
jpayne@1:   nodes=[]
jpayne@1:   for x in Final_list_passed:
jpayne@1:     if x[0].startswith("O-"):
jpayne@1:       nodes.append(x[0].split("___")[1].strip())
jpayne@1:     elif not x[0].startswith("fl"):
jpayne@1:       special_genes[x[0]]=x[2]#08172018, x[2] changed from x[-1]
jpayne@1:   #print "special_genes:",special_genes
jpayne@1:   c,d=Uniq(nodes)
jpayne@1:   #print "potential O antigen contig",c
jpayne@1:   final_O=[]
jpayne@1:   O_nodes_list=[]
jpayne@1:   for x in c:#c is the list for contigs
jpayne@1:     temp=0
jpayne@1:     for y in Final_list_passed:
jpayne@1:       if x in y[0] and y[0].startswith("O-"):
jpayne@1:         final_O.append(y)
jpayne@1:         break
jpayne@1:   ### O contig has the problem of two genes on same contig, so do additional test
jpayne@1:   potenial_new_gene=""
jpayne@1:   for x in final_O:
jpayne@1:     pointer=0 #for genes merged or not
jpayne@1:     #not consider O-1,3,19_not_in_3,10, too short compared with others
jpayne@1:     if "O-1,3,19_not_in_3,10" not in x[0] and int(x[0].split("__")[1].split("___")[0])*x[2]+850 <= 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@1:       pointer=x[0].split("___")[1].strip()#store the contig name
jpayne@1:       print(pointer)
jpayne@1:     if pointer!=0:#it has potential merge event
jpayne@1:       for y in Final_list:
jpayne@1:         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@1:           potenial_new_gene=y
jpayne@1:           #print(potenial_new_gene)
jpayne@1:           break
jpayne@1:   if potenial_new_gene!="":
jpayne@1:     print("two differnt genes in same contig, fix it for O antigen")
jpayne@1:     print(potenial_new_gene[:3])
jpayne@4:     pointer=0
jpayne@4:     for y in final_O:
jpayne@4:       if y[0].split("___")[-1]==potenial_new_gene[0].split("___")[-1]:
jpayne@4:         pointer=1
jpayne@4:     if pointer!=1:
jpayne@4:       final_O.append(potenial_new_gene)
jpayne@1:   ### end of the two genes on same contig test
jpayne@1:   final_O=sorted(final_O,key=lambda x: x[2], reverse=True)#sorted
jpayne@1:   if len(final_O)==0 or (len(final_O)==1 and "O-1,3,19_not_in_3,10" in final_O[0][0]):
jpayne@1:     #print "$$$No Otype, due to no hit"#may need to be changed
jpayne@1:     O_choice="-"
jpayne@1:   else:
jpayne@7:     highest_O_coverage=max([float(x[0].split("_cov_")[-1].split("_")[0]) for x in final_O if "O-1,3,19_not_in_3,10" not in x[0]])
jpayne@1:     O_list=[]
jpayne@1:     O_list_less_contamination=[]
jpayne@1:     for x in final_O:
jpayne@1:       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; to avoid contamination affect, use 0.15 of highest coverage as cut-off
jpayne@1:         O_list.append(x[0].split("__")[0])
jpayne@1:         O_nodes_list.append(x[0].split("___")[1])
jpayne@7:         if float(x[0].split("_cov_")[-1].split("_")[0])>highest_O_coverage*0.15:
jpayne@1:           O_list_less_contamination.append(x[0].split("__")[0])
jpayne@1:     ### special test for O9,46 and O3,10 family
jpayne@1:     if ("O-9,46_wbaV" in O_list or "O-9,46_wbaV-from-II-9,12:z29:1,5-SRR1346254" in O_list) and O_list_less_contamination[0].startswith("O-9,"):#not sure should use and float(O9_wbaV)/float(num_1) > 0.1
jpayne@1:       if "O-9,46_wzy" in O_list:#and float(O946_wzy)/float(num_1) > 0.1
jpayne@1:         O_choice="O-9,46"
jpayne@1:         #print "$$$Most possilble Otype:  O-9,46"
jpayne@1:       elif "O-9,46,27_partial_wzy" in O_list:#and float(O94627)/float(num_1) > 0.1
jpayne@1:         O_choice="O-9,46,27"
jpayne@1:         #print "$$$Most possilble Otype:  O-9,46,27"
jpayne@1:       else:
jpayne@1:         O_choice="O-9"#next, detect O9 vs O2?
jpayne@1:         O2=0
jpayne@1:         O9=0
jpayne@1:         for z in special_genes:
jpayne@1:           if "tyr-O-9" in z:
jpayne@1:             O9=special_genes[z]
jpayne@1:           elif "tyr-O-2" in z:
jpayne@1:             O2=special_genes[z]
jpayne@1:         if O2>O9:
jpayne@1:           O_choice="O-2"
jpayne@1:         elif O2<O9:
jpayne@1:           pass
jpayne@1:         else:
jpayne@1:           pass
jpayne@1:           #print "$$$No suitable one, because can't distinct it's O-9 or O-2, but O-9 has a more possibility."
jpayne@1:     elif ("O-3,10_wzx" in O_list) and ("O-9,46_wzy" in O_list) and (O_list[0].startswith("O-3,10") or O_list_less_contamination[0].startswith("O-9,46_wzy")):#and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1
jpayne@1:       if "O-3,10_not_in_1,3,19" in O_list:#and float(O310_no_1319)/float(num_1) > 0.1
jpayne@1:         O_choice="O-3,10"
jpayne@1:         #print "$$$Most possilble Otype:  O-3,10 (contain O-3,10_not_in_1,3,19)"
jpayne@1:       else:
jpayne@1:         O_choice="O-1,3,19"
jpayne@1:         #print "$$$Most possilble Otype:  O-1,3,19 (not contain O-3,10_not_in_1,3,19)"
jpayne@1:     ### end of special test for O9,46 and O3,10 family
jpayne@1:     else:
jpayne@1:       try: 
jpayne@1:         max_score=0
jpayne@1:         for x in final_O:
jpayne@7:           if x[2]>=max_score and float(x[0].split("_cov_")[-1].split("_")[0])>highest_O_coverage*0.15:#use x[2],08172018, the "coverage identity = cover_length * identity"; also meet coverage threshold
jpayne@1:             max_score=x[2]#change from x[-1] to x[2],08172018
jpayne@1:             O_choice=x[0].split("_")[0]
jpayne@1:         if O_choice=="O-1,3,19":
jpayne@1:           O_choice=final_O[1][0].split("_")[0]
jpayne@1:         #print "$$$Most possilble Otype: ",O_choice
jpayne@1:       except:
jpayne@1:         pass
jpayne@1:         #print "$$$No suitable Otype, or failure of mapping (please check the quality of raw reads)"
jpayne@1:   #print "O:",O_choice,O_nodes_list
jpayne@1:   Otypes=[]
jpayne@1:   for x in O_list:
jpayne@1:     if x!="O-1,3,19_not_in_3,10":
jpayne@1:       if "O-9,46_" not in x:
jpayne@1:         Otypes.append(x.split("_")[0])
jpayne@1:       else:
jpayne@1:         Otypes.append(x.split("-from")[0])#O-9,46_wbaV-from-II-9,12:z29:1,5-SRR1346254
jpayne@1:   #Otypes=[x.split("_")[0] for x in O_list if x!="O-1,3,19_not_in_3,10"]
jpayne@1:   Otypes_uniq,Otypes_fre=Uniq(Otypes)
jpayne@1:   contamination_O=""
jpayne@1:   if O_choice=="O-9,46,27" or O_choice=="O-3,10" or O_choice=="O-1,3,19":
jpayne@1:     if len(Otypes_uniq)>2:
jpayne@1:       contamination_O="potential contamination from O antigen signals"
jpayne@1:   else:
jpayne@1:     if len(Otypes_uniq)>1:
jpayne@1:       if O_choice=="O-4" and len(Otypes_uniq)==2 and "O-9,46,27" in Otypes_uniq: #for special 4,12,27 case such as Bredeney and Schwarzengrund
jpayne@1:         contamination_O=""
jpayne@1:       elif O_choice=="O-9,46" and len(Otypes_uniq)==2 and "O-9,46_wbaV" in Otypes_uniq and "O-9,46_wzy" in Otypes_uniq: #for special 4,12,27 case such as Bredeney and Schwarzengrund
jpayne@1:         contamination_O=""
jpayne@1:       else:
jpayne@1:         contamination_O="potential contamination from O antigen signals"
jpayne@4:   return O_choice,O_nodes_list,special_genes,final_O,contamination_O,Otypes_uniq
jpayne@1: ### End of SeqSero2 allele prediction and output
jpayne@1: 
jpayne@1: def get_input_files(make_dir,input_file,data_type,dirpath):
jpayne@1:   #tell input files from datatype
jpayne@1:   #"<int>: '1'(pair-end reads, interleaved),'2'(pair-end reads, seperated),'3'(single-end reads), '4'(assembly),'5'(nanopore fasta),'6'(nanopore fastq)"
jpayne@1:   for_fq=""
jpayne@1:   rev_fq=""
jpayne@1:   old_dir = os.getcwd()
jpayne@1:   os.chdir(make_dir)
jpayne@1:   if data_type=="1":
jpayne@1:     input_file=input_file[0].split("/")[-1]
jpayne@1:     if input_file.endswith(".sra"):
jpayne@1:       subprocess.check_output("fastq-dump --split-files "+input_file,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       for_fq=input_file.replace(".sra","_1.fastq")
jpayne@1:       rev_fq=input_file.replace(".sra","_2.fastq")
jpayne@1:     else:
jpayne@1:       core_id=input_file.split(".fastq")[0].split(".fq")[0]
jpayne@1:       for_fq=core_id+"_1.fastq"
jpayne@1:       rev_fq=core_id+"_2.fastq"
jpayne@1:       if input_file.endswith(".gz"):
jpayne@1:         subprocess.check_output("gzip -dc "+input_file+" | "+dirpath+"/deinterleave_fastq.sh "+for_fq+" "+rev_fq,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       else:
jpayne@1:         subprocess.check_output("cat "+input_file+" | "+dirpath+"/deinterleave_fastq.sh "+for_fq+" "+rev_fq,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   elif data_type=="2":
jpayne@1:     for_fq=input_file[0].split("/")[-1]
jpayne@1:     rev_fq=input_file[1].split("/")[-1]
jpayne@1:   elif data_type=="3":
jpayne@1:     input_file=input_file[0].split("/")[-1]
jpayne@1:     if input_file.endswith(".sra"):
jpayne@1:       subprocess.check_output("fastq-dump --split-files "+input_file,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       for_fq=input_file.replace(".sra","_1.fastq")
jpayne@1:     else:
jpayne@1:       for_fq=input_file
jpayne@1:   elif data_type in ["4","5","6"]:
jpayne@1:     for_fq=input_file[0].split("/")[-1]
jpayne@1:   os.chdir(old_dir)
jpayne@1:   return for_fq,rev_fq
jpayne@1: 
jpayne@4: def predict_O_and_H_types(Final_list,Final_list_passed,new_fasta):
jpayne@1:   #get O and H types from Final_list from blast parsing; allele mode
jpayne@4:   from Bio import SeqIO
jpayne@1:   fliC_choice="-"
jpayne@1:   fljB_choice="-"
jpayne@1:   fliC_contig="NA"
jpayne@1:   fljB_contig="NA"
jpayne@1:   fliC_region=set([0])
jpayne@1:   fljB_region=set([0,])
jpayne@7:   fliC_length=0 #can be changed to coverage in future; in 03292019, changed to ailgned length
jpayne@7:   fljB_length=0 #can be changed to coverage in future; in 03292019, changed to ailgned length
jpayne@1:   O_choice="-"#no need to decide O contig for now, should be only one
jpayne@4:   O_choice,O_nodes,special_gene_list,O_nodes_roles,contamination_O,Otypes_uniq=decide_O_type_and_get_special_genes(Final_list,Final_list_passed)#decide the O antigen type and also return special-gene-list for further identification
jpayne@1:   O_choice=O_choice.split("-")[-1].strip()
jpayne@1:   if (O_choice=="1,3,19" and len(O_nodes_roles)==1 and "1,3,19" in O_nodes_roles[0][0]) or O_choice=="":
jpayne@1:     O_choice="-"
jpayne@1:   H_contig_roles=decide_contig_roles_for_H_antigen(Final_list,Final_list_passed)#decide the H antigen contig is fliC or fljB
jpayne@7:   #add alignment locations, used for further selection, 03312019
jpayne@7:   for i in range(len(H_contig_roles)):
jpayne@7:     x=H_contig_roles[i]
jpayne@7:     for y in Final_list_passed:
jpayne@7:       if x[1] in y[0] and y[0].startswith(x[0]):
jpayne@7:         H_contig_roles[i]+=H_contig_roles[i]+(y[-1],)
jpayne@7:         break
jpayne@1:   log_file=open("SeqSero_log.txt","a")
jpayne@4:   extract_file=open("Extracted_antigen_alleles.fasta","a")
jpayne@4:   handle_fasta=list(SeqIO.parse(new_fasta,"fasta"))
jpayne@4: 
jpayne@4:   #print("O_contigs:")
jpayne@1:   log_file.write("O_contigs:\n")
jpayne@4:   extract_file.write("#Sequences with antigen signals (if the micro-assembled contig only covers the flanking region, it will not be used for contamination analysis)\n")
jpayne@4:   extract_file.write("#O_contigs:\n")
jpayne@1:   for x in O_nodes_roles:
jpayne@1:     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@4:       #print(x[0].split("___")[-1],x[0].split("__")[0],"blast score:",x[1],"identity%:",str(round(x[2]*100,2))+"%",str(min(x[-1]))+" to "+str(max(x[-1])))
jpayne@4:       log_file.write(x[0].split("___")[-1]+" "+x[0].split("__")[0]+"; "+"blast score: "+str(x[1])+" identity%: "+str(round(x[2]*100,2))+"%; alignment from "+str(min(x[-1]))+" to "+str(max(x[-1]))+" of antigen\n")
jpayne@4:       title=">"+x[0].split("___")[-1]+" "+x[0].split("__")[0]+"; "+"blast score: "+str(x[1])+" identity%: "+str(round(x[2]*100,2))+"%; alignment from "+str(min(x[-1]))+" to "+str(max(x[-1]))+" of antigen\n"      
jpayne@4:       seqs=""
jpayne@4:       for z in handle_fasta:
jpayne@4:         if x[0].split("___")[-1]==z.description:
jpayne@4:           seqs=str(z.seq)
jpayne@4:       extract_file.write(title+seqs+"\n")
jpayne@1:   if len(H_contig_roles)!=0:
jpayne@7:     highest_H_coverage=max([float(x[1].split("_cov_")[-1].split("_")[0]) for x in H_contig_roles]) #less than highest*0.1 would be regarded as contamination and noises, they will still be considered in contamination detection and logs, but not used as final serotype output
jpayne@1:   else:
jpayne@1:     highest_H_coverage=0
jpayne@1:   for x in H_contig_roles:
jpayne@1:     #if multiple choices, temporately select the one with longest length for now, will revise in further change
jpayne@7:     if "fliC" == x[0] and len(x[-1])>=fliC_length and x[1] not in O_nodes and float(x[1].split("_cov_")[-1].split("_")[0])>highest_H_coverage*0.13:#remember to avoid the effect of O-type contig, so should not in O_node list
jpayne@1:       fliC_contig=x[1]
jpayne@7:       fliC_length=len(x[-1])
jpayne@7:     elif "fljB" == x[0] and len(x[-1])>=fljB_length and x[1] not in O_nodes and float(x[1].split("_cov_")[-1].split("_")[0])>highest_H_coverage*0.13:
jpayne@1:       fljB_contig=x[1]
jpayne@7:       fljB_length=len(x[-1])
jpayne@1:   for x in Final_list_passed:
jpayne@1:     if fliC_choice=="-" and "fliC_" in x[0] and fliC_contig in x[0]:
jpayne@1:       fliC_choice=x[0].split("_")[1]
jpayne@1:     elif fljB_choice=="-" and "fljB_" in x[0] and fljB_contig in x[0]:
jpayne@1:       fljB_choice=x[0].split("_")[1]
jpayne@1:     elif fliC_choice!="-" and fljB_choice!="-":
jpayne@1:       break
jpayne@1:   #now remove contigs not in middle core part
jpayne@1:   first_allele="NA"
jpayne@1:   first_allele_percentage=0
jpayne@1:   for x in Final_list:
jpayne@1:     if x[0].startswith("fliC") or x[0].startswith("fljB"):
jpayne@1:       first_allele=x[0].split("__")[0] #used to filter those un-middle contigs
jpayne@1:       first_allele_percentage=x[2]
jpayne@1:       break 
jpayne@1:   additional_contigs=[]
jpayne@1:   for x in Final_list:
jpayne@1:     if first_allele in x[0]:
jpayne@1:       if (fliC_contig == x[0].split("___")[-1]): 
jpayne@1:         fliC_region=x[3]
jpayne@1:       elif fljB_contig!="NA" and (fljB_contig == x[0].split("___")[-1]):
jpayne@1:         fljB_region=x[3]
jpayne@1:       else:
jpayne@1:         if x[1]*1.1>int(x[0].split("___")[1].split("_")[3]):#loose threshold by multiplying 1.1
jpayne@1:           additional_contigs.append(x)
jpayne@1:         #else:
jpayne@1:           #print x[:3]
jpayne@1:   #we can just use the fljB region (or fliC depends on size), no matter set() or contain a large locations (without middle part); however, if none of them is fully assembled, use 500 and 1200 as conservative cut-off
jpayne@1:   if first_allele_percentage>0.9:
jpayne@1:     if len(fliC_region)>len(fljB_region) and (max(fljB_region)-min(fljB_region))>1000:
jpayne@1:       target_region=fljB_region|(fliC_region-set(range(min(fljB_region),max(fljB_region)))) #fljB_region|(fliC_region-set(range(min(fljB_region),max(fljB_region))))
jpayne@1:     elif len(fliC_region)<len(fljB_region) and (max(fliC_region)-min(fliC_region))>1000:
jpayne@1:       target_region=fliC_region|(fljB_region-set(range(min(fliC_region),max(fliC_region))))  #fljB_region|(fliC_region-set(range(min(fljB_region),max(fljB_region))))
jpayne@1:     else:
jpayne@1:       target_region=set()#doesn't do anything
jpayne@1:   else:
jpayne@1:     target_region=set()#doesn't do anything
jpayne@1:   #print(target_region)
jpayne@1:   #print(additional_contigs)
jpayne@1:   target_region2=set(list(range(0,525))+list(range(1200,1700)))#I found to use 500 to 1200 as special region would be best
jpayne@1:   target_region=target_region2|target_region
jpayne@1:   for x in additional_contigs:
jpayne@1:     removal=0
jpayne@1:     contig_length=int(x[0].split("___")[1].split("length_")[-1].split("_")[0])
jpayne@1:     if fljB_contig not in x[0] and fliC_contig not in x[0] and len(target_region&x[3])/float(len(x[3]))>0.65 and contig_length*0.5<len(x[3])<contig_length*1.5: #consider length and alignment length for now, but very loose,0.5 and 1.5 as cut-off
jpayne@1:       removal=1
jpayne@1:     else:
jpayne@1:       if first_allele_percentage > 0.9 and float(x[0].split("__")[1].split("___")[0])*x[2]/len(x[-1])>0.96:#if high similiarity with middle part of first allele (first allele >0.9, already cover middle part)
jpayne@1:         removal=1
jpayne@1:       else:
jpayne@1:         pass
jpayne@1:     if removal==1:
jpayne@1:       for y in H_contig_roles:
jpayne@1:         if y[1] in x[0]:
jpayne@1:           H_contig_roles.remove(y)
jpayne@1:     else:
jpayne@1:       pass
jpayne@1:       #print(x[:3],contig_length,len(target_region&x[3])/float(len(x[3])),contig_length*0.5,len(x[3]),contig_length*1.5)
jpayne@1:   #end of removing none-middle contigs
jpayne@4:   #print("H_contigs:")
jpayne@1:   log_file.write("H_contigs:\n")
jpayne@4:   extract_file.write("#H_contigs:\n")
jpayne@1:   H_contig_stat=[]
jpayne@1:   H1_cont_stat={}
jpayne@1:   H2_cont_stat={}
jpayne@1:   for i in range(len(H_contig_roles)):
jpayne@1:     x=H_contig_roles[i]
jpayne@1:     a=0
jpayne@1:     for y in Final_list_passed:
jpayne@1:       if x[1] in y[0] and y[0].startswith(x[0]):
jpayne@1:         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@1:           for y in Final_list_passed: #it's impossible to has the "first" and "last" allele as prediction, so re-do it
jpayne@1:             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@4:               #print(x[1],"can't_decide_fliC_or_fljB",y[0].split("_")[1],"blast_score:",y[1],"identity%:",str(round(y[2]*100,2))+"%",str(min(y[-1]))+" to "+str(max(y[-1])))
jpayne@4:               log_file.write(x[1]+" "+x[0]+" "+y[0].split("_")[1]+"; "+"blast score: "+str(y[1])+" identity%: "+str(round(y[2]*100,2))+"%; alignment from "+str(min(y[-1]))+" to "+str(max(y[-1]))+" of antigen\n")
jpayne@1:               H_contig_roles[i]="can't decide fliC or fljB, may be third phase"
jpayne@4:               title=">"+x[1]+" "+x[0]+" "+y[0].split("_")[1]+"; "+"blast score: "+str(y[1])+" identity%: "+str(round(y[2]*100,2))+"%; alignment from "+str(min(y[-1]))+" to "+str(max(y[-1]))+" of antiten\n"
jpayne@4:               seqs=""
jpayne@4:               for z in handle_fasta:
jpayne@4:                 if x[1]==z.description:
jpayne@4:                   seqs=str(z.seq)
jpayne@4:               extract_file.write(title+seqs+"\n")
jpayne@1:               break
jpayne@1:         else:
jpayne@4:           #print(x[1],x[0],y[0].split("_")[1],"blast_score:",y[1],"identity%:",str(round(y[2]*100,2))+"%",str(min(y[-1]))+" to "+str(max(y[-1])))
jpayne@4:           log_file.write(x[1]+" "+x[0]+" "+y[0].split("_")[1]+"; "+"blast score: "+str(y[1])+" identity%: "+str(round(y[2]*100,2))+"%; alignment from "+str(min(y[-1]))+" to "+str(max(y[-1]))+" of antigen\n")
jpayne@4:           title=">"+x[1]+" "+x[0]+" "+y[0].split("_")[1]+"; "+"blast score: "+str(y[1])+" identity%: "+str(round(y[2]*100,2))+"%; alignment from "+str(min(y[-1]))+" to "+str(max(y[-1]))+" of antigen\n"
jpayne@4:           seqs=""
jpayne@4:           for z in handle_fasta:
jpayne@4:             if x[1]==z.description:
jpayne@4:               seqs=str(z.seq)
jpayne@4:           extract_file.write(title+seqs+"\n")
jpayne@1:         if x[0]=="fliC":
jpayne@1:           if y[0].split("_")[1] not in H1_cont_stat:
jpayne@1:             H1_cont_stat[y[0].split("_")[1]]=y[2]
jpayne@1:           else:
jpayne@1:             H1_cont_stat[y[0].split("_")[1]]+=y[2]
jpayne@1:         if x[0]=="fljB":
jpayne@1:           if y[0].split("_")[1] not in H2_cont_stat:
jpayne@1:             H2_cont_stat[y[0].split("_")[1]]=y[2]
jpayne@1:           else:
jpayne@1:             H2_cont_stat[y[0].split("_")[1]]+=y[2]
jpayne@1:         break
jpayne@1:   #detect contaminations
jpayne@1:   #print(H1_cont_stat)
jpayne@1:   #print(H2_cont_stat)
jpayne@1:   H1_cont_stat_list=[x for x in H1_cont_stat if H1_cont_stat[x]>0.2]
jpayne@1:   H2_cont_stat_list=[x for x in H2_cont_stat if H2_cont_stat[x]>0.2]
jpayne@1:   contamination_H=""
jpayne@1:   if len(H1_cont_stat_list)>1 or len(H2_cont_stat_list)>1:
jpayne@1:     contamination_H="potential contamination from H antigen signals"
jpayne@1:   elif len(H2_cont_stat_list)==1 and fljB_contig=="NA":
jpayne@1:     contamination_H="potential contamination from H antigen signals, uncommon weak fljB signals detected"
jpayne@4:   #get additional antigens
jpayne@4:   """
jpayne@4:     if ("O-9,46_wbaV" in O_list or "O-9,46_wbaV-from-II-9,12:z29:1,5-SRR1346254" in O_list) and O_list_less_contamination[0].startswith("O-9,"):#not sure should use and float(O9_wbaV)/float(num_1) > 0.1
jpayne@4:       if "O-9,46_wzy" in O_list:#and float(O946_wzy)/float(num_1) > 0.1
jpayne@4:         O_choice="O-9,46"
jpayne@4:         #print "$$$Most possilble Otype:  O-9,46"
jpayne@4:       elif "O-9,46,27_partial_wzy" in O_list:#and float(O94627)/float(num_1) > 0.1
jpayne@4:         O_choice="O-9,46,27"
jpayne@4:         #print "$$$Most possilble Otype:  O-9,46,27"
jpayne@4:     elif ("O-3,10_wzx" in O_list) and ("O-9,46_wzy" in O_list) and (O_list[0].startswith("O-3,10") or O_list_less_contamination[0].startswith("O-9,46_wzy")):#and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1
jpayne@4:       if "O-3,10_not_in_1,3,19" in O_list:#and float(O310_no_1319)/float(num_1) > 0.1
jpayne@4:         O_choice="O-3,10"
jpayne@4:         #print "$$$Most possilble Otype:  O-3,10 (contain O-3,10_not_in_1,3,19)"
jpayne@4:       else:
jpayne@4:         O_choice="O-1,3,19"
jpayne@4:         #print "$$$Most possilble Otype:  O-1,3,19 (not contain O-3,10_not_in_1,3,19)"
jpayne@4:     ### end of special test for O9,46 and O3,10 family
jpayne@4: 
jpayne@4:     if O_choice=="O-9,46,27" or O_choice=="O-3,10" or O_choice=="O-1,3,19":
jpayne@4:     if len(Otypes_uniq)>2:
jpayne@4:       contamination_O="potential contamination from O antigen signals"
jpayne@4:   else:
jpayne@4:     if len(Otypes_uniq)>1:
jpayne@4:       if O_choice=="O-4" and len(Otypes_uniq)==2 and "O-9,46,27" in Otypes_uniq: #for special 4,12,27 case such as Bredeney and Schwarzengrund
jpayne@4:         contamination_O=""
jpayne@4:       elif O_choice=="O-9,46" and len(Otypes_uniq)==2 and "O-9,46_wbaV" in Otypes_uniq and "O-9,46_wzy" in Otypes_uniq: #for special 4,12,27 case such as Bredeney and Schwarzengrund
jpayne@4:         contamination_O=""
jpayne@4:   """
jpayne@4:   additonal_antigents=[]
jpayne@4:   #print(contamination_O)
jpayne@4:   #print(contamination_H)
jpayne@1:   log_file.write(contamination_O+"\n")
jpayne@1:   log_file.write(contamination_H+"\n")
jpayne@1:   log_file.close()
jpayne@4:   return O_choice,fliC_choice,fljB_choice,special_gene_list,contamination_O,contamination_H,Otypes_uniq,H1_cont_stat_list,H2_cont_stat_list
jpayne@1: 
jpayne@1: def get_input_K(input_file,lib_dict,data_type,k_size):
jpayne@1:   #kmer mode; get input_Ks from dict and data_type
jpayne@1:   kmers = []
jpayne@1:   for h in lib_dict:
jpayne@1:       kmers += lib_dict[h]
jpayne@1:   if data_type == '4':
jpayne@1:       input_Ks = target_multifasta_kmerizer(input_file, k_size, set(kmers))
jpayne@1:   elif data_type == '1' or data_type == '2' or data_type == '3':#set it for now, will change later
jpayne@1:       input_Ks = target_read_kmerizer(input_file, k_size, set(kmers))
jpayne@1:   elif data_type == '5':#minion_2d_fasta
jpayne@1:       input_Ks = minion_fasta_kmerizer(input_file, k_size, set(kmers))
jpayne@1:   if data_type == '6':#minion_2d_fastq
jpayne@1:       input_Ks = minion_fastq_kmerizer(input_file, k_size, set(kmers))
jpayne@1:   return input_Ks
jpayne@1: 
jpayne@1: def get_kmer_dict(lib_dict,input_Ks):
jpayne@1:   #kmer mode; get predicted types
jpayne@1:   O_dict = {}
jpayne@1:   H_dict = {}
jpayne@1:   Special_dict = {}
jpayne@1:   for h in lib_dict:
jpayne@1:       score = (len(lib_dict[h] & input_Ks) / len(lib_dict[h])) * 100
jpayne@1:       if score > 1:  # Arbitrary cut-off for similarity score very low but seems necessary to detect O-3,10 in some cases
jpayne@1:           if h.startswith('O-') and score > 25:
jpayne@1:               O_dict[h] = score
jpayne@1:           if h.startswith('fl') and score > 40:
jpayne@1:               H_dict[h] = score
jpayne@1:           if (h[:2] != 'fl') and (h[:2] != 'O-'):
jpayne@1:               Special_dict[h] = score
jpayne@1:   return O_dict,H_dict,Special_dict
jpayne@1: 
jpayne@1: def call_O_and_H_type(O_dict,H_dict,Special_dict,make_dir):
jpayne@1:   log_file=open("SeqSero_log.txt","a")
jpayne@1:   log_file.write("O_scores:\n")
jpayne@1:   #call O:
jpayne@1:   highest_O = '-'
jpayne@1:   if len(O_dict) == 0:
jpayne@1:       pass
jpayne@1:   else:
jpayne@1:       for x in O_dict:
jpayne@1:           log_file.write(x+"\t"+str(O_dict[x])+"\n")
jpayne@1:       if ('O-9,46_wbaV__1002' in O_dict and O_dict['O-9,46_wbaV__1002']>70) or ("O-9,46_wbaV-from-II-9,12:z29:1,5-SRR1346254__1002" in O_dict and O_dict['O-9,46_wbaV-from-II-9,12:z29:1,5-SRR1346254__1002']>70):  # not sure should use and float(O9_wbaV)/float(num_1) > 0.1
jpayne@1:           if 'O-9,46_wzy__1191' in O_dict:  # and float(O946_wzy)/float(num_1) > 0.1
jpayne@1:               highest_O = "O-9,46"
jpayne@1:           elif "O-9,46,27_partial_wzy__1019" in O_dict:  # and float(O94627)/float(num_1) > 0.1
jpayne@1:               highest_O = "O-9,46,27"
jpayne@1:           else:
jpayne@1:               highest_O = "O-9"  # next, detect O9 vs O2?
jpayne@1:               O2 = 0
jpayne@1:               O9 = 0
jpayne@1:               for z in Special_dict:
jpayne@1:                   if "tyr-O-9" in z:
jpayne@1:                       O9 = float(Special_dict[z])
jpayne@1:                   if "tyr-O-2" in z:
jpayne@1:                       O2 = float(Special_dict[z])
jpayne@1:               if O2 > O9:
jpayne@1:                   highest_O = "O-2"
jpayne@1:       elif ("O-3,10_wzx__1539" in O_dict) and (
jpayne@1:               "O-9,46_wzy__1191" in O_dict
jpayne@1:       ):  # and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1
jpayne@1:           if "O-3,10_not_in_1,3,19__1519" in O_dict:  # and float(O310_no_1319)/float(num_1) > 0.1
jpayne@1:               highest_O = "O-3,10"
jpayne@1:           else:
jpayne@1:               highest_O = "O-1,3,19"
jpayne@1:       ### end of special test for O9,46 and O3,10 family
jpayne@1:       else:
jpayne@1:           try:
jpayne@1:               max_score = 0
jpayne@1:               for x in O_dict:
jpayne@1:                   if float(O_dict[x]) >= max_score:
jpayne@1:                       max_score = float(O_dict[x])
jpayne@1:                       highest_O = x.split("_")[0]
jpayne@1:               if highest_O == "O-1,3,19":
jpayne@1:                   highest_O = '-'
jpayne@1:                   max_score = 0
jpayne@1:                   for x in O_dict:
jpayne@1:                       if x == 'O-1,3,19_not_in_3,10__130':
jpayne@1:                           pass
jpayne@1:                       else:
jpayne@1:                           if float(O_dict[x]) >= max_score:
jpayne@1:                               max_score = float(O_dict[x])
jpayne@1:                               highest_O = x.split("_")[0]
jpayne@1:           except:
jpayne@1:               pass
jpayne@1:   #call_fliC:
jpayne@1:   if len(H_dict)!=0:
jpayne@1:     highest_H_score_both_BC=H_dict[max(H_dict.keys(), key=(lambda k: H_dict[k]))] #used to detect whether fljB existed or not
jpayne@1:   else:
jpayne@1:     highest_H_score_both_BC=0
jpayne@1:   highest_fliC = '-'
jpayne@1:   highest_fliC_raw = '-'
jpayne@1:   highest_Score = 0
jpayne@1:   log_file.write("\nH_scores:\n")
jpayne@1:   for s in H_dict:
jpayne@1:       log_file.write(s+"\t"+str(H_dict[s])+"\n")
jpayne@1:       if s.startswith('fliC'):
jpayne@1:           if float(H_dict[s]) > highest_Score:
jpayne@1:               highest_fliC = s.split('_')[1]
jpayne@1:               highest_fliC_raw = s
jpayne@1:               highest_Score = float(H_dict[s])
jpayne@1:   #call_fljB
jpayne@1:   highest_fljB = '-'
jpayne@1:   highest_fljB_raw = '-'
jpayne@1:   highest_Score = 0
jpayne@1:   for s in H_dict:
jpayne@1:       if s.startswith('fljB'):
jpayne@1:           if float(H_dict[s]) > highest_Score and float(H_dict[s]) > highest_H_score_both_BC * 0.65: #fljB is special, so use highest_H_score_both_BC to give a general estimate of coverage, currently 0.65 seems pretty good; the reason use a high (0.65) is some fliC and fljB shared with each other
jpayne@1:               highest_fljB = s.split('_')[1]
jpayne@1:               highest_fljB_raw = s
jpayne@1:               highest_Score = float(H_dict[s])
jpayne@1:   log_file.write("\nSpecial_scores:\n")
jpayne@1:   for s in Special_dict:
jpayne@1:     log_file.write(s+"\t"+str(Special_dict[s])+"\n")
jpayne@1:   log_file.close()
jpayne@1:   return highest_O,highest_fliC,highest_fljB
jpayne@1: 
jpayne@1: def get_temp_file_names(for_fq,rev_fq):
jpayne@1:   #seqsero2 -a; get temp file names
jpayne@1:   sam=for_fq+".sam"
jpayne@1:   bam=for_fq+".bam"
jpayne@1:   sorted_bam=for_fq+"_sorted.bam"
jpayne@1:   mapped_fq1=for_fq+"_mapped.fq"
jpayne@1:   mapped_fq2=rev_fq+"_mapped.fq"
jpayne@1:   combined_fq=for_fq+"_combined.fq"
jpayne@1:   for_sai=for_fq+".sai"
jpayne@1:   rev_sai=rev_fq+".sai"
jpayne@1:   return sam,bam,sorted_bam,mapped_fq1,mapped_fq2,combined_fq,for_sai,rev_sai
jpayne@1: 
jpayne@1: def map_and_sort(threads,database,fnameA,fnameB,sam,bam,for_sai,rev_sai,sorted_bam,mapping_mode):
jpayne@1:   #seqsero2 -a; do mapping and sort
jpayne@1:   print("building database...")
jpayne@1:   subprocess.check_output("bwa index "+database,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   print("mapping...")
jpayne@1:   if mapping_mode=="mem":
jpayne@1:     subprocess.check_output("bwa mem -k 17 -t "+threads+" "+database+" "+fnameA+" "+fnameB+" > "+sam,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   elif mapping_mode=="sam":
jpayne@1:     if fnameB!="":
jpayne@1:       subprocess.check_output("bwa aln -t "+threads+" "+database+" "+fnameA+" > "+for_sai,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       subprocess.check_output("bwa aln -t "+threads+" "+database+" "+fnameB+" > "+rev_sai,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       subprocess.check_output("bwa sampe "+database+" "+for_sai+" "+ rev_sai+" "+fnameA+" "+fnameB+" > "+sam,shell=True, stderr=subprocess.STDOUT)
jpayne@1:     else:
jpayne@1:       subprocess.check_output("bwa aln -t "+threads+" "+database+" "+fnameA+" > "+for_sai,shell=True, stderr=subprocess.STDOUT)
jpayne@1:       subprocess.check_output("bwa samse "+database+" "+for_sai+" "+for_fq+" > "+sam, stderr=subprocess.STDOUT)
jpayne@1:   subprocess.check_output("samtools view -@ "+threads+" -F 4 -Sh "+sam+" > "+bam,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   ### check the version of samtools then use differnt commands
jpayne@1:   samtools_version=subprocess.Popen(["samtools"],stdout=subprocess.PIPE,stderr=subprocess.PIPE)
jpayne@1:   out, err = samtools_version.communicate()
jpayne@1:   version = str(err).split("ersion:")[1].strip().split(" ")[0].strip()
jpayne@1:   print("check samtools version:",version)
jpayne@1:   ### end of samtools version check and its analysis
jpayne@1:   if LooseVersion(version)<=LooseVersion("1.2"):
jpayne@1:     subprocess.check_output("samtools sort -@ "+threads+" -n "+bam+" "+fnameA+"_sorted",shell=True, stderr=subprocess.STDOUT)
jpayne@1:   else:
jpayne@1:     subprocess.check_output("samtools sort -@ "+threads+" -n "+bam+" >"+sorted_bam,shell=True, stderr=subprocess.STDOUT)
jpayne@1: 
jpayne@1: def extract_mapped_reads_and_do_assembly_and_blast(current_time,sorted_bam,combined_fq,mapped_fq1,mapped_fq2,threads,fnameA,fnameB,database,mapping_mode):
jpayne@1:   #seqsero2 -a; extract, assembly and blast
jpayne@1:   subprocess.check_output("bamToFastq -i "+sorted_bam+" -fq "+combined_fq,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   if fnameB!="":
jpayne@1:     subprocess.check_output("bamToFastq -i "+sorted_bam+" -fq "+mapped_fq1+" -fq2 "+mapped_fq2 ,shell=True, stderr=subprocess.STDOUT)#2> /dev/null if want no output
jpayne@1:   else:
jpayne@1:     pass
jpayne@1:   outdir=current_time+"_temp"
jpayne@1:   print("assembling...")
jpayne@1:   if int(threads)>4:
jpayne@1:     t="4"
jpayne@1:   else:
jpayne@1:     t=threads
jpayne@12:   new_fasta = None
jpayne@7:   if os.path.getsize(combined_fq)>100 and (fnameB=="" or os.path.getsize(mapped_fq1)>100):#if not, then it's "-:-:-"
jpayne@1:     if fnameB!="":
jpayne@1:       subprocess.check_output("spades.py --careful --pe1-s "+combined_fq+" --pe1-1 "+mapped_fq1+" --pe1-2 "+mapped_fq2+" -t "+t+" -o "+outdir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:     else:
jpayne@1:       subprocess.check_output("spades.py --careful --pe1-s "+combined_fq+" -t "+t+" -o "+outdir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:     new_fasta=fnameA+"_"+database+"_"+mapping_mode+".fasta"
jpayne@1:     subprocess.check_output("mv "+outdir+"/contigs.fasta "+new_fasta+ " 2> /dev/null",shell=True, stderr=subprocess.STDOUT)
jpayne@1:     #os.system("mv "+outdir+"/scaffolds.fasta "+new_fasta+ " 2> /dev/null") contigs.fasta
jpayne@1:     subprocess.check_output("rm -rf "+outdir+ " 2> /dev/null",shell=True, stderr=subprocess.STDOUT)
jpayne@1:     print("blasting...","\n")
jpayne@1:     xmlfile="blasted_output.xml"#fnameA+"-extracted_vs_"+database+"_"+mapping_mode+".xml"
jpayne@1:     subprocess.check_output('makeblastdb -in '+new_fasta+' -out '+new_fasta+'_db '+'-dbtype nucl',shell=True, stderr=subprocess.STDOUT) #temp.txt is to forbid the blast result interrupt the output of our program###1/27/2015
jpayne@1:     subprocess.check_output("blastn -query "+database+" -db "+new_fasta+"_db -out "+xmlfile+" -outfmt 5",shell=True, stderr=subprocess.STDOUT)###1/27/2015; 08272018, remove "-word_size 10"
jpayne@1:   else:
jpayne@1:     xmlfile="NA"
jpayne@4:   return xmlfile,new_fasta
jpayne@1: 
jpayne@1: def judge_subspecies(fnameA):
jpayne@1:   #seqsero2 -a; judge subspecies on just forward raw reads fastq
jpayne@12:   salmID_output=subprocess.check_output(os.path.join(dirpath, "SalmID/SalmID.py") + ' -i ' +fnameA, shell=True, stderr=subprocess.STDOUT)
jpayne@1:   #out, err = salmID_output.communicate()
jpayne@1:   #out=out.decode("utf-8")
jpayne@1:   out = salmID_output.decode("utf-8")
jpayne@1:   file=open("data_log.txt","a")
jpayne@1:   file.write(out)
jpayne@1:   file.close()
jpayne@1:   salm_species_scores=out.split("\n")[1].split("\t")[6:]
jpayne@1:   salm_species_results=out.split("\n")[0].split("\t")[6:]
jpayne@1:   max_score=0
jpayne@1:   max_score_index=1 #default is 1, means "I"
jpayne@1:   for i in range(len(salm_species_scores)):
jpayne@1:     if max_score<float(salm_species_scores[i]):
jpayne@1:       max_score=float(salm_species_scores[i])
jpayne@1:       max_score_index=i
jpayne@1:   prediction=salm_species_results[max_score_index].split(".")[1].strip().split(" ")[0]
jpayne@1:   if float(out.split("\n")[1].split("\t")[4]) > float(out.split("\n")[1].split("\t")[5]): #bongori and enterica compare
jpayne@1:     prediction="bongori" #if not, the prediction would always be enterica, since they are located in the later part
jpayne@1:   if max_score<10:
jpayne@1:     prediction="-"
jpayne@1:   return prediction
jpayne@1: 
jpayne@1: def judge_subspecies_Kmer(Special_dict):
jpayne@1:   #seqsero2 -k;
jpayne@1:   max_score=0
jpayne@1:   prediction="-" #default should be I
jpayne@1:   for x in Special_dict:
jpayne@1:     if "mer" in x:
jpayne@1:       if max_score<float(Special_dict[x]):
jpayne@1:         max_score=float(Special_dict[x])
jpayne@1:         prediction=x.split("_")[-1].strip()
jpayne@1:       if x.split("_")[-1].strip()=="bongori" and float(Special_dict[x])>95:#if bongori already, then no need to test enterica
jpayne@1:         prediction="bongori"
jpayne@1:         break
jpayne@1:   return prediction
jpayne@1: 
jpayne@1: def main():
jpayne@1:   #combine SeqSeroK and SeqSero2, also with SalmID
jpayne@1:   args = parse_args()
jpayne@1:   input_file = args.i
jpayne@1:   data_type = args.t
jpayne@1:   analysis_mode = args.m
jpayne@1:   mapping_mode=args.b
jpayne@1:   threads=args.p
jpayne@1:   make_dir=args.d
jpayne@1:   clean_mode=args.c
jpayne@1:   k_size=27 #will change for bug fixing
jpayne@1:   database="H_and_O_and_specific_genes.fasta"
jpayne@1:   
jpayne@1:   if len(sys.argv)==1:
jpayne@1:     subprocess.check_output(dirpath+"/SeqSero2_package.py -h",shell=True, stderr=subprocess.STDOUT)#change name of python file
jpayne@1:   else:
jpayne@1:     request_id = time.strftime("%m_%d_%Y_%H_%M_%S", time.localtime())
jpayne@1:     request_id += str(random.randint(1, 10000000))
jpayne@1:     if make_dir is None:
jpayne@1:       make_dir="SeqSero_result_"+request_id
jpayne@1:     if os.path.isdir(make_dir):
jpayne@1:       pass
jpayne@1:     else:
jpayne@1:       subprocess.check_output(["mkdir",make_dir])
jpayne@1:     #subprocess.check_output("cp "+dirpath+"/"+database+" "+" ".join(input_file)+" "+make_dir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:     subprocess.check_output("ln -srf "+dirpath+"/"+database+" "+" ".join(input_file)+" "+make_dir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:   ############################begin the real analysis 
jpayne@1:     if analysis_mode=="a":
jpayne@1:       if data_type in ["1","2","3"]:#use allele mode
jpayne@12:         H1_cont_stat_list = H2_cont_stat_list = []
jpayne@1:         for_fq,rev_fq=get_input_files(make_dir,input_file,data_type,dirpath)
jpayne@1:         os.chdir(make_dir)
jpayne@1:         ###add a function to tell input files
jpayne@1:         fnameA=for_fq.split("/")[-1]
jpayne@1:         fnameB=rev_fq.split("/")[-1]
jpayne@1:         current_time=time.strftime("%Y_%m_%d_%H_%M_%S", time.localtime())
jpayne@1:         sam,bam,sorted_bam,mapped_fq1,mapped_fq2,combined_fq,for_sai,rev_sai=get_temp_file_names(fnameA,fnameB) #get temp files id
jpayne@1:         map_and_sort(threads,database,fnameA,fnameB,sam,bam,for_sai,rev_sai,sorted_bam,mapping_mode) #do mapping and sort
jpayne@4:         xmlfile,new_fasta=extract_mapped_reads_and_do_assembly_and_blast(current_time,sorted_bam,combined_fq,mapped_fq1,mapped_fq2,threads,fnameA,fnameB,database,mapping_mode) #extract the mapped reads and do micro assembly and blast
jpayne@1:         if xmlfile=="NA":
jpayne@1:           O_choice,fliC_choice,fljB_choice,special_gene_list,contamination_O,contamination_H=("-","-","-",[],"","")
jpayne@1:         else:
jpayne@1:           Final_list=xml_parse_score_comparision_seqsero(xmlfile) #analyze xml and get parsed results
jpayne@1:           file=open("data_log.txt","a")
jpayne@1:           for x in Final_list:
jpayne@1:             file.write("\t".join(str(y) for y in x)+"\n")
jpayne@1:           file.close()
jpayne@7:           Final_list_passed=[x for x in Final_list if float(x[0].split("_cov_")[1].split("_")[0])>=0.9 and (x[1]>=int(x[0].split("__")[1]) or x[1]>=int(x[0].split("___")[1].split("_")[3]) or x[1]>1000)]
jpayne@4:           O_choice,fliC_choice,fljB_choice,special_gene_list,contamination_O,contamination_H,Otypes_uniq,H1_cont_stat_list,H2_cont_stat_list=predict_O_and_H_types(Final_list,Final_list_passed,new_fasta) #predict O, fliC and fljB
jpayne@1:         subspecies=judge_subspecies(fnameA) #predict subspecies
jpayne@1:         ###output
jpayne@1:         predict_form,predict_sero,star,star_line,claim=seqsero_from_formula_to_serotypes(O_choice,fliC_choice,fljB_choice,special_gene_list,subspecies)
jpayne@7:         claim="" #04132019, disable claim for new report requirement
jpayne@1:         contamination_report=""
jpayne@4:         H_list=["fliC_"+x for x in H1_cont_stat_list if len(x)>0]+["fljB_"+x for x in H2_cont_stat_list if len(x)>0]
jpayne@1:         if contamination_O!="" and contamination_H=="":
jpayne@4:           contamination_report="#Potential inter-serotype contamination detected from O antigen signals. All O-antigens detected:"+"\t".join(Otypes_uniq)+"."
jpayne@1:         elif contamination_O=="" and contamination_H!="":
jpayne@4:           contamination_report="#Potential inter-serotype contamination detected or potential thrid H phase from H antigen signals. All H-antigens detected:"+"\t".join(H_list)+"."
jpayne@1:         elif contamination_O!="" and contamination_H!="":
jpayne@4:           contamination_report="#Potential inter-serotype contamination detected from both O and H antigen signals.All O-antigens detected:"+"\t".join(Otypes_uniq)+". All H-antigens detected:"+"\t".join(H_list)+"."
jpayne@4:         if contamination_report!="":
jpayne@7:           #contamination_report="potential inter-serotype contamination detected (please refer below antigen signal report for details)." #above contamination_reports are for back-up and bug fixing #web-based mode need to be re-used, 04132019
jpayne@7:           contamination_report="Co-existence of multiple serotypes detected, indicating potential inter-serotype contamination. See 'Extracted_antigen_alleles.fasta' for detected serotype determinant alleles."
jpayne@7:         #claim="\n"+open("Extracted_antigen_alleles.fasta","r").read()#used to store H and O antigen sequeences #04132019, need to change if using web-version
jpayne@7:         if contamination_report+star_line+claim=="": #0413, new output style
jpayne@7:           note=""
jpayne@7:         else:
jpayne@7:           note="Note:"
jpayne@1:         if clean_mode:
jpayne@1:           subprocess.check_output("rm -rf ../"+make_dir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:           make_dir="none-output-directory due to '-c' flag"
jpayne@1:         else:
jpayne@1:           #new_file=open("Seqsero_result.txt","w")
jpayne@1:           if O_choice=="":
jpayne@1:             O_choice="-"
jpayne@1:           result = OrderedDict(
jpayne@1:           sample_name = input_file[0],
jpayne@1:           O_antigen_prediction = O_choice,
jpayne@1:           H1_antigen_prediction = fliC_choice,
jpayne@1:           H2_antigen_prediction = fljB_choice,
jpayne@1:           predicted_antigenic_profile = predict_form,
jpayne@1:           predicted_subspecies = subspecies,
jpayne@7:           predicted_serotype = "{}{}".format(predict_sero, star),
jpayne@1:           note=claim.replace('\n','')
jpayne@1:           )
jpayne@9:           result['*'] = star_line.replace('\t', ' ').replace('\n','')
jpayne@1:           with open("Seqsero_result.tsv","w") as new_file:
jpayne@1:          #new_file.write("Output_directory:"+make_dir+"\nInput files:\t"+input_file+"\n"+"O antigen prediction:\t"+O_choice+"\n"+"H1 antigen prediction(fliC):\t"+highest_fliC+"\n"+"H2 antigen prediction(fljB):\t"+highest_fljB+"\n"+"Predicted antigenic profile:\t"+predict_form+"\n"+"Predicted subspecies:\t"+subspecies+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+star+star_line+claim+"\n")#+##
jpayne@1:          #new_file.close()
jpayne@1:             wrtr = DictWriter(new_file, delimiter='\t', fieldnames=result.keys())
jpayne@1:             wrtr.writeheader()
jpayne@1:             wrtr.writerow(result)
jpayne@1: 
jpayne@1:           
jpayne@1:           # new_file.write("Output_directory:"+make_dir+"\nInput files:\t"+for_fq+" "+rev_fq+"\n"+"O antigen prediction:\t"+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 subspecies:\t"+subspecies+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+contamination_report+star+star_line+claim+"\n")#+##
jpayne@1:           # new_file.close()
jpayne@1:           #subprocess.check_output("cat Seqsero_result.txt",shell=True, stderr=subprocess.STDOUT)
jpayne@1:           #subprocess.call("rm H_and_O_and_specific_genes.fasta* *.sra *.bam *.sam *.fastq *.gz *.fq temp.txt *.xml "+fnameA+"*_db* 2> /dev/null",shell=True, stderr=subprocess.STDOUT)
jpayne@1:           subprocess.call("rm H_and_O_and_specific_genes.fasta* *.sra *.bam *.sam *.fastq *.gz *.fq temp.txt "+fnameA+"*_db* 2> /dev/null",shell=True, stderr=subprocess.STDOUT)
jpayne@1:         print("Output_directory:"+make_dir+"\nInput files:\t"+for_fq+" "+rev_fq+"\n"+"O antigen prediction:\t"+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 subspecies:\t"+subspecies+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+contamination_report+star+star_line+claim+"\n")#+##
jpayne@1:       else:
jpayne@1:         print("Allele modes only support raw reads datatype, i.e. '-t 1 or 2 or 3'; please use '-m k'")
jpayne@1:     elif analysis_mode=="k":
jpayne@1:       ex_dir = os.path.dirname(os.path.realpath(__file__))
jpayne@1:       #output_mode = args.mode
jpayne@1:       for_fq,rev_fq=get_input_files(make_dir,input_file,data_type,dirpath)
jpayne@1:       input_file = for_fq #-k will just use forward because not all reads were used
jpayne@1:       os.chdir(make_dir)
jpayne@1:       f = open(ex_dir + '/antigens.pickle', 'rb')
jpayne@1:       lib_dict = pickle.load(f)
jpayne@1:       f.close
jpayne@1:       input_Ks = get_input_K(input_file,lib_dict,data_type,k_size)
jpayne@1:       O_dict,H_dict,Special_dict = get_kmer_dict(lib_dict,input_Ks)
jpayne@1:       highest_O,highest_fliC,highest_fljB = call_O_and_H_type(O_dict,H_dict,Special_dict,make_dir)
jpayne@1:       subspecies = judge_subspecies_Kmer(Special_dict)
jpayne@1:       if subspecies=="IIb" or subspecies=="IIa":
jpayne@1:         subspecies="II"
jpayne@1:       predict_form,predict_sero,star,star_line,claim = seqsero_from_formula_to_serotypes(
jpayne@1:           highest_O.split('-')[1], highest_fliC, highest_fljB, Special_dict,subspecies)
jpayne@7:       claim="" #no claim any more based on new output requirement
jpayne@7:       if star_line+claim=="": #0413, new output style
jpayne@7:         note=""
jpayne@7:       else:
jpayne@7:         note="Note:"
jpayne@1:       if clean_mode:
jpayne@1:         subprocess.check_output("rm -rf ../"+make_dir,shell=True, stderr=subprocess.STDOUT)
jpayne@1:         make_dir="none-output-directory due to '-c' flag"
jpayne@1:       else:
jpayne@1:         if highest_O.split('-')[-1]=="":
jpayne@1:           O_choice="-"
jpayne@1:         else:
jpayne@1:           O_choice=highest_O.split('-')[-1]
jpayne@1:         #print("Output_directory:"+make_dir+"\tInput_file:"+input_file+"\tPredicted subpecies:"+subspecies + '\tPredicted antigenic profile:' + predict_form + '\tPredicted serotype(s):' + predict_sero)
jpayne@1:         result = OrderedDict(
jpayne@1:           sample_name = input_file,
jpayne@1:           O_antigen_prediction = O_choice,
jpayne@1:           H1_antigen_prediction = highest_fliC,
jpayne@1:           H2_antigen_prediction = highest_fljB,
jpayne@1:           predicted_antigenic_profile = predict_form,
jpayne@1:           predicted_subspecies = subspecies,
jpayne@7:           predicted_serotype = "{}{}".format(predict_sero, star),
jpayne@1:           note=claim.replace('\n','')
jpayne@1:         )
jpayne@9:         result['*'] = star_line.replace('\t', ' ').replace('\n','')
jpayne@1:         with open("Seqsero_result.tsv","w") as new_file:
jpayne@1:          #new_file.write("Output_directory:"+make_dir+"\nInput files:\t"+input_file+"\n"+"O antigen prediction:\t"+O_choice+"\n"+"H1 antigen prediction(fliC):\t"+highest_fliC+"\n"+"H2 antigen prediction(fljB):\t"+highest_fljB+"\n"+"Predicted antigenic profile:\t"+predict_form+"\n"+"Predicted subspecies:\t"+subspecies+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+star+star_line+claim+"\n")#+##
jpayne@1:          #new_file.close()
jpayne@1:          wrtr = DictWriter(new_file, delimiter='\t', fieldnames=result.keys())
jpayne@1:          wrtr.writeheader()
jpayne@1:          wrtr.writerow(result)
jpayne@1:         subprocess.call("rm *.fasta* *.fastq *.gz *.fq temp.txt *.sra 2> /dev/null",shell=True, stderr=subprocess.STDOUT)
jpayne@1:         print("Output_directory:"+make_dir+"\nInput files:\t"+input_file+"\n"+"O antigen prediction:\t"+O_choice+"\n"+"H1 antigen prediction(fliC):\t"+highest_fliC+"\n"+"H2 antigen prediction(fljB):\t"+highest_fljB+"\n"+"Predicted antigenic profile:\t"+predict_form+"\n"+"Predicted subspecies:\t"+subspecies+"\n"+"Predicted serotype(s):\t"+predict_sero+star+"\n"+star+star_line+claim+"\n")#+##
jpayne@1:         return 0
jpayne@1: 
jpayne@1: if __name__ == '__main__':
jpayne@1:   try:
jpayne@1:     quit(main())
jpayne@1:   except subprocess.CalledProcessError as e:
jpayne@1:     print(str(e.output))
jpayne@1:     print(e)
jpayne@1:     quit(e.returncode)