Samtools mpileup - why does it call an "N" when there ARE correct alignments?

[TabeaK]

Hi all!

Happy 2013!

I am having a bit of an issue with generating consensus sequences from an Alignment BAM file.

The alignment was done with Stampy. It is a mix of single end and paired end data. Looking at the BAM in IGV shows me I have good coverage across the board and reasonable mapping qualites.

I then do the routine samtools mpileup (using -u -f -r paramters; also tried the -A); then bcftools view (- c -g ); followed by conversion to fastq and ultimately fasta.

In my consensus fasta however; I have an inordinate amount of Ns; even in positions; where I can clearly see correct alignment/mapping.

Does anyone know when mpileup is calling an N and what I might be doing wrong?

Funnily enough; I have aligned the exact same data with BWA and done the same process; even though BWA maps less reads properly; with lower mapping quality and coverage; calling the consensus only spits out Ns if there is either no coverage at a given position or there is already an N in the reference.

So, why is this not working from my stampy-generated BAM?
Argh!

Any other recommendations for tools to generate consensus sequences?

[TabeaK]

Bumpety-bump.

[TabeaK]

Ha! I think this is the problem: http://seqanswers.com/forums/showthread.php?t=17486

I also get only 'N' in the ref column...no I just need to figure out what the hell is going on.

[TabeaK]

Okay, problem solved...

This little thing:

11.2 Parsing NCBI fasta files

NCBI fasta files use a >gi|nnn|ref|xxx identifier. By default
Stampy parses this and uses only the "xxx" part. Use --noparseNCBI
to switch off this behaviour and use the full NCBI identifier.

from the stampy manual broke my neck. In the BAM file created with Stampy; the identifiers were parse. While creating the fasta index; however samtools; DOES NOT parse. Ergo the identifers did not match.

I had to do some rather elaborate messing around with the reference fasta to get the strings to match...*sigh*
All good in the end....just in case anyone runs into this problem in the future.