samtools does not report positions with zero coverage (with some
exceptions).
There is no coverage in the first region:
$st view test1.bam 2:186689185-186689205 |wc -l
0
and there are five reads in the second:
$st view test1.bam 2:179511757-179511777 |wc -l
5
If you look at the alignments, you'll see that also the second region is
skipped, the cigars are like 78M264827N23M
$st tview test1.bam $ref -p 2:179511757
petr
On Wed, 2014-10-29 at 17:07 -0400, jbr950 wrote:
> Hi,
> I am finding I can't rely on mpileup to report what is happening at
> each position (described in a bed file) in a .bam file. Some
> positions are omitted and I can't predict the reason this occurs.
> Example below:
>
>
> # First, make a small bed file:
> cat smallBed.bed
> 2 186689185 186689205
> 2 179511757 179511777
>
>
> # Then, get a subset of a larger bam file:
> ${samtools_htslib_path}/samtools view -b $bamFile
> 2:179511657-186689305 > test1.bam
>
>
>
> # Now, run mpileup:
> ${samtools_htslib_path}/samtools mpileup -B -d10000000 -l smallBed.bed
> -t DPR,INFO/DPR,DP4 -uf ref/${refFile} test1.bam |
> ${bcftools_htslib_path}/bcftools view - > result1.vcf
>
>
>
> # Check a position in first interval:
> grep 186689186 result1.vcf
>
>
> # nothing reported.
>
>
> # Check a position in the second interval:
> grep 179511758 result1.vcf
> 2 179511758 . T <X> 0 .
> DP=0;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0;DPR=0,0
> PL:DP4:DPR 0,0,0:0,0,0,0:0,0
>
>
> Locus reported with a depth of 0.
>
>
> I repeated the call to mpileup without the -B or the -d parameters to
> see if this had something to do with it, but the result was the same.
>
>
> I can repeat this experiment using a second bam file, and see that in
> the second case, both positions are reported, even though both have a
> depth of 0 just like in the first sample.
>
>
>
> ${samtools_htslib_path}/samtools view -b $bamFile2
> 2:179511657-186689305 > test2.bam
>
> ${samtools_htslib_path}/samtools mpileup -B -d10000000 -l smallBed.bed
> -t DPR,INFO/DPR,DP4 -uf ref/${refFile} test2.bam |
> ${bcftools_htslib_path}/bcftools view - > result2.vcf
>
> grep 186689186 result2.vcf
> 2 186689186 . C <X> 0 .
> DP=0;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0;DPR=0,0
> PL:DP4:DPR 0,0,0:0,0,0,0:0,0
>
>
> grep 179511758 result2.vcf
> 2 179511758 . T <X> 0 .
> DP=0;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0;DPR=0,0
> PL:DP4:DPR 0,0,0:0,0,0,0:0,0
>
>
>
>
>
>
> Because I am not comfortable drawing conclusions about positions that
> are invisible to me, it would be great to learn what is happening at
> the omitted locus, and others like it. In earlier experiments using
> Samtools 0.1.18, I found some positions of importance to us with
> plenty of coverage were omitted, and I presume this may still be
> occurring in the latest samtools. Or perhaps I'm just missing some
> key info in the documentation that would explain this phenomenon?
>
>
> I am pasting a link to the test files, in case this is helpful for
> anyone interested in troubleshooting the software (I am using GRCh37
> for the ref).
>
>
> https://drive.google.com/folderview?id=0B_o0arNq2sU0eTFYMTZUV2xzQTg&usp=sharing
>
>
>
> Thanks for any assistance.
>
>
> Best,
> Jonathan
>
>
> ------------------------------------------------------------------------------
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/samtools-help
--
The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.
------------------------------------------------------------------------------
_______________________________________________
Samtools-help mailing list
Samtools-help@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/samtools-help
