bcftools mpileup issue calling indels

[tdunham19]

Hello,

I have been recently working to develop a pipeline that performs iterative alignments, variant calling, and consensus generation for virus sequences. I have been using bcftools mpileup to perform the variant calling and have noticed insertions/deletions that shouldn't be there.

I downloaded publicly available sequence data and the consensus sequence generated for that data. I utilized the consensus sequence as the reference. The final consensus sequence has insertions and deletions that shouldn't be there and I presume this is due to settings within mpileup. I am using the following code within a nextflow module of my pipeline:

bcftools mpileup -f ${reference} -Ou -o ${meta.id}.vcf ${input}

Would anybody have any suggestions on other methods or ways that I can edit this code that will help me accurately call indels?

Thank you!
Tillie

[jkbonfield]

There are newer indel callers available. Please explore bcftools mpileup --indels-cns for an alternative consensus-oriented caller. We can't guarantee every correct call in the default tool is found by the new one, but overall the benchmarks are highly positive. Note this also also automatically enabled if you specify machine specific profiles. See bcftools mpileup -X list for details. The default caller is tuned around Illumina as that's what the original method was optimised for, but it can be substantial benefits to specifying the instrument type when using another technology.

Edit: for faster turnaround, you may also wish to explore samtools consensus. It's obviously not as robust at indel realignment and variant calling as bcftools, but it may be sufficient for your needs and it goes from BAM/CRAM to new consensus in a single step.

[tdunham19]

@jkbonfield Thank you for the information, this has been incredibly helpful! Would you still recommend bcftools mpileup --indels-cns for viral genomes which are not diploid? The bcftools manual page specifically says that this has not been tested on non-diploid genomes.

For some more background on my desired outcome, I am currently in need of robust variant and indel calling for both Illumina and Nanopore platforms. Additionally, I have already implemented bcftools consensus following bcftools mpileup within my workflow.

[jkbonfield]

I cannot make any recommendations. While it's been designed with diploid genomes in mind, I cannot say whether it is still better than the default caller or not on non-diploid data. All I can suggest is that you try both and use bcftools isec to insect the two call sets and eyeball a few cases yourself to see what the differences are and which you think is more accurate on your data.

Note the diploid nature of it doesn't mean it explicitly won't work in other scenarios, rather it's just the number of paths it considers.

Edit: The method works by identifying a candidate indel to evaluate. There may be many alt alleles present and it will evaluate all of them (not just 2). However it produces a pair of consensus sequences confirming each candidate allele. This is the "diploid" part of it. It then realigns all the sequence fragments to both consensus sequences per candidate indel so it can assign reads to indels in order to get AD counts and to work out which candidate indels actually have data confirming them and which are just caused by misalignment artifacts. If you have a triploid say, then it may create poor alignments if neither of the two diploid consensuses are a close match, however it's still better than the original method which didn't generate a pair of consensus sequences at all (just 1). So while not perfect, it would expect it to still be an improvement. However as I say you need to manually verify this yourself.

[tdunham19]

@jkbonfield Thank you so much for the advice and the information on how the diploid part of this works. This has been very helpful and I appreciate your time!

[jkbonfield]

Do you have any feedback on whether the alternative indel callers are working better for you?

[tdunham19]

@jkbonfield I have not had the opportunity to implement them within my pipeline yet but if/when I do I will let you know if they helped solve my issue! Thank you again!