Work around colons in reference names (GRCh38 HLA).

This extends the callback function used by the hts iterator code to
the ref:start-end range parser, permitting the parser to validate
whether the whole thing matches a reference name.

This works around parsing conflicts when given range queries like
"HLA-DRB1*12:17", which is either the whole of "HLA-DRB1*12:17" (this
exists in GRCh38) or base 17 onwards of ""HLA-DRB1*12" (which does
not).

Note there are still some undecided questions here.  Do we want to
handle the special names used in the iterator at this level to?  Eg
"*" meaning unmapped data and "." meaning whole file?

Author: jkbonfield

hts.c

@@ -2507,15 +2507,39 @@ long long hts_parse_decimal(const char *str, char **strend, int flags)
     return (sign == '+')? n : -n;
 }
 
-const char *hts_parse_reg(const char *s, int *beg, int *end)
+/*
+ * A variant of hts_parse_reg which is reference-id aware.  It uses
+ * the iterator name2id callbacks to validate the region tokenisation works.
+ *
+ * This is necessary due to GRCh38 HLA additions which have reference names
+ * like "HLA-DRB1*12:17".
+ *
+ * On success the end of the reference is returned (colon or end of string).
+ * On failure NULL is returned, and if tid/getid are supplied *tid will be -1.
+ */
+const char *hts_parse_reg2(const char *s, int *tid, int *beg, int *end,
+                           hts_name2id_f getid, void *hdr)
 {
-    char *hyphen;
+    // FIXME: do we need to permit tid=-1 for reference "*" to indicate unmapped
+    // reads, and strictly have NULL as failure test?
+    int tid_, s_len = strlen(s); // int is sufficient given beg/end types
+    if (!tid) tid = &tid_;       // simplifies code below
+
     const char *colon = strrchr(s, ':');
     if (colon == NULL) {
         *beg = 0; *end = INT_MAX;
-        return s + strlen(s);
+        *tid = getid ? getid(hdr, s) : 0;
+        return *tid >= 0 ? s + s_len : NULL;
     }
 
+    // Has a colon, but check whole name first
+    if (getid) {
+        *beg = 0; *end = INT_MAX;
+        if ((*tid = getid(hdr, s)) >= 0)
+            return s + s_len;
+    }
+
+    char *hyphen;
     *beg = hts_parse_decimal(colon+1, &hyphen, HTS_PARSE_THOUSANDS_SEP) - 1;
     if (*beg < 0) *beg = 0;
 
@@ -2524,37 +2548,42 @@ const char *hts_parse_reg(const char *s, int *beg, int *end)
     else return NULL;
 
     if (*beg >= *end) return NULL;
+    if (getid) {
+        kstring_t ks = { 0, 0, NULL };
+        kputsn(s, colon-s, &ks); // convert to nul terminated string
+        if (!ks.s) {
+            *tid = -1;
+            return NULL;
+        }
+        *tid = getid(hdr, ks.s);
+        free(ks.s);
+    } else {
+        *tid = 0;
+    }
     return colon;
 }
 
+const char *hts_parse_reg(const char *s, int *beg, int *end)
+{
+    return hts_parse_reg2(s, NULL, beg, end, NULL, NULL);
+}
+
 hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f getid, void *hdr, hts_itr_query_func *itr_query, hts_readrec_func *readrec)
 {
     int tid, beg, end;
-    const char *q;
 
     if (strcmp(reg, ".") == 0)
         return itr_query(idx, HTS_IDX_START, 0, 0, readrec);
     else if (strcmp(reg, "*") == 0)
         return itr_query(idx, HTS_IDX_NOCOOR, 0, 0, readrec);
 
-    q = hts_parse_reg(reg, &beg, &end);
-    if (q) {
-        char tmp_a[1024], *tmp = tmp_a;
-        if (q - reg + 1 > 1024)
-            if (!(tmp = malloc(q - reg + 1)))
-                return NULL;
-        strncpy(tmp, reg, q - reg);
-        tmp[q - reg] = 0;
-        tid = getid(hdr, tmp);
-        if (tmp != tmp_a)
-            free(tmp);
-    }
-    else {
-        // not parsable as a region, but possibly a sequence named "foo:a"
-        tid = getid(hdr, reg);
+    if ((tid = getid(hdr, reg)) >= 0) {
         beg = 0; end = INT_MAX;
+        return itr_query(idx, tid, beg, end, readrec);
     }
 
+    hts_parse_reg2(reg, &tid, &beg, &end, getid, hdr);
+
     if (tid < 0) return NULL;
     return itr_query(idx, tid, beg, end, readrec);
 }

htslib/hts.h

@@ -663,13 +663,27 @@ long long hts_parse_decimal(const char *str, char **strend, int flags);
     @return  Pointer to the colon or '\0' after the reference sequence name,
              or NULL if @a str could not be parsed.
 */
+
+typedef int (*hts_name2id_f)(void*, const char*);
+typedef const char *(*hts_id2name_f)(void*, int);
 const char *hts_parse_reg(const char *str, int *beg, int *end);
 
+/// Parse a "CHR:START-END"-style region string
+/** @param str   String to be parsed
+    @param tid   Set on return (if not NULL) to be reference index (-1 if invalid)
+    @param beg   Set on return to the 0-based start of the region
+    @param end   Set on return to the 1-based end of the region
+    @param getid Function pointer.  Called if not NULL to set tid.
+    @param hdr   Caller data passed to getid.
+    @return  Pointer to the colon or '\0' after the reference sequence name,
+             or NULL if @a str could not be parsed.
+*/
+const char *hts_parse_reg2(const char *str, int *tid, int *beg, int *end,
+                           hts_name2id_f getid, void *hdr);
+
     hts_itr_t *hts_itr_query(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec);
     void hts_itr_destroy(hts_itr_t *iter);
 
-    typedef int (*hts_name2id_f)(void*, const char*);
-    typedef const char *(*hts_id2name_f)(void*, int);
     typedef hts_itr_t *hts_itr_query_func(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec);
 
     hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f getid, void *hdr, hts_itr_query_func *itr_query, hts_readrec_func *readrec);

htslib/sam.h

@@ -275,6 +275,7 @@ typedef struct {
     int bam_hdr_write(BGZF *fp, const bam_hdr_t *h) HTS_RESULT_USED;
     void bam_hdr_destroy(bam_hdr_t *h);
     int bam_name2id(bam_hdr_t *h, const char *ref);
+    const char *sam_parse_reg(bam_hdr_t *h, const char *s, int *tid, int *beg, int *end);
     bam_hdr_t* bam_hdr_dup(const bam_hdr_t *h0);
 
     bam1_t *bam_init1(void);

sam.c

@@ -285,6 +285,11 @@ int bam_name2id(bam_hdr_t *h, const char *ref)
     return k == kh_end(d)? -1 : kh_val(d, k);
 }
 
+const char *sam_parse_reg(bam_hdr_t *h, const char *s, int *tid, int *beg, int *end) {
+    // FIXME: do we need to also use cram_name2id?
+    return hts_parse_reg2(s, tid, beg, end, (hts_name2id_f)bam_name2id, h);
+}
+
 /*************************
  *** BAM alignment I/O ***
  *************************/