SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to extract a regional consensus sequence of reads from BAM file jeffseq Bioinformatics 8 01-29-2015 01:05 PM
Consensus FASTA from BAM files mixter Bioinformatics 14 05-19-2014 10:07 AM
How to extract base coverage from bam/bai-files Coru S General 9 01-27-2014 01:07 PM
Merge bam files to create a consensus? bwubb Bioinformatics 3 03-20-2012 12:40 PM
Getting pileup consensus from BAM files using Bio::DB::Sam ragowthaman Bioinformatics 2 08-03-2010 09:21 AM

Reply
 
Thread Tools
Old 03-20-2014, 09:37 AM   #1
gr8guns
Junior Member
 
Location: Delaware

Join Date: Mar 2014
Posts: 5
Default Extract genotype consensus from BAM files at specific locations

Hi all,

I'm very new to this forum, and haven't looked around much to see if this has been answered already.

But is there an easy way to extract the genotype (consensus) for a sample (across all reads) at a very specific location (not a region), given its BAM file? I know this could be very simple & possible with samtools, but my search/research hasn't been great.

Example: The base at chr1 pos 45800 for that sample from its BAM file?
And also based on certain coverage depth & quality (Like say similar to Q30 & DP3 in VCF language).

To give you the context, the problem is we have analyzed SNPs for a large set of samples, using bowtie/samtools/vcftools, reporting the SNP calls in VCF format. Obviously, this captures only variants (wrt reference) individually. But, when comparing these samples across the population to define haplotype groups, using the VCF files... looking at the overlapping SNP loci, we do not know if a missing value for a sample means the genotype is same as reference or there hasn't been any coverage at all.

Hence now, I retroactively want to fetch these missing bases for each of the samples at all the SNP loci called, for more accurate haplotype groupings.

I only have the sorted BAM file (with index) and VCF file for each sample & reference of course. Hope someone can help.

Thanks so much in advance
gr8guns is offline   Reply With Quote
Old 03-20-2014, 12:09 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

You can force samtools mpileup to give you an entry for every single point in your reference.

Also, you know that if you run all your .bams together in mpileup, it will take every point where one sample has a discrepancy, and tell you what all the samples look like at that point.
swbarnes2 is offline   Reply With Quote
Old 03-20-2014, 12:38 PM   #3
gr8guns
Junior Member
 
Location: Delaware

Join Date: Mar 2014
Posts: 5
Default

Quote:
Originally Posted by swbarnes2 View Post
You can force samtools mpileup to give you an entry for every single point in your reference.

Also, you know that if you run all your .bams together in mpileup, it will take every point where one sample has a discrepancy, and tell you what all the samples look like at that point.
Thanks for your reply swbarnes2! That sounds like exactly what I would need. The calls for a set of loci across all samples & I have the files setup exactly that way. But can you delve into a bit more detail? How can I accomplish that?

I've been looking up a lot all morning into samtools mpileup, but still couldn't coax it for my liking, although I surely know there is a solution out there! Can anyone help?

So, I've used the -r option, and can extract a specific location precisely, although it points to REFERENCE BASE, while I need the consensus from the SAMPLE READS...

Quote:
samtools mpileup -I -Bf ref.fa -r chr1:45600-45600 -Q 30 -d 100 in.bam
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-sample depth to 8000
chr1 45600 T 19 CcCCccCccCC$cCcccccc IIIIIIIIIIIIIIIIIII
Although I'm not sure if the filters work either. Thanks so much in advance guys.

And yes, I have about 100 BAM files & I'll looking up about a million loci

Last edited by gr8guns; 03-20-2014 at 12:42 PM.
gr8guns is offline   Reply With Quote
Old 03-20-2014, 12:47 PM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,814
Default

Here is a previous thread on the actual commands: http://seqanswers.com/forums/showthread.php?t=28281

Original reference for the mpileup: http://samtools.sourceforge.net/mpileup.shtml
GenoMax is offline   Reply With Quote
Old 03-20-2014, 03:49 PM   #5
gr8guns
Junior Member
 
Location: Delaware

Join Date: Mar 2014
Posts: 5
Default

Hi folks,

I really appreciate you bothering to reply & pointing to the commands, but my situation is a bit different & they don't work. Let me explain.

I need to retrieve just ONE base from the sample BAM (variant or not!) at a particular location. Not a long consensus sequence (using vcfutils.pl vcf2fq), nor the SNPs only (using vcfutils.pl varFilter -D100 > var.flt.vcf)

Extracting the consensus seq & then taking out just those loci could be an option, but would be time-consuming I fear.

For example, I have a list of loci in a file... "SNP_loci.txt"

Chr1 100017386
Chr1 100019654
Chr1 100019657
Chr1 100019756
Chr1 100020740
Chr1 100022994
Chr1 100023027
Chr1 100030383
Chr1 100032933
Chr1 100033225
...

And I have a bunch of samples-specific BAM files (along with their indexed files)...
s1.bam
s2.bam
s3.bam
s4.bam
s5.bam
s6.bam
...

And a reference file
ref.fa

So, I need to extract the calls at those locations from each of the bam files.

Alternately, using something like this..
Quote:
samtools mpileup -I -uf ref.fa s1.bam | bcftools view -bvcg - > s1.raw.bcf
bcftools view s1.raw.bcf | vcfutils.pl varFilter -D100 > s1.flt.vcf
Would give me VARIANTS only in VCF file i.e., when REF != ALT, correct?? But I basically want to know the ALT, even if its same as REF or even absent (based on quality Q>=30 & DP>=3)... hence checking back retroactively this way... for a particular set of locations.

I'm thinking on these lines, but can't put it together...
Are there other tools possibly too??

Quote:
samtools mpileup -I -Bf ref.fa -l SNP_loci.txt -Q 30 -d 100 -D -S s1.bam | bcftools view -cg - | vcfutils.pl
Thanks again for your interest & time!
gr8guns is offline   Reply With Quote
Old 03-21-2014, 11:09 AM   #6
gr8guns
Junior Member
 
Location: Delaware

Join Date: Mar 2014
Posts: 5
Red face

bump... any suggestions guys?
gr8guns is offline   Reply With Quote
Old 03-21-2014, 12:07 PM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,814
Default

Quote:
Originally Posted by gr8guns View Post
So, I need to extract the calls at those locations from each of the bam files.
If that is all you need then why not use the view option from samtools. This post explains how to do that: http://www.biostars.org/p/47945/
GenoMax is offline   Reply With Quote
Old 03-24-2014, 08:04 AM   #8
gr8guns
Junior Member
 
Location: Delaware

Join Date: Mar 2014
Posts: 5
Default

Quote:
Originally Posted by GenoMax View Post
If that is all you need then why not use the view option from samtools. This post explains how to do that: http://www.biostars.org/p/47945/
Again thanks for your response genoMax! But that wouldn't work. I need the ONE consensus genotype (major allele) across all reads for that sample... something that mpileup would help with? Not just extracting the base for all reads.

Example:
At position Chr1:45680, what is the consensus genotype for sample-1, given its BAM file?
It could be same as reference, a variant, a het or no consensus at all.
gr8guns is offline   Reply With Quote
Old 06-19-2015, 09:04 AM   #9
adumitri
Member
 
Location: Cambridge, MA

Join Date: Jan 2010
Posts: 27
Default

Hi gr8guns,

Did you figure out how to use samtools to properly address your question? If so, it would be great if you could post your solution.
adumitri 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 09:18 PM.


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