Extract genotype consensus from BAM files at specific locations

[gr8guns]

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

[swbarnes2]

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.

[gr8guns]

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...

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

[GenoMax]

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

[gr8guns]

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..

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??

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]

bump... any suggestions guys?

[GenoMax]

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/

[gr8guns]

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.

[adumitri]

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.