differences in coverage BEDtools - samtools mpileup

KNS · 05-17-2011, 01:40 AM

Hi,
I have been using samtools mpileup and BEDtools coverageBed under default settings to count coverages at SNP sites identified with samtools mpileup. Both programs used the same paired BAM files produced by the Bioscope reference mapper as input. Samtools mpileup should report coverages under the DP and (quality checked) DP4 values in the vcf-fileoutput, while coverageBed reports all features in the file covering identified sites. I noted that BEDtools counts are much higher than those reported by samtools.

What am I getting wrong? Or if not, based on what criteria are reads counted? Any suggestions much appreciated!

samtools 0.1.16 (r963:234):

samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

coverageBed (v2.12.0):
samtools view -b /path/to/folder/A_paired_sorted.bam | ./coverageBed -abam /path/to/folder/A_paired_sorted.bam -b /path/to/folder/SNPs_zerobased.bed > /path/to/folder/counts_A_s.txt

fileformat SNPs_zerobased.bed file(includes only sites which interest me):
scaffold_1 7354 7355

example vcf line (4 separate samples sequenced):
scaffold_1 7355 . A T 257 . DP=15;AF1=0.8585;CI95=0.625,0.875;DP4=0,1,0,12;MQ=34;FQ=-7.6;PV4=1,0.34,0.058,0.41 GT:PL

P:SP:GQ 1/1:36,6,0:2:0:11 1/1:82,9,0:3:32668:13 0/1:81,0,28:4:8249856:24 1/1:101,12,0:4:0:16

separate counts in coverageBeds for each sample
scaffold_1 7354 7355 4 1 1 1.0000000
scaffold_1 7354 7355 5 1 1 1.0000000
scaffold_1 7354 7355 6 1 1 1.0000000
scaffold_1 7354 7355 6 1 1 1.0000000

counts from samtools mpileup: DP=15 (DP4=13)
counts from coverageBed: 21

KNS · 05-17-2011, 08:39 PM

can anyone help?

aslihan · 12-01-2011, 10:41 AM

Did you find a solution to your problem? I have also similar questions.

Thanks

KNS · 12-01-2011, 10:55 PM

Alishan,
two ways to get at counts etc.
1) Use one of the old samtools versions - 0.1.12a seemed ok.
2) Even better: Use GATK

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
samtools mpileup and bcftools command lines to compare individual differences?	aslihan	Bioinformatics	1	11-17-2011 03:20 PM
exome sequencing - big coverage differences	fadista	Bioinformatics	1	10-20-2011 08:35 AM
exome sequencing - big coverage differences	fadista	Illumina/Solexa	0	10-17-2011 02:44 AM
samtools pileup vs mpileup for computing per-base coverage	SamH	Bioinformatics	5	08-10-2011 03:14 PM
coverage depth, bedtools and RNA-Seq experiment	PFS	Bioinformatics	2	03-18-2011 08:27 AM

05-17-2011, 01:40 AM	#1
KNS Member Location: CH Join Date: Aug 2010 Posts: 17	differences in coverage BEDtools - samtools mpileup Hi, I have been using samtools mpileup and BEDtools coverageBed under default settings to count coverages at SNP sites identified with samtools mpileup. Both programs used the same paired BAM files produced by the Bioscope reference mapper as input. Samtools mpileup should report coverages under the DP and (quality checked) DP4 values in the vcf-fileoutput, while coverageBed reports all features in the file covering identified sites. I noted that BEDtools counts are much higher than those reported by samtools. What am I getting wrong? Or if not, based on what criteria are reads counted? Any suggestions much appreciated! samtools 0.1.16 (r963:234): samtools mpileup -uf ref.fa aln1.bam aln2.bam \| bcftools view -bvcg - > var.raw.bcf bcftools view var.raw.bcf \| vcfutils.pl varFilter -D100 > var.flt.vcf coverageBed (v2.12.0): samtools view -b /path/to/folder/A_paired_sorted.bam \| ./coverageBed -abam /path/to/folder/A_paired_sorted.bam -b /path/to/folder/SNPs_zerobased.bed > /path/to/folder/counts_A_s.txt fileformat SNPs_zerobased.bed file(includes only sites which interest me): scaffold_1 7354 7355 example vcf line (4 separate samples sequenced): scaffold_1 7355 . A T 257 . DP=15;AF1=0.8585;CI95=0.625,0.875;DP4=0,1,0,12;MQ=34;FQ=-7.6;PV4=1,0.34,0.058,0.41 GT:PLP:SP:GQ 1/1:36,6,0:2:0:11 1/1:82,9,0:3:32668:13 0/1:81,0,28:4:8249856:24 1/1:101,12,0:4:0:16 separate counts in coverageBeds for each sample scaffold_1 7354 7355 4 1 1 1.0000000 scaffold_1 7354 7355 5 1 1 1.0000000 scaffold_1 7354 7355 6 1 1 1.0000000 scaffold_1 7354 7355 6 1 1 1.0000000 counts from samtools mpileup: DP=15 (DP4=13) counts from coverageBed: 21

05-17-2011, 08:39 PM	#2
KNS Member Location: CH Join Date: Aug 2010 Posts: 17	can anyone help?

12-01-2011, 10:41 AM	#3
aslihan Member Location: USA Join Date: Jun 2011 Posts: 23	Did you find a solution to your problem? I have also similar questions. Thanks

12-01-2011, 10:55 PM	#4
KNS Member Location: CH Join Date: Aug 2010 Posts: 17	Alishan, two ways to get at counts etc. 1) Use one of the old samtools versions - 0.1.12a seemed ok. 2) Even better: Use GATK