Switch to regexp for matching Illumina read names (FASTQ UMI).

The default regexp is still the 8th component:

    "^[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:([^:#]+)"

This is changed using the fastq_umi_regex option.

Author: jkbonfield

hts.c

@@ -1212,6 +1212,10 @@ int hts_opt_add(hts_opt **opts, const char *c_arg) {
         strcmp(o->arg, "FASTQ_UMI") == 0)
         o->opt = FASTQ_OPT_UMI, o->val.s = val;
 
+    else if (strcmp(o->arg, "fastq_umi_regex") == 0 ||
+        strcmp(o->arg, "FASTQ_UMI_REGEX") == 0)
+        o->opt = FASTQ_OPT_UMI_REGEX, o->val.s = val;
+
     else {
         hts_log_error("Unknown option '%s'", o->arg);
         free(o->arg);
@@ -1255,6 +1259,7 @@ int hts_opt_apply(htsFile *fp, hts_opt *opts) {
             case FASTQ_OPT_AUX:
             case FASTQ_OPT_BARCODE:
             case FASTQ_OPT_UMI:
+            case FASTQ_OPT_UMI_REGEX:
                 if (hts_set_opt(fp,  opts->opt,  opts->val.s) != 0)
                     return -1;
                 break;
@@ -1847,6 +1852,7 @@ int hts_set_opt(htsFile *fp, enum hts_fmt_option opt, ...) {
 
     case FASTQ_OPT_BARCODE:
     case FASTQ_OPT_UMI:
+    case FASTQ_OPT_UMI_REGEX:
         if (fp->format.format == fastq_format ||
             fp->format.format == fasta_format) {
             va_start(args, opt);

htslib/hts.h

@@ -372,6 +372,10 @@ enum hts_fmt_option {
     // On write, it queries the tags in turn and copies the first found
     // to the read name suffix, converting any non-alpha to "+".
     FASTQ_OPT_UMI,
+
+    // Regex to use for matching read name.
+    // Def: "^[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:([^:#]+)"
+    FASTQ_OPT_UMI_REGEX,
 };
 
 // Profile options for encoding; primarily used at present in CRAM

sam.c

@@ -36,6 +36,7 @@ DEALINGS IN THE SOFTWARE.  */
 #include <signal.h>
 #include <inttypes.h>
 #include <unistd.h>
+#include <regex.h>
 
 #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
 #include "fuzz_settings.h"
@@ -3715,6 +3716,7 @@ typedef struct {
     khash_t(tag) *tags; // which aux tags to use (if empty, use all).
     char nprefix;
     int sra_names;
+    regex_t regex;
 } fastq_state;
 
 // Initialise fastq state.
@@ -3725,6 +3727,12 @@ static fastq_state *fastq_state_init(int name_char) {
         return NULL;
     strcpy(x->BC, "BC");
     x->nprefix = name_char;
+    // Default Illumina naming convention
+    char *re = "^[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:[^:]+:([^:#]+)";
+    if (regcomp(&x->regex, re, REG_EXTENDED) != 0) {
+        free(x);
+        return NULL;
+    }
 
     return x;
 }
@@ -3737,6 +3745,7 @@ void fastq_state_destroy(htsFile *fp) {
         ks_free(&x->name);
         ks_free(&x->seq);
         ks_free(&x->qual);
+        regfree(&x->regex);
         free(fp->state);
     }
 }
@@ -3798,6 +3807,7 @@ int fastq_state_set(samFile *fp, enum hts_fmt_option opt, ...) {
     }
 
     case FASTQ_OPT_UMI: {
+        // UMI tag: an empty string disables UMI by setting x->UMI[0] to \0\0\0
         va_start(args, opt);
         char *bc = va_arg(args, char *), *bc_orig = bc;
         va_end(args);
@@ -3817,6 +3827,7 @@ int fastq_state_set(samFile *fp, enum hts_fmt_option opt, ...) {
                 break;
             }
             bc+=(*bc==',');
+            x->UMI[ntags][2] = 0;
         }
         for (; ntags < UMI_TAGS; ntags++)
             x->UMI[ntags][0] = x->UMI[ntags][1] = x->UMI[ntags][2] = 0;
@@ -3828,6 +3839,19 @@ int fastq_state_set(samFile *fp, enum hts_fmt_option opt, ...) {
         break;
     }
 
+    case FASTQ_OPT_UMI_REGEX: {
+        va_start(args, opt);
+        char *re = va_arg(args, char *);
+        va_end(args);
+
+        regfree(&x->regex);
+        if (regcomp(&x->regex, re, REG_EXTENDED) != 0) {
+            hts_log_error("Regular expression '%s' is not supported", re);
+            return -1;
+        }
+        break;
+    }
+
     case FASTQ_OPT_RNUM:
         x->rnum = 1;
         break;
@@ -3944,30 +3968,23 @@ static int fastq_parse1(htsFile *fp, bam1_t *b) {
     char UMI_seq[256]; // maximum length in spec
     size_t UMI_len = 0;
     if (x->UMI[0][0]) {
-        // 7th colon
-        char *cp = x->name.s;
-        int nc = 0;
-        for (nc = 0; nc < 7; nc++, cp++) {
-            if ((cp = strchr(cp, ':')) == NULL)
-                break;
-        }
-        if (cp) {
-            x->name.l = cp - x->name.s - 1;
-
-            // UMI ends at string end or '#num'.
-            char *UMI_start = cp;
-            while (*cp && *cp != '#') {
-                UMI_len++;
-                cp++;
-            }
+        regmatch_t match[3];
+        if (regexec(&x->regex, x->name.s, 2, match, 0) == 0
+            && match[0].rm_so >= 0     // whole regex
+            && match[1].rm_so >= 0) {  // bracketted UMI component
+            UMI_len = match[1].rm_eo - match[1].rm_so;
             if (UMI_len > 255) {
                 hts_log_error("SAM read name is too long");
                 return -2;
             }
-            memcpy(UMI_seq, UMI_start, UMI_len);
+            memcpy(UMI_seq, x->name.s + match[1].rm_so, UMI_len);
             UMI_seq[UMI_len++] = 0;
 
             // Move any trailing #num earlier in the name
+            x->name.l = match[1].rm_so;
+            if (x->name.l > 0 && x->name.s[x->name.l-1] == ':')
+                x->name.l--; // remove colon too
+            char *cp = x->name.s + match[1].rm_eo;
             while (*cp)
                 x->name.s[x->name.l++] = *cp++;
             x->name.s[x->name.l] = 0;