Fhyujk

This is related to the question I asked here. Consider a vcf file that contains duplicate variants, but where the duplicates aren't simply the same thing in the same notation but instead one is a subset of the other. For example:

##fileformat=VCFv4.1

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1

chr12   529516  .   C   T   .   PASS    .   GT  0/1

These the two variants are actually the same. They result in exactly the same genotype. Changing AACAC to AATAC at position 529514 just means change C to T at position 529516.

Is there any tool that can detect such duplicates and remove them? I tried vcfuniq from vcflib, but that doesn't seem to recognize this as a duplicate. I think it only looks at the 1st 4 fields and only considers duplicates those variants with exactly the same values in the 1st 4 fields:

$ ./bin/vcfuniq test.vcf

##fileformat=VCFv4.1

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1

chr12   529516  .   C   T   .   PASS    .   GT  0/1

However, as explained in the linked question, EBI's vcf_validator considers this invalid. And it doesn't really make sense to have these duplicates in any case, so is there any way I can detect and remove them? Preferably an existing tool, but I am open to scripting solutions as well.

This is further complicated by cases like this one:

##fileformat=VCFv4.1

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529514  529514  AACAC   AAT,AATAC       0.00    .       .     GT    0/1

chr12   529516  529516  C       T       .       PASS    .     GT    0/1

Unfortunatley, this one won't be caught by the approach in Daniel's clever script:

$ cat test2.vcf | foo.py

##fileformat=VCFv4.1

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529514  529514  AACAC   AAT,AATAC   0.00    .   .   GT  0/1

chr12   529516  529516  C   T   .   PASS    .   GT  0/1

It turns out that bcftools can do this (tested on bcftools-1.8), if you give it the reference genome to test against:

$ bcftools norm -d none -f hg19.fa test.vcf

##fileformat=VCFv4.1

##FILTER=<ID=PASS,Description="All filters passed">

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

##bcftools_normVersion=1.8+htslib-1.8

##bcftools_normCommand=norm -d none -f hg19.fa test.vcf; Date=Wed Feb 27 16:08:44 2019

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529516  .       C       T       .       PASS    .       GT      0/1

Lines   total/split/realigned/skipped:  2/0/1/0

For the more complex case of the multi-allelic variant in the second VCF example from the question, you can run it through bcftools twice. Once using norm to left-align and split multi-allelic variants, and then again t remove the duplicates:

$ bcftools norm -m -any -NO z -O v -o - ~/test2.vcf |

    bcftools norm -d none -f hg19.fa

Lines   total/split/realigned/skipped:  2/1/0/0

##fileformat=VCFv4.1

##FILTER=<ID=PASS,Description="All filters passed">

##reference=foo

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##contig=<ID=chr12>

##bcftools_normVersion=1.8+htslib-1.8

##bcftools_normCommand=norm -m -any -NO z -O v -o - test2.vcf; Date=Wed Feb 27 18:18:32 2019

##bcftools_normCommand=norm -d none -f hg19.fa -; Date=Wed Feb 27 18:18:32 2019

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1

chr12   529516  529514  CAC     T       0       .       .       GT      0/1

chr12   529516  529514  C       T       0       .       .       GT      0/0

Lines   total/split/realigned/skipped:  3/0/2/0

I'm no expert with VCF (few can say they are!) but I have worked a lot with VCF data in the last few years, both tools to consume and produce VCF. I've never seen variants encoded in this fashion, and it seems to be non-canonical. Typically:

• Single nucleotide variants (SNVs) are encoded with a single base as the REF allele and a single base as the ALT allele.


• In the case of insertions or deletions, the shorter of the REF and ALT alleles will be a single base, the base preceding the inserted/deleted sequence. Thus the first base of the REF and ALT alleles is always the same.


• In the rarer case of two or more consecutive substitutions forming a multinucleotide variant (MNV) the REF and ALT alleles will have the same length.

Using multi-bp strings of the same length to encode SNVs is unnecessary and, as you've pointed out, problematic. This makes me think its a bug or a "feature" of the variant predictor that produced the VCF.

In this case, I'd write a small script that would check for variants where the REF and ALT alleles have the same length. If the base is the same for REF and ALT in any position, drop it, and adjust the position accordingly.

The script below will convert these funky SNVs to the canonical representation, and will also work on MNVs. Standard tools should then work to remove the duplicates.

#!/usr/bin/env python3





def canonicalize(instream):

    for line in instream:

        if not line.startswith('#'):

            values = line.split('t')

            pos = int(values[1])

            ref, alt = values[3:5]

            if len(ref) > 1 and len(ref) == len(alt):

                # How many bp to trim off the end

                for n, (r, a) in enumerate(zip(ref[::-1], alt[::-1])):

                    if r != a:

                        revoffset = -1 * n

                        break



                # How many bp to trim off the front

                for n, (r, a) in enumerate(zip(ref, alt)):

                    if r != a:

                        offset = n

                        values[1] = str(pos + offset)

                        values[3] = ref[offset:revoffset]

                        values[4] = alt[offset:revoffset]

                        break

                line = 't'.join(values)

        yield line





if __name__ == '__main__':

    import sys

    for line in canonicalize(sys.stdin):

        print(line, end='')