problem with samtools and SNP calling

[lpn]

I am trying to do some SNP calling with samtools following this:

http://samtools.sourceforge.net/mpileup.shtml

The issue is that I don't get any SNPs or indels. To check what happens, I modify this:

samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf

and split it into 2 parts:

samtools mpileup -uf hg19.fa file.bam > file.tmp
and
bcftools view -bvcg file.tmp > var.raw.bcf

There is a lot of data in file.tmp, but var.raw.bcf seems too small (3k). Then, as suggested, I do:

bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

the resulting var.flt.vcf is just headers.

I guess I am missing something obvious but can't figure it now. Any ideas what I am doing wrong?

Edit: or any ideas for finding SNPs based on tophat .bam output?

[peromhc]

try letting samtools use anomalous pairs with the -A flag.

[lpn]

try letting samtools use anomalous pairs with the -A flag.

Same result

Any ideas for finding SNPs based on tophat .bam output?

[swbarnes2]

Okay, the obvious things to try... stop filtering with vcfutils. Just look at the raw vcf.

Second, make a pileup file, and look at it, and see if it looks right. Is coverage right? Does mpileup have the right reference sequences?

[lpn]

Okay, the obvious things to try... stop filtering with vcfutils. Just look at the raw vcf.

I looked at the var.raw.bcf; seems too small (3k). The raw vcf can't be large either.

Second, make a pileup file, and look at it, and see if it looks right. Is coverage right? Does mpileup have the right reference sequences?

How do I look at it? seems like a binary file, but of a good size, several GB.

seems like something happens in:
bcftools view -bvcg file.tmp > var.raw.bcf

where file.tmp is the pileup file

[vv85]

Use the following to command to create a readable vcf file:

bcftools view -vcg file.tmp > var.raw.vcf

Since you are using the -v option it will only output variants. You could also try creating a vcf only using -cg and check for coverage on all the other regions.

Edit: Of course , that last option will generate a huge file.

[swbarnes2]

mpileup pileup files are perfectly human readable.

mpileup -f reference.fasta > pileup.txt

[lpn]

Use the following to command to create a readable vcf file:

bcftools view -vcg file.tmp > var.raw.vcf

Same results

Since you are using the -v option it will only output variants. You could also try creating a vcf only using -cg and check for coverage on all the other regions.

Edit: Of course , that last option will generate a huge file.

Yes, I stopped it; it produces something.

[vv85]

Totally true. I just started from the first step you mentioned. The file.tmp file was generated with -u so it's not human readable.

[lpn]

mpileup pileup files are perfectly human readable.

yes, converted that to readable; seems like a pileup file.

[vv85]

do as swbarnes2 suggested and check coverage and reference names on the pileup file.
also what aligment parameters did you use?

[lpn]

do as swbarnes2 suggested and check coverage and reference names on the pileup file.
also what aligment parameters did you use?

In the pileup file I get things like:
chr2 218208 G 30 ....................,......... +(+++)++++!#+++!+)+++!*****(#(

In the reference the names are like chr2, so they are the same.

The tophat parameters are:
-g 1 -r 150 --coverage-search --no-novel-indels --microexon-search --solexa1.3-quals

(I know, these are not optimal, but should not affect SNPs, I think)


and the from the tophat output, one chromosome is chosen as:

samtools index accepted_hits.bam
samtools view -b accepted_hits.bam chr2 > accepted_hits.chr2.bam

samtools sort accepted_hits.chr2.bam file.bam

[chadn737]

This is RNA-seq data (you are using Tophat)? How many reads?

[lpn]

This is RNA-seq data (you are using Tophat)? How many reads?

Yes, RNAseq. GAII: 30-40M reads, PE 2x100.
There should be at least some SNPs detected with such data.

[vv85]

another program for variant search is GATK but I don't think that is the problem. By the way, a quick inspection of the pileup file shows variants? are there any known variants that you could directly check in the pileup file?