Why are there all these As in my consensus sequence: An issue with samtools pileup

[dagarfield]

Hey folks,

I've run into an interesting problem. I've aligned a bunch of 454 reads from one species against the reference genome of a closely related species. The goal is to construct, the the extent our low-coverage allows, a map of the differences between species A and species B.

I've been using the samtools pileup tool followed by the samtools.pl pileup2fq to pull out the consensus sequence as suggested by the SAMtools FAQ (http://tinyurl.com/64r9mdf).

Frustratingly, the program seems to do something odd even with high confidence indels: Even when all of the reads show a deletion, some portion of the reference sequence is transformed into As. An example is shown below.



For reasons that are totally unclear to me, the program changes the reference from GTT to AAA starting at base 172292 even though the reference GTT is the only information. Arg.

Any thoughts?


On a related note, pileup2fq skips all insertions. Does anyone know of a method to more directly transform a pileup into a consensus fasta or fastq?

Many thanks!!

[dagarfield]

Oh, I get to respond to my own post.
tview seems to always produce the As, as does pileup.
If you use mpileup, things are fine.

[swbarnes2]

I wish there were a way to get pileup2fq to deal with indels properly. I don't know of any automated way to do that either.

There is a known bug with pileup producing A's in the consensus when it shouldn't. IT's mentioned in a wikia somewhere.

[dagarfield]

Yeah, I found the reference on the SAMTOOLS FAQ. By posting in the forum over at SAMTOOLS, I learned that this has been dealt with in mpileup. There appear to be no plans to fix pileup itself any time soon.

Unfortunately, mpileup doesn't appear to produce a consensus sequence (no -c option), so that may be the end of the utility of pileup2fq in any event.