samtools mpileup VCF output

Question

[igor]

On the samtools github wiki, there is an example for variant calling via pileup:

samtools mpileup -uf ref.fa aln.bam | bcftools call -mv -Oz > calls.vcf.gz

I've seen a similar example in other places. It's piping pileup output to bcftools. This made sense before. However, newer versions of samtools (1.0+) have a --VCF flag. Why is it still necessary to pipe to bcftools? Wasn't bcftools there just to convert to VCF format?

Answer 1

[Kevin Blighe]

Apologies for a late response, but I have recently been looking at samtools mpileup. I found that the following works very well (here, I use multiple samples over which variants are called):

samtools mpileup --redo-BAQ --min-BQ 30 --per-sample-mF --output-tags DP,AD -f human_g1k_v37.fasta --BCF BAM1.bam BAM2.bam BAM3.bam BAM4.bam | bcftools call --multiallelic-caller --variants-only -Ov -Ob > Output.bcf

...then normalise:

bcftools norm -m-any Output.bcf | bcftools norm -Ov --check-ref w -f human_g1k_v37.fasta > Output.vcf

I then do further filtering using the vcfutils perl script that comes bundled with BCFtools:

bcftools view Output.vcf | vcfutils.pl varFilter -d 18 -w 1 -W 3 -a 1 -1 0.05 -2 0.05 -3 0.05 -4 0.05 -e 0.05 -p > Output.Filt.vcf

Answer 2

[Govardhan Anande]

Pipe just redirecting output of mpileup to input for bcftools.

Hope this helps.

Cheers

[igor]

Yes, that is what is happening. But why pipe to bcftools instead of just using the --VCF option and save the extra step?

[Shane McCarthy]

Because bcftools call is actually doing the variant calling. mpileup is only outputting genotype likelihoods (in VCF or BCF format), which are then used by bcftools to call variants (in VCF of BCF format). Does that make sense?

[iraun]

You could try to run samtools mpileup with --VCF argument and see how the "vcf" output looks like. I'm telling you this because some time ago I had the same doubt, and after doing the check, I answered myself :-)