I am interested in generating consensus sequences from a bwa alignment. I generated the consensus sequences as follows:
samtools mpileup -uf reference.fasta sorted.bam -d 8000 | bcftools call -mv -Oz -o sorted.vcf.gz tabix sorted.vcf.gz cat reference.fasta | bcftools consensus sorted.vcf.gz > consensus.fasta
However, I noticed that even for sequences in my reference where none of the reads matched, I got a potential "consensus" sequence back. I assume this is because there were no entries in the VCF file that indicated SNPs, which could mean a perfect alignment.
However, this is misleading in my case, because if nothing mapped to a sequence in the reference set then I don't want to be confused with a consensus sequence. On the other hand, if a lot of reads mapped perfectly to a sequence in the reference set, then of course I'd like to see the correct consensus sequence.
How would you recommend I modify this approach?
Is it better to filter the bam somehow to only keep things that have a certain mapping depth? Or should I filter the VCF?
Thanks in advance!