Using truvari to merge short variant calls with SVs

[JosephLalli]

Hi there! I'm trying to use truvari collapse, but not as intended. While this code appears to be working, I wanted to check in and make sure that there aren't any glaring problems that I am missing.

I am trying to combine Pangenie calls with deepvariant calls from the same sample. The issue I am encountering is that sometimes, deepvariant will call several low-coverage SNPs and short indels in the middle of a Pangenie deletion call. Eg:

Pangenie:    ACT------------------TAAG
Reference:   ACTCGTCGTCAAGTACTAGGATAAG
Deepvariant: ACTC----TCAAGTTCTAG-ATAAG
                ^          ^    ^

These variants are generally of low quality and low coverage. However, Google suggests retaining all variants, even low quality ones, to maximize sensitivity. Meanwhile, I have a pangenie deletion (let's say a 100bp deletion) that is a confident call, and is too large to be called by Deepvariant.

I'd like to identify Deepvariant calls that are inside pangenie calls, and remove them where they are contradictory.

To do so, I am running this code. Are there any unintended side effects of these settings that I should be aware of?

bcftools merge -Ou -m none $SV_vcf_filtered $filtered_shortread_calls | bcftools norm -cs -m- --fasta $ref_fasta -Ou - | bcftools sort -Oz - > $merged_SV_shortread_calls  && bcftools index -f -t $merged_SV_shortread_calls

truvari collapse -i $merged_SV_shortread_calls -o $synthesized_calls -c $truvari_variant_collapse_notes -f $ref_fasta \
    --sizemin 0 \   # allows comparison of SNPs and SVs
    --passonly \    # Ignore ref calls
    --intra \       # Produce a final vcf that contains the combined variant callset for the sample.
    --pctsize 0 \   # Ignore size differences
    --typeignore \  # Ignore type (SV/Indel/SNP/) differences
    --refdist 0 \   # Do not compare variants that do not overlap.
    --pctseq 0      # Do compare variants that overlap in a completely incompatible way.

Thanks,
Joe

[ACEnglish]

Hello,

That isn't going to give you the results you're expecting. There might be ways to get truvari collapse to behave closer to what you're trying to accomplish, but even then you'll have unintended side effects. The best way to do this would be to make a new script. Truvari has a python package that can be leveraged for creating new tools (docs). So you could take advantage of that to make the coding easier.

A rough sketch of what you'd want to build is below. Note that I haven't tested this so please use with caution.

import sys
import truvari
import pysam

def variant_is_within_larger_variant(small, larges):
    """
    Check if small variant is within a large variant's span
    """
    for large in larges:
        start, end = truvari.entry_boundaries(large)
        if truvari.entry_within(small, start, end):
            return True
    return False

if __name__ == '__main__':
    vcf_fn1 = sys.argv[1]
    vcf_fn2 = sys.argv[2]

    vcf1 = pysam.VariantFile(vcf_fn1)
    vcf2 = pysam.VariantFile(vcf_fn2)

    # If you have a bed file, use regions to subset the input vcfs
    regions = truvari.build_region_tree(vcf1, vcf2, bed_filename)

    # Only needed when there are potentially overlapping bed entries
    truvari.merge_region_tree_overlaps(regions)

    # Make iterators that will only return variants in regions of interest
    # This also helps zip vcfs that have different chromosome orders
    vcf1_iter = truvari.region_filter(vcf1, regions)
    vcf2_iter = truvari.region_filter(vcf2, regions)

    # Matchers control what variants are filtered from VCFs as well as thresholds for matching
    # For this, you probably(?) won't need to worry about thresholds.
    matcher = truvari.Matcher()
    matcher.params.sizemin = 0
    matcher.params.sizefilt = 0
    matcher.params.passonly = True

    # Make the chunker which will return a list of dictionaries
    chunks = truvari.chunker(matcher, ('small', vcf1_iter), ('large', vcf2_iter))

    # Make an output VCF
    # Note that this assumes the two headers are compatible such that vcf1 and vcf2 entries
    # can be written to the same file. If they aren't (e.g. vcf2 has an INFO header not in vcf1)
    # you'll need to make a custom header
    out = pysam.VariantFile("/dev/stdout", header=vcf1.header)

    for chunk, m_id in chunks:
        # At this point, a chunk holds all small and large variants within matcher.params.chunksize
        # of one another
        for small in chunk['small']:
            # Return true if this is something that should be filtered
            if not variant_is_within_larger_variant(small, chunk['large']):
                out.write(small)
        for large in chunk['large']:
            out.write(large)

[JosephLalli]

Oh wow, that's wonderful! Thank you.

I've also got a bedtools based approach, but this might be cleaner/faster/avoid I/O time.

One last thought/question: integrating short variant calls with SV calls seems like an important problem to me, but many people I've spoken to about the problem seem to have not considered it before. Am I putting a lot of effort into something that is just not that important? Producing a consolidated, internally consistent representation of all forms of genetic variation seems vitally important for identifying causal variants in disease and QTLs, but is almost never done.

[ACEnglish]

I think its likely useful, but unexplored. Some 'cutting edge' work is still trying to figure out how to even find all classes of variants efficiently and barely stepping into joint evaluation.