This is a test version of Biostars. For the public version, visit https://www.biostars.org.
samtools mpileup VCF output

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?

samtools
• 12,253 views • link
updated by Kevin Blighe Kevin Blighe • written by igor igor

2 answers

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
• 0 views • link
written by Kevin Blighe Kevin Blighe

Pipe just redirecting output of mpileup to input for bcftools.

Hope this helps.

Cheers

• 1 views • link
written by Govardhan Anande Govardhan Anande

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

• 1 views • link
written by igor igor

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?

• 1 views • link
written by Shane McCarthy Shane McCarthy

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 :-)

• 1 views • link
written by iraun iraun

Log in to answer this question.

More posts like this