TopHat samtools snp calling

[JasonR]

Hey all,

I am trying to call SNPs on tophat output and am missing known SNPs with good coverage. I hoping to get advice on tuning samtools/bcftools for RNA-seq data or recommendations for alternative pipelines.

I have directional, polyA-selected, paired-end 101bp libraries sequenced on an Illumina HiSeq machine. I t used tophat (v1.1.4) to map trimmed reads and samtools (version 0.1.11 (r851)) to sort the accepted.hits.bam file and used mpileup/bcftools (I also tried pileup) to call variants. I followed the suggestions in the samtools manual (samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf), except I played with the -D parameter a bit.

Using these parameters I call around 22,000 snps. However, I miss around half of a set of known snps with genes.

As an example, a T to C snp, which I know is in my sample from Sanger on the DNA, didn't show up in this analysis. Even if I call pileup with only -c -f, I can't find this position, which doesn't make any sense to me because I was under the impression that omitting -v would provide information on all covered positions. Using tview, I know that I have 119 reads that cover the position in the bam file and 38 are ref, 79 are C and 2 are G.

This maybe as good as it gets, but I would appreciate any suggestions on altering samtools parameters or alternative methods. Thanks!

[chenjy]

Hey all,

I am trying to call SNPs on tophat output and am missing known SNPs with good coverage. I hoping to get advice on tuning samtools/bcftools for RNA-seq data or recommendations for alternative pipelines.

I have directional, polyA-selected, paired-end 101bp libraries sequenced on an Illumina HiSeq machine. I t used tophat (v1.1.4) to map trimmed reads and samtools (version 0.1.11 (r851)) to sort the accepted.hits.bam file and used mpileup/bcftools (I also tried pileup) to call variants. I followed the suggestions in the samtools manual (samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf), except I played with the -D parameter a bit.

Using these parameters I call around 22,000 snps. However, I miss around half of a set of known snps with genes.

As an example, a T to C snp, which I know is in my sample from Sanger on the DNA, didn't show up in this analysis. Even if I call pileup with only -c -f, I can't find this position, which doesn't make any sense to me because I was under the impression that omitting -v would provide information on all covered positions. Using tview, I know that I have 119 reads that cover the position in the bam file and 38 are ref, 79 are C and 2 are G.

This maybe as good as it gets, but I would appreciate any suggestions on altering samtools parameters or alternative methods. Thanks!

hi,

do you have any method bypassing this problems now?

thank you. looking forward to your reply.

[m_elena_bioinfo]

Hi!
I've the same problem: known deletion and SNP in a sample that I can't detected by variant calling in RNASeq alignment from Tophat.
Have you fix the problem?
Thanx a lot

[gennadyM]

I've just started to do the same thing - call SNPs from RNA-Seq alignments - and ran into the same issue, so I think this still might be useful:

Try using '-A' flag with samtools mpileup - this flag is well hidden in the man page, because it doesn't start from a new line, but goes after '-6' flag. But you can also see it by simply typing 'samtools mpileup' without any arguments. Also, the following post might be useful (and the first reply to it mentions -A): http://massgenomics.org/2012/03/5-th...s-mpileup.html

[swbarnes2]

Try adding -B to the pileup command.

Sometimes, the BAQ "corrections" eat real SNPs. -B will turn them off.

[gennadyM]

I tried to use either -B or -E (it's also discussed in the link I mentioned), but neither helped. the interval of ~55 bases was missing from the samtools mpileup output, while in fact there are tens of RNA-Seq reads aligned there, with good base and mapping qualities. Even those positions at the flanks of that interval also counted a small fraction of reads (1 or 2 out of many). The SNP I was looking for is a known mutation in a gene (exon), and it's clearly seen in the data (e.g., IGV browser). I was just testing if I see something "obvious" and I didn't ... until I've put '-A'. I think you're right though that -B might give additional improvements if the realignment to transcriptome is not going to work well.