<p>Hey all :]</p>

<p>I use samtools&#39;s depth, and occasionally samtool&#39;s pileup commands to calculate coverage of my reads to the genome before binning for coverage. I&#39;m pretty sure everyone else does too ;)</p>

<p>One very common problem is that people find &quot;samtools depth&quot; and &quot;samtools mpileup&quot; don&#39;t match up - and it&#39;s commonly attributed to filtering poor quality reads, duplicates, etc (mpileup doing more filtering).</p>

<p>But there&#39;s a ton of other questions I just can&#39;t get the answers too from the sametools docs, namely:</p>

<p>Does depth/pileup count the region between a paired reads, or just the reads themselves?<br />
Does samtools depth/pileup count singletons in paired end sequencing? Does it extrapolate based on average read length to fillout the whole fragment?<br />
If a read maps to multiple locations, is it counted multiple times?</p>

<p>I&#39;m working with ChIP-Seq data, so i might have to correct for peak-shift. What do you guys think?<br />
<br />
Thanks! :)</p>

> Does depth/pileup count the region between a paired reads, or just the reads themselves?

samtools depth 'just' runs a multi-mpileup and count the number of bases under each base of the reference. smal l'qlen', deletions and base quality below 'baseQ' are ignored. Here is the core code of 'depth':

    if (!(b->core.flag&BAM_FUNMAP)) {
        if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
        else if (aux->min_len && bam_cigar2qlen(&b->core, bam1_cigar(b)) < aux->min_len) b->core.flag |= BAM_FUNMAP;

(...)

    while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { 
        if (pos < beg || pos >= end) continue; // out of range; skip
        if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue; 
        fputs(h->target_name[tid], stdout); printf("\t%d", pos+1); 
        for (i = 0; i < n; ++i) { // base level filters have to go here
            int j, m = 0;
            for (j = 0; j < n_plp[i]; ++j) {
                const bam_pileup1_t *p = plp[i] + j;
                if (p->is_del || p->is_refskip) ++m;
                else if (bam1_qual(p->b)[p->qpos] < baseQ) ++m; // low
            }
            printf("\t%d", n_plp[i] - m); // this the depth to output
        }
    }