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?

  • $\begingroup$ What does your VCF look like? Many tools produce VCFs that include the type of variant. $\endgroup$
    – terdon
    Jun 16, 2017 at 10:59
  • $\begingroup$ Safe to assume there's not a type annotation you can just grep out, though running an annotation tool first could be one answer. $\endgroup$
    – blmoore
    Jun 16, 2017 at 11:09
  • $\begingroup$ Oh well, things are never that easy. Still, could you please edit your question and give an example of the kind of multi-allelic record you're thinking of? Ideally, give us a minimal example of your vcf file and include a couple of entries you want to keep and a couple you want to skip and show the desired output. That way we can use your example to test our solutions and we can be sure we're giving you what you need. $\endgroup$
    – terdon
    Jun 16, 2017 at 11:14

5 Answers 5


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:

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:

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.


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

  • $\begingroup$ Very nice, I had only ever used bcftools view to uncompress and then parse the output into either the header, or the data records. $\endgroup$
    – mdperry
    Sep 22, 2018 at 22:00

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

  • 1
    $\begingroup$ bcftools norm will also split multi-allelic records with the -m - option $\endgroup$
    – gringer
    Jun 16, 2017 at 12:24

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



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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.