samtools mpileup consensus calls A instead of a gap

[StefanP]

Hello,

I'm using samtools with the following command to call the consensus of my mapping:

samtools mpileup -uf Ref SEQ | bcftools view -cg - | vcfutils.pl vcf2fq > SEQ.fq

The problem is that I have sequences that have 1 and some have 2 copies of a short repeat - the reference has 2. If I use the reference with 2 copies I end up with 2 copies in the consensus since it sometimes maps the reads to the first and sometimes to the second copy of the repeat. So I exchanged one copy in the reference with N's., but end up with A's instead of the gap in the consensus.

Example:

Ref AANNNTATTATTA
Read1 AA---TATTATTA
Read2 AA---TATTATTA

it calls AAAAATATTATTA instead of AATATTATTA.

I also get the following error if I use the reference with N's instead of the second copy (that I don't get if I use the reference with 2 copies):

Unknown command "vcf2fq".


I've had the same problem with the previous pileup function!

Any idea how to change that?

Cheers,
Stefan

[maubp]

I wonder if this is an issue with treating N in the reference as A, or down to the fact that the reads have an A before the gap? Do you have some other examples with different bases used instead of the N's?

[StefanP]

Hi,

Replacing "N" by "-" seams to help. It looks like it's really reading N's as A's.

Thanks!

Stefan

[maubp]

So if you replaced the "N" with "C" or "T" or "G" do you get those letters instead of "A"? That would pretty much confirm the nature of this bug / design limitation, and give a good clue about where in the code to start looking.

Reminds me of this tweet from Vince Buffalo,
@vsbuffalo: Any time you only handle A, T, C, or G, the members of IUPAC kill an innocent kitten. We've all done it. #bioinfo

https://twitter.com/vsbuffalo/status/258852545634648064

[StefanP]

Hi,

Yes, I checked and it calls the reference in case the position is not covered by any read.

I should probably report that as a bug.

Thanks!

Stefan

[StefanP]

It seams they changed it a little bit since the pileup function, which also interprets "-" as "A", but didn't fully resolve the problem.

Cheers,
Stefan

[StefanP]

It seams to be a problem only when positions are "covered" by gaps in the read. Since it's calling N in the consensus if the position is not covered by a single read.

Stefan