Update to the depth filtering for mpileup.

We can now count the depth anywhere along the read rather than only at
the left end.  It's not buffering reads up so it doesn't know what's
about to come, so for any individual column it'll still favour reads
that end near a column than start near a column, but by counting at
the right hand end we avoid the boom & bust approach to before where
the coverage could plumet to zero immediately following a large spike.

With iter->depth_pos = 1.0 the coverage ends up being a minimum of
maxcnt (assuming the true coverage is higher), rather than a maximum
of maxcnt.  So we go over-depth, but avoid the drop outs.

Author: jkbonfield

htslib/sam.h

@@ -1650,6 +1650,17 @@ typedef struct __bam_mplp_t *bam_mplp_t;
     HTSLIB_EXPORT
     void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt);
 
+    /**
+     *  bam_plp_set_cntpos() - indicate which end of reach to apply maxcnt.
+     *  @plp:    The bam_plp_t initialised using bam_plp_init.
+     *  @end:    A fraction along the read from 0.0 (start) to 1.0 (end).
+     *           Defaults to 0.0 which makes maxcnt a true upper-bound value.
+     *           Using 1.0 makes maxcnt a lower-bound, with upper-bound
+     *           varying based on length and coordinate distribution.
+     */
+    HTSLIB_EXPORT
+    void bam_plp_set_cntpos(bam_plp_t iter, double end);
+
     HTSLIB_EXPORT
     void bam_plp_reset(bam_plp_t iter);
 
@@ -1712,6 +1723,18 @@ typedef struct __bam_mplp_t *bam_mplp_t;
     HTSLIB_EXPORT
     void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt);
 
+    /**
+     *  bam_mplp_set_cntpos() - indicate which end of reach to apply maxcnt.
+     *  @mplp:   The bam_mplp_t initialised using bam_mplp_init.
+     *  @end:    A fraction along the read from 0.0 (start) to 1.0 (end).
+     *           Defaults to 0.0 which makes maxcnt a true upper-bound value.
+     *           Using 1.0 makes maxcnt a lower-bound, with upper-bound
+     *           varying based on length and coordinate distribution.
+     */
+    HTSLIB_EXPORT
+    void bam_mplp_set_cntpos(bam_mplp_t iter, double end);
+
+
     HTSLIB_EXPORT
     int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp);
 

sam.c

@@ -4139,6 +4139,8 @@ int bam_plp_insertion(const bam_pileup1_t *p, kstring_t *ins, int *del_len) {
 KHASH_MAP_INIT_STR(olap_hash, lbnode_t *)
 typedef khash_t(olap_hash) olap_hash_t;
 
+#define MAX_AHEAD 4096
+
 struct __bam_plp_t {
     mempool_t *mp;
     lbnode_t *head, *tail;
@@ -4153,6 +4155,10 @@ struct __bam_plp_t {
     void *data;
     olap_hash_t *overlaps;
 
+    uint32_t depth[MAX_AHEAD];
+    uint64_t last_depth_pos;
+    double depth_pos_fract;
+
     // For notification of creation and destruction events
     // and associated client-owned pointer.
     int (*plp_construct)(void *data, const bam1_t *b, bam_pileup_cd *cd);
@@ -4172,6 +4178,9 @@ bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data)
         iter->data = data;
         iter->b = bam_init1();
     }
+    memset(iter->depth, 0, MAX_AHEAD * sizeof(*iter->depth));
+    iter->depth_pos_fract = 0;  // fraction between 0 = left and 1 = right
+    iter->last_depth_pos = 0;
     return iter;
 }
 
@@ -4494,21 +4503,60 @@ const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_
 
 int bam_plp_push(bam_plp_t iter, const bam1_t *b)
 {
+    uint64_t endpos = 0;
+
     if (iter->error) return -1;
     if (b) {
         if (b->core.tid < 0) { overlap_remove(iter, b); return 0; }
         // Skip only unmapped reads here, any additional filtering must be done in iter->func
         if (b->core.flag & BAM_FUNMAP) { overlap_remove(iter, b); return 0; }
-        if (iter->tid == b->core.tid && iter->pos == b->core.pos && iter->mp->cnt > iter->maxcnt)
-        {
+        if (iter->depth_pos_fract == 0
+            && iter->tid == b->core.tid
+            && iter->pos == b->core.pos
+            && iter->mp->cnt > iter->maxcnt) {
+            // Removal from left-end depth; mp->cnt
             overlap_remove(iter, b);
             return 0;
+        } else if (iter->depth_pos_fract > 0
+                   && iter->tid == b->core.tid
+                   && iter->pos == b->core.pos) {
+            // Removal from depth somewhere else along read
+            endpos = bam_endpos(b);
+            uint32_t len = (uint32_t)(endpos - b->core.pos - 1);
+            len = len * iter->depth_pos_fract + 0.4999;
+            if (b->core.pos + len < iter->last_depth_pos
+                && iter->depth[(b->core.pos+len)%MAX_AHEAD] > iter->maxcnt) {
+                // NB this means read longer than MAX_AHEAD (after downsizing
+                // by depth_pos fraction) will be retained as by definition
+                // it'll have zero coverage that far out.
+                overlap_remove(iter, b);
+                return 0;
+            }
+
+            // Zero newly observed iter depth elemtns.
+            uint64_t end_clipped = endpos, p;
+            if (end_clipped > iter->pos + MAX_AHEAD)
+                end_clipped = iter->pos + MAX_AHEAD;
+            if (iter->last_depth_pos < end_clipped) {
+                //iter->last_depth_pos = end_clipped;
+                if (iter->last_depth_pos < end_clipped-MAX_AHEAD)
+                    iter->last_depth_pos = end_clipped-MAX_AHEAD;
+                for (p = iter->last_depth_pos; p < end_clipped; p++)
+                    iter->depth[p % MAX_AHEAD] = 0;
+                iter->last_depth_pos = end_clipped;
+            }
+
+            // Increment depth
+            for (p = b->core.pos; p < end_clipped; p++)
+                iter->depth[p % MAX_AHEAD]++;
         }
         if (bam_copy1(&iter->tail->b, b) == NULL)
             return -1;
         iter->tail->b.id = iter->id++;
         iter->tail->beg = b->core.pos;
-        iter->tail->end = bam_endpos(b);
+        if (!endpos)
+            endpos = bam_endpos(b);
+        iter->tail->end = endpos;
         iter->tail->s = g_cstate_null; iter->tail->s.end = iter->tail->end - 1; // initialize cstate_t
         if (b->core.tid < iter->max_tid) {
             hts_log_error("The input is not sorted (chromosomes out of order)");
@@ -4602,6 +4650,11 @@ void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt)
     iter->maxcnt = maxcnt;
 }
 
+void bam_plp_set_cntpos(bam_plp_t iter, double end)
+{
+    iter->depth_pos_fract = end;
+}
+
 /************************
  *** Mpileup iterator ***
  ************************/
@@ -4651,6 +4704,13 @@ void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt)
         iter->iter[i]->maxcnt = maxcnt;
 }
 
+void bam_mplp_set_cntpos(bam_mplp_t iter, double end)
+{
+    int i;
+    for (i = 0; i < iter->n; ++i)
+        iter->iter[i]->depth_pos_fract = end;
+}
+
 void bam_mplp_destroy(bam_mplp_t iter)
 {
     int i;