How to generate a consensus sequence from a multi-reference BAM file?

Question

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!

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.