Pileup / extract information from BAM/SAM files

liu_xt005 · 05-01-2012, 09:00 AM

I am aiming to extract the full information at each site: read depth, # ref reads, # variant reads (for all non-reference alleles), strand information, etc. An old thread started by nilmot13 suggested "genomeCoverageBed", but the output seems to contain only simple read depth.

My problem started from an observation that SAMtools and GATK generate very different REF/ALT read depth at some site from the same BAM file, while BAMview and IGV both tend to support GATK counts at the site.

clk · 09-28-2012, 07:45 AM

Hi,
looking for an answer to the same question I came across your post... Did you ever find an answer to this question? The discrepancies I see between IGV and samtools mpileup are huge in my data, which is RNAseq data...

Thanks!

liu_xt005 · 10-01-2012, 12:56 PM

Sorry I did not find a good solution.
But my study showed that SAMtools tend to keep longer tails for indels than GATK and others. For example, SAMtools gives TAAAA:TAAA (REF:ALT), while GATK gives TA:T.

clk · 10-01-2012, 01:07 PM

Thanks for responding. I finally found a solution to this, so I'll post it here in case is useful for others.

I found out that samtools filters reads before including them in the pileup; it reads the flag field in the bam file and discards reads that
a) are not paired
b) not properly mapped
c) mate is not mapped
d) alignment is not primary
e) reads fail quality control of vendor
f) is marked as PCR duplicates.

If the filters (a) and (c) are not desired, you can use the parameter -A.
In addition, samtools performs realignment unless the parameter -B is used, and discards low quality reads unless -Q0 is used. Finally, it stops at a certain number of reads unless the -d parameter is invoqued.

I really needed a good quantification of the reads at each position, so I needed to make sure I could trust the pileup (or vcf) files generated by samtools. So I wrote a small script that parses the bam file, reading the flag field, and removes specific reads from the alignment. In this way, I was finally able to produce an alignment that gave me the exact same counts with IGV and samtools pileup.

It would be really great if all these little details were more clear in the documentation, but in the end, the filtering criteria used by samtools was adequate for my needs. Except for the "anomalous read pairs" parameter (-A), which not very appropriate for RNA-seq data.

Hope that helps somebody!

zyxue · 01-19-2015, 10:14 AM

Hi clk, where did you find the information? Did you analyze the source code for samtools mpileup? Thanks!

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
0 depth in pileup files from BAM files ?	jjlaisnoopy	Bioinformatics	3	04-22-2012 08:38 AM
Pileup from BAM or SAM	nilmot13	General	7	02-03-2012 03:54 AM
Extract perfectly mapped reads from SAM/BAM file	Graham Etherington	Bioinformatics	2	07-21-2011 08:27 AM
NEw to Chip-seq and have .bam/.sam/.bam.bai files... then what?	NGS newbie	Bioinformatics	11	05-25-2011 08:48 AM
Getting pileup consensus from BAM files using Bio::DB::Sam	ragowthaman	Bioinformatics	2	08-03-2010 09:21 AM

05-01-2012, 09:00 AM	#1
liu_xt005 Member Location: Iowa City, IA Join Date: Jun 2011 Posts: 24	Pileup / extract information from BAM/SAM files I am aiming to extract the full information at each site: read depth, # ref reads, # variant reads (for all non-reference alleles), strand information, etc. An old thread started by nilmot13 suggested "genomeCoverageBed", but the output seems to contain only simple read depth. My problem started from an observation that SAMtools and GATK generate very different REF/ALT read depth at some site from the same BAM file, while BAMview and IGV both tend to support GATK counts at the site.

09-28-2012, 07:45 AM	#2
clk Junior Member Location: Montreal Join Date: Nov 2010 Posts: 4	Hi, looking for an answer to the same question I came across your post... Did you ever find an answer to this question? The discrepancies I see between IGV and samtools mpileup are huge in my data, which is RNAseq data... Thanks!

10-01-2012, 12:56 PM	#3
liu_xt005 Member Location: Iowa City, IA Join Date: Jun 2011 Posts: 24	Sorry I did not find a good solution. But my study showed that SAMtools tend to keep longer tails for indels than GATK and others. For example, SAMtools gives TAAAA:TAAA (REF:ALT), while GATK gives TA:T.

10-01-2012, 01:07 PM	#4
clk Junior Member Location: Montreal Join Date: Nov 2010 Posts: 4	Thanks for responding. I finally found a solution to this, so I'll post it here in case is useful for others. I found out that samtools filters reads before including them in the pileup; it reads the flag field in the bam file and discards reads that a) are not paired b) not properly mapped c) mate is not mapped d) alignment is not primary e) reads fail quality control of vendor f) is marked as PCR duplicates. If the filters (a) and (c) are not desired, you can use the parameter -A. In addition, samtools performs realignment unless the parameter -B is used, and discards low quality reads unless -Q0 is used. Finally, it stops at a certain number of reads unless the -d parameter is invoqued. I really needed a good quantification of the reads at each position, so I needed to make sure I could trust the pileup (or vcf) files generated by samtools. So I wrote a small script that parses the bam file, reading the flag field, and removes specific reads from the alignment. In this way, I was finally able to produce an alignment that gave me the exact same counts with IGV and samtools pileup. It would be really great if all these little details were more clear in the documentation, but in the end, the filtering criteria used by samtools was adequate for my needs. Except for the "anomalous read pairs" parameter (-A), which not very appropriate for RNA-seq data. Hope that helps somebody!

01-19-2015, 10:14 AM	#5
zyxue Junior Member Location: Canada Join Date: May 2014 Posts: 3	Hi clk, where did you find the information? Did you analyze the source code for samtools mpileup? Thanks!