Major rewrite of hts_parse_region.

The old function (hts_parse_reg) has been put back, meaning the new
API can now have the getid func pointer as a mandatory requirement.

Added tests.

Tweaked sam_parse_region to have the flags parameter.  This is still
required as bcftools and samtools use different parameters (parsing
style is tool based rather than file format based).

Author: jkbonfield

Makefile

@@ -78,7 +78,8 @@ BUILT_TEST_PROGRAMS = \
 	test/test-vcf-api \
 	test/test-vcf-sweep \
 	test/test-bcf-sr \
-	test/test-bcf-translate
+	test/test-bcf-translate \
+	test/test-parse-reg
 
 BUILT_THRASH_PROGRAMS = \
 	test/thrash_threads1 \
@@ -358,6 +359,7 @@ check test: $(BUILT_PROGRAMS) $(BUILT_TEST_PROGRAMS)
 	test/fieldarith test/fieldarith.sam
 	test/hfile
 	test/test_bgzf test/bgziptest.txt
+	test/test-parse-reg -t test/colons.bam
 	cd test/tabix && ./test-tabix.sh tabix.tst
 	REF_PATH=: test/sam test/ce.fa test/faidx.fa
 	test/test-regidx
@@ -384,6 +386,9 @@ test/test_realn: test/test_realn.o libhts.a
 test/test-regidx: test/test-regidx.o libhts.a
 	$(CC) $(LDFLAGS) -o $@ test/test-regidx.o libhts.a $(LIBS) -lpthread
 
+test/test-parse-reg: test/test-parse-reg.o libhts.a
+	$(CC) $(LDFLAGS) -o $@ test/test-parse-reg.o libhts.a $(LIBS) -lpthread
+
 test/test_view: test/test_view.o libhts.a
 	$(CC) $(LDFLAGS) -o $@ test/test_view.o libhts.a $(LIBS) -lpthread
 

hts.c

@@ -2490,6 +2490,12 @@ long long hts_parse_decimal(const char *str, char **strend, int flags)
         if (esign == '-') e = -e;
     }
 
+    switch (*s) {
+    case 'k': case 'K': e += 3; s++; break;
+    case 'm': case 'M': e += 6; s++; break;
+    case 'g': case 'G': e += 9; s++; break;
+    }
+
     e -= decimals;
     while (e > 0) n *= 10, e--;
     while (e < 0) lost += n % 10, n /= 10, e++;
@@ -2501,45 +2507,90 @@ long long hts_parse_decimal(const char *str, char **strend, int flags)
     if (strend) {
         *strend = (char *)s;
     } else if (*s) {
-        hts_log_warning("Ignoring unknown characters after %.*s[%s]", (int)(s - str), str, s);
+        if ((flags & HTS_PARSE_THOUSANDS_SEP) || (!(flags & HTS_PARSE_THOUSANDS_SEP) && *s != ','))
+            hts_log_warning("Ignoring unknown characters after %.*s[%s]", (int)(s - str), str, s);
     }
 
     return (sign == '+')? n : -n;
 }
 
+static void *hts_memrchr(const void *s, int c, size_t n) {
+    size_t i;
+    unsigned char *u = (unsigned char *)s;
+    for (i = n; i > 0; i--) {
+        if (u[i-1] == c)
+            return u+i-1;
+    }
+
+    return NULL;
+}
+
 /*
  * 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".
  *
- * getid is optional and may be passed in as NULL.  If given it is used to
- * validate the reference name exists and is unambiguously parseable.  If not
- * given the best guess will be made but no has guarantees in validity.
+ * All parameters are mandatory.
+ *
+ * To work around ambiguous parsing issues, eg both "chr1" and "chr1:100-200"
+ * are reference names, we may quote using curly braces.
+ * Thus "{chr1}:100-200" and "{chr1:100-200}" disambiguate the above example.
+ *
+ * Flags are used to control how parsing works, and can be one of the below.
+ *
+ * HTS_PARSE_LIST:
+ *     If present, the region is assmed to be a comma separated list and
+ *     position parsing will not contain commas (this implicitly
+ *     clears HTS_PARSE_THOUSANDS_SEP in the call to hts_parse_decimal).
+ *     On success the return pointer will be the start of the next region, ie
+ *     the character after the comma.  (If *ret != '\0' then the caller can
+ *     assume another region is present in the list.)
+ *
+ *     If not set then positions may contain commas.  In this case the return
+ *     value should point to the end of the string, or NULL on failure.
  *
- * To work around these issues quoting is also permitted via {ref}:start-end.
- * In this case, the return value will point to '}' and not the end of the
- * reference (but this is a useful indication that it started with '{').
+ * HTS_PARSE_ONE_COORD:
+ *     If present, X:100 is treated as the single base pair region X:100-100.
+ *     In this case X:-100 is shorthand for X:1-100 and X:100- is X:100-<end>.
+ *     (This is the standard bcftools region convention.)
+ *
+ *     When not set X:100 is considered to be X:100-<end> where <end> is
+ *     the end of chromosome X (set to INT_MAX here).  X:100- and X:-100 are
+ *     invalid.
+ *     (This is the standard samtools region convention.)
+ *
+ * Note the supplied string expects 1 based inclusive coordinates, but the
+ * returned coordinates start from 0 and are half open, so pos0 is valid
+ * for use in e.g. "for (pos0 = beg; pos0 < end; pos0++) {...}"
  *
  * On success the end of the reference is returned (colon or end of string)
  *            beg/end will be set, plus tid if getid has been supplied.
  * On failure NULL is returned.
  */
-const char *hts_parse_reg2(const char *s, int *tid, int *beg, int *end,
-                           hts_name2id_f getid, void *hdr)
+const char *hts_parse_region(const char *s, int *tid, int *beg, int *end,
+                             hts_name2id_f getid, void *hdr, int flags)
 {
-    // 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
+    if (!s || !tid || !beg || !end || !getid)
+        return NULL;
 
-    const char *colon = NULL;
+    int s_len = strlen(s); // int is sufficient given beg/end types
+    kstring_t ks = { 0, 0, NULL };
+
+    const char *colon = NULL, *comma = NULL;
     int quoted = 0;
 
+    if (flags & HTS_PARSE_LIST)
+        flags &= ~HTS_PARSE_THOUSANDS_SEP;
+    else
+        flags |= HTS_PARSE_THOUSANDS_SEP;
+
+    const char *s_end = s + s_len;
+
     // Braced quoting of references is permitted to resolve ambiguities.
     if (*s == '{') {
-        const char *close = strrchr(s, '}');
+        const char *close = memchr(s, '}', s_len);
         if (!close) {
             hts_log_error("Mismatching braces in \"%s\"", s);
             return NULL;
@@ -2549,36 +2600,56 @@ const char *hts_parse_reg2(const char *s, int *tid, int *beg, int *end,
         if (close[1] == ':')
             colon = close+1;
         quoted = 1; // number of trailing characters to trim
+
+        // Truncate to this item only, if appropriate.
+        if (flags & HTS_PARSE_LIST) {
+            comma = strchr(close, ',');
+            if (comma) {
+                s_len = comma-s;
+                s_end = comma+1;
+            }
+        }
     } else {
-        colon = strrchr(s, ':');
+        // Truncate to this item only, if appropriate.
+        if (flags & HTS_PARSE_LIST) {
+            comma = strchr(s, ',');
+            if (comma) {
+                s_len = comma-s;
+                s_end = comma+1;
+            }
+        }
+
+        colon = hts_memrchr(s, ':', s_len);
     }
 
+    // No colon is simplest case; just check and return.
     if (colon == NULL) {
         *beg = 0; *end = INT_MAX;
-        if (getid) {
-            kstring_t ks = { 0, 0, NULL };
-            kputsn(s, s_len-quoted, &ks); // convert to nul terminated string
-            if (!ks.s) {
-                *tid = -1;
-                return NULL;
-            }
-
-            *tid = getid(hdr, ks.s);
-            free(ks.s);
-        } else {
-            *tid = 0;
+        kputsn(s, s_len-quoted, &ks); // convert to nul terminated string
+        if (!ks.s) {
+            *tid = -1;
+            return NULL;
         }
-        return *tid >= 0 ? s + s_len : NULL;
+
+        *tid = getid(hdr, ks.s);
+        free(ks.s);
+
+        return *tid >= 0 ? s_end : NULL;
     }
 
     // Has a colon, but check whole name first.
-    if (!quoted && getid) {
+    if (!quoted) {
         *beg = 0; *end = INT_MAX;
-        if ((*tid = getid(hdr, s)) >= 0) {
+        kputsn(s, s_len, &ks); // convert to nul terminated string
+        if (!ks.s) {
+            *tid = -1;
+            return NULL;
+        }
+        if ((*tid = getid(hdr, ks.s)) >= 0) {
             // Entire name matches, but also check this isn't
             // ambiguous.  eg we have ref chr1 and ref chr1:100-200
             // both present.
-            kstring_t ks = { 0, 0, NULL };
+            ks.l = 0;
             kputsn(s, colon-s, &ks); // convert to nul terminated string
             if (!ks.s) {
                 *tid = -1;
@@ -2594,39 +2665,81 @@ const char *hts_parse_reg2(const char *s, int *tid, int *beg, int *end,
             }
             free(ks.s);
 
-            return s + s_len;
+            return s_end;
         }
     }
 
-    char *hyphen;
-    *beg = hts_parse_decimal(colon+1, &hyphen, HTS_PARSE_THOUSANDS_SEP) - 1;
-    if (*beg < 0) *beg = 0;
-
-    if (*hyphen == '\0') *end = INT_MAX;
-    else if (*hyphen == '-') *end = hts_parse_decimal(hyphen+1, NULL, HTS_PARSE_THOUSANDS_SEP);
-    else return NULL;
+    // Quoted, or unquoted and whole string isn't a name.
+    // Check the pre-colon part is valid.
+    ks.l = 0;
+    kputsn(s, colon-s-quoted, &ks); // convert to nul terminated string
+    if (!ks.s) {
+        *tid = -1;
+        return NULL;
+    }
+    *tid = getid(hdr, ks.s);
+    free(ks.s);
+    if (*tid < 0)
+        return NULL;
 
-    if (*beg >= *end) return NULL;
-    if (getid) {
-        kstring_t ks = { 0, 0, NULL };
-        kputsn(s, colon-s-quoted, &ks); // convert to nul terminated string
-        if (!ks.s) {
-            *tid = -1;
+    // Finally parse the post-colon coordinates
+    char *hyphen;
+    *beg = hts_parse_decimal(colon+1, &hyphen, flags) - 1;
+    if (*beg < 0) {
+        if (isdigit(*hyphen) || *hyphen == '\0' || *hyphen == ',') {
+            // interpret chr:-100 as chr:1-100
+            *end = *beg==-1 ? INT_MAX : -(*beg+1);
+            *beg = 0;
+            return s_end;
+        } else if (*hyphen == '-') {
+            *beg = 0;
+        } else {
+            hts_log_error("Unexpected string \"%s\" after region", hyphen);
             return NULL;
         }
-        *tid = getid(hdr, ks.s);
-        free(ks.s);
-        if (*tid < 0)
+    }
+
+    if (*hyphen == '\0' || ((flags & HTS_PARSE_LIST) && *hyphen == ',')) {
+        *end = flags & HTS_PARSE_ONE_COORD ? *beg+1 : INT_MAX;
+    } else if (*hyphen == '-') {
+        *end = hts_parse_decimal(hyphen+1, &hyphen, flags);
+        if (*hyphen != '\0' && *hyphen != ',') {
+            hts_log_error("Unexpected string \"%s\" after region", hyphen);
             return NULL;
+        }
     } else {
-        *tid = 0;
+        hts_log_error("Unexpected string \"%s\" after region", hyphen);
+        return NULL;
     }
-    return colon;
+
+    if (*end == 0)
+        *end = INT_MAX; // interpret chr:100- as chr:100-<end>
+
+    if (*beg >= *end) return NULL;
+
+    return s_end;
 }
 
+// Next release we should mark this as deprecated?
+// Use hts_parse_region above instead.
 const char *hts_parse_reg(const char *s, int *beg, int *end)
 {
-    return hts_parse_reg2(s, NULL, beg, end, NULL, NULL);
+    char *hyphen;
+    const char *colon = strrchr(s, ':');
+    if (colon == NULL) {
+        *beg = 0; *end = INT_MAX;
+        return s + strlen(s);
+    }
+
+    *beg = hts_parse_decimal(colon+1, &hyphen, HTS_PARSE_THOUSANDS_SEP) - 1;
+    if (*beg < 0) *beg = 0;
+
+    if (*hyphen == '\0') *end = INT_MAX;
+    else if (*hyphen == '-') *end = hts_parse_decimal(hyphen+1, NULL, HTS_PARSE_THOUSANDS_SEP);
+    else return NULL;
+
+    if (*beg >= *end) return NULL;
+    return colon;
 }
 
 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)
@@ -2638,7 +2751,7 @@ hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f g
     else if (strcmp(reg, "*") == 0)
         return itr_query(idx, HTS_IDX_NOCOOR, 0, 0, readrec);
 
-    if (!hts_parse_reg2(reg, &tid, &beg, &end, getid, hdr))
+    if (!hts_parse_region(reg, &tid, &beg, &end, getid, hdr, HTS_PARSE_THOUSANDS_SEP))
         return NULL;
 
     return itr_query(idx, tid, beg, end, readrec);

htslib/hts.h

@@ -641,6 +641,8 @@ int hts_idx_set_meta(hts_idx_t *idx, uint32_t l_meta, uint8_t *meta, int is_copy
 
 
 #define HTS_PARSE_THOUSANDS_SEP 1  ///< Ignore ',' separators within numbers
+#define HTS_PARSE_ONE_COORD     2  ///< chr:pos means chr:pos-pos and not chr:pos-end
+#define HTS_PARSE_LIST          4  ///< Expect a comma separated list of regions. (Disables HTS_PARSE_THOUSANDS_SEP)
 
 /// Parse a numeric string
 /** The number may be expressed in scientific notation, and optionally may
@@ -666,20 +668,31 @@ long long hts_parse_decimal(const char *str, char **strend, int flags);
 
 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);
 
+//typedef struct {
+//    const char *chr_beg, *chr_end;
+//    int tid;
+//    int start, end; // half open, zero based
+//} hts_parsed_region;
+//const char *hts_parse_region_(const char *str, hts_parsed_region *reg,
+//                              hts_name2id_f getid, void *hdr, int flags);
+
+
 /// 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.
+    @param flags Bitwise HTS_PARSE_* flags listed above.
     @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);
+const char *hts_parse_region(const char *str, int *tid, int *beg, int *end,
+                             hts_name2id_f getid, void *hdr, int flags);
 
     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);

htslib/sam.h

@@ -275,7 +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);
+    const char *sam_parse_region(bam_hdr_t *h, const char *s, int *tid, int *beg, int *end, int flags);
     bam_hdr_t* bam_hdr_dup(const bam_hdr_t *h0);
 
     bam1_t *bam_init1(void);

sam.c

@@ -285,9 +285,8 @@ 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);
+const char *sam_parse_region(bam_hdr_t *h, const char *s, int *tid, int *beg, int *end, int flags) {
+    return hts_parse_region(s, tid, beg, end, (hts_name2id_f)bam_name2id, h, flags);
 }
 
 /*************************