Mercurial > repos > jpayne > seqsero_v2

annotate SeqSero2/SalmID/SalmID.py @ 14:496f5f2b5e75

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author estrain
date Tue, 14 Feb 2023 12:11:27 -0500
parents 4ac593d4b40f
children
rev   line source
jpayne@12 1 #!/usr/bin/env python3
jpayne@12 2
jpayne@12 3
jpayne@12 4 import gzip
jpayne@12 5 import io
jpayne@12 6 import pickle
jpayne@12 7 import os
jpayne@12 8 import sys
jpayne@12 9
jpayne@12 10 from argparse import ArgumentParser
jpayne@12 11 try:
jpayne@12 12 from version import SalmID_version
jpayne@12 13 except:
jpayne@12 14 SalmID_version = "version unknown"
jpayne@12 15
jpayne@12 16
jpayne@12 17 def reverse_complement(sequence):
jpayne@12 18 complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N': 'N', 'M': 'K', 'R': 'Y', 'W': 'W',
jpayne@12 19 'S': 'S', 'Y': 'R', 'K': 'M', 'V': 'B', 'H': 'D', 'D': 'H', 'B': 'V'}
jpayne@12 20 return "".join(complement[base] for base in reversed(sequence))
jpayne@12 21
jpayne@12 22
jpayne@12 23 def parse_args():
jpayne@12 24 "Parse the input arguments, use '-h' for help."
jpayne@12 25 parser = ArgumentParser(description='SalmID - rapid Kmer based Salmonella identifier from sequence data')
jpayne@12 26 # inputs
jpayne@12 27 parser.add_argument('-v','--version', action='version', version='%(prog)s ' + SalmID_version)
jpayne@12 28 parser.add_argument(
jpayne@12 29 '-i','--input_file', type=str, required=False, default= 'None', metavar = 'your_fastqgz',
jpayne@12 30 help='Single fastq.gz file input, include path to file if file is not in same directory ')
jpayne@12 31 parser.add_argument(
jpayne@12 32 '-e', '--extension', type=str, required=False, default= '.fastq.gz', metavar = 'file_extension',
jpayne@12 33 help='File extension, if specified without "--input_dir", SalmID will attempt to ID all files\n' +
jpayne@12 34 ' with this extension in current directory, otherwise files in input directory')
jpayne@12 35
jpayne@12 36 parser.add_argument(
jpayne@12 37 '-d','--input_dir', type=str, required=False, default='.', metavar = 'directory',
jpayne@12 38 help='Directory which contains data for identification, when not specified files in current directory will be analyzed.')
jpayne@12 39 parser.add_argument(
jpayne@12 40 '-r', '--report', type=str, required=False, default='percentage', metavar = 'percentage, coverage or taxonomy',
jpayne@12 41 help='Report either percentage ("percentage") of clade specific kmers recovered, average kmer-coverage ("cov"), or '
jpayne@12 42 'taxonomy (taxonomic species ID, plus observed mean k-mer coverages and expected coverage).')
jpayne@12 43 parser.add_argument(
jpayne@12 44 '-m', '--mode', type=str, required=False, default='quick', metavar = 'quick or thorough',
jpayne@12 45 help='Quick [quick] or thorough [thorough] mode')
jpayne@12 46 if len(sys.argv)==1:
jpayne@12 47 parser.print_help(sys.stderr)
jpayne@12 48 sys.exit(1)
jpayne@12 49 return parser.parse_args()
jpayne@12 50
jpayne@12 51 def get_av_read_length(file):
jpayne@12 52 i = 1
jpayne@12 53 n_reads = 0
jpayne@12 54 total_length = 0
jpayne@12 55 if file.endswith(".gz"):
jpayne@12 56 file_content=io.BufferedReader(gzip.open(file))
jpayne@12 57 else:
jpayne@12 58 file_content=open(file,"r").readlines()
jpayne@12 59 for line in file_content:
jpayne@12 60 if i % 4 == 2:
jpayne@12 61 total_length += len(line.strip())
jpayne@12 62 n_reads +=1
jpayne@12 63 i += 1
jpayne@12 64 if n_reads == 100:
jpayne@12 65 break
jpayne@12 66 return total_length/100
jpayne@12 67
jpayne@12 68
jpayne@12 69 def createKmerDict_reads(list_of_strings, kmer):
jpayne@12 70 kmer_table = {}
jpayne@12 71 for string in list_of_strings:
jpayne@12 72 sequence = string.strip('\n')
jpayne@12 73 for i in range(len(sequence)-kmer+1):
jpayne@12 74 new_mer =sequence[i:i+kmer]
jpayne@12 75 new_mer_rc = reverse_complement(new_mer)
jpayne@12 76 if new_mer in kmer_table:
jpayne@12 77 kmer_table[new_mer.upper()] += 1
jpayne@12 78 else:
jpayne@12 79 kmer_table[new_mer.upper()] = 1
jpayne@12 80 if new_mer_rc in kmer_table:
jpayne@12 81 kmer_table[new_mer_rc.upper()] += 1
jpayne@12 82 else:
jpayne@12 83 kmer_table[new_mer_rc.upper()] = 1
jpayne@12 84 return kmer_table
jpayne@12 85
jpayne@12 86
jpayne@12 87 def target_read_kmerizer_multi(file, k, kmerDict_1, kmerDict_2, mode):
jpayne@12 88 mean_1 = None
jpayne@12 89 mean_2 = None
jpayne@12 90 i = 1
jpayne@12 91 n_reads_1 = 0
jpayne@12 92 n_reads_2 = 0
jpayne@12 93 total_coverage_1 = 0
jpayne@12 94 total_coverage_2 = 0
jpayne@12 95 reads_1 = []
jpayne@12 96 reads_2 = []
jpayne@12 97 total_reads = 0
jpayne@12 98 if file.endswith(".gz"):
jpayne@12 99 file_content=io.BufferedReader(gzip.open(file))
jpayne@12 100 else:
jpayne@12 101 file_content=open(file,"r").readlines()
jpayne@12 102 for line in file_content:
jpayne@12 103 start = int((len(line) - k) // 2)
jpayne@12 104 if i % 4 == 2:
jpayne@12 105 total_reads += 1
jpayne@12 106 if file.endswith(".gz"):
jpayne@12 107 s1 = line[start:k + start].decode()
jpayne@12 108 line=line.decode()
jpayne@12 109 else:
jpayne@12 110 s1 = line[start:k + start]
jpayne@12 111 if s1 in kmerDict_1:
jpayne@12 112 n_reads_1 += 1
jpayne@12 113 total_coverage_1 += len(line)
jpayne@12 114 reads_1.append(line)
jpayne@12 115 if s1 in kmerDict_2:
jpayne@12 116 n_reads_2 += 1
jpayne@12 117 total_coverage_2 += len(line)
jpayne@12 118 reads_2.append(line)
jpayne@12 119 i += 1
jpayne@12 120 if mode == 'quick':
jpayne@12 121 if total_coverage_2 >= 800000:
jpayne@12 122 break
jpayne@12 123
jpayne@12 124 if len(reads_1) == 0:
jpayne@12 125 kmer_Dict1 = {}
jpayne@12 126 else:
jpayne@12 127 kmer_Dict1 = createKmerDict_reads(reads_1, k)
jpayne@12 128 mers_1 = set([key for key in kmer_Dict1])
jpayne@12 129 mean_1 = sum([kmer_Dict1[key] for key in kmer_Dict1])/len(mers_1)
jpayne@12 130 if len(reads_2) == 0:
jpayne@12 131 kmer_Dict2 = {}
jpayne@12 132 else:
jpayne@12 133 kmer_Dict2 = createKmerDict_reads(reads_2, k)
jpayne@12 134 mers_2 = set([key for key in kmer_Dict2])
jpayne@12 135 mean_2 = sum([kmer_Dict2[key] for key in kmer_Dict2])/len(mers_2)
jpayne@12 136 return kmer_Dict1, kmer_Dict2, mean_1, mean_2, total_reads
jpayne@12 137
jpayne@12 138 def mean_cov_selected_kmers(iterable, kmer_dict, clade_specific_kmers):
jpayne@12 139 '''
jpayne@12 140 Given an iterable (list, set, dictrionary) returns mean coverage for the kmers in iterable
jpayne@12 141 :param iterable: set, list or dictionary containing kmers
jpayne@12 142 :param kmer_dict: dictionary with kmers as keys, kmer-frequency as value
jpayne@12 143 :param clade_specific_kmers: list, dict or set of clade specific kmers
jpayne@12 144 :return: mean frequency as float
jpayne@12 145 '''
jpayne@12 146 if len(iterable) == 0:
jpayne@12 147 return 0
jpayne@12 148 return sum([kmer_dict[value] for value in iterable])/len(clade_specific_kmers)
jpayne@12 149
jpayne@12 150 def kmer_lists(query_fastq_gz, k,
jpayne@12 151 allmers,allmers_rpoB,
jpayne@12 152 uniqmers_bongori,
jpayne@12 153 uniqmers_I,
jpayne@12 154 uniqmers_IIa,
jpayne@12 155 uniqmers_IIb,
jpayne@12 156 uniqmers_IIIa,
jpayne@12 157 uniqmers_IIIb,
jpayne@12 158 uniqmers_IV,
jpayne@12 159 uniqmers_VI,
jpayne@12 160 uniqmers_VII,
jpayne@12 161 uniqmers_VIII,
jpayne@12 162 uniqmers_bongori_rpoB,
jpayne@12 163 uniqmers_S_enterica_rpoB,
jpayne@12 164 uniqmers_Escherichia_rpoB,
jpayne@12 165 uniqmers_Listeria_ss_rpoB,
jpayne@12 166 uniqmers_Lmono_rpoB,
jpayne@12 167 mode):
jpayne@12 168 dict_invA, dict_rpoB, mean_invA, mean_rpoB , total_reads = target_read_kmerizer_multi(query_fastq_gz, k, allmers,
jpayne@12 169 allmers_rpoB, mode)
jpayne@12 170 target_mers_invA = set([key for key in dict_invA])
jpayne@12 171 target_mers_rpoB = set([key for key in dict_rpoB])
jpayne@12 172 if target_mers_invA == 0:
jpayne@12 173 print('No reads found matching invA, no Salmonella in sample?')
jpayne@12 174 else:
jpayne@12 175 p_bongori = (len(uniqmers_bongori & target_mers_invA) / len(uniqmers_bongori)) * 100
jpayne@12 176 p_I = (len(uniqmers_I & target_mers_invA) / len(uniqmers_I)) * 100
jpayne@12 177 p_IIa = (len(uniqmers_IIa & target_mers_invA) / len(uniqmers_IIa)) * 100
jpayne@12 178 p_IIb = (len(uniqmers_IIb & target_mers_invA) / len(uniqmers_IIb)) * 100
jpayne@12 179 p_IIIa = (len(uniqmers_IIIa & target_mers_invA) / len(uniqmers_IIIa)) * 100
jpayne@12 180 p_IIIb = (len(uniqmers_IIIb & target_mers_invA) / len(uniqmers_IIIb)) * 100
jpayne@12 181 p_VI = (len(uniqmers_VI & target_mers_invA) / len(uniqmers_VI)) * 100
jpayne@12 182 p_IV = (len(uniqmers_IV & target_mers_invA) / len(uniqmers_IV)) * 100
jpayne@12 183 p_VII = (len(uniqmers_VII & target_mers_invA) / len(uniqmers_VII)) * 100
jpayne@12 184 p_VIII = (len(uniqmers_VIII & target_mers_invA) / len(uniqmers_VIII)) * 100
jpayne@12 185 p_bongori_rpoB = (len(uniqmers_bongori_rpoB & target_mers_rpoB) / len(uniqmers_bongori_rpoB)) * 100
jpayne@12 186 p_Senterica = (len(uniqmers_S_enterica_rpoB & target_mers_rpoB) / len(uniqmers_S_enterica_rpoB)) * 100
jpayne@12 187 p_Escherichia = (len(uniqmers_Escherichia_rpoB & target_mers_rpoB) / len(uniqmers_Escherichia_rpoB)) * 100
jpayne@12 188 p_Listeria_ss = (len(uniqmers_Listeria_ss_rpoB & target_mers_rpoB) / len(uniqmers_Listeria_ss_rpoB)) * 100
jpayne@12 189 p_Lmono = (len(uniqmers_Lmono_rpoB & target_mers_rpoB) / len(uniqmers_Lmono_rpoB)) * 100
jpayne@12 190 bongori_invA_cov = mean_cov_selected_kmers(uniqmers_bongori & target_mers_invA, dict_invA, uniqmers_bongori)
jpayne@12 191 I_invA_cov = mean_cov_selected_kmers(uniqmers_I & target_mers_invA, dict_invA, uniqmers_I)
jpayne@12 192 IIa_invA_cov = mean_cov_selected_kmers(uniqmers_IIa & target_mers_invA, dict_invA, uniqmers_IIa)
jpayne@12 193 IIb_invA_cov = mean_cov_selected_kmers(uniqmers_IIb & target_mers_invA, dict_invA, uniqmers_IIb)
jpayne@12 194 IIIa_invA_cov = mean_cov_selected_kmers(uniqmers_IIIa & target_mers_invA, dict_invA, uniqmers_IIIa)
jpayne@12 195 IIIb_invA_cov = mean_cov_selected_kmers(uniqmers_IIIb & target_mers_invA, dict_invA, uniqmers_IIIb)
jpayne@12 196 IV_invA_cov = mean_cov_selected_kmers(uniqmers_IV & target_mers_invA, dict_invA, uniqmers_IV)
jpayne@12 197 VI_invA_cov = mean_cov_selected_kmers(uniqmers_VI & target_mers_invA, dict_invA, uniqmers_VI)
jpayne@12 198 VII_invA_cov = mean_cov_selected_kmers(uniqmers_VII & target_mers_invA, dict_invA, uniqmers_VII)
jpayne@12 199 VIII_invA_cov = mean_cov_selected_kmers(uniqmers_VIII & target_mers_invA, dict_invA, uniqmers_VIII)
jpayne@12 200 S_enterica_rpoB_cov = mean_cov_selected_kmers((uniqmers_S_enterica_rpoB & target_mers_rpoB), dict_rpoB,
jpayne@12 201 uniqmers_S_enterica_rpoB)
jpayne@12 202 S_bongori_rpoB_cov = mean_cov_selected_kmers((uniqmers_bongori_rpoB & target_mers_rpoB), dict_rpoB,
jpayne@12 203 uniqmers_bongori_rpoB)
jpayne@12 204 Escherichia_rpoB_cov = mean_cov_selected_kmers((uniqmers_Escherichia_rpoB & target_mers_rpoB), dict_rpoB,
jpayne@12 205 uniqmers_Escherichia_rpoB)
jpayne@12 206 Listeria_ss_rpoB_cov = mean_cov_selected_kmers((uniqmers_Listeria_ss_rpoB & target_mers_rpoB), dict_rpoB,
jpayne@12 207 uniqmers_Listeria_ss_rpoB)
jpayne@12 208 Lmono_rpoB_cov = mean_cov_selected_kmers((uniqmers_Lmono_rpoB & target_mers_rpoB), dict_rpoB,
jpayne@12 209 uniqmers_Lmono_rpoB)
jpayne@12 210 coverages = [Listeria_ss_rpoB_cov, Lmono_rpoB_cov, Escherichia_rpoB_cov, S_bongori_rpoB_cov,
jpayne@12 211 S_enterica_rpoB_cov, bongori_invA_cov, I_invA_cov, IIa_invA_cov, IIb_invA_cov,
jpayne@12 212 IIIa_invA_cov, IIIb_invA_cov, IV_invA_cov, VI_invA_cov, VII_invA_cov, VIII_invA_cov]
jpayne@12 213 locus_scores = [p_Listeria_ss, p_Lmono, p_Escherichia, p_bongori_rpoB, p_Senterica, p_bongori,
jpayne@12 214 p_I, p_IIa,p_IIb, p_IIIa, p_IIIb, p_IV, p_VI, p_VII, p_VIII]
jpayne@12 215 return locus_scores, coverages, total_reads
jpayne@12 216
jpayne@12 217 def report_taxon(locus_covs, average_read_length, number_of_reads):
jpayne@12 218 list_taxa = [ 'Listeria ss', 'Listeria monocytogenes', 'Escherichia sp.',
jpayne@12 219 'Salmonella bongori (rpoB)', 'Salmonella enterica (rpoB)',
jpayne@12 220 'Salmonella bongori (invA)', 'S. enterica subsp. enterica (invA)',
jpayne@12 221 'S. enterica subsp. salamae (invA: clade a)','S. enterica subsp. salamae (invA: clade b)',
jpayne@12 222 'S. enterica subsp. arizonae (invA)', 'S. enterica subsp. diarizonae (invA)',
jpayne@12 223 'S. enterica subsp. houtenae (invA)', 'S. enterica subsp. indica (invA)',
jpayne@12 224 'S. enterica subsp. VII (invA)', 'S. enterica subsp. salamae (invA: clade VIII)']
jpayne@12 225 if sum(locus_covs) < 1:
jpayne@12 226 rpoB = ('No rpoB matches!', 0)
jpayne@12 227 invA = ('No invA matches!', 0)
jpayne@12 228 return rpoB, invA, 0.0
jpayne@12 229 else:
jpayne@12 230 # given list of scores get taxon
jpayne@12 231 if sum(locus_covs[0:5]) > 0:
jpayne@12 232 best_rpoB = max(range(len(locus_covs[1:5])), key=lambda x: locus_covs[1:5][x])+1
jpayne@12 233 all_rpoB = max(range(len(locus_covs[0:5])), key=lambda x: locus_covs[0:5][x])
jpayne@12 234 if (locus_covs[best_rpoB] != 0) & (all_rpoB == 0):
jpayne@12 235 rpoB = (list_taxa[best_rpoB], locus_covs[best_rpoB])
jpayne@12 236 elif (all_rpoB == 0) & (round(sum(locus_covs[1:5]),1) < 1):
jpayne@12 237 rpoB = (list_taxa[0], locus_covs[0])
jpayne@12 238 else:
jpayne@12 239 rpoB = (list_taxa[best_rpoB], locus_covs[best_rpoB])
jpayne@12 240 else:
jpayne@12 241 rpoB = ('No rpoB matches!', 0)
jpayne@12 242 if sum(locus_covs[5:]) > 0:
jpayne@12 243 best_invA = max(range(len(locus_covs[5:])), key=lambda x: locus_covs[5:][x])+5
jpayne@12 244 invA = (list_taxa[best_invA], locus_covs[best_invA])
jpayne@12 245 else:
jpayne@12 246 invA = ('No invA matches!', 0)
jpayne@12 247 if 'Listeria' in rpoB[0]:
jpayne@12 248 return rpoB, invA, (average_read_length * number_of_reads) / 3000000
jpayne@12 249 else:
jpayne@12 250 return rpoB, invA, (average_read_length * number_of_reads) / 5000000
jpayne@12 251
jpayne@12 252
jpayne@12 253
jpayne@12 254 def main():
jpayne@12 255 ex_dir = os.path.dirname(os.path.realpath(__file__))
jpayne@12 256 args = parse_args()
jpayne@12 257 input_file = args.input_file
jpayne@12 258 if input_file != 'None':
jpayne@12 259 files = [input_file]
jpayne@12 260 else:
jpayne@12 261 extension = args.extension
jpayne@12 262 inputdir = args.input_dir
jpayne@12 263 files = [inputdir + '/'+ f for f in os.listdir(inputdir) if f.endswith(extension)]
jpayne@12 264 report = args.report
jpayne@12 265 mode = args.mode
jpayne@12 266 f_invA = open(ex_dir + "/invA_mers_dict", "rb")
jpayne@12 267 sets_dict_invA = pickle.load(f_invA)
jpayne@12 268 f_invA.close()
jpayne@12 269 allmers = sets_dict_invA['allmers']
jpayne@12 270 uniqmers_I = sets_dict_invA['uniqmers_I']
jpayne@12 271 uniqmers_IIa = sets_dict_invA['uniqmers_IIa']
jpayne@12 272 uniqmers_IIb = sets_dict_invA['uniqmers_IIb']
jpayne@12 273 uniqmers_IIIa = sets_dict_invA['uniqmers_IIIa']
jpayne@12 274 uniqmers_IIIb = sets_dict_invA['uniqmers_IIIb']
jpayne@12 275 uniqmers_IV = sets_dict_invA['uniqmers_IV']
jpayne@12 276 uniqmers_VI = sets_dict_invA['uniqmers_VI']
jpayne@12 277 uniqmers_VII = sets_dict_invA['uniqmers_VII']
jpayne@12 278 uniqmers_VIII = sets_dict_invA['uniqmers_VIII']
jpayne@12 279 uniqmers_bongori = sets_dict_invA['uniqmers_bongori']
jpayne@12 280
jpayne@12 281 f = open(ex_dir + "/rpoB_mers_dict", "rb")
jpayne@12 282 sets_dict = pickle.load(f)
jpayne@12 283 f.close()
jpayne@12 284
jpayne@12 285 allmers_rpoB = sets_dict['allmers']
jpayne@12 286 uniqmers_bongori_rpoB = sets_dict['uniqmers_bongori']
jpayne@12 287 uniqmers_S_enterica_rpoB = sets_dict['uniqmers_S_enterica']
jpayne@12 288 uniqmers_Escherichia_rpoB = sets_dict['uniqmers_Escherichia']
jpayne@12 289 uniqmers_Listeria_ss_rpoB = sets_dict['uniqmers_Listeria_ss']
jpayne@12 290 uniqmers_Lmono_rpoB = sets_dict['uniqmers_L_mono']
jpayne@12 291 #todo: run kmer_lists() once, create list of tuples containing data to be used fro different reports
jpayne@12 292 if report == 'taxonomy':
jpayne@12 293 print('file\trpoB\tinvA\texpected coverage')
jpayne@12 294 for f in files:
jpayne@12 295 locus_scores, coverages, reads = kmer_lists(f, 27,
jpayne@12 296 allmers, allmers_rpoB,
jpayne@12 297 uniqmers_bongori,
jpayne@12 298 uniqmers_I,
jpayne@12 299 uniqmers_IIa,
jpayne@12 300 uniqmers_IIb,
jpayne@12 301 uniqmers_IIIa,
jpayne@12 302 uniqmers_IIIb,
jpayne@12 303 uniqmers_IV,
jpayne@12 304 uniqmers_VI,
jpayne@12 305 uniqmers_VII,
jpayne@12 306 uniqmers_VIII,
jpayne@12 307 uniqmers_bongori_rpoB,
jpayne@12 308 uniqmers_S_enterica_rpoB,
jpayne@12 309 uniqmers_Escherichia_rpoB,
jpayne@12 310 uniqmers_Listeria_ss_rpoB,
jpayne@12 311 uniqmers_Lmono_rpoB,
jpayne@12 312 mode)
jpayne@12 313 pretty_covs = [round(cov, 1) for cov in coverages]
jpayne@12 314 report = report_taxon(pretty_covs, get_av_read_length(f), reads)
jpayne@12 315 print(f.split('/')[-1] + '\t' + report[0][0] + '[' + str(report[0][1]) + ']' + '\t' + report[1][0] +
jpayne@12 316 '[' + str(report[1][1]) + ']' +
jpayne@12 317 '\t' + str(round(report[2], 1)))
jpayne@12 318 else:
jpayne@12 319 print(
jpayne@12 320 'file\tListeria sensu stricto (rpoB)\tL. monocytogenes (rpoB)\tEscherichia spp. (rpoB)\tS. bongori (rpoB)\tS. enterica' +
jpayne@12 321 '(rpoB)\tS. bongori (invA)\tsubsp. I (invA)\tsubsp. II (clade a: invA)\tsubsp. II' +
jpayne@12 322 ' (clade b: invA)\tsubsp. IIIa (invA)\tsubsp. IIIb (invA)\tsubsp.IV (invA)\tsubsp. VI (invA)\tsubsp. VII (invA)' +
jpayne@12 323 '\tsubsp. II (clade VIII : invA)')
jpayne@12 324 if report == 'percentage':
jpayne@12 325 for f in files:
jpayne@12 326 locus_scores, coverages , reads = kmer_lists( f, 27,
jpayne@12 327 allmers,allmers_rpoB,
jpayne@12 328 uniqmers_bongori,
jpayne@12 329 uniqmers_I,
jpayne@12 330 uniqmers_IIa,
jpayne@12 331 uniqmers_IIb,
jpayne@12 332 uniqmers_IIIa,
jpayne@12 333 uniqmers_IIIb,
jpayne@12 334 uniqmers_IV,
jpayne@12 335 uniqmers_VI,
jpayne@12 336 uniqmers_VII,
jpayne@12 337 uniqmers_VIII,
jpayne@12 338 uniqmers_bongori_rpoB,
jpayne@12 339 uniqmers_S_enterica_rpoB,
jpayne@12 340 uniqmers_Escherichia_rpoB,
jpayne@12 341 uniqmers_Listeria_ss_rpoB,
jpayne@12 342 uniqmers_Lmono_rpoB,
jpayne@12 343 mode)
jpayne@12 344 pretty_scores = [str(round(score)) for score in locus_scores]
jpayne@12 345 print(f.split('/')[-1] +'\t' + '\t'.join(pretty_scores))
jpayne@12 346 else:
jpayne@12 347 for f in files:
jpayne@12 348 locus_scores, coverages , reads = kmer_lists( f, 27,
jpayne@12 349 allmers,allmers_rpoB,
jpayne@12 350 uniqmers_bongori,
jpayne@12 351 uniqmers_I,
jpayne@12 352 uniqmers_IIa,
jpayne@12 353 uniqmers_IIb,
jpayne@12 354 uniqmers_IIIa,
jpayne@12 355 uniqmers_IIIb,
jpayne@12 356 uniqmers_IV,
jpayne@12 357 uniqmers_VI,
jpayne@12 358 uniqmers_VII,
jpayne@12 359 uniqmers_VIII,
jpayne@12 360 uniqmers_bongori_rpoB,
jpayne@12 361 uniqmers_S_enterica_rpoB,
jpayne@12 362 uniqmers_Escherichia_rpoB,
jpayne@12 363 uniqmers_Listeria_ss_rpoB,
jpayne@12 364 uniqmers_Lmono_rpoB,
jpayne@12 365 mode)
jpayne@12 366 pretty_covs = [str(round(cov, 1)) for cov in coverages]
jpayne@12 367 print(f.split('/')[-1] + '\t' + '\t'.join(pretty_covs))
jpayne@12 368
jpayne@12 369 if __name__ == '__main__':
jpayne@12 370 main()
jpayne@12 371