-
Notifications
You must be signed in to change notification settings - Fork 0
/
VCFgenie.py
executable file
·2019 lines (1674 loc) · 102 KB
/
VCFgenie.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#! /usr/bin/env python
"""
Purpose: Dynamically filter within-sample (pooled sequencing) variants to control false-positive count
Author : Chase W. Nelson <[email protected]>
Cite : https://github.com/chasewnelson/SARS-CoV-2-ORF3d
Date : 2021-12-04
Details: The script applies a dynamic (site-by-site) filter using the binomial
distribution and the user-provided parameters. Must be called in a working
directory with the target VCF files, automatically detected using their
'.vcf' extension. All VCF files in the directory will be used.
The user-provided parameter num_samples may correspond to the number of VCF
files in the directory if none others have been analyzed in the same
analysis; alternatively, the parameter may exceed the number of VCF files
if other samples have been analyzed as part of the same analysis.
The user may specify which FORMAT keys to use for read depth (coverage)
(DP) and allele frequency (AF). ALT allele counts (AC) will be inferred.
Multiallelic records (more than one ALT allele) will be printed across
multiple lines, one per ALT.
Each VCF file must contain the called SNPs for exactly 1 pooled sequencing
sample, and therefore must contain exactly 8-10 columns:
(1) CHROM
(2) POS
(3) ID
(4) REF
(5) ALT
(6) QUAL
(7) FILTER
(8) INFO
(9) FORMAT [OPTIONAL]
(10) <sample_name> [OPTIONAL]
The binomial filter uses the scipy function:
binom.cdf(x, n, p)
where:
n = total reads
p = error rate / 3
x = number of reads of allele
FP cutoff = (1 - binom.cdf(x - 1, n, p)) * seq_len * num_samples
Note that x - 1 is used because we subtract the cumulative probability
BEFORE the given number of variants.
The equivalent function in base R is: pbinom(x, n, p)
Definitions:
-REF = reference sequence against which variants were called
-ALT = alternative allele
-AF = allele frequency of alternative allele
-MAF = minor allele frequency
-DP = depth (coverage)
-fixedREF = site fixed for the REF allele (100% reference)
-fixedALT = site fixed for a particular ALT allele (100% non-reference)
-fail = allele fails criteria (all-inclusive)
-failDP = read depth (coverage) fails min_DP
-failZeroAC = allele fails when 0 reads support it
-failAC = allele count (whether REF or ALT) fails min_AC
-failMinAF = allele frequency (whether REF or ALT) fails min_AF
-failMaxAF = allele frequency (whether REF or ALT) fails max_AF
-failINFO: one or more of the user-provided INFO_rules fails
-failsample: one or more of the user-provided sample_rules fails
-failp = allele fails p_cutoff
-pass = allele passes all criteria
TODO:
-whether INFO or INDIVIDUAL SAMPLES examined
EXAMPLES:
$ VCFgenie.py --error_rate=0.0001 --VCF_files=example.vcf
$ VCFgenie.py --error_rate=0.0001 --p_cutoff=0.000001 --AC_key=FAO --AC_key_new=FAO --AF_key=AF --AF_key_new=AF --DP_key=FDP --overwrite_INFO --VCF_files=example.vcf
Returns:
1) summary statistics about the variants that were processed (STDOUT)
2) one *_filtered.vcf file for each *.vcf file in the working directory (to --out_dir)
"""
import argparse
import numpy as np
import operator as op
import os
import re
import sys
import vcf
from collections import defaultdict
from scipy.stats import binom
from typing import Any, Dict, List, NamedTuple, Optional, TextIO
usage = """# -----------------------------------------------------------------------------
VCFgenie.py - Dynamically filter within-sample (pooled sequencing) variants to control false positive calls
# -----------------------------------------------------------------------------
For DOCUMENTATION, run:
$ VCFgenie.py --help
$ pydoc ./VCFgenie.py
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
EXAMPLE:
$ VCFgenie.py --error_rate=0.0001 --VCF_files=example.vcf
# -----------------------------------------------------------------------------
"""
class Args(NamedTuple):
""" Command-line arguments """
# REQUIRED
VCF_files: List[str]
error_rate: float
# seq_len: int
# num_samples: int
# FP_cutoff: float
# OPTIONAL
out_dir: str
p_cutoff: float
AC_key: str
AC_key_new: str
AF_key: str
AF_key_new: str
DP_key: str
PVR_key_new: str
PVA_key_new: str
min_AC: float
min_AF: float
max_AF: float
min_DP: float
INFO_rules: str
sample_rules: str
overwrite_INFO: bool
# -----------------------------------------------------------------------------
def get_args() -> Args:
""" Get command-line arguments """
parser = argparse.ArgumentParser(
description='Dynamically filter variants to control false-positive count | EXAMPLE: VCFgenie.py --error_rate=0.0001 --VCF_files=example.vcf',
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
# Rename "optional" arguments
parser._optionals.title = 'Named arguments'
# -------------------------------------------------------------------------
# REQUIRED
parser.add_argument('-i',
'--VCF_files',
metavar='FILE',
help='input, one or more variant call format (VCF) files [REQUIRED]',
required=True, # 'required' won't work for positionals
nargs='+',
# type=argparse.FileType('r')) # was 'rt' but want .gz as option?
type=str)
parser.add_argument('-e',
'--error_rate',
metavar='float',
help='sequencing error rate for single nucleotide substitutions per base sequenced [REQUIRED]',
required=True,
type=float)
# help='error rate per site expected as the result of sample preparation, library preparation, '
# 'and sequencing (currently assumes all nucleotide errors equally probable) [REQUIRED]',
# parser.add_argument('-L',
# '--seq_len',
# metavar='int',
# help='length of reference sequence (e.g., contig, chromosome, or genome) in nucleotides '
# '[REQUIRED]',
# required=True,
# type=int)
# parser.add_argument('-n',
# '--num_samples',
# metavar='int',
# help='number of samples (VCF files) in full analysis [REQUIRED]',
# required=True,
# type=int)
# parser.add_argument('-f',
# '--FP_cutoff',
# metavar='float',
# help='analysis-wide false-positive (FP) count cutoff [REQUIRED]',
# required=True,
# type=float)
# -------------------------------------------------------------------------
# OPTIONAL - WITH DEFAULT
parser.add_argument('-o',
'--out_dir',
metavar='DIR',
help='output directory name [OPTIONAL]',
required=False,
type=str,
default='VCFgenie_output')
parser.add_argument('-p',
'--p_cutoff',
metavar='float',
help='p-value cutoff [OPTIONAL]',
required=False,
type=float,
default=1.)
parser.add_argument('-c',
'--AC_key',
metavar='str',
help='FORMAT key to use for obtaining ALT allele count [OPTIONAL]',
required=False,
type=str,
default='AC') # 'FAO' for Ion Torrent
parser.add_argument('-C',
'--AC_key_new',
metavar='str',
help='FORMAT key to use for new (filtered) ALT allele count [OPTIONAL]',
required=False,
type=str,
default='NAC')
parser.add_argument('-a',
'--AF_key',
metavar='str',
help='FORMAT key to use for ALT allele frequency [OPTIONAL]',
required=False,
type=str,
default='AF') # 'AF' remains for Ion Torrent
parser.add_argument('-A',
'--AF_key_new',
metavar='str',
help='FORMAT key to use for new (filtered) ALT allele frequency [OPTIONAL]',
required=False,
type=str,
default='NAF')
parser.add_argument('-d',
'--DP_key',
metavar='str',
help='FORMAT key to use for read depth (coverage) [OPTIONAL]',
required=False,
type=str,
default='DP') # 'FDP' for Ion Torrent
parser.add_argument('-v',
'--PVR_key_new',
metavar='str',
help='INFO and FORMAT key to use for p-value for REF allele [OPTIONAL]',
required=False,
type=str,
default='PVR')
parser.add_argument('-V',
'--PVA_key_new',
metavar='str',
help='INFO and FORMAT key to use for p-values for ALT allele(s) [OPTIONAL]',
required=False,
type=str,
default='PVA')
parser.add_argument('-t',
'--min_AC',
metavar='int',
help='allele count cutoff (min reads allowed to support allele, REF or ALT) [OPTIONAL]',
required=False,
type=float, # OK to be a float
default=0)
parser.add_argument('-q',
'--min_AF',
metavar='float',
help='minor allele frequency cutoff (min allowed, REF or ALT) [OPTIONAL]',
required=False,
type=float,
default=0)
parser.add_argument('-Q',
'--max_AF',
metavar='float',
help='major allele frequency cutoff (max allowed, REF or ALT) [OPTIONAL]',
required=False,
type=float,
default=1)
parser.add_argument('-D',
'--min_DP',
metavar='int',
help='read depth (coverage) cutoff (min allowed, REF or ALT) [OPTIONAL]',
required=False,
type=float, # OK to be a float
default=1)
parser.add_argument('-I',
'--INFO_rules',
metavar='str',
help='custom rules for filtering based on the INFO column, formatted as a comma-'
'separated string in the case of multiple rules; must use keys present in the '
'INFO column, and must use the comparison operators "=="/"!="/"<"/"<="/">="/">" '
'(e.g., --INFO_rules="STB>0.5,STB<0.9") [OPTIONAL]',
required=False,
type=str,
default=None)
parser.add_argument('-s',
'--sample_rules',
metavar='str',
help='custom rules for filtering based on the sample column, formatted as a comma-'
'separated string in the case of multiple rules; must use keys present in the '
'FORMAT column, and must use the comparison operators "=="/"!="/"<"/"<="/">="/">" '
'(e.g., --sample_rules="FSRF>=5,FSRR>=5,FSAF>=5,FSAR>=5") [OPTIONAL]',
required=False,
type=str,
default=None)
parser.add_argument('-O',
'--overwrite_INFO',
help='for any data keys matching between INFO and FORMAT, OVERWRITE INFO data with <sample> '
'data [OPTIONAL]',
action='store_true')
args = parser.parse_args()
# -------------------------------------------------------------------------
# VALIDATE arguments
# REQUIRED
# VCF files EXIST
if bad_files := [file for file in args.VCF_files if not os.path.isfile(file)]:
parser.error(f'\n### ERROR: VCF file(s) do not exist: {",".join(bad_files)}')
# VCF files have RECOGNIZABLE HEADER
re_start_pound = re.compile(r'^#')
re_header_start = re.compile(r'^#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t')
for this_VCF_file in args.VCF_files:
missing_header = True
this_VCF_fh = open(this_VCF_file)
for line in this_VCF_fh:
if not re_start_pound.match(line):
break
if re_header_start.match(line):
missing_header = False
break
if missing_header:
parser.error(f'\n### ERROR: VCF_file="{this_VCF_file}" does not contain a recognizable header. '
'Wrong format?\n')
# VCF files RECOGNIZABLE to pyvcf module
# TODO validate AF_key and DP_key here?
for this_VCF_file in args.VCF_files:
try:
this_vcf_Reader = vcf.Reader(filename=this_VCF_file) # opens without error if exists, even if not VCF
except AttributeError:
parser.error(f'\n### ERROR: VCF_file="{this_VCF_file}" not readable by pyvcf. VCF file wrong format?\n')
try:
next(this_vcf_Reader)
except ValueError as e:
parser.error(f'### ERROR: ValueError: {e}\n'
f'### Records in VCF_file="{this_VCF_file}" not readable by pyvcf.\n'
'### VCF file wrong format?\n')
except StopIteration as e:
parser.error(f'\n### ERROR: StopIteration: {e}\n'
f'### VCF_file="{this_VCF_file}" not iterable by pyvcf.\n'
'### VCF may not contain any records (rows; FAOvariants).\n')
# Appropriate numeric values
if args.error_rate < 0 or args.error_rate >= 1 or not args.error_rate:
parser.error(f'\n### ERROR: error_rate "{args.error_rate}" must be >= 0 and < 1')
# if args.seq_len <= 0 or not args.seq_len:
# parser.error(f'\n### ERROR: seq_len "{args.seq_len}" must be > 0 nucleotides')
#
# if args.num_samples <= 0 or not args.num_samples:
# parser.error(f'\n### ERROR: num_samples "{args.num_samples}" must be > 0')
# if args.FP_cutoff <= 0 or not args.FP_cutoff:
# parser.error(f'\n### ERROR: FP_cutoff "{args.FP_cutoff}" must be > 0')
# OPTIONAL
if os.path.isdir(args.out_dir):
# sys.exit(f'\n### ERROR: out_dir "{args.out_dir}" already exists')
parser.error(f'\n### ERROR: out_dir="{args.out_dir}" already exists')
else:
os.makedirs(args.out_dir)
if args.p_cutoff < 0 or args.p_cutoff > 1 or not args.p_cutoff:
parser.error(f'\n### ERROR: p_cutoff "{args.p_cutoff}" must be between 0 and 1 (inclusive)')
# if args.AC_key == args.AC_key_new:
# # parser.error(f'AC_key="{args.AC_key}" may not match AC_key_new="{args.AC_key_new}"')
# print(f'\n### WARNING: AC_key="{args.AC_key_new}" already exists and will be OVERWRITTEN')
#
# if args.AF_key == args.AF_key_new:
# # parser.error(f'AF_key="{args.AF_key}" may not match AF_key_new="{args.AF_key_new}"')
# print(f'\n### WARNING: AF_key="{args.AF_key_new}" already exists and will be OVERWRITTEN')
if args.min_AC < 0:
parser.error(f'\n### ERROR: min_AC "{args.min_AC}" must be >= 0')
if args.min_AF < 0:
parser.error(f'\n### ERROR: min_AF "{args.min_AF}" must be >= 0')
if args.max_AF > 1:
parser.error(f'\n### ERROR: min_AF "{args.max_AF}" must be <= 1')
if args.min_DP < 1:
parser.error(f'\n### ERROR: min_DP "{args.min_DP}" must be >= 1')
# INFO_rules and sample_rules will be validated first thing in main() when they're regexed
return Args(VCF_files=args.VCF_files,
error_rate=args.error_rate,
# seq_len=args.seq_len,
# num_samples=args.num_samples,
# FP_cutoff=args.FP_cutoff,
out_dir=args.out_dir,
p_cutoff=args.p_cutoff,
AC_key=args.AC_key,
AC_key_new=args.AC_key_new,
AF_key=args.AF_key,
AF_key_new=args.AF_key_new,
DP_key=args.DP_key,
PVR_key_new=args.PVR_key_new,
PVA_key_new=args.PVA_key_new,
min_AC=args.min_AC,
min_AF=args.min_AF,
max_AF=args.max_AF,
min_DP=args.min_DP,
INFO_rules=args.INFO_rules,
sample_rules=args.sample_rules,
overwrite_INFO=args.overwrite_INFO)
# -----------------------------------------------------------------------------
def main() -> None:
""" The best criticism of the bad is the practice of the better """
# -------------------------------------------------------------------------
# GATHER arguments
args = get_args()
VCF_files = args.VCF_files
error_rate = args.error_rate
# seq_len = args.seq_len
# num_samples = args.num_samples
# FP_cutoff = args.FP_cutoff
out_dir = args.out_dir # optional
p_cutoff = args.p_cutoff # optional
AC_key = args.AC_key # optional
AC_key_new = args.AC_key_new # optional
AF_key = args.AF_key # optional
AF_key_new = args.AF_key_new # optional
DP_key = args.DP_key # optional
PVR_key_new = args.PVR_key_new # optional
PVA_key_new = args.PVA_key_new # optional
min_AC = args.min_AC # optional
min_AF = args.min_AF # optional
max_AF = args.max_AF # optional
min_DP = args.min_DP # optional
INFO_rules = args.INFO_rules # optional
sample_rules = args.sample_rules # optional
overwrite_INFO = args.overwrite_INFO # optional
# -------------------------------------------------------------------------
# INITIALIZE OUTPUT AND LOG
print(usage)
print('# -----------------------------------------------------------------------------')
# print('LOG')
print(f'LOG:command="{" ".join(sys.argv)}"')
print(f'LOG:cwd="{os.getcwd()}"')
print(f'LOG:error_rate="{error_rate}"')
# print(f'LOG:seq_len="{seq_len}"')
# print(f'LOG:num_samples="{num_samples}"')
# print(f'LOG:FP_cutoff="{FP_cutoff}"')
print(f'LOG:out_dir="{out_dir}"')
print(f'LOG:p_cutoff="{p_cutoff}"')
print(f'LOG:AC_key="{AC_key}"')
print(f'LOG:AC_key_new="{AC_key_new}"')
print(f'LOG:AF_key="{AF_key}"')
print(f'LOG:AF_key_new="{AF_key_new}"')
print(f'LOG:DP_key="{DP_key}"')
print(f'LOG:PVR_key_new="{PVR_key_new}"')
print(f'LOG:PVA_key_new="{PVA_key_new}"')
print(f'LOG:min_AC="{min_AC}"')
print(f'LOG:min_AF="{min_AF}"')
print(f'LOG:max_AF="{max_AF}"')
print(f'LOG:min_DP="{min_DP}"')
print(f'LOG:INFO_rules="{INFO_rules}"')
print(f'LOG:sample_rules="{sample_rules}"')
print(f'LOG:overwrite_INFO="{overwrite_INFO}"', flush=True)
if overwrite_INFO:
print('### WARNING: INFO will be overwritten with <sample> data when possible!')
# WARN about unintended overwriting due to the specified new keys
if AC_key == AC_key_new:
print(f'### WARNING: AC_key_new="{AC_key_new}" already exists as AC_key="{AC_key_new}" and will be OVERWRITTEN')
if AF_key == AF_key_new:
print(f'### WARNING: AF_key_new="{AF_key_new}" already exists as AF_key="{AF_key_new}" and will be OVERWRITTEN')
# TODO: check that ANY of the keys being used aren't already present? Perhaps overkill
# -------------------------------------------------------------------------
# REGEX & TUPLES
# match user-supplied rules
re_rule = re.compile(r'(\w+)([!=<>]+)([.\d]+)')
# VCF file regex
# NOTE: DP and AF are sufficient because the allele count, so inconsistently coded across VCFs, can simply be inferred
# VCF_line_pattern = r"([\w\_\.]+)\t(\d+)\t([^\s]+)\t([acgtuACGTU]+)\t([acgtuACGTU]+)\t(\d+)\t(\w+)\t([\w\=\;\.\,]+)"
# re_VCF_INFO = re.compile(VCF_INFO_pattern)
re_VCF_INFO = re.compile(r'(\w+)=([\w\d/,.\-<>]+)')
# VCF_SB_pattern = r"SB=\d+" # SB=0
# re_VCF_SB = re.compile(VCF_SB_pattern)
# VCF_DP4_grppattern = r"DP4=(\d+),(\d+),(\d+),(\d+)" # DP4=27,235,0,4
# VCF_DP4_grpregex = re.compile(VCF_DP4_pattern)
# VCF_DP_pattern = r"DP=\d+" # DP=266
# re_VCF_DP = re.compile(VCF_DP_pattern)
# VCF_AF_pattern = r"AF=[\d\.]+" # AF=0.015038
# re_VCF_AF = re.compile(r'AF=[\d.,]+')
# Data regex
# re_datum = re.compile(r'[\d/.\-]+')
# re_AF_datum = re.compile(r'([\s\t;])AF=[\d.NA]+')
# re_DP_datum = re.compile(r'([\s\t;])DP=[\d.NA]+')
# re_CSV_datum = re.compile(r'[,\d\w.\-]+')
# VCF metadata block regex
# re_metadata_block = re.compile(r'^##(\w+)=.+$')
# ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
# re_VCF_FORMAT_metadata = re.compile(r'^##(\w+)=<ID=(\w+),Number=(\w+),Type=(\w+),Description="(.+)">$')
re_VCF_INFO_metadata = re.compile(r'^##(\w+)=<ID=(\w+),Number=(\w+),Type=(\w+),Description="(.+)".*>$')
# TODO update regex for recording Source and Version:
# ##INFO=<ID=ID,Number=number,Type=type,Description="description",Source="source",Version="version">
# non-float, non-int regex
# re_non_int = re.compile(r'\D') # not a digit or dot or minus sign
# re_non_float = re.compile(r'^\d^.^-') # not a digit or dot or minus sign
operator_choices = ('==', '!=', '<', '<=', '>=', '>')
decision_choices = ('pass',
'fail', 'failZeroAC', 'failDP', 'failAC', 'failMinAF', 'failMaxAF', 'failINFO', 'failsample',
'failp',
'fixedREF', 'fixedALT',) # 'failFP',
decision_choices_FAIL = ('failZeroAC', 'failDP', 'failAC', 'failMinAF', 'failMaxAF', 'failINFO', 'failsample', 'failp') # 'failFP'
# -------------------------------------------------------------------------
# PARSE INFO_rules and sample_rules into LISTS of 3-TUPLES using re_rule
# rule_FILTER_metadata: List[str] = []
INFO_rule_FILTER_dict: Dict[str, List[str]] = defaultdict(list)
sample_rule_FILTER_dict: Dict[str, List[str]] = defaultdict(list)
# INFO rules
INFO_rule_lt = []
if INFO_rules is not None:
INFO_rule_list = INFO_rules.split(',')
for i, INFO_rule in enumerate(INFO_rule_list):
# key, operator, value = re_rule.match(INFO_rule).groups()
INFO_rule_match = re_rule.match(INFO_rule)
if INFO_rule_match is not None:
INFO_rule_match_groups = INFO_rule_match.groups()
key: str = INFO_rule_match_groups[0]
operator: str = INFO_rule_match_groups[1]
value: Any = INFO_rule_match_groups[2]
if operator in operator_choices:
print(f'LOG:INFO_rule_{i}:key="{key}",operator="{operator}",value="{value}"')
INFO_rule_lt.append((key, operator, value))
INFO_rule_FILTER_dict[key].append(INFO_rule)
# rule_FILTER_metadata.append(f'##FILTER=<ID={},Description="User-provided INFO rule: {INFO_rule}">')
else:
sys.exit(f'\n### ERROR: INFO_rule="{INFO_rule}" does not use an acceptable operator {operator_choices}')
else:
sys.exit(f'\n### ERROR: INFO_rule="{INFO_rule}" does not use an acceptable format or operator {operator_choices}')
# Loop cleanup
del INFO_rule_match_groups
del key
del operator
del value
# print(f'INFO_rule_lt={INFO_rule_lt}')
# SAMPLE rules
sample_rule_lt = []
if sample_rules is not None:
sample_rule_list = sample_rules.split(',')
for i, sample_rule in enumerate(sample_rule_list):
# key, operator, value = re_rule.match(sample_rule).groups()
# print(f'LOG:sample_rule_{i}:key="{key}",operator="{operator}",value="{value}"')
# sample_rule_lt.append((key, operator, value))
# key, operator, value = re_rule.match(sample_rule).groups()
sample_rule_match = re_rule.match(sample_rule)
if sample_rule_match is not None:
key, operator, value = sample_rule_match.groups()
if operator in operator_choices:
print(f'LOG:sample_rule_{i}:key="{key}",operator="{operator}",value="{value}"')
sample_rule_lt.append((key, operator, value))
sample_rule_FILTER_dict[key].append(sample_rule)
# rule_FILTER_metadata.append(f'##FILTER=<ID=sample_rule_{i},Description="User-provided sample rule: {sample_rule}">')
else:
sys.exit(
f'\n### ERROR: sample_rule="{sample_rule}" does not use an acceptable operator {operator_choices}')
else:
sys.exit(f'\n### ERROR: sample_rule="{sample_rule}" does not use an acceptable format or operator {operator_choices}')
# print(f'sample_rule_lt={sample_rule_lt}')
# VERIFY there are no duplicate keys; may have to deal with this later
if not set(INFO_rule_FILTER_dict.keys()).isdisjoint(set(sample_rule_FILTER_dict.keys())):
sys.exit(f'\n### ERROR: overlapping key(s) utilized for --INFO_rules ({set(INFO_rule_FILTER_dict.keys())}) and '
f'--sample_rules ({set(sample_rule_FILTER_dict.keys())}). Contact author for update')
print()
# NEW LINES TO PRINT TO METADATA
new_FILTER_lines: List[str] = []
for key in INFO_rule_FILTER_dict:
new_FILTER_lines.append(f'##FILTER=<ID={key},Description="User-provided INFO rule(s): {",".join(INFO_rule_FILTER_dict[key])}">')
for key in sample_rule_FILTER_dict:
new_FILTER_lines.append(f'##FILTER=<ID={key},Description="User-provided sample rule(s): {",".join(sample_rule_FILTER_dict[key])}">')
# -------------------------------------------------------------------------
# INITIALIZE descriptions of summary statistics types
summaryKey_to_description = {
'fixedREF': 'site fixed for the REF allele (100% reference)',
'fixedALT': 'site fixed for a particular ALT allele (100% non-reference)',
'fail': 'allele fails criteria (all-inclusive)',
'failZeroAC': 'allele fails when 0 reads support it (REF or ALT)',
'failDP': 'read depth (coverage) fails --min_DP',
'failAC': 'allele count (REF or ALT) fails --min_AC',
'failMinAF': 'allele frequency (REF or ALT) fails --min_AF',
'failMaxAF': 'allele frequency (REF or ALT) fails --max_AF',
'failINFO': 'fails one or more of the user-provided --INFO_rules',
'failsample': 'fails one or more of the user-provided --sample_rules',
# 'failFP': 'fails --FP_cutoff',
'failp': 'fails --p_cutoff',
# 'failToFixedALT': 'failing minor allele matched REF; converted to fixed ALT allele',
'pass': 'allele passes all criteria'
}
# # -------------------------------------------------------------------------
# # Gather VCF files in working directory
# VCF_filenames = sorted([name for name in os.listdir('.') if os.path.isfile(name) and name.endswith('.vcf')])
# print('VCF files to examine, gathered from working directory: ' + str(VCF_filenames))
# print('Will write output to files with names of the form: <VCF_name>_filtered.vcf')
# -------------------------------------------------------------------------
# Name VCF files to examine
print('# -----------------------------------------------------------------------------')
# re_file_ext = re.compile(r'.\w+$')
print('VCF files to process (output files will have names of the form "<VCF_root_name>_filtered.vcf"): '
'[IN_FILE_NAME] -> [OUT_FILE_NAME]')
# Store input -> output file names
VCF_to_out_file_dict: Dict[str, str] = defaultdict(str)
VCF_files = sorted(VCF_files)
# for this_VCF_file in sorted(VCF_files):
for this_VCF_file in VCF_files:
# Create out_file name
this_inFile_root, this_inFile_ext = os.path.splitext(this_VCF_file) # takes str, includes '.'
# Input VCF file must end with '.vcf' extension and its implied out_file must not already exist
if this_inFile_ext == '.vcf':
this_outFile_name = this_inFile_root + '_filtered.vcf'
if not os.path.exists(this_outFile_name):
VCF_to_out_file_dict[this_VCF_file] = this_outFile_name
print(f'{this_VCF_file} -> {this_outFile_name}')
else:
sys.exit(f'\n### ERROR: output VCF file {this_VCF_file} already exists\n')
# elif this_inFile_ext == '.gz':
# this_outFile_name = this_inFile_root + '.gz_filtered.vcf'
#
# if not os.path.exists(this_outFile_name):
# VCF_to_out_file_dict[this_VCF_file] = this_outFile_name
# print(f'{this_VCF_file} -> {this_outFile_name}')
# else:
# sys.exit(f'\n### ERROR: output VCF file {this_VCF_file} already exists\n')
else:
# sys.exit(f'\n### ERROR: input VCF file {this_VCF_file} must end with the extension ".vcf" or ".gz"\n')
sys.exit(f'\n### ERROR: input VCF file {this_VCF_file} must end with the extension ".vcf"\n')
print()
# -------------------------------------------------------------------------
# Prepare summary statistics
sample_n = 0
total_record_n = 0 # lines/rows
total_ALLELE_n = 0
total_ALLELE_pass_n = 0
total_ALLELE_fail_n = 0
total_REF_n = 0
total_REF_pass_n = 0
total_REF_fail_n = 0
total_ALT_n = 0
total_ALT_pass_n = 0
total_ALT_fail_n = 0
total_MAJOR_n = 0
total_MAJOR_pass_n = 0
total_MAJOR_fail_n = 0
total_MINOR_n = 0
total_MINOR_pass_n = 0
total_MINOR_fail_n = 0
# count_dd: Dict[str, Dict[str, int]] = defaultdict(dict)
# AF_ddl: Dict[str, Dict[str, float]] = defaultdict(dict)
# AC_ddl: Dict[str, Dict[str, int]] = defaultdict(dict)
# DP_ddl: Dict[str, Dict[str, int]] = defaultdict(dict)
ALLELE_n_dict: Dict[str, int] = defaultdict(int)
ALLELE_AF_dl: Dict[str, List[float]] = defaultdict(list)
ALLELE_AC_dl: Dict[str, List[int]] = defaultdict(list)
ALLELE_DP_dl: Dict[str, List[int]] = defaultdict(list)
REF_n_dict: Dict[str, int] = defaultdict(int)
REF_AF_dl: Dict[str, List[float]] = defaultdict(list)
REF_AC_dl: Dict[str, List[int]] = defaultdict(list)
REF_DP_dl: Dict[str, List[int]] = defaultdict(list)
ALT_n_dict: Dict[str, int] = defaultdict(int)
ALT_AF_dl: Dict[str, List[float]] = defaultdict(list)
ALT_AC_dl: Dict[str, List[int]] = defaultdict(list)
ALT_DP_dl: Dict[str, List[int]] = defaultdict(list)
MAJOR_n_dict: Dict[str, int] = defaultdict(int)
MAJOR_AF_dl: Dict[str, List[float]] = defaultdict(list)
MAJOR_AC_dl: Dict[str, List[int]] = defaultdict(list)
MAJOR_DP_dl: Dict[str, List[int]] = defaultdict(list)
MINOR_n_dict: Dict[str, int] = defaultdict(int)
MINOR_AF_dl: Dict[str, List[float]] = defaultdict(list)
MINOR_AC_dl: Dict[str, List[int]] = defaultdict(list)
MINOR_DP_dl: Dict[str, List[int]] = defaultdict(list)
# -------------------------------------------------------------------------
# LOOP VCF files to EXTRACT SNP DATA
for this_VCF_file in VCF_files:
this_VCF_fh = open(this_VCF_file, 'rt')
# if this_VCF_file.endswith('_filtered.vcf'):
# sys.exit(f'\n### ERROR: VCF file {this_VCF_file} may not already end with "_filtered.vcf", because '
# 'this script will give it that suffix\n')
# # Check file exists - already did this in arguments
# if os.path.isfile(this_VCF_file):
print('# -----------------------------------------------------------------------------')
print(f'Analyzing file={this_VCF_file}')
# ---------------------------------------------------------------------
# OPEN the VCF file object
# vcf_file_object = vcf.Reader(fsock=this_VCF_fh) # CAN be .gz!, CAN be fh
# this_vcf_reader = vcf.Reader(filename=this_VCF_file)
# OPEN output file for writing
# this_outFile_hdl = open(VCF_to_out_file_dict[this_VCF_fh.name], "w")
this_out_file = os.path.join(out_dir, VCF_to_out_file_dict[this_VCF_file])
this_out_fh = open(this_out_file, "wt")
# this_vcf_writer = vcf.Writer(this_outFile_hdl, this_vcf_reader)
# ---------------------------------------------------------------------
# VCF DOCUMENTATION: http://samtools.github.io/hts-specs/VCFv4.2.pdf
# INFO
# ##INFO=<ID=ID,Number=number,Type=type,Description="description",Source="source",Version="version">
# 'Number': an Integer that describes the number of values that can be included with the INFO field
# 1=always a single number
# 2=always a pair of numbers
# A=the field has one value per alternate allele ALT # TODO this is what we need for collapsing MULTIALLELICS back into one record
# R=the field has one value for each possible allele (including the reference REF)
# G=the field has one value for each possible genotype (more relevant to the FORMAT tags)
# .=the number of possible values varies, is unknown, or is unbounded
# 0=for ‘Flag’ type, for which the INFO field does not contain a Value entry
# 'Description': must be surround by "double quotes"
# Types for INFO fields are: Integer, Float, Flag, Character, and String
# ##INFO=<ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
# FILTER
# ## FILTER=<ID=ID,Description="description">
# ---------------------------------------------------------------------
# Prepare additional metadata entries (##), to print for each file directly above the header (#)
# new_FILTER_lines = f'##FILTER=<ID=seq_len,Description="Reference sequence genome length (nucleotides): {seq_len}">\n' \
# f'##FILTER=<ID=error_rate,Description="Sequencing error rate per site (assumes all nucleotides equally probable): {error_rate}">\n' \
# f'##FILTER=<ID=num_samples,Description="Number of samples (VCF files) in full analysis: {num_samples}">\n' \
# f'##FILTER=<ID=FP_cutoff,Description="Analysis-wide false discovery rate (FP) cutoff: {FP_cutoff}">\n' \
# f'##FILTER=<ID=min_DP,Description="Read depth (coverage) cutoff (min allowed): {min_DP}">\n' \
# f'##FILTER=<ID=min_AC,Description="Allele count cutoff (min reads allowed to support minor allele): {min_AC}">\n' \
# f'##FILTER=<ID=min_AF,Description="Minor allele frequency cutoff (min allowed): {min_AF}">'
# new_INFO_lines = f'##INFO=<ID=MULTIALLELIC,Number=0,Type=Flag,Description="Indicates whether a site is multiallelic (more than one ALT allele)\">'
# new_metadata_lines = f'##FILTER=<ID=error_rate,Description="Sequencing error rate per site (assumes all nucleotides equally probable): {error_rate}">\n' \
# f'##FILTER=<ID=seq_len,Description="Reference sequence genome length (nucleotides): {seq_len}">\n' \
# f'##FILTER=<ID=num_samples,Description="Number of samples (VCF files) in full analysis: {num_samples}">\n' \
# f'##FILTER=<ID=FP_cutoff,Description="Analysis-wide false-positive (FP) count cutoff: {FP_cutoff}">\n' \
# f'##FILTER=<ID=min_DP,Description="Read depth (coverage) cutoff (min allowed): {min_DP}">\n' \
# f'##FILTER=<ID=min_AC,Description="Allele count cutoff (min reads allowed to support allele): {min_AC}">\n' \
# f'##FILTER=<ID=min_AF,Description="Minor allele frequency cutoff (min allowed): {min_AF}">\n' \
# f'##FILTER=<ID=max_AF,Description="Major allele frequency cutoff (min allowed): {max_AF}">\n'
new_metadata_lines = f'##FILTER=<ID=error_rate,Description="VCFgenie, sequencing error rate for single nucleotide substitutions per base sequenced: {error_rate}">\n' \
f'##FILTER=<ID=p_cutoff,Description="VCFgenie, p-value cutoff: {p_cutoff}">\n' \
f'##FILTER=<ID=min_DP,Description="VCFgenie, read depth (coverage) cutoff (min allowed): {min_DP}">\n' \
f'##FILTER=<ID=min_AC,Description="VCFgenie, allele count cutoff (min reads allowed to support allele): {min_AC}">\n' \
f'##FILTER=<ID=min_AF,Description="VCFgenie, minor allele frequency cutoff (min allowed): {min_AF}">\n' \
f'##FILTER=<ID=max_AF,Description="VCFgenie, major allele frequency cutoff (min allowed): {max_AF}">\n'
for new_FILTER_line in new_FILTER_lines:
new_metadata_lines += f'{new_FILTER_line}\n'
new_metadata_lines += '##INFO=<ID=DECISION,Number=.,Type=String,Description="VCFgenie, decision(s) made to yield STATUS of PASS and/or FAIL">\n' \
'##INFO=<ID=STATUS,Number=A,Type=String,Description="VCFgenie, whether the ALT allele(s) is PASS or FAIL">\n' \
'##INFO=<ID=MULTIALLELIC,Number=0,Type=Flag,Description="VCFgenie, site is multiallelic (more than one ALT allele)">\n' \
'##INFO=<ID=REF_FAIL,Number=0,Type=Flag,Description="VCFgenie, REF allele failed to meet the criteria">\n' \
f'##INFO=<ID={PVR_key_new},Number=1,Type=Float,Description="VCFgenie, p-value for REF allele">\n'\
f'##INFO=<ID={PVA_key_new},Number=A,Type=Float,Description="VCFgenie, p-value for ALT allele(s)">\n' \
f'##FORMAT=<ID={PVR_key_new},Number=1,Type=Float,Description="VCFgenie, p-value for REF allele">\n' \
f'##FORMAT=<ID={PVA_key_new},Number=A,Type=Float,Description="VCFgenie, p-value for ALT allele(s)">'
# Others may be written if encountered just below for AC_key_new, AF_key_new, or DP_key_new
# TODO: here, add p_value
# Questions:
# 1) does p-value get written to BOTH format AND INFO?
# 2) should it be REF,ALTs or separate for REF and then ALTs?
# Keep track of numbers of records/samples/pass
this_record_n = 0
# this_variant_n = 0
# this_variant_pass_n = 0
this_ALLELE_n = 0
this_ALLELE_pass_n = 0
this_ALLELE_fail_n = 0
this_REF_n = 0
this_REF_pass_n = 0
this_REF_fail_n = 0
this_ALT_n = 0
this_ALT_pass_n = 0
this_ALT_fail_n = 0
this_MAJOR_n = 0
this_MAJOR_pass_n = 0
this_MAJOR_fail_n = 0
this_MINOR_n = 0
this_MINOR_pass_n = 0
this_MINOR_fail_n = 0
sample_n += 1
this_sample = ''
INFO_metadata_dd: Dict[str, dict] = {}
FORMAT_metadata_dd: Dict[str, dict] = {}
seen_AC_key_new_FLAG = False
new_AC_FORMAT_line = f'##FORMAT=<ID={AC_key_new},Number=A,Type=Integer,Description="VCFgenie, ALT allele count after processing">'
seen_AF_key_new_FLAG = False
new_AF_FORMAT_line = f'##FORMAT=<ID={AF_key_new},Number=A,Type=Float,Description="VCFgenie, ALT allele frequency after processing">'
for line in this_VCF_fh:
# chomp newline
line = line.rstrip()
if line.startswith("##"):
# SAVE metadata in the INFO_metadata_dd or FORMAT_metadata_dd
# r'^##(\w+)=<ID=(\w+),Number=(\w+),Type=(\w+),Description="(.+)"$'
# r'^##( 0 )=<ID=( 1 ),Number=( 2 ),Type=( 3 ),Description="( 4)"$'
if match := re_VCF_INFO_metadata.match(line):
if line.startswith("##INFO="):
# print(f'INFO_line={line}')
grps = match.groups() # 0-based for groups BUT 1-based if using .sub - hate Python
ID = grps[1]
INFO_metadata_dd[ID] = {}
INFO_metadata_dd[ID]['Number'] = grps[2]
INFO_metadata_dd[ID]['Type'] = grps[3]
INFO_metadata_dd[ID]['Description'] = grps[4]
this_out_fh.write(line + "\n")
elif line.startswith("##FORMAT="):
# print(f'FORMAT_line={line}')
grps = match.groups() # 0 is whole match; groups are 1-based
ID = grps[1]
FORMAT_metadata_dd[ID] = {}
FORMAT_metadata_dd[ID]['Number'] = grps[2]
FORMAT_metadata_dd[ID]['Type'] = grps[3]
FORMAT_metadata_dd[ID]['Description'] = grps[4]
if ID == AC_key_new:
print(f'### WARNING: FORMAT key="{ID}" already exists and will be OVERWRITTEN')
# WRITE *new* definition
this_out_fh.write(f'{new_AC_FORMAT_line}\n')
elif ID == AF_key_new:
print(f'### WARNING: FORMAT key="{ID}" already exists and will be OVERWRITTEN')
# WRITE *new* definition
this_out_fh.write(f'{new_AF_FORMAT_line}\n')
else:
# WRITE existing metadata headers
this_out_fh.write(line + "\n")
else:
# WRITE existing metadata headers
this_out_fh.write(line + "\n")
elif line.startswith("#"):
# WRITE new AC, AF FORMAT metadata if they were not already seen and replaced
if not seen_AC_key_new_FLAG:
this_out_fh.write(f'{new_AC_FORMAT_line}\n')
if not seen_AF_key_new_FLAG:
this_out_fh.write(f'{new_AF_FORMAT_line}\n')
# WRITE new metadata
this_out_fh.write(new_metadata_lines + "\n")
# PRINT gathered metadata information
# print('INFO_metadata_dd:')
# pprint(INFO_metadata_dd)
# print(f'FORMAT_metadata_dd:')
# pprint(FORMAT_metadata_dd)
# TEST that the keys given at top of script are present in FORMAT
needed_keys_set = {AC_key, AF_key, DP_key}
FORMAT_keys_set = set(FORMAT_metadata_dd.keys())
if not needed_keys_set.issubset(FORMAT_keys_set):
sys.exit(f'\n### ERROR: one or more keys="{needed_keys_set}" are not present in FORMAT="{FORMAT_keys_set}"')
# WRITE header line
this_out_fh.write(line + "\n")
# RECORD sample name
this_sample = line.split('\t')[9]
print(f'Analyzing sample={this_sample}')
else:
this_record_n += 1
total_record_n += 1 # all records examined, no matter how they're categorized
line_list = line.split("\t")
if len(line_list) != 10:
sys.exit('\n# ERROR: each VCF file must contain the called SNPs for exactly 1\n'
'# pooled sequencing sample, and therefore must contain exactly 10 columns:\n'
'\t(1) CHROM\n'
'\t(2) POS\n'
'\t(3) ID\n'
'\t(4) REF\n'
'\t(5) ALT\n'
'\t(6) QUAL\n'
'\t(7) FILTER\n'
'\t(8) INFO\n'
'\t(9) FORMAT [OPTIONAL]\n'
'\t(10) <sample_name>\n')
# -------------------------------------------------------------
# EXTRACT DATA for this SNP
this_CHROM = line_list[0]
this_POS = int(line_list[1])
this_ID = line_list[2]
this_REF = str(line_list[3])
this_ALT_rec = str(line_list[4])
this_QUAL = line_list[5]
this_FILTER = line_list[6] # formatted as a semicolon-separated list, so not possible to CSV by variant
this_INFO_rec = line_list[7]
this_FORMAT_rec = line_list[8]
this_sample_rec = line_list[9]
# -------------------------------------------------------------
# PROCESS DATA
# ==> ALT <==
# this_ALT_rec_comma_n = this_ALT_rec.count(',')
this_ALT_list = this_ALT_rec.split(',')
# this_ALT_n = len(this_ALT_list) # ALREADY USE this_ALT_n for another purpose