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

<pre><code>$ bcftools norm -d none -f hg19.fa test.vcf
##fileformat=VCFv4.1
##FILTER=&lt;ID=PASS,Description="All filters passed"&gt;
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
##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
</code></pre>

<p>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 <code>norm</code> to left-align and split multi-allelic variants, and then again to remove the duplicates:</p>

<pre><code>$ 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=&lt;ID=PASS,Description="All filters passed"&gt;
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
##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
</code></pre>

CC BY-SA 4.0
2019-03-01T15:18:18.850Z
Great information to have in the back pocket. Nice solution!
d_kennetz
CC BY-SA 4.0
2019-03-01T15:54:05.533Z
Wish I could +1 this again.
Daniel Standage
CC BY-SA 4.0
2019-02-27T16:10:34.370Z
2019-03-01T15:21:03.897Z
2019-03-01T15:21:03.897Z
terdon
<p>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:</p>

<ul>
<li>Single nucleotide variants (SNVs) are encoded with a single base as the REF allele and a single base as the ALT allele.</li>
<li>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.</li>
<li>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.</li>
</ul>

<p>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.</p>

<p>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.</p>

<p>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.</p>

<pre class="lang-py prettyprint-override"><code>#!/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) &gt; 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='')
</code></pre>

<p><strong>UPDATE</strong>: Upon further reflection, the more complicated example you listed actually makes sense. In the reference we have the sequence <code>AACAC</code>, and the alternate alleles represent two variations on this: deletion of the final two bp (in the first case), and a point mutation of the middle <code>C</code> to <code>T</code> (in both cases). Usually, only a single bp precedes the definition of an indel, so I would've encoded this complex variant as <code>ref=ACAC alt=AT,ATAC</code>.</p>

<p>So the SNV is "implied by"/"encoded in"/"redundant with" the complex variant, but it's not strictly a duplicate. I'm curious whether the VCF validator complains about these cases as well?</p>

CC BY-SA 4.0
2019-02-27T16:27:33.343Z
It is indeed atypical, but this question was prompted because I actually encountered this in the wild. I didn't generate the vcf, it was sent to me, but the header suggests it was produced by `freebayes` and the merge of two separate files. So this was probably an artefact of the merging. Unfortunately, I need to deal with VCF files that are given to me by clients, so while I can insist that they conform to the standards, this [_does_ conform](https://bioinformatics.stackexchange.com/q/7105/298) (AFAIK), so I needed a way of fixing it.
terdon
CC BY-SA 4.0
2019-02-27T18:20:17.170Z
Nice one! Unfortunately (and sorry, I should have made this clear) the actual file I had contains multi-allelic variants where only one was a dupe (see updated question). Your approach won't catch those, but bcftools does (see my answer). This works great for singles though.
terdon
CC BY-SA 4.0
2019-02-27T18:34:28.907Z
Dude, that is some gnarly VCF.
Daniel Standage
CC BY-SA 4.0
2019-02-27T18:36:27.977Z
Welcome to my world :(
terdon
CC BY-SA 4.0
2019-03-01T15:16:19.610Z
Still a cool solution @DanielStandage!
d_kennetz
CC BY-SA 4.0
2019-02-27T16:15:59.773Z
2019-02-27T20:05:10.100Z
2019-02-27T20:05:10.100Z
Daniel Standage
<p>This is related to the question I asked <a href="https://bioinformatics.stackexchange.com/q/7105/298">here</a>. 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:</p>

<pre><code>##fileformat=VCFv4.1
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
#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
</code></pre>

<p>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.</p>

<p>Is there any tool that can detect such duplicates and remove them? I tried <a href="https://github.com/vcflib/vcflib#vcfuniq" rel="nofollow noreferrer"><code>vcfuniq</code></a> from <code>vcflib</code>, 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:</p>

<pre><code>$ ./bin/vcfuniq test.vcf
##fileformat=VCFv4.1
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
#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
</code></pre>

<p>However, as explained in the <a href="https://bioinformatics.stackexchange.com/q/7105/298">linked question</a>, EBI's <code>vcf_validator</code> 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.</p>

<hr>

<p>This is further complicated by cases like this one:</p>

<pre><code>##fileformat=VCFv4.1
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
#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
</code></pre>

<p>Unfortunately, this one won't be caught by the approach in <a href="https://bioinformatics.stackexchange.com/a/7128/298">Daniel's clever script</a>:</p>

<pre><code>$ cat test2.vcf | foo.py
##fileformat=VCFv4.1
##reference=foo
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description="Genotype"&gt;
##contig=&lt;ID=chr12&gt;
#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
</code></pre>

<vcf><variation>
How can I remove (non-trivial) duplicates from a VCF file?
# How can I remove (non-trivial) duplicates from a VCF file?

(Asked by: username_0 on 2019-02-27)

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=
##contig=
#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=
##contig=
#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=
##contig=
#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

Unfortunately, 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=
##contig=
#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

**Tags:** vcf, variation

---

## 2 answers

### Answer by username_0 on 2019-02-27. Score: 3 | Accepted answer

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=
##reference=foo
##FORMAT=
##contig=
##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 to 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=
##reference=foo
##FORMAT=
##contig=
##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

> **Comments:**
> **username_1**: Great information to have in the back pocket. Nice solution!  
> **username_2**: Wish I could +1 this again.  

---

### Answer by username_2 on 2019-02-27. Score: 2

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='')

UPDATE: Upon further reflection, the more complicated example you listed actually makes sense. In the reference we have the sequence AACAC, and the alternate alleles represent two variations on this: deletion of the final two bp (in the first case), and a point mutation of the middle C to T (in both cases). Usually, only a single bp precedes the definition of an indel, so I would've encoded this complex variant as ref=ACAC alt=AT,ATAC.
So the SNV is "implied by"/"encoded in"/"redundant with" the complex variant, but it's not strictly a duplicate. I'm curious whether the VCF validator complains about these cases as well?

> **Comments:**
> **username_0**: It is indeed atypical, but this question was prompted because I actually encountered this in the wild. I didn't generate the vcf, it was sent to me, but the header suggests it was produced by `freebayes` and the merge of two separate files. So this was probably an artefact of the merging. Unfortunately, I need to deal with VCF files that are given to me by clients, so while I can insist that they conform to the standards, this [_does_ conform](https://bioinformatics.stackexchange.com/q/7105/298) (AFAIK), so I needed a way of fixing it.  
> **username_0**: Nice one! Unfortunately (and sorry, I should have made this clear) the actual file I had contains multi-allelic variants where only one was a dupe (see updated question). Your approach won't catch those, but bcftools does (see my answer). This works great for singles though.  
> **username_2**: Dude, that is some gnarly VCF.  
> **username_0**: Welcome to my world :(  
> **username_1**: Still a cool solution @DanielStandage!  

---

