How can I extract only insertions from a VCF file?

Question

I'm looking to subset a standard VCF file to generate one which only includes insertions (i.e. not indels).

I can get part of the way there with:

bcftools view -v indels <vcf> | awk '{if(length($4) == 1) print}'

However this wouldn't catch an insertion that was part of a multi-allelic record with an insertion and deletion, where the reference length would be greater than 1 bp. Potentially then one way to go is a chain of decomposition, normalisation and then the same reference length filtering — surely there's a better way in one of the many VCF manipulation utils?

Answer 1

A one-liner:

zcat my.vcf.gz |
  perl -ane '$x=0;for $y (split(",",$F[4])){$x=1 if length($y)>length($F[3])}print if /^#/||$x'

or equivalently

zcat my.vcf.gz |
  perl -ane '$x=0;map{$x=1 if length>length($F[3])}split(",",$F[4]);print if /^#/||$x'

For simple VCF operations, I generally recommend to write a script. This may be faster than those using heavy libraries. With a script, you only parse fields you care about; most libraries unnecessarily parse every field.

---

On a related note, I recommend not to decompose multi-allelic sites unless necessary. Decomposing is tricky, makes VCF harder to parse and to understand and may be a potential source of errors. Here is an example:

#CHROM POS ID REF  ALT               QUAL FILTER  INFO  FORMAT  S1  S2  S3  S4
11     101 .  GCGT G,GCGA,GTGA,CCGT  199  PASS    .     GT      0/1 1/2 2/3 2/4

vt decompose+normalize produces the following VCF:

#CHROM POS ID REF  ALT QUAL FILTER INFO FORMAT  S1   S2   S3   S4
11     101 .  GCGT G   199  PASS   .    GT      0/1  1/.  ./.  ./.
11     101 .  G    C   199  PASS   .    GT      0/.  ./.  ./.  ./1
11     102 .  CGT  TGA 199  PASS   .    GT      0/.  ./.  ./1  ./.
11     104 .  T    A   199  PASS   .    GT      0/.  ./1  1/.  1/.

In theory, you can reconstruct the original VCF from this output. However, it is very challenging for a program to do that. When you compute allele frequency line-by-line, this VCF will give you wrong results. bcftools norm -m- replaces "." with "0". You can get a correct ALT allele frequency from the bcftools output, but a wrong REF allele frequency. Furthermore, vt is also imperfect in that "CGT=>TGA" is not decomposed.

My preferred output is:

#CHROM POS    ID     REF    ALT    QUAL   FILTER INFO   FORMAT S1     S2     S3     S4
11     101    .      GCGT   G,<M>  0      .      .      GT     0/1    1/2    2/2    2/2
11     101    .      G      C,<M>  0      .      .      GT     0/2    2/0    0/0    0/1
11     102    .      C      T,<M>  0      .      .      GT     0/2    2/0    0/1    0/0
11     104    .      T      A,<M>  0      .      .      GT     0/2    2/1    1/1    1/0

Here we use a symbolic allele <M> to represent "another ALT allele". You can calculate the allele frequency by looking at one line, and won't confuse other ALT alleles with REF. bgt can produce such a VCF indirectly. However, it discards all INFO, so is not a practical solution, either.

In summary, it is very difficult to decompose multi-allelic sites. When you get decomposing wrong, your downstream analyses may be inaccurate. Decomposition should be used with caution.

Answer 2

Just to highlight that all the steps can be done within bcftools capabilities, and since I can't just comment on @blmoore 's answer:

bcftools view --types indels <vcf> |
  bcftools norm -m - |
  bcftools filter --include 'strlen(REF)<strlen(ALT)' |
  bcftools view -H

more bult-in functions for 'bcftools filter' here

Answer 3

One method is to decompose multi-allelic records so that they're represented as one-allele, one-record using vt or similar:

bcftools view -v indels <vcf> |
  vt decompose - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

Some decomposed alleles will come with excess reference padding. To left-shift and trim these, add a normalize step (NB matching these back to your input VCF becomes non-trivial):

bcftools view -v indels <vcf> |
  vt decompose - |
  vt normalize -r <reference.fasta - |
  awk '{if(length($4)==1) print}'

edit: As gringer suggests, this can also be done without vt:

bcftools view -Ou -v indels <vcf> |
  bcftools norm -Ou -Nm - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

To also include complex alleles, use view -V snps (!snvs) instead of -v indels

Answer 4

Using vcffilterjs

• get the length of the REF;
• loop over the ALT, ignore the symbolic
• accept the variant if it's an insertion , eq: len(ALT)>len(REF)

.

java -jar dist/vcffilterjs.jar -e 'function accept(vc){var a=vc.getAlleles();var lenRef=a.get(0).length();for(i=1;i<a.size();++i) {var alt=a.get(i);if(alt.isSymbolic()) continue;var lenAlt=alt.length(); if(lenRef<lenAlt) return true; } return false; }accept(variant);' input.vcf

.

Answer 5

If you're doing set operations, you could use vcf2bed:

$ vcf2bed --insertions < in.vcf > out.bed

Based on the VCF v4.2 specification, --snvs, --insertions, and --deletions are options available to filter input. In each case, the length of the reference and alternate alleles is used to determine which type of variant is being handled.