I am trying to generate a consensus sequence from a BAM file that was generated by mapping reads to a reference FASTA containing multiple sequences.

Usually, I generate consensus sequences from BAM files using samtools and bcftools:

samtools mpileup -vf reference.fasta filename.sorted.bam | bcftools call -m -O z - > filename.vcf.gz

bcftools index filename.vcf.gz

bcftools consensus -f reference.fasta filename.vcf.gz > filename.consensus.fasta

However, this doesn't work when the reference FASTA contains more than 1 sequence.

Is there a way to generate a consensus sequence from a multi-reference BAM file?

It would also be awesome if the method filled in gaps in coverage with Ns or Xs in the generated consensus sequence but this isn't essential.

Thank you for your help!


1 Answer 1


This post on biostars has several suggestions.

THe one that seems most straightforward, starting where you are, is to use vcftools:

cat ref.fa | vcf-consensus file.vcf.gz > out.fa

Based on that I think that you should also be able to solve this problem using alfred and GATK.

If all of these suffer the same problem of needing a single sequence (which I very much doubt that they do), you could split the bam / ref files into multiple files, run whatever consensus workflow, and then combine them again. But that doesn't seem like it should be necessary.

See also this post and this post.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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