SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools mpileup using -l or -r vplagnol Genomic Resequencing 3 05-08-2013 10:23 PM
Many Ns in samtools mpileup mikheyev Bioinformatics 5 03-22-2012 09:56 AM
samtools mpileup hv4 Bioinformatics 2 10-29-2011 01:21 PM
Why are there all these As in my consensus sequence: An issue with samtools pileup dagarfield Bioinformatics 3 02-03-2011 10:14 AM
Samtools pileup coverage issue colindaven Bioinformatics 2 06-11-2010 07:49 AM

Reply
 
Thread Tools
Old 01-08-2012, 10:25 PM   #1
ashkot
Member
 
Location: Cupertino, CA

Join Date: Nov 2011
Posts: 59
Default SAMTools mpileup issue

Hi there,
I am following the following workflow:

Sample ID from 1000 Genomes project: HG00096

First, indexing the bam file using samtool index
Second, sorting the bam file using samtools sort eg2.bam eg2.sorted
Third, generating a bcf file using samtools mpileup -uf reference_genome_file eg2.sorted.bam | bcftools view -bvcg - > eg2.raw.bcf

The reference human genome file is the one obtained from 1kGenomes, Human_37_*.fasta

Lastly, converting my bcf file to vcf using bcftools view *.bcf | perl /usr/share/samtools/vcfutils.pl varFilter -D100 > *.vcf

I am interested in knowing genotypes for rs numbered dbSNPs and other variants for which I have HG37 genomic coordinates.

When i follow the above workflow, i am not getting 0/0 genotypes and at the same time some of the genomic locations also do not show up. This is why my first question was if absence of variation should be considered a 0/0 genotype?

Thanks,
Ashwin
ashkot is offline   Reply With Quote
Old 01-09-2012, 08:30 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

'bcftools view -bvcg' will only show you the variant locations. So you won't get sites where all the samples in the .bam are homozygous reference. 'bcftools view -bcg' should do the trick.
swbarnes2 is offline   Reply With Quote
Old 01-09-2012, 08:54 AM   #3
jflowers
Member
 
Location: New York, NY

Join Date: Oct 2011
Posts: 41
Default

bcftools view -bvcg will give a variants-only vcf. However, whether this is a good way of getting genotypes is a bit unclear, but may not be appropriate (depending on your application).

In this post (http://biostar.stackexchange.com/que...quality-scores), lh3 suggests that mpileup / bcftools is not really meant for calling genotypes, but instead making inferences about the data when genotypes are unknown. In fact, lh3 suggests that calling genotypes is a "flaw" of bcftools (see also the original paper, Bioinformatics. 27(21) 2987-2993).

This seems like a subtle, but important point.
jflowers is offline   Reply With Quote
Old 01-10-2012, 11:10 AM   #4
ashkot
Member
 
Location: Cupertino, CA

Join Date: Nov 2011
Posts: 59
Default

i think this is a really good point. what i am after is to identify the genotype at all genomic locations but it seems that mplieup is really meat to find variations.

also, i deleted the -v switch and i still cannot get certain genomic locations to show up in my VCF file. should i assume the absence of that location to be homozygous normal at that location?
ashkot is offline   Reply With Quote
Old 01-10-2012, 12:31 PM   #5
jflowers
Member
 
Location: New York, NY

Join Date: Oct 2011
Posts: 41
Default

I have some of the same questions as you ashkot when it comes to getting a consensus genome, [especially in light of the "flaw" concerning genotypes inferred by bcftools mentioned by lh3 I mentioned in my previous post on this thread]. So, I'm not really sure I can be of much help.

The topic of generating a consensus has also been discussed here:
http://seqanswers.com/forums/showthread.php?t=9818

I would defer to comments made in that thread since lh3 developed the methods in question. But, what's not clear to me is whether the "flaws" with using bcftools to get genotypes in the vcf (mentioned previously in this thread) lead to bad genotypes being propagated throughout the consensus. Or does the method of calling genotypes in the consensus using vcfutils.pl vcf2fq differ from the method used by bcftools to call genotypes in the vcf. Can someone more knowledgable chime in here?

Anyhow, generating a consensus is described here:
http://samtools.sourceforge.net/mpileup.shtml

As far as missing sites in your vcf file, have you looked at the sites that are missing using samtools tview to see if there are reads spanning those sites?
jflowers is offline   Reply With Quote
Old 01-11-2012, 06:02 AM   #6
aggp11
Member
 
Location: Wisconsin

Join Date: Jun 2011
Posts: 87
Default

Quote:
Originally Posted by ashkot View Post
Hi there,

First, indexing the bam file using samtool index
Second, sorting the bam file using samtools sort eg2.bam eg2.sorted
Third, generating a bcf file using samtools mpileup -uf reference_genome_file eg2.sorted.bam | bcftools view -bvcg - > eg2.raw.bcf
Hi Ashwin,

I don't know if this would matter, but I have always sorted the bam file first and then created it's index.

The thing is, based on the steps that you have mentioned, your index file would be based on your original bam file, and since you are using the sorted bam file in your mpileup, this might be causing some issues.

Again this is just a thought, but you might want to try it once.

Thanks,
Praful
aggp11 is offline   Reply With Quote
Old 01-11-2012, 08:42 AM   #7
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Quote:
Originally Posted by aggp11 View Post
Hi Ashwin,
The thing is, based on the steps that you have mentioned, your index file would be based on your original bam file, and since you are using the sorted bam file in your mpileup, this might be causing some issues.

Again this is just a thought, but you might want to try it once.

Thanks,
Praful
You are right, but mpileup doesn't look at the index, so that's not the issue. The issue is that mpileup really isn't designed to do what the poster wants to do with it, so it's working correctly, but not in the off-label way that the poster hoped it would.

I'm not sure what tool would be more appropriate. What I think will work is making a .bed file of your genotype loci, and using intersectBed with that and your .bam. That'll filter out all the reads that don't cover your genotype loci, which I think will make the vcf smaller, and easier to manage. Then make the all-points vcf by leaving out the -v, then use intersectBed again with that .bed against the vcf. That should leave only the mpileup entries that hit your genotype loci, and since it's all points, the 0/0 entries shoud be there.
swbarnes2 is offline   Reply With Quote
Old 01-11-2012, 09:10 AM   #8
aggp11
Member
 
Location: Wisconsin

Join Date: Jun 2011
Posts: 87
Default

Quote:
Originally Posted by swbarnes2 View Post
You are right, but mpileup doesn't look at the index, so that's not the issue. The issue is that mpileup really isn't designed to do what the poster wants to do with it, so it's working correctly, but not in the off-label way that the poster hoped it would.

I'm not sure what tool would be more appropriate. What I think will work is making a .bed file of your genotype loci, and using intersectBed with that and your .bam. That'll filter out all the reads that don't cover your genotype loci, which I think will make the vcf smaller, and easier to manage. Then make the all-points vcf by leaving out the -v, then use intersectBed again with that .bed against the vcf. That should leave only the mpileup entries that hit your genotype loci, and since it's all points, the 0/0 entries shoud be there.
Thanks for the information on mpileup. And I like your idea of using intersectBed here.
aggp11 is offline   Reply With Quote
Old 01-11-2012, 10:06 AM   #9
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Oh, and I'll add if you get rid of the -v, and a locus doesn't show up, I'm pretty srue that's an indicator that you have no reads there, or at least no reads of any quality.
swbarnes2 is offline   Reply With Quote
Old 01-12-2012, 04:22 PM   #10
ashkot
Member
 
Location: Cupertino, CA

Join Date: Nov 2011
Posts: 59
Default

hi there, thanks for you input, i am going to give intersetBed a try and see what i get. I'll keep everybody informed about the outcome.

thanks again.

ashwin
ashkot is offline   Reply With Quote
Reply

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 02:54 AM.


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