[lpn]

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?

How do I find variants in the pileup file? (sorry I am new to the SNP business).
Don't have SNP array or gDNA sequencing data for these samples unfortunately.

[vv85]

this is one of the lines in your pileup file:
chr2 218208 G 30 ....................,......... +(+++)++++!#+++!+)+++!*****(#(

the 5th column with the dots and commas is the one of interest. Those are the "reads" that map to that position. What you are looking for is letters instead of dots and commas, those are "mismatches" to the reference and in theory your SNPs (of course , you have to have plenty of mismatches to the reference for it to be considered a variant).

[chadn737]

30-40 million total (i.e. 15-20 million pairs) or 30-40 million pairs? So you are running this just on chr2? How many reads are in your input bam file? Are you sure you have enough coverage?

[lpn]

this is one of the lines in your pileup file:
chr2 218208 G 30 ....................,......... +(+++)++++!#+++!+)+++!*****(#(

the 5th column with the dots and commas is the one of interest. Those are the "reads" that map to that position. What you are looking for is letters instead of dots and commas, those are "mismatches" to the reference and in theory your SNPs (of course , you have to have plenty of mismatches to the reference for it to be considered a variant).

Ok, thanks, will look for letters there.

[lpn]

30-40 million total (i.e. 15-20 million pairs) or 30-40 million pairs? So you are running this just on chr2? How many reads are in your input bam file? Are you sure you have enough coverage?

30-40M pairs. I was running this on the whole transcriptome, but it gave me segmentation fault (shared node on cluster, probably not enough resources), so I decided to run it per chromosome for several chromosomes.
Typically I would have 2/3 mapped reads or so. The coverage should be enough, some genes have coverage of 100s and there should be at least one SNP in such genes. It is highly unlikely that the SNPs are in low coverage areas. Again -- this is RNAseq.

[swbarnes2]

this is one of the lines in your pileup file:
chr2 218208 G 30 ....................,......... +(+++)++++!#+++!+)+++!*****(#(

the 5th column with the dots and commas is the one of interest. Those are the "reads" that map to that position. What you are looking for is letters instead of dots and commas, those are "mismatches" to the reference and in theory your SNPs (of course , you have to have plenty of mismatches to the reference for it to be considered a variant).

Those reads are rather poor quality. Maybe double-check the orignal fastq, to see that that the quality really is that bad.

[kstamm]

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

Similar problem here.

I have tophat generated BAM files and samtools mpileup makes an empty vcf, just headers.

I think the problem is that the reference genome chromosome names are not matching. When I use a different reference genome, mpileup finds no matching reads and happily makes the empty VCF. When I use the reference I used with tophat, the bowtie NCBI37 has chromosome names like
@SQ SN:gi|224384756|gb|CM000675.1| LN:115169878

and after much digging it seems that samtools mpileup fails to process a chr name with a | or a . character. They say the solution is to convert chr names down to something simple.

I don't have a good way to do that in a BAM file so we're not fixed yet.

Does that sound like the same problem?