Mercurial > repos > estrain > seqsero2_v1_0_1

annotate SalmID.py @ 16:81fce62ea8c1

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