Edit distance seems higher than it should #23

mcolpus · 2024-03-25T10:43:02Z

mcolpus
Mar 25, 2024

I'm not quite sure what the edit distance outputs mean, as they seem to suggest more edits than needed, particularly for indels.

I ran vcfdist with -d flag with two vcfs which each had a single snp:

vcf 1:

NC_000962.3	1977	.	A	G	25.50	PASS	P	GT:GQ:DP:AD:AF	1:25:78:0,75:0.9615

vcf2:

NC_000962.3	1977	.	A	G	36.28	PASS	F	GT:GQ:DP:AD:AF	1:36:78:0,75:0.9615

But then distance summary claims to have 2 indel edits

VAR_TYPE	MIN_QUAL	EDIT_DIST	DISTINCT_EDITS	ED_QSCORE	DE_QSCORE	ALN_QSCORE
ALL	4	2	2	1.760913	1.760913	1.760913
ALL	4	2	2	1.760913	1.760913	1.760913
ALL	61	3	3	0.000000	0.000000	0.000000
SNP	4	0	0	100.000000	100.000000	0.000000
SNP	4	0	0	100.000000	100.000000	0.000000
SNP	61	1	1	0.000000	0.000000	0.000000
INDEL	4	2	2	0.000000	0.000000	0.000000
INDEL	4	2	2	0.000000	0.000000	0.000000
INDEL	61	2	2	0.000000	0.000000	0.000000

And the edits file shows:

CONTIG	START	HAP	TYPE	SIZE	SUPERCLUSTER	MIN_QUAL	MAX_QUAL
NC_000962.3	1977	0	DEL	1	0	0	26
NC_000962.3	1976	0	DEL	1	0	26	62
NC_000962.3	1977	0	SNP	1	0	26	62
NC_000962.3	1977	1	DEL	1	0	0	62

surely we should only expect a distance for QUAL between 26 and 36?

I thought it might be a Genotype problem. But changing the GT to 1/1 doesn't remove the issue, but edits changes to:

CONTIG	START	HAP	TYPE	SIZE	SUPERCLUSTER	MIN_QUAL	MAX_QUAL
NC_000962.3	1977	0	DEL	1	0	0	26
NC_000962.3	1976	0	DEL	1	0	26	62
NC_000962.3	1977	0	SNP	1	0	26	62
NC_000962.3	1977	1	DEL	1	0	0	26
NC_000962.3	1976	1	DEL	1	0	26	62
NC_000962.3	1977	1	SNP	1	0	26	62

Answered by TimD1

Mar 25, 2024

This was also a bug, caused by an off-by-one error (the truth string included an additional reference base, causing spurious INDEL edits); it is fixed in v2.5.1.

Please note that the quality threshold is only used to filter query variants.

View full answer

TimD1 · 2024-03-25T16:31:15Z

TimD1
Mar 25, 2024
Maintainer

This was also a bug, caused by an off-by-one error (the truth string included an additional reference base, causing spurious INDEL edits); it is fixed in v2.5.1.

Please note that the quality threshold is only used to filter query variants.

0 replies

mcolpus · 2024-03-26T14:35:41Z

mcolpus
Mar 26, 2024
Author

Thanks, that works now. But it now seems to have a problem with taking exponentially longer to run with larger vcfs.
For example trying to run it on a subset of my vcfs with about 100 lines in the truth and alt, takes a minute to run (see output below), whilst increasing to about 140 lines then takes 5min to run.

[INFO  vcfdist 14:40:34] Command: 'vcfdist small_alt.vcf small_truth.vcf ../../h37rv_20231215.fa -p small_test/ --min-qual 4 -d -v 2 -t 12 -sv 20 -rt -rq'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [0/8] Loading reference FASTA '../../h37rv_20231215.fa'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [Q 0/8] Parsing QUERY VCF 'small_alt.vcf'
[WARN  vcfdist 14:40:34] 'PS' tag not defined in QUERY VCF header, assuming one phase set per contig
[INFO  vcfdist 14:40:34] 1 variants of low quality (<4) in QUERY VCF, skipped
[INFO  vcfdist 14:40:34]   Genotypes:
[INFO  vcfdist 14:40:34]       1: 108
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Variant types:
[INFO  vcfdist 14:40:34]     Haplotype 1
[INFO  vcfdist 14:40:34]       REF: 0
[INFO  vcfdist 14:40:34]       SNP: 99
[INFO  vcfdist 14:40:34]       INS: 4
[INFO  vcfdist 14:40:34]       DEL: 5
[INFO  vcfdist 14:40:34]       CPX: 0
[INFO  vcfdist 14:40:34]     Haplotype 2
[INFO  vcfdist 14:40:34]       REF: 0
[INFO  vcfdist 14:40:34]       SNP: 0
[INFO  vcfdist 14:40:34]       INS: 0
[INFO  vcfdist 14:40:34]       DEL: 0
[INFO  vcfdist 14:40:34]       CPX: 0
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Contigs:
[INFO  vcfdist 14:40:34]     [ 0] NC_000962.3: 108 | 0 variants
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   QUERY VCF overview:
[INFO  vcfdist 14:40:34]     TOTAL: 109
[INFO  vcfdist 14:40:34]     KEPT : 108
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [T 0/8] Parsing TRUTH VCF 'small_truth.vcf'
[WARN  vcfdist 14:40:34] 'PS' tag not defined in TRUTH VCF header, assuming one phase set per contig
[INFO  vcfdist 14:40:34]   Genotypes:
[INFO  vcfdist 14:40:34]       1: 106
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Variant types:
[INFO  vcfdist 14:40:34]     Haplotype 1
[INFO  vcfdist 14:40:34]       REF: 0
[INFO  vcfdist 14:40:34]       SNP: 96
[INFO  vcfdist 14:40:34]       INS: 6
[INFO  vcfdist 14:40:34]       DEL: 4
[INFO  vcfdist 14:40:34]       CPX: 0
[INFO  vcfdist 14:40:34]     Haplotype 2
[INFO  vcfdist 14:40:34]       REF: 0
[INFO  vcfdist 14:40:34]       SNP: 0
[INFO  vcfdist 14:40:34]       INS: 0
[INFO  vcfdist 14:40:34]       DEL: 0
[INFO  vcfdist 14:40:34]       CPX: 0
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Contigs:
[INFO  vcfdist 14:40:34]     [ 0] NC_000962.3: 106 | 0 variants
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   TRUTH VCF overview:
[INFO  vcfdist 14:40:34]     TOTAL: 106
[INFO  vcfdist 14:40:34]     KEPT : 106
[INFO  vcfdist 14:40:34]   Writing original query VCF to './small_test/orig-query.vcf'
[INFO  vcfdist 14:40:34]   Writing original truth VCF to './small_test/orig-truth.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Checking contigs:
[INFO  vcfdist 14:40:34]     All contig checks passed!
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [Q 1/8] Wavefront clustering QUERY VCF 'small_alt.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [Q 2/8] Realigning QUERY VCF 'small_alt.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [Q 3/8] Wavefront reclustering QUERY VCF 'small_alt.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [T 1/8] Wavefront clustering TRUTH VCF 'small_truth.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [T 2/8] Realigning TRUTH VCF 'small_truth.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [T 3/8] Wavefront reclustering TRUTH VCF 'small_truth.vcf'
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [4/8] Superclustering TRUTH and QUERY variants
[INFO  vcfdist 14:40:34]            Total superclusters: 100
[INFO  vcfdist 14:40:34]   Largest supercluster (bases): 28
[INFO  vcfdist 14:40:34]   Largest supercluster  (vars): 6
[INFO  vcfdist 14:40:34]   Average supercluster (bases): 3.790
[INFO  vcfdist 14:40:34]   Average supercluster  (vars): 2.140
[INFO  vcfdist 14:40:34]               QUERY phase sets: 1
[INFO  vcfdist 14:40:34]         QUERY phase block NG50: 4411532
[INFO  vcfdist 14:40:34]               TRUTH phase sets: 1
[INFO  vcfdist 14:40:34]         TRUTH phase block NG50: 4411532
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34]   Sorting superclusters by size
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [5/8] Calculating precision and recall
[INFO  vcfdist 14:40:34]   Superclusters using   0.000 to   5.333 GB RAM each ( 12 threads):      100
[INFO  vcfdist 14:40:34]   Superclusters using   5.333 to  10.667 GB RAM each (  6 threads):        0
[INFO  vcfdist 14:40:34]   Superclusters using  10.667 to  21.333 GB RAM each (  3 threads):        0
[INFO  vcfdist 14:40:34]   Superclusters using  21.333 to  64.000 GB RAM each (  1 threads):        0
[INFO  vcfdist 14:40:34]     done with precision-recall
[INFO  vcfdist 14:40:34]  
[INFO  vcfdist 14:40:34] [6/8] Calculating edit distance metrics
[INFO  vcfdist 14:40:34]   Contigs:
[INFO  vcfdist 14:41:27]     [ 0] NC_000962.3: 131516
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Total edit distance: 131516
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Writing distance results to './small_test/distance.tsv'
[INFO  vcfdist 14:41:27]   Writing distance summary to './small_test/distance-summary.tsv'
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] ALIGNMENT DISTANCE SUMMARY
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       EDIT_DIST       DISTINCT_EDITS  ED_QSCORE       DE_QSCORE       ALN_QSCORE
[INFO  vcfdist 14:41:27] ALL    NONE Q >= 4     131512          83              0.001024        2.239367        0.004714
[INFO  vcfdist 14:41:27] ALL    BEST Q >= 10    131511          82              0.001057        2.292010        0.004813
[INFO  vcfdist 14:41:27] ALL    REF  Q >= 61    131543          139             0.000000        0.000000        0.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       EDIT_DIST       DISTINCT_EDITS  ED_QSCORE       DE_QSCORE
[INFO  vcfdist 14:41:27] SNP    NONE Q >= 4     5               5               10.413927       10.413927
[INFO  vcfdist 14:41:27] SNP    BEST Q >= 10    4               4               11.383027       11.383027
[INFO  vcfdist 14:41:27] SNP    REF  Q >= 61    55              55              0.000000        0.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       EDIT_DIST       DISTINCT_EDITS  ED_QSCORE       DE_QSCORE
[INFO  vcfdist 14:41:27] INS    NONE Q >= 4     131452          76              0.000000        0.330003
[INFO  vcfdist 14:41:27] INS    BEST Q >= 4     131452          76              0.000000        0.330003
[INFO  vcfdist 14:41:27] INS    REF  Q >= 61    131447          82              0.000000        0.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       EDIT_DIST       DISTINCT_EDITS  ED_QSCORE       DE_QSCORE
[INFO  vcfdist 14:41:27] DEL    NONE Q >= 4     55              2               0.000000        0.000000
[INFO  vcfdist 14:41:27] DEL    BEST Q >= 20    37              1               0.445821        3.010300
[INFO  vcfdist 14:41:27] DEL    REF  Q >= 61    41              2               0.000000        0.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       EDIT_DIST       DISTINCT_EDITS  ED_QSCORE       DE_QSCORE
[INFO  vcfdist 14:41:27] INDEL  NONE Q >= 4     131507          78              0.000000        0.321847
[INFO  vcfdist 14:41:27] INDEL  BEST Q >= 4     131507          78              0.000000        0.321847
[INFO  vcfdist 14:41:27] INDEL  REF  Q >= 61    131488          84              0.000000        0.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Writing edit results to './small_test/edits.tsv'
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] [7/8] Phasing superclusters
[INFO  vcfdist 14:41:27]   Contigs:
[INFO  vcfdist 14:41:27]     [ 0] NC_000962.3: 0 switch errors, 0 flip errors, 1 phase blocks
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]              Total  phase blocks: 1
[INFO  vcfdist 14:41:27]              Total switch errors: 0
[INFO  vcfdist 14:41:27]              Total   flip errors: 0
[INFO  vcfdist 14:41:27]   Supercluster switch error rate: 0.000000%
[INFO  vcfdist 14:41:27]     Supercluster flip error rate: 0.000000%
[INFO  vcfdist 14:41:27]                Phase block  NG50: 0
[INFO  vcfdist 14:41:27]   (switch)     Phase block NGC50: 0
[INFO  vcfdist 14:41:27]   (switchflip) Phase block NGC50: 0
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Writing phasing summary to './small_test/phasing-summary.tsv'
[INFO  vcfdist 14:41:27]   Writing switchflip results to './small_test/switchflips.tsv'
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] [8/8] Writing results
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Writing precision-recall results to './small_test/precision-recall.tsv'
[INFO  vcfdist 14:41:27]   Writing precision-recall summary to './small_test/precision-recall-summary.tsv'
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] PRECISION-RECALL SUMMARY
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] TYPE   THRESHOLD       TRUTH_TP        QUERY_TP        TRUTH_FN        QUERY_FP        PREC            RECALL          F1_SCORE        F1_QSCORE
[INFO  vcfdist 14:41:27] SNP    NONE Q >= 4     97              97              0               3               0.970000        1.000000        0.984772        18.173460
[INFO  vcfdist 14:41:27] SNP    BEST Q >= 10    97              97              0               1               0.989796        1.000000        0.994872        22.900347
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] INDEL  NONE Q >= 4     7               7               1               1               0.875000        0.875000        0.875000        9.030900
[INFO  vcfdist 14:41:27] INDEL  BEST Q >= 4     7               7               1               1               0.875000        0.875000        0.875000        9.030900
[INFO  vcfdist 14:41:27]  
[WARN  vcfdist 14:41:27] No QUERY SV variants pass all filters.
[INFO  vcfdist 14:41:27] SV     NONE Q >= 4     0               0               1               0               1.000000        0.000000        0.000000        0.000000
[WARN  vcfdist 14:41:27] No QUERY SV variants pass all filters.
[WARN  vcfdist 14:41:27] No TRUTH SV variants pass all filters.
[INFO  vcfdist 14:41:27] SV     BEST Q >= 0     0               0               0               0               1.000000        1.000000        1.000000        100.000000
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] ALL    NONE Q >= 4     104             104             2               4               0.962963        0.981132        0.971963        15.522618
[INFO  vcfdist 14:41:27] ALL    BEST Q >= 10    104             104             2               2               0.981132        0.981132        0.981132        17.242762
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27]   Writing phasing results to './small_test/phase-blocks.tsv'
[INFO  vcfdist 14:41:27]   Writing superclustering results to './small_test/superclusters.tsv'
[INFO  vcfdist 14:41:27]   Writing query variant results to './small_test/query.tsv'
[INFO  vcfdist 14:41:27]   Writing truth variant results to './small_test/truth.tsv'
[INFO  vcfdist 14:41:27]   Writing GA4GH-compatible summary VCF to './small_test/summary.vcf'
[INFO  vcfdist 14:41:27]  
[INFO  vcfdist 14:41:27] Timers:
[INFO  vcfdist 14:41:27]   [0] reading         :    0.014s
[INFO  vcfdist 14:41:27]   [1] clustering      :    0.005s
[INFO  vcfdist 14:41:27]   [2] realigning      :    0.001s
[INFO  vcfdist 14:41:27]   [3] reclustering    :    0.004s
[INFO  vcfdist 14:41:27]   [4] superclustering :    0.000s
[INFO  vcfdist 14:41:27]   [5] precision/recall:    0.002s
[INFO  vcfdist 14:41:27]   [6] edit distance   :   53.394s
[INFO  vcfdist 14:41:27]   [7] phasing         :    0.000s
[INFO  vcfdist 14:41:27]   [8] writing         :    0.004s
[INFO  vcfdist 14:41:27]   [9] total           :   53.424s

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Edit distance seems higher than it should #23

{{title}}

Replies: 2 comments

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Select a reply

Edit distance seems higher than it should #23

mcolpus Mar 25, 2024

Replies: 2 comments

TimD1 Mar 25, 2024 Maintainer

mcolpus Mar 26, 2024 Author

mcolpus
Mar 25, 2024

TimD1
Mar 25, 2024
Maintainer

mcolpus
Mar 26, 2024
Author