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!