Question: How to do consensus calling from bam files
gravatar for AW
3.8 years ago by
United Kingdom
AW350 wrote:


I want to extract consensus sequences for a number of bam files with the intent to align sequences across individuals and build trees.

I am following the advice here via vcf2fq: samtools mpileup -uf ref.fa aln.bam | bcftools call -c | vcf2fq > cns.fq

via bcftools consensus: samtools mpileup -uf ref.fa aln.bam | bcftools call -mv -Oz -o calls.vcf.gz tabix calls.vcf.gz cat ref.fa | bcftools consensus calls.vcf.gz > cns.fa

But I have some questions I cannot find the answer to:

  1. How does bcftools consensus work differently from vcf2fq?

  2. How does either deal with missing data? For example if a site is not covered by reads in the bam file, ideally I would like the consensus sequence to contain an N. Is this the case or is the allele in the reference sequence used?

  3. Does the consensus callers filter the SNPs? And if not, how do you do this? is this achieved by running bcftools filter after bcftools calls?


ADD COMMENTlink modified 3.8 years ago by Brice Sarver3.5k • written 3.8 years ago by AW350

I would say that inserting N might not be the right choice. If a region is not covered by reads it is probably not present - perhaps you might insert the reference sequence for that but definitely not Ns.

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by Istvan Albert ♦♦ 84k
gravatar for Brice Sarver
3.8 years ago by
Brice Sarver3.5k
United States
Brice Sarver3.5k wrote:

Some consensus generators do strange things, especially if you're using the results for downstream analysis and not, say, as a crude sample-specific reference. Especially if you're doing phylogenetic analysis, I'd recommend a slightly different approach:

Get a set of high-quality, filtered SNPs based on whatever criteria you think is good for your system. Then, use code or the GATK's FastaAlternateReferenceMaker to inject that variation back into the reference. For regions that you are uncertain about, such as positions with low coverage or low quality, you can mask them using bedtools' maskfasta.

ADD COMMENTlink written 3.8 years ago by Brice Sarver3.5k

This is sound like the a better approach - we've often tried to call consensus from mpileup outputs but that only works for the simplest of cases in general the situation is more complicated.

ADD REPLYlink written 3.8 years ago by Istvan Albert ♦♦ 84k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1534 users visited in the last hour