Mpileup consensus Ns

[kirstyn]

Hi

I have been generating consensus sequences from a BWA-aligned BAM file using the following script:
samtools mpileup -Euf Reference.fasta input.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

In Tablet view the BAM file shows decent coverage across the genome and my reads have been quality filtered. My issue is that the consensus sequence has Ns even in positions where I can clearly see reads that map to the reference. At first I thought this was to do with a default read depth cutoff, but even accounting for this the numbers don't add up!

Does anyone know why mpileup is calling an N at certain positions?

[auratus]

16S and 23S consensus replaced with N

I have the same issue. In particular with 16S and 13S in bacteria. I have plenty of coverage in this area but I got "N" on my consensus sequence.
I used this:

samtools mpileup -u foo.bam | bcftools view -cg - | vcfutils.pl vcf2fq > foo.fq

Any answer?

[syfo]

Same here.
Ns in the consensus fasta while the region seems to be well covered. I guess it is about some quality filtering but I do not manage to get the right parameters.
Anyone with a proposition please?

Command line:

samtools mpileup -u infile.bam | bcftools view -cg - | vcfutils.pl vcf2fq | seqtk fq2fa - > outfile.fa

[scami]

May be indexing.....

Hi Guys
the N in the consensus may me due to a missing indexed reference file. Try the following:

samtools index alignmentFile.bam

and then try the mpileup command again.

Hope it helps

[syfo]

Hi, thanks for your help.

Well, in my case I precisely do not want samtools to use the reference sequence for the calling, I am looking for the consensus from the reads only.

Also, most of the bases are correctly called but a couple of positions are missing and end up in Ns although base information seems available in the alignments. I am suspecting base quality to be an issue but I don't manage to get a decent set of parameters for the command line.

Even with something like this
samtools mpileup -RAB -d10000000 -Q0 -q0 -u | bcftools view -g -p1 - | vcfutils.pl vcf2fq -d0 -Q0 | seqtk fq2fa -

some bases are discarded.

I believe vcfutils.pl is involved because if I replace the "AF1=0" tag by "AF1=1" in the VCF (after the bcftools view command) N are replaced by a proper base information.

Still, even small and "obvious" deletions in the reference sequence are not handled the way I want because the consensus includes "N" where it should be uninterrupted, but that is another issue I guess.
Note: the bam files are already sorted and indexed.