Have run into an interesting issue using samtools consensus with very large amounts of input data. The general command we use for our consensus generation is:

samtools consensus --threads 8 \ -f fasta -m simple -a -q -c 0.95 -d 600 --show-del yes --show-ins yes \ -o out.fasta input.bam

On a recent consensus generation of data from a large sequencing run we found our sequence was only 40% correct. This same library had been run on lower throughput machines and returned 100% identity leading us to question the accuracy of consensus generation. After some experimentation we discovered if we removed -q we no longer had the issue.

To further determine why this was happening we subsampled our bam file to various average target depths. The original input had ~8,300,000x coverage on average across all sites. When we subsampled to 4,000,000 average across sties our consensus was correct. However, if we increased the coverage to 8,000,000 we began to see failures.

We do not currently know why this may be occurring but are curious if it may be an overflow occurring somewhere in the code?