Mercurial > repos > degregory > sam_refiner

comparison SAM_Refiner_Gal.py @ 0:22e3f843f1db draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded main file
author degregory
date Mon, 13 Sep 2021 14:12:21 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:22e3f843f1db
1 #!/bin/env python3
2 # Writen by Devon Gregory with assistance from Christopher Bottoms
3 # University of Missouri
4 # Distributed under GNU GENERAL PUBLIC LICENSE v3
5 import os
6 import sys
7 import argparse
8 import itertools
9 import time
10 from multiprocessing import Process
11 from pathlib import Path
12
13 # process for collecing command line arguments
14 def arg_parser():
15
16 parser = argparse.ArgumentParser(
17 description='process Sam files for variant information'
18 )
19
20 parser.add_argument(
21 '--samp',
22 type=str,
23 default='Sample',
24 help='Sample Name'
25 )
26 parser.add_argument(
27 '-r', '-reference',
28 type=argparse.FileType('r'),
29 dest='ref',
30 help='reference fasta'
31 )
32 parser.add_argument(
33 '-S', '--Sam_file',
34 type=argparse.FileType('r'),
35 dest='Sam_file',
36 help='.sam file'
37 )
38 parser.add_argument(
39 '-o1', '--out_file1',
40 dest='outfile1',
41 help='output file'
42 )
43 parser.add_argument(
44 '-o2', '--out_file2',
45 dest='outfile2',
46 help='output file'
47 )
48 parser.add_argument(
49 '-o3', '--out_file3',
50 dest='outfile3',
51 help='output file'
52 )
53 parser.add_argument(
54 '-o4', '--out_file4',
55 dest='outfile4',
56 help='output file'
57 )
58 parser.add_argument(
59 '-o5', '--out_file5',
60 dest='outfile5',
61 help='output file'
62 )
63 parser.add_argument(
64 '-o6', '--out_file6',
65 dest='outfile6',
66 help='output file'
67 )
68 parser.add_argument(
69 '-o7', '--out_file7',
70 dest='outfile7',
71 help='output file'
72 )
73 parser.add_argument(
74 '--use_count',
75 type=int,
76 default=1,
77 choices=[0, 1],
78 help='Enable/Disable (1/0) use of counts in sequence IDs, default enabled (--use_count 1)'
79 )
80 parser.add_argument(
81 '--min_count',
82 type=int,
83 default=10,
84 help='Minimum observations required to be included in sample reports; >= 1 occurance count; < 1 %% observed (.1 = 10%%), (default: .001)'
85 )
86 parser.add_argument(
87 '--min_samp_abund',
88 type=float,
89 default=0.001,
90 help='Minimum abundance required for inclusion in sample reports; %% observed (.1 = 10%%), (default: .001)'
91 )
92 parser.add_argument(
93 '--ntabund',
94 type=float,
95 default=.001,
96 help='Minimum abundance relative to total reads required for a position to be reported in the nt call output; must be non-negative and < 1, %% observed (.1 = 10%%), (default: .001)'
97 )
98 parser.add_argument(
99 '--max_dist',
100 type=int,
101 default=40,
102 help='Maximum number of variances from the reference a sequence can have to be consider in covars processing (default: 40)'
103 )
104 parser.add_argument(
105 '--max_covar',
106 type=int,
107 default=8,
108 help='Maximum number of variances from the reference to be reported in covars (default: 8)'
109 )
110 parser.add_argument(
111 '--AAreport',
112 type=int,
113 default=1,
114 choices=[0, 1],
115 help='Enable/Disable (1/0) amino acid reporting, default enabled (--AAreport 1)'
116 )
117 parser.add_argument(
118 '--AAcodonasMNP',
119 type=int,
120 default=1,
121 choices=[0, 1],
122 help='Enable/Disable (1/0) reporting multiple nt changes in a single codon as one polymorphism, default enabled (--AAcodonasMNP 1), requires AAreport enabled'
123 )
124 parser.add_argument(
125 '--alpha',
126 type=float,
127 default=1.2,
128 help='Modifier for chim_rm chimera checking, default 1.2. Higher = more sensitive, more false chimeras removed; lower = less sensitive, fewer chimeras removed'
129 )
130 parser.add_argument(
131 '--foldab',
132 type=float,
133 default=1.8,
134 help='Threshold for potential parent / chimera abundance ratio for chim_rm; default is 1.8'
135 )
136 parser.add_argument(
137 '--redist',
138 type=int,
139 default=1,
140 choices=[0, 1],
141 help='Enable/Disable (1/0) redistribution of chimera counts for chim_rm, default enabled (--redist 1)'
142 )
143 parser.add_argument(
144 '--beta',
145 type=float,
146 default=1,
147 help='Modifier for covar pass checking, default 1. Higher = more sensitive, more failed checks; lower = less sensitive, fewer failed checks'
148 )
149 parser.add_argument(
150 '--autopass',
151 type=float,
152 default=.3,
153 help='threshold for a sequence to automatically pass the covar pass checking'
154 )
155 parser.add_argument(
156 '--nt_call',
157 type=int,
158 default=1,
159 choices=[0, 1],
160 help='Enable/Disable (1/0) nt_call output, default enabled (--nt_call 1)'
161 )
162 parser.add_argument(
163 '--ntvar',
164 type=int,
165 default=0,
166 choices=[0, 1],
167 help='Enable/Disable (1/0) nt_var output, default enabled (--nt_call 1)'
168 )
169 parser.add_argument(
170 '--indel',
171 type=int,
172 default=1,
173 choices=[0, 1],
174 help='Enable/Disable (1/0) indel output, default enabled (--indel 1)'
175 )
176 parser.add_argument(
177 '--seq',
178 type=int,
179 default=1,
180 choices=[0, 1],
181 help='Enable/Disable (1/0) unique seq output, default enabled (--seq 1)'
182 )
183 parser.add_argument(
184 '--covar', type=int,
185 default=1,
186 choices=[0, 1],
187 help='Enable/Disable (1/0) covar output, default enabled (--covar 1)'
188 )
189 parser.add_argument(
190 '--chim_rm',
191 type=int,
192 default=1,
193 choices=[0, 1],
194 help='Enable/Disable (1/0) chim_rm output, default enabled (--chim_rm 1)'
195 )
196 parser.add_argument(
197 '--deconv',
198 type=int,
199 default=1,
200 choices=[0, 1],
201 help='Enable/Disable (1/0) covar deconv, default enabled (--deconv 1)'
202 )
203 parser.add_argument(
204 '--wgs',
205 type=int,
206 default=0,
207 choices=[0, 1],
208 help='Enable/Disable (1/0) covar deconv, default enabled (--deconv 1)'
209 )
210
211 args = parser.parse_args()
212
213 # checking for proper range of some parameters and consistency/compatibility
214 if args.wgs == 1:
215 if args.deconv == 1 or args.chim_rm == 1:
216 args.deconv = 0
217 args.chim_rm = 0
218 print('WGS mode enabled, disabling chimera removal methods')
219
220 return(args)
221
222 def get_ref(ref): # get the reference ID and sequence from the FASTA file. Will only get the first.
223
224 n=0
225 refID = ''
226 refseq = ''
227 if ref:
228 for line in ref:
229 if line.startswith('>'):
230 n+=1
231 if n > 1:
232 break
233 refID = line[1:].strip("\n\r")
234 elif n == 1:
235 refseq = refseq + line.strip("\n\r")
236 refseq = refseq.upper()
237
238
239 return(refID, refseq)
240
241 def AAcall(codon): # amino acid / codon dictionary to return encoded AAs
242 AAdict = {
243 'TTT' : 'F',
244 'TTC' : 'F',
245 'TTA' : 'L',
246 'TTG' : 'L',
247 'TCT' : 'S',
248 'TCC' : 'S',
249 'TCA' : 'S',
250 'TCG' : 'S',
251 'TAT' : 'Y',
252 'TAC' : 'Y',
253 'TAA' : '*',
254 'TAG' : '*',
255 'TGT' : 'C',
256 'TGC' : 'C',
257 'TGA' : '*',
258 'TGG' : 'W',
259
260 'CTT' : 'L',
261 'CTC' : 'L',
262 'CTA' : 'L',
263 'CTG' : 'L',
264 'CCT' : 'P',
265 'CCC' : 'P',
266 'CCA' : 'P',
267 'CCG' : 'P',
268 'CAT' : 'H',
269 'CAC' : 'H',
270 'CAA' : 'Q',
271 'CAG' : 'Q',
272 'CGT' : 'R',
273 'CGC' : 'R',
274 'CGA' : 'R',
275 'CGG' : 'R',
276
277 'ATT' : 'I',
278 'ATC' : 'I',
279 'ATA' : 'I',
280 'ATG' : 'M',
281 'ACT' : 'T',
282 'ACC' : 'T',
283 'ACA' : 'T',
284 'ACG' : 'T',
285 'AAT' : 'N',
286 'AAC' : 'N',
287 'AAA' : 'K',
288 'AAG' : 'K',
289 'AGT' : 'S',
290 'AGC' : 'S',
291 'AGA' : 'R',
292 'AGG' : 'R',
293
294 'GTT' : 'V',
295 'GTC' : 'V',
296 'GTA' : 'V',
297 'GTG' : 'V',
298 'GCT' : 'A',
299 'GCC' : 'A',
300 'GCA' : 'A',
301 'GCG' : 'A',
302 'GAT' : 'D',
303 'GAC' : 'D',
304 'GAA' : 'E',
305 'GAG' : 'E',
306 'GGT' : 'G',
307 'GGC' : 'G',
308 'GGA' : 'G',
309 'GGG' : 'G'
310 }
311 AA = '?'
312 if codon in AAdict:
313 AA = AAdict[codon]
314
315 return(AA)
316
317 def singletCodon(ntPOS, nt, ref): # process to return the AA and protein seq. position based on the reference and provided nt seq position and nt
318 AAPOS = (ntPOS-1)//3
319 AAmod = (ntPOS-1)%3
320 codon = ""
321 try:
322 if AAmod == 0:
323 codon = nt+ref[1][ntPOS]+ref[1][ntPOS+1]
324 elif AAmod == 1:
325 codon = ref[1][ntPOS-2]+nt+ref[1][ntPOS]
326 elif AAmod == 2:
327 codon = ref[1][ntPOS-3]+ref[1][ntPOS-2]+nt
328 except:
329 codon = "XXX"
330
331 return(AAPOS+1, AAcall(codon))
332
333 def getCombos(qlist, clen): # returns combinations of single polymorphisms in a sequence
334 combos = []
335 if (clen == 0 or clen > len(qlist)):
336 clen = len(qlist)
337 for N in range(1, clen+1):
338 for comb in itertools.combinations(qlist, N):
339 combos.append(' '.join(comb))
340 return(combos)
341
342 def SAMparse(args, ref, refprot, file): # process SAM files
343 covar_array = []
344 seq_array = []
345
346 nt_call_dict_dict = {}
347 indel_dict = {}
348 seq_species = {}
349 sam_read_count = 0
350 sam_line_count = 0
351 firstPOS = 0
352 lastPOS = 0
353 coverage = {}
354 for line in args.Sam_file:
355 if not line.startswith('@'): # ignore header lines
356 splitline = line.split("\t")
357 if ref[0].upper().startswith(splitline[2].upper()): # check map ID matches referecne ID
358 if int(splitline[4]) > 0: # Check mapping score is positive
359
360 abund_count=1
361 if args.use_count == 1: # get the unique sequence counts
362 if '-' in splitline[0] and '=' in splitline[0]:
363 eq_split = splitline[0].split('=')
364 dash_split = splitline[0].split('-')
365 if len(eq_split[-1]) > len(dash_split[-1]):
366 abund_count=int(dash_split[-1])
367 else:
368 abund_count=int(eq_split[-1])
369
370 elif '-' in splitline[0]:
371 try:
372 abund_count=int(splitline[0].split('-')[-1])
373 except:
374 pass
375
376 elif '=' in splitline[0]:
377 try:
378 abund_count=int(splitline[0].split('=')[-1])
379 except:
380 pass
381
382
383 sam_read_count += abund_count
384
385 CIGAR = splitline[5]
386 POS = int(splitline[3])
387 if firstPOS == 0:
388 firstPOS = POS
389 elif POS < firstPOS:
390 firstPOS = POS
391
392 readID = splitline[0]
393 query_seq = splitline[9].upper()
394 run_length = 0
395 query_seq_parsed = ''
396 query_pos = 0
397 q_pars_pos = 0
398 mutations = []
399
400 for C in CIGAR: # process sequence based on standard CIGAR line
401 if C == 'M' or C == 'I' or C == 'D' or C == 'S' or C == 'H':
402 if C == 'S':
403 query_pos = query_pos + run_length
404
405 if C == 'I':
406 if query_pos > 0:
407 # add insertion to dict
408 iPOS = q_pars_pos+POS
409
410 iSeq = query_seq[query_pos: query_pos+run_length]
411 istring = str(iPOS)+'-insert'+iSeq
412
413 try:
414 indel_dict[istring]
415 except:
416 indel_dict[istring] = abund_count
417 else:
418 indel_dict[istring] += abund_count
419
420 if args.AAreport == 1 and (run_length % 3 == 0):
421
422 iProt = ''
423 if iPOS % 3 == 1:
424 for x in range(0, (run_length//3)):
425 AA = AAcall(iSeq[x*3]+iSeq[x*3+1]+iSeq[x*3+2])
426 iProt = iProt + AA
427 mutations.append(istring + '(' + str((iPOS//3)+1) + iProt + ')')
428 elif iPOS % 3 == 2:
429 if query_pos > 0:
430 ipSeq = query_seq[query_pos-1:query_pos+run_length+2]
431 else:
432 ipSeq = "XXX"+query_seq[query_pos+2:query_pos+run_length+2]
433 for x in range(0, (run_length//3)+1):
434 AA = AAcall(ipSeq[x*3]+ipSeq[x*3+1]+ipSeq[x*3+2])
435 iProt = iProt + AA
436 mutations.append(istring + '(' + str((iPOS//3)+1) + iProt + ')')
437 else:
438 if query_pos > 1:
439 ipSeq = query_seq[query_pos-2:query_pos+run_length+1]
440 else:
441 ipSeq = "XXX"+query_seq[query_pos+1:query_pos+run_length+1]
442
443 for x in range(0, (run_length//3)+1):
444 AA = AAcall(ipSeq[x*3]+ipSeq[x*3+1]+ipSeq[x*3+2])
445 iProt = iProt + AA
446 mutations.append(istring + '(' + str((iPOS//3)+1) + iProt + ')')
447 else:
448 mutations.append(istring)
449
450 query_pos = query_pos + run_length
451
452 elif C == 'D':
453 for X in range(0, run_length):
454 query_seq_parsed += '-'
455
456 delstring = str(q_pars_pos+POS)+'-'+str(q_pars_pos+POS+run_length-1)+'Del'
457
458 if args.AAreport == 1 and (run_length % 3 == 0) and not ((q_pars_pos+POS) % 3 == 1 ):
459 if (q_pars_pos+POS) % 3 == 2:
460 newcodon = query_seq[query_pos-1:query_pos+2]
461 newAArefpos = (q_pars_pos+POS) // 3
462 mutations.append(delstring + '(' + refprot[newAArefpos] + str(newAArefpos+1) + AAcall(newcodon) + ')')
463 else:
464 newcodon = query_seq[query_pos-2:query_pos+1]
465 newAArefpos = (q_pars_pos+POS) // 3
466 mutations.append(delstring + '(' + refprot[newAArefpos] + str(newAArefpos+1) + AAcall(newcodon) + ')')
467 else:
468 mutations.append(delstring)
469
470 if args.nt_call == 1:
471 for N in range(q_pars_pos+POS, q_pars_pos+POS+int(run_length)):
472 try:
473 nt_call_dict_dict[N]
474 except:
475 nt_call_dict_dict[N] = {'A' : 0,
476 'T' : 0,
477 'C' : 0,
478 'G' : 0,
479 '-' : 0}
480 nt_call_dict_dict[N]['-'] = abund_count
481 else:
482 nt_call_dict_dict[N]['-'] += abund_count
483
484 try:
485 indel_dict[str(q_pars_pos+POS)+'-'+str(q_pars_pos+POS+run_length-1)+'Del']
486 except:
487 indel_dict[str(q_pars_pos+POS)+'-'+str(q_pars_pos+POS+run_length-1)+'Del'] = int(abund_count)
488 else:
489 indel_dict[str(q_pars_pos+POS)+'-'+str(q_pars_pos+POS+run_length-1)+'Del'] += int(abund_count)
490
491 q_pars_pos = q_pars_pos + run_length
492
493 elif C == 'M':
494 offset = q_pars_pos-query_pos
495 refPOS = POS+offset
496
497 for ntPOS in range(query_pos, query_pos+run_length):
498 if query_seq[ntPOS] == 'A' or query_seq[ntPOS] == 'T' or query_seq[ntPOS] == 'C' or query_seq[ntPOS] == 'G':
499 if query_seq[ntPOS] != ref[1][refPOS+ntPOS-1]:
500 if args.AAreport == 1 and args.AAcodonasMNP == 0:
501 AAinfo = singletCodon(refPOS+ntPOS, query_seq[ntPOS], ref)
502 mutations.append(str(refPOS+ntPOS)+query_seq[ntPOS]+'('+refprot[AAinfo[0]-1]+str(AAinfo[0])+AAinfo[1]+')')
503 else:
504 mutations.append(str(refPOS+ntPOS)+query_seq[ntPOS])
505 if args.nt_call == 1:
506 try:
507 nt_call_dict_dict[refPOS+ntPOS]
508 except:
509 nt_call_dict_dict[refPOS+ntPOS] = {'A' : 0,
510 'T' : 0,
511 'C' : 0,
512 'G' : 0,
513 '-' : 0}
514 nt_call_dict_dict[refPOS+ntPOS][query_seq[ntPOS]] = abund_count
515 else:
516 nt_call_dict_dict[refPOS+ntPOS][query_seq[ntPOS]] += abund_count
517
518
519 q_pars_pos = q_pars_pos + run_length
520 query_pos = query_pos + run_length
521
522 run_length = 0
523
524
525 else:
526 run_length = (10 * run_length) + int(C)
527 # END CIGAR PARSE
528
529
530
531 if len(mutations) == 0: # record reference counts
532 if args.wgs == 0:
533 try:
534 seq_species['Reference'] += abund_count
535 except:
536 seq_species['Reference'] = abund_count
537 else:
538 try:
539 seq_species[str(POS)+' Ref '+str(POS+q_pars_pos)] += abund_count
540 except:
541 seq_species[str(POS)+' Ref '+str(POS+q_pars_pos)] = abund_count
542
543 else: # record variants and counts
544 if args.AAreport == 1 and args.AAcodonasMNP == 1:
545 codonchecked = []
546 codon = ''
547 skip = 0
548 MNP = ''
549 for i in range(0, len(mutations)): # checking for MNP
550 if 'Del' in mutations[i] or 'insert' in mutations[i]:
551 codonchecked.append(mutations[i])
552 elif skip > 0:
553 skip -= 1
554 else:
555 mut1POS = int(''.join([c for c in mutations[i] if c.isdigit()]))
556 try:
557 mutations[i+1]
558 except:
559 AAinfo = singletCodon(mut1POS, mutations[i][-1], ref)
560 codonchecked.append(mutations[i]+'('+refprot[AAinfo[0]-1]+str(AAinfo[0])+AAinfo[1]+')')
561 else:
562
563 if mut1POS % 3 == 0:
564 AAinfo = singletCodon(mut1POS, mutations[i][-1], ref)
565 codonchecked.append(mutations[i]+'('+refprot[AAinfo[0]-1]+str(AAinfo[0])+AAinfo[1]+')')
566 else:
567 mut2POS = int(''.join([c for c in mutations[i+1] if c.isdigit()]))
568 if mut2POS - mut1POS < 3:
569 AAPOS = mut1POS // 3
570 if mut1POS % 3 == 1:
571 if mut2POS % 3 == 0:
572 codon = mutations[i][-1]+ref[1][mut1POS]+mutations[i+1][-1]
573 MNP = mutations[i][-1]+'r'+mutations[i+1][-1]
574 skip = 1
575 else:
576 try:
577 mutations[i+2]
578 except:
579 codon = mutations[i][-1]+mutations[i+1][-1]+ref[1][mut2POS]
580 MNP = mutations[i][-1]+mutations[i+1][-1]+'r'
581 skip = 1
582 else:
583 mut3POS = int(''.join([c for c in mutations[i+2] if c.isdigit()]))
584 if mut2POS == mut3POS -1:
585 codon = mutations[i][-1]+mutations[i+1][-1]+mutations[i+2][-1]
586 MNP = mutations[i][-1]+mutations[i+1][-1]+mutations[i+2][-1]
587 skip = 2
588 else:
589 codon = mutations[i][-1]+mutations[i+1][-1]+ref[1][mut2POS]
590 MNP = mutations[i][-1]+mutations[i+1][-1]+'r'
591 skip = 1
592 codonchecked.append(str(mut1POS)+MNP+'('+refprot[AAPOS]+str(AAPOS+1)+AAcall(codon)+')')
593 elif mut2POS - mut1POS == 1:
594 codon = ref[1][mut1POS-2]+mutations[i][-1]+mutations[i+1][-1]
595 MNP = mutations[i][-1]+mutations[i+1][-1]
596 skip = 1
597 codonchecked.append(str(mut1POS)+MNP+'('+refprot[AAPOS]+str(AAPOS+1)+AAcall(codon)+')')
598 else:
599 AAinfo = singletCodon(mut1POS, mutations[i][-1], ref)
600 codonchecked.append(mutations[i]+'('+refprot[AAinfo[0]-1]+str(AAinfo[0])+AAinfo[1]+')')
601
602
603 else:
604 AAinfo = singletCodon(mut1POS, mutations[i][-1], ref)
605 codonchecked.append(mutations[i]+'('+refprot[AAinfo[0]-1]+str(AAinfo[0])+AAinfo[1]+')')
606 mutations = " ".join(codonchecked)
607 else:
608 mutations = " ".join(mutations)
609
610 if args.wgs == 0:
611 try:
612 seq_species[mutations] += abund_count
613 except:
614 seq_species[mutations] = abund_count
615 else:
616 try:
617 seq_species[str(POS)+' '+mutations+' '+str(POS+q_pars_pos)] += abund_count
618 except:
619 seq_species[str(POS)+' '+mutations+' '+str(POS+q_pars_pos)] = abund_count
620
621 if lastPOS < POS+q_pars_pos:
622 lastPOS = POS+q_pars_pos
623 for i in range(POS, POS+q_pars_pos): # update coverage
624 try:
625 coverage[i] += abund_count
626 except:
627 coverage[i] = abund_count
628
629 args.Sam_file.close()
630 # END SAM LINES
631
632 total_mapped_reads = sam_read_count
633
634 if args.seq == 1: # output the sequence
635 seq_fh = open(args.outfile1, "w")
636 seq_fh.write(args.samp+"("+str(sam_read_count)+")\n")
637 seq_fh.write("Unique Sequence\tCount\tAbundance\n")
638
639 sorted_seq = sorted(seq_species, key=seq_species.__getitem__, reverse=True)
640 for key in sorted_seq:
641 if seq_species[key] >= args.min_count:
642 if (seq_species[key] / sam_read_count >= args.min_samp_abund) and args.wgs == 0:
643 seq_fh.write(f"{key}\t{seq_species[key]}\t{(seq_species[key]/sam_read_count):.3f}\n")
644 seq_array.append(f"{key}\t{seq_species[key]}\t{(seq_species[key]/sam_read_count):.3f}")
645 elif args.wgs == 1:
646 splitseqs = key.split()
647 cov = []
648 for x in range(int(splitseqs[0]), int(splitseqs[-1])):
649 cov.append(coverage[x])
650 min_cov = min(cov)
651 if (seq_species[key]/min_cov >= args.min_samp_abund):
652 seq_fh.write(f"{key}\t{seq_species[key]}\t{(seq_species[key]/min_cov):.3f}\n")
653 else:
654 break
655
656 seq_fh.close()
657 # END SEQ OUT
658
659 if args.indel == 1: # and len(indel_dict) > 0: # output indels, if there are any
660 sorted_indels = sorted(indel_dict, key=indel_dict.__getitem__, reverse=True)
661 indels_to_write = []
662 for key in sorted_indels:
663 if indel_dict[key] >= args.min_count:
664 if indel_dict[key] / sam_read_count >= args.min_samp_abund and args.wgs == 0:
665 indels_to_write.append(f"{key}\t{indel_dict[key]}\t{(indel_dict[key]/sam_read_count):.3f}\n")
666 elif args.wgs == 1:
667 indelPOS = ''
668 for c in key:
669 if c.isdigit():
670 indelPOS += c
671 else:
672 break
673 indelPOS = int(indelPOS)
674 if indel_dict[key] / coverage[indelPOS] >= args.min_samp_abund:
675 indels_to_write.append(f"{key}\t{indel_dict[key]}\t{(indel_dict[key] / coverage[indelPOS]):.3f}\n")
676 else:
677 break
678 # if len(indels_to_write) > 0:
679 indel_fh = open(args.outfile2, "w")
680 indel_fh.write(args.samp+"("+str(sam_read_count)+")\n")
681 indel_fh.write("Indel\tCount\tAbundance\n")
682 for indel_entry in indels_to_write:
683 indel_fh.write(indel_entry)
684 indel_fh.close()
685 # END INDEL OUT
686
687 if args.nt_call == 1: # out put nt calls
688 ntcall_fh = open(args.outfile3, "w")
689 ntcall_fh.write(args.samp+"("+str(sam_read_count)+")\n")
690 if args.ntvar == 1:
691 ntcallv_fh = open(args.outfile7, "w")
692 ntcallv_fh.write(args.samp+"("+str(sam_read_count)+")\n")
693 sorted_POS = sorted(nt_call_dict_dict)
694 if args.AAreport == 1:
695 ntcall_fh.write("Position\tref NT\tAA POS\tref AA\tA\tT\tC\tG\t-\tTotal\tPrimary NT\tCounts\tAbundance\tPrimary Seq AA\tsingle nt AA\tSecondary NT\tCounts\tAbundance\tAA\tTertiary NT\tCounts\tAbundance\tAA\n")
696 if args.ntvar == 1:
697 ntcallv_fh.write("Position\tref NT\tAA POS\tref AA\tA\tT\tC\tG\t-\tTotal\tPrimary NT\tCounts\tAbundance\tPrimary Seq AA\tsingle nt AA\tSecondary NT\tCounts\tAbundance\tAA\tTertiary NT\tCounts\tAbundance\tAA\n")
698 for POS in sorted_POS:
699 try:
700 total = coverage[POS]
701 except:
702 total = 0
703 if total >= (sam_read_count * args.ntabund):
704 AAinfo = singletCodon(POS, ref[1][POS-1], ref)
705 POS_calls = {}
706 for key in nt_call_dict_dict[POS]:
707 POS_calls[key] = nt_call_dict_dict[POS][key]
708 sorted_calls = sorted(POS_calls, key=POS_calls.__getitem__, reverse=True)
709
710 ntcall_fh.write(str(POS)+"\t"+ref[1][POS-1]+"\t"+str(AAinfo[0])+"\t"+AAinfo[1])
711 ntcall_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
712 ntcall_fh.write("\t"+str(total)+"\t"+sorted_calls[0]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[0]]))
713 ntcall_fh.write(f"\t{(nt_call_dict_dict[POS][sorted_calls[0]]/total):.3f}")
714 if sorted_calls[0] != ref[1][POS-1]:
715 if args.ntvar == 1:
716 ntcallv_fh.write(str(POS)+"\t"+ref[1][POS-1]+"\t"+str(AAinfo[0])+"\t"+AAinfo[1])
717 ntcallv_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
718 ntcallv_fh.write("\t"+str(total)+"\t"+sorted_calls[0]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[0]]))
719 ntcallv_fh.write(f"\t{(nt_call_dict_dict[POS][sorted_calls[0]]/total):.3f}")
720
721 mod = (POS)%3
722
723 if mod == 0:
724 try:
725 codon = ref[1][POS-3]+ref[1][POS-2]+sorted_calls[0]
726 except:
727 codon = 'NNN'
728 elif mod == 2:
729 try:
730 codon = ref[1][POS-2]+sorted_calls[0]+ref[1][POS]
731 except:
732 codon = 'NNN'
733 elif mod == 1:
734 try:
735 codon = sorted_calls[0]+ref[1][POS]+ref[1][POS+1]
736 except:
737 codon = 'NNN'
738 ntcall_fh.write("\t"+AAcall(codon)+"\t"+singletCodon(POS, sorted_calls[0], ref)[1])
739 if args.ntvar == 1:
740 ntcallv_fh.write("\t"+AAcall(codon)+"\t"+singletCodon(POS, sorted_calls[0], ref)[1])
741 if (nt_call_dict_dict[POS][sorted_calls[1]] >= args.min_count) and (nt_call_dict_dict[POS][sorted_calls[1]] /total >= args.min_samp_abund):
742 if sorted_calls[1] != ref[1][POS-1] and args.ntvar == 1:
743 ntcallv_fh.write(f"\t{sorted_calls[1]}\t{nt_call_dict_dict[POS][sorted_calls[1]]}\t{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}"+"\t"+singletCodon(POS, sorted_calls[1], ref)[1])
744 ntcall_fh.write(f"\t{sorted_calls[1]}\t{nt_call_dict_dict[POS][sorted_calls[1]]}\t{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}"+"\t"+singletCodon(POS, sorted_calls[1], ref)[1])
745 if nt_call_dict_dict[POS][sorted_calls[2]] >= args.min_count:
746 if nt_call_dict_dict[POS][sorted_calls[2]] / total >= args.min_samp_abund:
747 if sorted_calls[2] != ref[1][POS-1] and args.ntvar == 1:
748 ntcallv_fh.write(f"\t{sorted_calls[2]}\t{nt_call_dict_dict[POS][sorted_calls[2]]}\t{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}\t{singletCodon(POS, sorted_calls[2], ref)[1]}")
749 ntcall_fh.write(f"\t{sorted_calls[2]}\t{nt_call_dict_dict[POS][sorted_calls[2]]}\t{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}\t{singletCodon(POS, sorted_calls[2], ref)[1]}")
750
751 if args.ntvar == 1:
752 ntcallv_fh.write("\n")
753 elif (nt_call_dict_dict[POS][sorted_calls[1]] >= args.min_count) and (nt_call_dict_dict[POS][sorted_calls[1]] / total >= args.min_samp_abund):
754 if args.ntvar == 1:
755 ntcallv_fh.write(str(POS)+"\t"+ref[1][POS-1]+"\t"+str(AAinfo[0])+"\t"+AAinfo[1])
756 ntcallv_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
757 ntcallv_fh.write("\t"+str(total)+"\t\t")
758 ntcallv_fh.write(f"\t")
759 ntcallv_fh.write("\t\t")
760 ntcallv_fh.write(f"\t{sorted_calls[1]}\t{nt_call_dict_dict[POS][sorted_calls[1]]}\t{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}"+"\t"+singletCodon(POS, sorted_calls[1], ref)[1])
761 ntcall_fh.write("\t\t")
762 ntcall_fh.write(f"\t{sorted_calls[1]}\t{nt_call_dict_dict[POS][sorted_calls[1]]}\t{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}"+"\t"+singletCodon(POS, sorted_calls[1], ref)[1])
763
764 if (nt_call_dict_dict[POS][sorted_calls[2]] >= args.min_count) and (nt_call_dict_dict[POS][sorted_calls[2]] /total >= args.min_samp_abund):
765 ntcall_fh.write(f"\t{sorted_calls[2]}\t{nt_call_dict_dict[POS][sorted_calls[2]]}\t{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}\t{singletCodon(POS, sorted_calls[2], ref)[1]}")
766 if sorted_calls[2] != ref[1][POS-1] and args.ntvar == 1:
767 ntcallv_fh.write(f"\t{sorted_calls[2]}\t{nt_call_dict_dict[POS][sorted_calls[2]]}\t{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}\t{singletCodon(POS, sorted_calls[2], ref)[1]}")
768 if args.ntvar == 1:
769 ntcallv_fh.write("\n")
770
771 ntcall_fh.write("\n")
772 else:
773 ntcall_fh.write("Position\tref NT\tA\tT\tC\tG\t-\tTotal\tPrimary NT\tCounts\tAbundance\tSecondary NT\tCounts\tAbundance\tTertiary NT\tCounts\tAbundance\n")
774 if args.ntvar == 1:
775 ntcallv_fh.write("Position\tref NT\tA\tT\tC\tG\t-\tTotal\tPrimary NT\tCounts\tAbundance\tSecondary NT\tCounts\tAbundance\tTertiary NT\tCounts\tAbundance\n")
776
777 for POS in sorted_POS:
778 try:
779 total = coverage[POS] # sum(nt_call_dict_dict[POS].values())
780 except:
781 total = 0
782 if total >= (sam_read_count * args.ntabund):
783 POS_calls = {}
784 for key in nt_call_dict_dict[POS]:
785 POS_calls[key] = nt_call_dict_dict[POS][key]
786 sorted_calls = sorted(POS_calls, key=POS_calls.__getitem__, reverse=True)
787
788 ntcall_fh.write(str(POS)+"\t"+ref[1][POS-1])
789 ntcall_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
790 ntcall_fh.write("\t"+str(total)+"\t"+sorted_calls[0]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[0]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[0]]/total):.3f}")
791
792 if sorted_calls[0] != ref[1][POS-1]:
793 if args.ntvar == 1:
794 ntcallv_fh.write(str(POS)+"\t"+ref[1][POS-1])
795 ntcallv_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
796 ntcallv_fh.write("\t"+str(total)+"\t"+sorted_calls[0]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[0]]))
797 ntcallv_fh.write(f"\t{(nt_call_dict_dict[POS][sorted_calls[0]]/total):.3f}")
798 if (nt_call_dict_dict[POS][sorted_calls[1]] > args.min_count) and (nt_call_dict_dict[POS][sorted_calls[1]] /total > args.min_samp_abund):
799 ntcall_fh.write("\t"+sorted_calls[1]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[1]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}")
800 if sorted_calls[1] != ref[1][POS-1] and args.ntvar == 1:
801 ntcallv_fh.write("\t"+sorted_calls[1]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[1]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}")
802 if (nt_call_dict_dict[POS][sorted_calls[2]] > args.min_count) and (nt_call_dict_dict[POS][sorted_calls[2]] /total > args.min_samp_abund):
803 ntcall_fh.write("\t"+sorted_calls[2]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[2]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}")
804 if sorted_calls[2] != ref[1][POS-1] and args.ntvar == 1:
805 ntcallv_fh.write("\t"+sorted_calls[2]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[2]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}")
806 if args.ntvar == 1:
807 ntcallv_fh.write("\n")
808
809 elif (nt_call_dict_dict[POS][sorted_calls[1]] > args.min_count) and (nt_call_dict_dict[POS][sorted_calls[1]] /total > args.min_samp_abund):
810 if args.ntvar == 1:
811 ntcallv_fh.write(str(POS)+"\t"+ref[1][POS-1])
812 ntcallv_fh.write("\t"+str(nt_call_dict_dict[POS]['A'])+"\t"+str(nt_call_dict_dict[POS]['T'])+"\t"+str(nt_call_dict_dict[POS]['C'])+"\t"+str(nt_call_dict_dict[POS]['G'])+"\t"+str(nt_call_dict_dict[POS]['-']))
813 ntcallv_fh.write("\t"+str(total)+"\t\t")
814 ntcallv_fh.write(f"\t")
815 ntcallv_fh.write("\t"+sorted_calls[1]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[1]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}")
816 ntcall_fh.write("\t"+sorted_calls[1]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[1]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[1]]/total):.3f}")
817 if (nt_call_dict_dict[POS][sorted_calls[2]] > args.min_count) and(nt_call_dict_dict[POS][sorted_calls[2]] /total > args.min_samp_abund):
818 ntcall_fh.write("\t"+sorted_calls[2]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[2]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}")
819 if sorted_calls[2] != ref[1][POS-1] and args.ntvar == 1:
820 ntcallv_fh.write("\t"+sorted_calls[2]+"\t"+str(nt_call_dict_dict[POS][sorted_calls[2]])+"\t"+f"{(nt_call_dict_dict[POS][sorted_calls[2]]/total):.3f}")
821 if args.ntvar == 1:
822 ntcallv_fh.write("\n")
823
824 ntcall_fh.write("\n")
825
826 ntcall_fh.close()
827 if args.ntvar == 1:
828 ntcallv_fh.close()
829 # END NT CALL OUT
830 if args.covar == 1: # output covariants
831 testtrack = 0
832 combinations = {}
833 for sequence in seq_species:
834 if args.wgs == 0:
835 singles = sequence.split()
836 else:
837 singles = (sequence.split())[1:-1]
838 if len(singles) <= args.max_dist and singles[0] != 'Ref':
839 for combo in getCombos(singles, args.max_covar):
840 if not combo in combinations:
841 combinations[combo] = seq_species[sequence]
842 else:
843 combinations[combo] += seq_species[sequence]
844
845 covar_fh = open(args.outfile4, "w")
846 covar_fh.write(args.samp+"("+str(sam_read_count)+")\n")
847 covar_fh.write("Co-Variants\tCount\tAbundance\n")
848 sortedcombos = sorted(combinations, key=combinations.__getitem__, reverse=True)
849 for key in sortedcombos:
850 if (combinations[key] >= args.min_count):
851 if (combinations[key] / sam_read_count >= args.min_samp_abund) and args.wgs == 0:
852 covar_fh.write(key+"\t"+str(combinations[key])+"\t"+f"{(combinations[key]/sam_read_count):.3f}\n")
853 covar_array.append(key+"\t"+str(combinations[key])+"\t"+f"{(combinations[key]/sam_read_count):.3f}")
854 elif args.wgs == 1:
855 coveragepercent = 0
856 splitcombos = key.split()
857 if len(splitcombos) == 1:
858 coveragePOS = ''
859 for c in key:
860 if c.isdigit():
861 coveragePOS += c
862 else:
863 break
864 coveragepercent = combinations[key] / coverage[int(coveragePOS)]
865 else:
866 startcovPOS = ''
867 for c in splitcombos[0]:
868 if c.isdigit():
869 startcovPOS += c
870 else:
871 break
872 endcovPOS = ''
873 for c in splitcombos[-1]:
874 if c.isdigit():
875 endcovPOS += c
876 else:
877 break
878 coveragevals = []
879 for i in range(int(startcovPOS), int(endcovPOS)+1):
880 coveragevals.append(coverage[i])
881 mincov = min(coverval for coverval in coveragevals)
882 coveragepercent = combinations[key] / mincov
883 if coveragepercent >= args.min_samp_abund:
884 covar_fh.write(f"{key}\t{combinations[key]}\t{coveragepercent:.3f}\n")
885
886 covar_fh.close()
887
888
889 # END COVAR OUT
890
891 return(covar_array, seq_array, sam_read_count)
892
893 def cvdeconv(args, covardict, seqdict): # covar deconvolution process
894
895 passedseqs = {}
896 preservedseqs = {}
897 for seq in seqdict: # pass check for actual : expected abundance
898 if seq != 'total' and seq != 'singles':
899
900 splitseq = seq.split(' ')
901 abund = 1
902 for sing in splitseq:
903 try:
904 abund = abund * (covardict[sing] / covardict['total'])
905 except:
906 abund = abund * (seqdict[seq] / seqdict['total'])
907
908 try:
909 covarabund = covardict[seq]/covardict['total']
910 except:
911 covarabund = seqdict[seq]/seqdict['total']
912 covardict[seq] = seqdict[seq]
913
914 if covarabund >= args.autopass:
915 preservedseqs[seq] = max(1, args.beta, (covarabund / abund))
916 elif covarabund >= abund * args.beta:
917 passedseqs[seq] = covarabund / abund
918 elif len(seq.split(' ')) == 1:
919 passedseqs[seq] = max(1, args.beta)
920
921 if args.min_samp_abund < 1:
922 min_count = args.min_samp_abund * covardict['total']
923 else:
924 min_count = args.min_samp_abund
925
926 # sort passed covars
927 lensortedpassed = sorted(passedseqs, key = lambda key : len(key.split(' ')), reverse=True)
928 ratiolensortedpassed = sorted(lensortedpassed, key = lambda key : passedseqs[key], reverse = True)
929 sortedsingles = sorted(covardict['singles'], key = covardict['singles'].__getitem__)
930 deconved = {}
931 for seq in ratiolensortedpassed: # reassign counts
932 singles = seq.split(' ')
933 first = 0
934 rem_count = 0
935 for sing in sortedsingles:
936 if sing in singles:
937 if covardict['singles'][sing] > 0:
938 if first == 0:
939 first = 1
940 rem_count = covardict['singles'][sing]
941 covardict['singles'][sing] = 0
942 deconved[seq] = rem_count
943 else:
944 covardict['singles'][sing] = covardict['singles'][sing] - rem_count
945 else:
946 break
947 sortedsingles = sorted(covardict['singles'], key = covardict['singles'].__getitem__)
948
949 sortedpreserved = sorted(preservedseqs, key = lambda key : covardict[key])
950
951 for seq in sortedpreserved:
952 singles = seq.split(' ')
953 first = 0
954 rem_count = 0
955 for sing in sortedsingles:
956 if sing in singles:
957 if covardict['singles'][sing] > 0:
958 if first == 0:
959 first = 1
960 rem_count = covardict['singles'][sing]
961 covardict['singles'][sing] = 0
962 deconved[seq] = rem_count
963 else:
964 covardict['singles'][sing] = covardict['singles'][sing] - rem_count
965 else:
966 break
967 sortedsingles = sorted(covardict['singles'], key = covardict['singles'].__getitem__)
968
969
970 newtotal = sum(deconved.values())
971 fh_deconv = open(args.outfile6, "w")
972 fh_deconv.write(f"{args.samp}({covardict['total']}) | ({newtotal})\tCount\tAbundance\n")
973 sorted_deconved = sorted(deconved, key = deconved.__getitem__, reverse = True)
974 for seq in sorted_deconved: # write deconv
975 if deconved[seq] >= min_count:
976 fh_deconv.write(f"{seq}\t{deconved[seq]}\t{(deconved[seq]/newtotal):.3f}\n")
977 fh_deconv.close()
978 # END COVAR DECONV OUT
979
980 return()
981
982 def dechim(args, seqs): # processing sequence dictionary to remove chimeras
983 total = seqs['total']
984 del seqs['total']
985 sorted_seqs = sorted(seqs, key=seqs.__getitem__) # sort sequences by abundance, least to greatest
986 chimeras = []
987 for seq in sorted_seqs:
988 pot_chim = ['Reference']+seq.split()+['Reference']
989 chim_halves = []
990 for i in range(0, len(pot_chim)-1): # create dimera halves
991 chim_halves.append([pot_chim[:i+1], pot_chim[i+1:]])
992 parent_pairs = []
993
994 for dimera in chim_halves: # check for potential parents
995
996 pot_left = []
997 pot_rt = []
998 lft_len = len(dimera[0])
999 rt_len = len(dimera[1])
1000 for pseq in sorted_seqs:
1001 if not seq == pseq:
1002 if (seqs[pseq] >= (seqs[seq] * args.foldab)):
1003 pot_par = ['Reference']+pseq.split()+['Reference']
1004 if dimera[0] == pot_par[:lft_len]:
1005 pot_left.append(pot_par)
1006 if ((len(pot_par) >= rt_len) and (dimera[1] == pot_par[(len(pot_par)-rt_len):])):
1007 pot_rt.append(pot_par)
1008
1009 if (len(pot_left) > 0 and len(pot_rt) > 0 ):
1010 for left_par in pot_left: # check potential parents' pairing
1011 for rt_par in pot_rt:
1012 if left_par != rt_par:
1013 left_break = left_par[lft_len]
1014 rt_break = rt_par[(len(rt_par)-rt_len)-1]
1015 if left_break == 'Reference' or rt_break == 'Reference':
1016 parent_pairs.append([' '.join(left_par[1:-1]), ' '.join(rt_par[1:-1])])
1017 else:
1018 left_break_POS = ''
1019 for c in left_break:
1020 if c.isdigit():
1021 left_break_POS += c
1022 else:
1023 break
1024
1025 rt_break_POS = ''
1026 for c in rt_break:
1027 if c.isdigit():
1028 rt_break_POS += c
1029 else:
1030 break
1031
1032 if int(left_break_POS) > int(rt_break_POS):
1033 parent_pairs.append([' '.join(left_par[1:-1]), ' '.join(rt_par[1:-1])])
1034
1035 par_tot_abund = 0
1036 pair_probs = []
1037 for parents in parent_pairs: # calc 'expected' abundance
1038 pair_prob = (seqs[parents[0]] / total) * (seqs[parents[1]] / total)
1039 par_tot_abund += pair_prob
1040 pair_probs.append(pair_prob)
1041
1042 recomb_count = par_tot_abund * total
1043
1044 if not seqs[seq] >= recomb_count * args.alpha: # chimera check
1045 redist_count = float(seqs[seq])
1046 chimeras.append(seq)
1047 seqs[seq] = 0
1048 if args.redist == 1: # redist counts of chimera
1049 toadd = {}
1050 for i in range(0, len(parent_pairs)):
1051 counts_to_redist = (redist_count * (pair_probs[i]/par_tot_abund))/2
1052 seqs[parent_pairs[i][0]] += counts_to_redist
1053 seqs[parent_pairs[i][1]] += counts_to_redist
1054
1055
1056
1057 for chim in chimeras: # remove chimeras
1058 del seqs[chim]
1059
1060 total = sum(seqs.values())
1061
1062
1063 # total = sum(seqs.values)
1064 seqs['total'] = total
1065
1066 return(seqs)
1067
1068 def chimrm(args, seqs): # chimera removed process
1069
1070 pre_len = len(seqs)
1071 inf_loop_shield = 0
1072 while True: # send sequences for chimera removal while chimeras are still found
1073 dechim(args, seqs)
1074 post_len = len(seqs)
1075 inf_loop_shield += 1
1076 if post_len >= pre_len:
1077 break
1078 if inf_loop_shield > 100:
1079 break
1080 pre_len = len(seqs)
1081
1082 total = seqs['total']
1083 del seqs['total']
1084 if args.min_samp_abund < 1:
1085 min_count = args.min_samp_abund * total
1086 else:
1087 min_count = args.min_samp_abund
1088 sorted_seqs = sorted(seqs, key=seqs.__getitem__, reverse=True)
1089 fh_dechime = open(args.outfile5, 'w')
1090 fh_dechime.write(f"{args.samp}({int(total)})\n")
1091 fh_dechime.write("Sequences\tCount\tAbundance\n")
1092 for seq in seqs: # write chim_rm seqs
1093 abund = seqs[seq]/total
1094 if seqs[seq] >= min_count:
1095 fh_dechime.write(f"{seq}\t{round(seqs[seq])}\t{abund:.3f}\n")
1096
1097 fh_dechime.close()
1098 return()
1099
1100 def main():
1101
1102 args = arg_parser() # getting command line arguments
1103
1104 covar_array = []
1105 seq_array = []
1106 total_mapped_reads = 0
1107
1108 if args.ref:
1109 ref = get_ref(args.ref) # get the reference ID and sequence from the FASTA file
1110 args.ref.close()
1111 if ref[1] == '':
1112 print('Reference not recognized as a Fasta format, skipping SAM parsing')
1113 else:
1114 refprot = ''
1115 if args.AAreport == 1: # make an Amino Acid sequence based on the reference sequence
1116 for x in range(0, (len(ref[1])-1)//3):
1117 AA = AAcall(ref[1][x*3]+ref[1][x*3+1]+ref[1][x*3+2])
1118 refprot = refprot + AA
1119 if (len(ref[1])-1)%3 != 0:
1120 refprot = refprot + '?'
1121
1122 covar_array, seq_array, total_mapped_reads = SAMparse(args, ref, refprot, args.Sam_file)
1123
1124 else:
1125 print('No reference provided, skipping SAM parsing')
1126
1127 in_covars = {'total' : total_mapped_reads,
1128 'singles' : {}
1129 }
1130 in_seqs = {'total' : total_mapped_reads}
1131 # Begin chimera removal if enabled
1132 if args.chim_rm == 1 or args.deconv == 1:
1133 if args.deconv == 1:
1134 for line in covar_array:
1135 lineparts = line.split("\t")
1136 in_covars[lineparts[0]] = int(lineparts[1])
1137 if len(lineparts[0].split(' ')) == 1:
1138 in_covars['singles'][lineparts[0]] = int(lineparts[1])
1139
1140 for line in seq_array:
1141 lineparts = line.split("\t")
1142 in_seqs[lineparts[0]] = float(lineparts[1])
1143
1144 if args.deconv == 1:
1145 cvdeconv(args, in_covars, in_seqs)
1146
1147 if args.chim_rm == 1:
1148 chimrm(args, in_seqs)
1149
1150 if __name__ == '__main__':
1151 main()