SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
TopHat samtools snp calling JasonR Bioinformatics 5 05-13-2013 09:43 AM
Samtools SNP calling vidhya Bioinformatics 3 04-07-2011 07:17 AM
Samtools: SNP calling from only coding regions newbietonextgen Bioinformatics 5 12-10-2010 02:42 PM
SAMTOOLS SNP calling question harrb Bioinformatics 2 12-10-2010 07:37 AM
SAMtools and SNP calling Jan Bioinformatics 2 09-16-2010 02:01 PM

Reply
 
Thread Tools
Old 03-01-2012, 08:30 AM   #1
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default problem with samtools and SNP calling

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?

Last edited by lpn; 03-01-2012 at 09:36 AM.
lpn is offline   Reply With Quote
Old 03-01-2012, 09:09 AM   #2
peromhc
Senior Member
 
Location: Durham, NH

Join Date: Sep 2009
Posts: 108
Default

try letting samtools use anomalous pairs with the -A flag.
peromhc is offline   Reply With Quote
Old 03-01-2012, 09:35 AM   #3
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by peromhc View Post
try letting samtools use anomalous pairs with the -A flag.
Same result

Any ideas for finding SNPs based on tophat .bam output?
lpn is offline   Reply With Quote
Old 03-01-2012, 09:50 AM   #4
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

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?
swbarnes2 is offline   Reply With Quote
Old 03-01-2012, 09:54 AM   #5
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by swbarnes2 View Post
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.

Quote:
Originally Posted by swbarnes2 View Post
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

Last edited by lpn; 03-01-2012 at 09:57 AM.
lpn is offline   Reply With Quote
Old 03-01-2012, 10:03 AM   #6
vv85
Member
 
Location: Mexico

Join Date: Feb 2011
Posts: 17
Default

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.

Last edited by vv85; 03-01-2012 at 10:10 AM.
vv85 is offline   Reply With Quote
Old 03-01-2012, 10:08 AM   #7
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

mpileup pileup files are perfectly human readable.

Quote:
mpileup -f reference.fasta > pileup.txt
swbarnes2 is offline   Reply With Quote
Old 03-01-2012, 10:15 AM   #8
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by vv85 View Post
Use the following to command to create a readable vcf file:

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

Quote:
Originally Posted by vv85 View Post
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.
lpn is offline   Reply With Quote
Old 03-01-2012, 10:16 AM   #9
vv85
Member
 
Location: Mexico

Join Date: Feb 2011
Posts: 17
Default

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.
vv85 is offline   Reply With Quote
Old 03-01-2012, 10:23 AM   #10
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by swbarnes2 View Post
mpileup pileup files are perfectly human readable.
yes, converted that to readable; seems like a pileup file.
lpn is offline   Reply With Quote
Old 03-01-2012, 11:23 AM   #11
vv85
Member
 
Location: Mexico

Join Date: Feb 2011
Posts: 17
Default

do as swbarnes2 suggested and check coverage and reference names on the pileup file.
also what aligment parameters did you use?
vv85 is offline   Reply With Quote
Old 03-01-2012, 11:36 AM   #12
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by vv85 View Post
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

Last edited by lpn; 03-01-2012 at 11:41 AM.
lpn is offline   Reply With Quote
Old 03-01-2012, 11:52 AM   #13
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

This is RNA-seq data (you are using Tophat)? How many reads?
chadn737 is offline   Reply With Quote
Old 03-01-2012, 11:56 AM   #14
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by chadn737 View Post
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.
lpn is offline   Reply With Quote
Old 03-01-2012, 11:56 AM   #15
vv85
Member
 
Location: Mexico

Join Date: Feb 2011
Posts: 17
Default

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?
vv85 is offline   Reply With Quote
Old 03-01-2012, 11:59 AM   #16
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by vv85 View Post
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.
lpn is offline   Reply With Quote
Old 03-01-2012, 12:12 PM   #17
vv85
Member
 
Location: Mexico

Join Date: Feb 2011
Posts: 17
Default

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).
vv85 is offline   Reply With Quote
Old 03-01-2012, 12:19 PM   #18
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

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?
chadn737 is offline   Reply With Quote
Old 03-01-2012, 12:46 PM   #19
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by vv85 View Post
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 is offline   Reply With Quote
Old 03-01-2012, 12:49 PM   #20
lpn
Member
 
Location: west coast

Join Date: May 2011
Posts: 17
Default

Quote:
Originally Posted by chadn737 View Post
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.
lpn is offline   Reply With Quote
Reply

Tags
bam, samtools, snp, tophat

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:31 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO