Truvari bench filtered SV counts

[jonperdomo]

Hi all,
I am learning Truvari bench to use it for benchmarking SV callers. I set some filtering parameters, but I want to determine how many were filtered for each criteria. Is this possible to determine this, such as from the debug output?

For example, for a Sniffles2 dataset I use --passonly which accounts for 2735 filtered SVs, but Truvari filters 2748, and 2749 seem to be filtered for the final statistics (call_cnt):

Parameters:

truvari bench -b $benchmark -c  $input_vcf -o "$output_dir --reference $ref_file  --pctsim 0 --refdist 1000 --sizemax 1000000 --passonly

Here are example Truvari reported counts:

[INFO] Zipped 63032 variants Counter({'base': 35662, 'comp': 27370})
[INFO] 27424 chunks of 63032 variants Counter({'base': 35662, 'comp': 24622, '__filtered': 2748})
[DEBUG] Comparing chunk defaultdict(<class 'list'>, {'__filtered': [<pysam.libcbcf.VariantRecord object at 0x7fffd958dfd0>], 'comp': [<pysam.libcbcf.VariantRecord object at 0x7fffd959afd0>]})
[DEBUG] All FP chunk -> {ret}
[INFO] Stats: {
    "base cnt": 35662,
    "call cnt": 24621,

Thank you!

[ACEnglish]

The logs do not have information about why variants were filtered. You'd need to make a custom script that replicates the filtering code around here and have your own counter.

There does appear to be a variant that goes missing since Counter['comp'] - Stats['call cnt'] == 1. However, I'd need to have the base/comp vcf files (and --includebed if that was also used) in order to trace what happened to that variant.

[jonperdomo]

Okay I can try that, thank you. I did not filter by bed file, and here are the VCFs I used:

Base:
benchmark_NA12878_TR_INCL.vcf.gz

Comp:
whole_genome.vcf.gz

[jonperdomo]

I actually used the following sizemax value for the posted result:

--sizemax 100000000 --pctsim 0 --passonly --refdist 1000

[ACEnglish]

Okay, so there's a bunch of running notes pasted below of what I looked for. But the main points are:

1. --pctsim is a very old (v3?) name for the parameter. So I'm going to assume you should update to the new v5.2 version I released this morning.
2. Without know which version you're using, I can't say for sure what variants are being dropped to account for the discrepancy in comp counts. But in the latest version, I figured out that a discrepancy was being caused by BNDs falling through the cracks. After fixing the problem, the comp counts in the log match the output counts. I'll push the changes to develop in a minute after I run the functional tests.

Update: Thinking about it, I think the 1 dropped variant was the BND with a size (somehow) above --sizefilt that causing the original discrepancy. The comp variants between --sizefilt to --sizemin are only counted if they're TP, so there's no guarantee that the logs' comp count will always equal the output's comp count. Still, for BNDs, which were only recently added, the latest changes have the desired behavior.

Begin Notes:

$ truvari bench -b benchmark_NA12878_TR_INCL.vcf.gz -c whole_genome.vcf.gz -o test1 --passonly --refdist 1000 --pctseq 0 --sizemax 1000000
2025-02-16 16:57:29,808 [INFO] 27050 chunks of 63032 variants Counter({'base': 35662, 'comp': 24632, '__filtered': 2738})
$ grep -w 'cnt\|FN\|FP\|TP' test1/summary.json | grep -v gt
    "TP-base": 18393,
    "TP-comp": 18393,
    "FP": 6226,
    "FN": 17269,
    "base cnt": 35662,
    "comp cnt": 24619,

So 24632 - 24619 = 13 variants missing

# 6226 + 18393 == 24619 match
# 17269 + 18393 == 35662 match

$ zgrep -vc "#" test1/*.vcf.gz
test1/fn.vcf.gz:17269   # matches 
test1/fp.vcf.gz:6226    # matches
test1/tp-base.vcf.gz:18393  # matches
test1/tp-comp.vcf.gz:18393  # matches

$ bcftools concat -a test1/fp.vcf.gz test1/tp-comp.vcf.gz | bcftools sort -O z -o analyzed.vcf.gz
$ zgrep -vc "#" analyzed.vcf.gz
24619 # matches

$ truvari consistency -o table.txt whole_genome.vcf.gz analyzed.vcf.gz
#
# Total 27370 calls across 2 VCFs
#
#File	NumCalls
whole_genome.vcf.gz	27370
analyzed.vcf.gz	24619
#
# Summary of consistency
#
#VCFs	Calls	Pct
2	24619	89.95%
1	2751	10.05%
#
# Breakdown of VCFs' consistency
#
#Group	Total	TotalPct	PctOfFileCalls
11	24619	89.95%	89.95% 100.00%
10	2751	10.05%	10.05% 0%

2751 variants not analyzed by truvari 

Question 1 - how many variants are on chromosomes that are not in the baseline's set of contigs

>>> import truvari
base>>> base = truvari.VariantFile("benchmark_NA12878_TR_INCL.vcf.gz")
>>> comp = truvari.VariantFile("whole_genome.vcf.gz")
>>> b_ctg = set(base.header.contigs)
>>> c_ctg = set(comp.header.contigs)
>>> len(b_ctg); len(c_ctg)
195
194
>>> c_ctg - b_ctg
set() # None
>>> b_ctg - c_ctg
{'chrEBV'}

So nothing useful here

Question 2 - get comp variants not in analyzed.vcf.gz, check their filter/size/type

import pandas as pd
import truvari

d = pd.read_csv("table.txt", sep='\t')
view = d[d['FLAG'] == 2] # unique to second file

vcf = truvari.VariantFile("whole_genome.vcf.gz")
out = truvari.VariantFile("not_analyzed.vcf", 'w', header=vcf.header)
l = list(view['ID'])

for entry in vcf:
    if entry.id in l:
        out.write(entry)
        print(entry.id, ','.join(entry.filter.keys()), entry.var_size(), entry.var_type().name, sep='\t')```

Most of them fail because of their filter is not PASS
The rest:

Above sizemax
Sniffles2.INV.C418S22	PASS	1349923	INV
Sniffles2.INV.BC6DS0	PASS	129371498	INV
Sniffles2.INV.4893S11	PASS	1116309	INV

BNDs
Sniffles2.BND.87D3S37	PASS	33	BND
Sniffles2.BND.6C4CS2	PASS	17	BND
Sniffles2.BND.6C4BS2	PASS	17	BND
Sniffles2.BND.4BA5S13	PASS	17	BND
Sniffles2.BND.8893S37	PASS	16	BND
Sniffles2.BND.787ES1	PASS	16	BND
Sniffles2.BND.761ASBF	PASS	16	BND
Sniffles2.BND.4E48S5	PASS	16	BND
Sniffles2.BND.BF66S0	PASS	15	BND
Sniffles2.BND.7579SBF	PASS	15	BND
Sniffles2.BND.6CBBS4	PASS	15	BND
Sniffles2.BND.EFS12	PASS	13	BND
Sniffles2.BND.BE1FS0	PASS	13	BND

There's our 13 missing variants

Question 3 - I confirmed that these BNDs make it through filtering, so now we have to figure out where they're going.

... I'm going to guess it has something to do with sizemin/sizefilt. We can check that by setting sizemin to 0

Yeah... 

truvari/bench.py line 388
            elif match.comp.var_size() >= self.m_params.sizemin or match.comp.is_bnd():
I added the `is_bnd` because they don't follow the size rules. If the call made it this far, it needs to be output

Rerunning truvari

2025-02-16 17:34:31,184 [INFO] 27050 chunks of 63032 variants Counter({'base': 35662, 'comp': 24632, '__filtered': 2738})
    "base cnt": 35662,
    "comp cnt": 24632,

Fixed

[jonperdomo]

Okay I was using v3.5.0 but I updated to the latest version and the issue seems resolved now thank you!

[jonperdomo]

I am still getting filtered BNDs after the update on a different dataset, attached are the bench (NA12878), comp (whole_genome.vcf.gz), and the not_analyzed.vcf containing filtered records with PASS if you have time to look into it, thank you!

truvari-bnd.zip

Sniffles2.BND.BE1FS0    PASS    13      BND
Sniffles2.BND.BF66S0    PASS    15      BND
Sniffles2.BND.787ES1    PASS    16      BND
Sniffles2.BND.6C4BS2    PASS    17      BND
Sniffles2.BND.6C4CS2    PASS    17      BND
Sniffles2.BND.6CBBS4    PASS    15      BND
Sniffles2.BND.4E48S5    PASS    16      BND
Sniffles2.BND.EFS12     PASS    13      BND
Sniffles2.BND.4BA5S13   PASS    17      BND
Sniffles2.BND.87D3S37   PASS    33      BND
Sniffles2.BND.8893S37   PASS    16      BND
Sniffles2.BND.7579SBF   PASS    15      BND
Sniffles2.BND.761ASBF   PASS    16      BND

[ACEnglish]

You sure?

unzip truvari-bnd.zip
bgzip not_analyzed.vcf
tabix not_analyzed.vcf.gz
tabix NA12878.filtered.PASS.vcf.gz

truvari bench -b NA12878.filtered.PASS.vcf.gz  -c not_analyzed.vcf.gz -o test1/ --passonly
#     "comp cnt": 14,      "FP": 14,

tabix whole_genome.vcf.gz 
truvari bench -b NA12878.filtered.PASS.vcf.gz -c whole_genome.vcf.gz -o test2/ --sizemax 100000
truvari consistency not_analyzed.vcf.gz test2/fp.vcf.gz

Result:

#
# Total 9747 calls across 2 VCFs
#
#File	NumCalls
not_analyzed.vcf.gz	14
test2/fp.vcf.gz	9746
#
# Summary of consistency
#
#VCFs	Calls	Pct
2	13	0.13%
1	9734	99.87%
#
# Breakdown of VCFs' consistency
#
#Group	Total	TotalPct	PctOfFileCalls
01	9733	99.86%	0% 99.87%
11	13	0.13%	92.86% 0.13%
10	1	0.01%	7.14% 0%

With the 1 not_analyzed.vcf.gz variant being the 129,371,498bp INV which was below my sizemax.

[jonperdomo]

Okay I have the same result now, my conda environment had selected a different version

[jonperdomo]

I also found that if the base VCF specifies contigs in the header but with no SV calls, the contig's SV calls in the comp VCF are also filtered out. Not sure if these should be included as false positives. Otherwise my results are consistent now thank you!

[ACEnglish]

I'm unable to observe this behavior.

truvari bench -b NA12878.filtered.PASS.vcf.gz -c whole_genome.vcf.gz -o full
# Chromosomes in comparison VCF which aren't in the tp-comp/fp VCFs
comm -3 <(bcftools query -f "%CHROM\n" whole_genome.vcf.gz | sort | uniq) <(bcftools concat -a full/tp-comp.vcf.gz full/fp.vcf.gz | bcftools query -f "%CHROM\n" | sort | uniq)
# Returns zero lines

# Chromosomes in baseline VCF which aren't in the tp-comp/fp VCFs (thus satisfying the condition described)
comm -3 <(bcftools query -f "%CHROM\n" NA12878.filtered.PASS.vcf.gz | sort | uniq) <(bcftools concat -a full/tp-comp.vcf.gz full/fp.vcf.gz | bcftools query -f "%CHROM\n" | sort | uniq) | wc -l
# 77 chromosomes which are inside the baseline header but without baseline calls (e.g. chrY) and have comparison VCF entries

The only thing close to this behavior would be comparison calls on contigs which aren't defined in the baseline header being ignored, which is expected. https://github.com/ACEnglish/truvari/wiki/bench#include-bed--vcf-header-contigs

[jonperdomo]

Okay that is the behavior I see, it was in other VCF files. I tested this by updating the script, and all SVs are accounted for with the size and contig filters:

def getVCFContigsWithSVCalls(vcf_file):
    vcf = truvari.VariantFile(vcf_file)
    contigs = set()
    for entry in vcf:
        contigs.add(entry.chrom)
    vcf.close()
    return contigs

def main(cons_table, vcf_file, out_file, base_file):
    d = pd.read_csv(cons_table, sep='\t')
    view = d[d['FLAG'] == 2]  # unique to second file

    basefile_contigs = getVCFContigsWithSVCalls(base_file)
    logging.info("Contigs in base file:")
    for contig in basefile_contigs:
        logging.info(contig)

    vcf = truvari.VariantFile(vcf_file)
    out = truvari.VariantFile(out_file, 'w', header=vcf.header)
    logging.basicConfig(filename='filter.log', level=logging.INFO)
    l = list(view['ID'])

    PASS_filter_counts = 0
    SIZE_filter_lower_counts = 0
    SIZE_filter_upper_counts = 0
    PASS_counts = 0
    NO_BASE_CONTIG_counts = 0
    for entry in vcf:
        if entry.id in l:
            # Continue if PASS is not the filter
            if 'PASS' not in entry.filter.keys():
                PASS_filter_counts += 1
                continue

            # Continue if abs(SVLEN) < 50 or > 100000000
            if abs(entry.var_size()) < 50:
                SIZE_filter_lower_counts += 1
                continue

            if abs(entry.var_size()) > 100000000:
                SIZE_filter_upper_counts += 1
                continue

            # Continue if the contig is not in the base file
            if entry.chrom not in basefile_contigs:
                NO_BASE_CONTIG_counts += 1
                continue

            # Update the PASS counts
            PASS_counts += 1
            out.write(entry)
            print(entry.id, ','.join(entry.filter.keys()), entry.var_size(), entry.var_type().name, sep='\t')

    print("PASS filter counts:", PASS_filter_counts)
    print("SIZE filter lower counts:", SIZE_filter_lower_counts)
    print("SIZE filter upper counts:", SIZE_filter_upper_counts)
    print("Contig not in base file counts:", NO_BASE_CONTIG_counts)
    print("non-PASS filter counts:", PASS_counts)
    print("Total entries:", PASS_filter_counts + SIZE_filter_lower_counts + SIZE_filter_upper_counts + NO_BASE_CONTIG_counts + PASS_counts)

    out.close()
    vcf.close()