7
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ What does your VCF look like? Many tools produce VCFs that include the type of variant. $\endgroup$ – terdon Jun 16 '17 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 '17 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 '17 at 11:14
4
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

.

$\endgroup$
2
$\begingroup$

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

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

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

$\endgroup$
  • $\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 '18 at 22:00
0
$\begingroup$

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.

$\endgroup$

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.