Pileup shows variants - Not called as SNPs

[BertieWooster]

Hello,

I notice that with samtools pileup/mpileup, there are several positions which have variants that can be classified as SNPs. However, they fail to show up when using the subsequent SNP program in its most basic form (i.e. without any filters or flags). What can be causing these 'calls' to drop off?

Thanks,
BertieW

[gaffa]

Maybe you could post some examples?

[swbarnes2]

First guess, the quality of the alternate letters is too low. Post the line of the mpileup that you think is a putative SNP, and see if that's the case.

[hrajasim]

I have a similar situation with several SNPs. I filtered out reads that align to this SNP position to investigate the mapping quality and base qualities. All reads aligning to the position where I observe this SNP on IGV have a mapping quality of 37 which is pretty good. All of these reads match perfectly without mismatch to the reference and a majority of these reads (>40) differ at this position compared to reference. But Samtools pileup does not call it a SNP. Is it just the limitation of samtools or is there something I could try to overcome this? The SNP I observe is at position 33388713

328MB44S:8:42:10743:11542#0 0 1 33388674 37 70M * 0 0 CCTTAGGTCACACTGCAGAACAGAGTGATATTGCAAAAAACCAAATTATCAGTCTACGTGCAACCACCAT gcghhhffdhhghfhhhhhhhhhfhehfhhhghghhgcfhghghhhhghhhhgdfhhhhhhffheahhhf XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:39G30

[nilshomer]

Looks like your base qualities are in Illumina's phred+64 format (see http://en.wikipedia.org/wiki/FASTQ_format). They should in Sanger format. Since the base qualities are out of range, all bets are off.

[hrajasim]

Thanks Nils. I am redoing the workflow with phred+33 qualities now. However, I am thinking what difference should I expect? It looks like with phred+64, I was overestimating the base quality and hence I was seeing lot many SNPs and indels than expected. With the corrected base qualities, this may go down. I wonder how I can expect this fix to enable the detection of the SNP that samtools had missed earlier. Any thoughts?

[nilshomer]

Post your results when you are done.

[SWP]

BWA aln has a -I option if your fastq files are in Illumina format rather than Sanger. If this option is used are quality scores changed to phred+33 in the output?

[Papillon]

I have a similar question and I hope someone can help me with this. If you need more info, just let me know.

I aligned my 2 x 100 bp reads with BWA - I and applied realignment with GATK in a later stage and looking at the results in IGV, everything looks quite good. There are some problems though.

The bam files of two family members show an AC-deletion at chr1:165,378,729 and that one is being called. At position 165,378,733 and 165,378,735 I find a C > T variants that are supported by 62 out of 62 and 58 out of 68 reads in one patient and by 36 out of 42 and 22 out of 42 reads in the second patient. Seems like clear snp's to me and the base call quality is definitely in SANGER format and between 30 and 40.

Does samtools allow snp's near indels or could the problem be somewhere else? (or maybe just because it is Monday)

Thank you!

EDIT: I know that both the deletion and the C>T snp are actual known variants (dbSNP132), although I can imagine it to be a bit complex, because they are in a repetitive region/simple repeats.

[JimC]

Similar Problem

I have a similar problem and don't want to start a separate thread unless necessary.

Illumina HiSeq paired-end genomic DNA for Arabidopsis.
Reads aligned using BWA (aln and sampe modules).
Alignment is sorted and indexed.

I used Samtools mpileup:
samtools mpileup -uf TAIR10.fa arab1.sorted.bam | bcftools view -c - > abar1.var.raw.vcf

When I compare a particular region: Chr5:4,764,955 between the output of mpileup & bcftools with the samtools tview -- I don't see the same results. I will post here if possible:

from the vcf file output:
Chr5 4763954 . T . 45 . DP=65;AF1=4.92e-07;CI95=1.5,0;DP4=1,4,0,0;MQ=54 PL 0
Chr5 4763955 . C . 3.67 . DP=64;AF1=1;CI95=0.5,1;DP4=0,0,0,3;MQ=60 PL 40
Chr5 4763956 . C . 42 . DP=64;AF1=4.837e-07;CI95=1.5,0;DP4=1,3,0,0;MQ=60 PL 0

The tview "view" is attached. I'm not sure why, with 64 reads, ALL with a T at position 4763955, that it is not called as a SNP. I'm trying to assess why? Is it the gap nearby?

Please note that this "gap" inserted into the reference is caused by a single low-quality read (not shown in the png file). How can I get BWA or mpileup to not insert such low quality inserts ??

Many Thanks !

Jim

[Papillon]

Do'h, do I feel stupid!

I am also using SAMTools' varFilter to loose the variants which have too low coverage. But it actually does much more and this is why I am loosing these SNP's.