Many Ns in samtools mpileup

[mikheyev]

Hi,

I am trying to wrangle samtools into giving me a consensus, using reads mapped with bwa to a reference. After I generate a BAM, the results look pretty good using samtools tview:


So, there is good coverage, etc.

Now, when I try to generat a consensus, for example, by

samtools  pileup -6Aug sorted.bam |bcftools view -cg - | vcfutils.pl vcf2fq

I get a whole lot of Ns in my consensus

@gi|159883886|emb|CU329670.1|
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnn
+
!!!!@@@@@@@@@@CEEEHHHHHHHHHHHHHHNKKNKQQQWWWWWWccux{{{{{~~~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~~~~~~~~~~~~~~

I swear, I tried every combination of samtools and bcftools I could think of. Adding a reference sequence also doesn't make any difference. Any ideas would be most very much appreciated.

Sasha

[macrowave]

You didn't provide the reference sequences when calling pileup? If that's the case, all the ref bases are Ns and you probably wouldn't get the right consensus.

[mikheyev]

Unfortunately, having a reference doesn't help:

samtools  mpileup -6Augf sequence.fasta -r "gi|159883886|emb|CU329670.1|:1-200"  sorted.bam |bcftools view -cg - | vcfutils.pl vcf2fq

@gi|159883886|emb|CU329670.1|
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnn
+
!!!!@@@@@@@@@@CEEEHHHHHHHHHHHHHHNKKNKQQQWWWWWWccux{{{{{~~~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~~~~~~~~~~~~~~

[Tally]

I'm not sure about your pileup issue, but see the thread "samtools tview reference sequence" regarding the string of N's in your tview. For me, this was caused by having a colon character in the reference sequence fasta header.

[mikheyev]

I think I figured out what it was -- I had an old from a previous bam file hanging around. I re-indexed the file and everything was OK.

[cwisch88]

For me, this was caused by having a colon character in the reference sequence fasta header.

Wow, this sentence has solved a problem I have been having all morning.Thanks so much!