Add UMI support to FASTQ input/output.

The fastq_umi option (FASTQ_OPT_UMI enum) is used to enable UMI
parsing in read names.

When reading FASTQ it converts the 8th Illumina field to an aux tag
(default to RX).  We may need to amend this if people require it to
work on earlier Illumina naming systems, but for now we target the
current software.

When writing FASTQ we hunt for a series of tags and choose the first
one found.  RX is the usual one, but users may wish to also use OX if
they potentially have error corrected data in RX and want to
regenerate fastq from the original uncorrect tag.

Complexities arrive when dealing with /1 or /2 and #num multi-plexing
strings, meaning we have to use temporary buffers rather than simply
truncating the read names.

Note we convert dual-index UMIs of SEQ+SEQ to RX:Z:SEQ-SEQ as per the
SAMtags recommendation.  However it's a bit wild west out there and
not everyone does this.  (I've seen 10X outputs with CR using
underscores for example.)  This is a best efforts approach to reduce
the likelihood of incompatibility with downstream pipelines.

Author: jkbonfield

hts.c

@@ -1208,6 +1208,14 @@ int hts_opt_add(hts_opt **opts, const char *c_arg) {
         strcmp(o->arg, "FASTQ_NAME2") == 0)
         o->opt = FASTQ_OPT_NAME2, o->val.i = 1;
 
+    else if (strcmp(o->arg, "fastq_umi") == 0 ||
+        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);
@@ -1250,6 +1258,8 @@ int hts_opt_apply(htsFile *fp, hts_opt *opts) {
             case HTS_OPT_FILTER:
             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;
@@ -1841,6 +1851,8 @@ int hts_set_opt(htsFile *fp, enum hts_fmt_option opt, ...) {
         return 0;
 
     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

@@ -366,6 +366,16 @@ enum hts_fmt_option {
     // name to the second field and insert a constructed <run>.<number>
     // name in its place.
     FASTQ_OPT_NAME2,
+
+    // Process the UMI tag.  Tag or Tag,tag,tag...
+    // On read, this converts the last read-name element (Illumina) to the tag.
+    // 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"
@@ -3701,6 +3702,7 @@ int sam_set_threads(htsFile *fp, int nthreads) {
     return 0;
 }
 
+#define UMI_TAGS 5
 typedef struct {
     kstring_t name;
     kstring_t comment; // NB: pointer into name, do not free
@@ -3710,9 +3712,11 @@ typedef struct {
     int aux;
     int rnum;
     char BC[3];         // aux tag ID for barcode
+    char UMI[UMI_TAGS][3]; // aux tag list for UMIs.
     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.
@@ -3723,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;
 }
@@ -3735,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);
     }
 }
@@ -3795,6 +3806,52 @@ int fastq_state_set(samFile *fp, enum hts_fmt_option opt, ...) {
         break;
     }
 
+    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);
+        if (!bc || strcmp(bc, "1") == 0)
+            bc = "RX";
+        int ntags = 0, err = 0;
+        for (ntags = 0; *bc && ntags < UMI_TAGS; ntags++) {
+            if (!isalpha(bc[0]) || !isalnum_c(bc[1])) {
+                err = 1;
+                break;
+            }
+
+            strncpy(x->UMI[ntags], bc, 3);
+            bc += 2;
+            if (*bc && *bc != ',') {
+                err = 1;
+                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;
+
+
+        if (err)
+            hts_log_warning("Bad UMI tag list '%s'", bc_orig);
+
+        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;
@@ -3907,6 +3964,43 @@ static int fastq_parse1(htsFile *fp, bam1_t *b) {
         x->name.s[x->name.l-=2] = 0;
     }
 
+    // Strip Illumina formatted UMI off read-name
+    char UMI_seq[256]; // maximum length in spec
+    size_t UMI_len = 0;
+    if (x->UMI[0][0]) {
+        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;
+            }
+
+            // The SAMTags spec recommends (but not requires) separating
+            // barcodes with hyphen ('-').
+            size_t i;
+            for (i = 0; i < UMI_len; i++)
+                UMI_seq[i] = isalpha_c(x->name.s[i+match[1].rm_so])
+                    ? x->name.s[i+match[1].rm_so]
+                    : '-';
+
+            // Move any trailing #num earlier in the name
+            if (UMI_len) {
+                UMI_seq[UMI_len++] = 0;
+
+                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;
+            }
+        }
+    }
+
     // Convert to BAM
     ret = bam_set1(b,
                    x->name.s + x->name.l - name, name,
@@ -3918,6 +4012,12 @@ static int fastq_parse1(htsFile *fp, bam1_t *b) {
                    0);
     if (ret < 0) return -2;
 
+    // Add UMI tag if removed from read-name above
+    if (UMI_len) {
+        if (bam_aux_append(b, x->UMI[0], 'Z', UMI_len, (uint8_t *)UMI_seq) < 0)
+            ret = -2;
+    }
+
     // Identify Illumina CASAVA strings.
     // <read>:<is_filtered>:<control_bits>:<barcode_sequence>
     char *barcode = NULL;
@@ -4260,6 +4360,39 @@ int fastq_format1(fastq_state *x, const bam1_t *b, kstring_t *str)
     if (kputc(x->nprefix, str) == EOF || kputs(bam_get_qname(b), str) == EOF)
         return -1;
 
+    // UMI tag
+    if (x && *x->UMI[0]) {
+        // Temporary copy of '#num' if present
+        char plex[256];
+        size_t len = str->l;
+        while (len && str->s[len] != ':' && str->s[len] != '#')
+            len--;
+
+        if (str->s[len] == '#' && str->l - len < 255) {
+            memcpy(plex, &str->s[len], str->l - len);
+            plex[str->l - len] = 0;
+            str->l = len;
+        } else {
+            *plex = 0;
+        }
+
+        uint8_t *bc = NULL;
+        int n;
+        for (n = 0; !bc && n < UMI_TAGS; n++)
+            bc = bam_aux_get(b, x->UMI[n]);
+        if (bc && *bc == 'Z') {
+            int err = kputc(':', str) < 0;
+            // Replace any non-alpha with '+'
+            while (*++bc)
+                err |= kputc(isalpha_c(*bc) ? toupper_c(*bc) : '+', str) < 0;
+            if (err)
+                return -1;
+        }
+
+        if (*plex && kputs(plex, str) < 0)
+            return -1;
+    }
+
     // /1 or /2 suffix
     if (x && x->rnum && (flag & BAM_FPAIRED)) {
         int r12 = flag & (BAM_FREAD1 | BAM_FREAD2);

test/fastq/UMI.fq

@@ -0,0 +1,28 @@
+@HS25:09827:FID:2:1201:1505:59794
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59795:A
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59796:ACG+TTTTTA
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59797:TACG+TTTTTA/1
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59797:TACG+TTTTTA/2
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59798:TTA#49/1
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+@HS25:09827:FID:2:1201:1505:59798:TTA#49/2
+CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT
++
+CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI

test/fastq/UMI.sam

@@ -0,0 +1,7 @@
+HS25:09827:FID:2:1201:1505:59794	4	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI
+HS25:09827:FID:2:1201:1505:59795	4	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:A
+HS25:09827:FID:2:1201:1505:59796	4	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:ACG-TTTTTA
+HS25:09827:FID:2:1201:1505:59797	77	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:TACG-TTTTTA
+HS25:09827:FID:2:1201:1505:59797	141	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:TACG-TTTTTA
+HS25:09827:FID:2:1201:1505:59798#49	77	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:TTA
+HS25:09827:FID:2:1201:1505:59798#49	141	*	0	0	*	*	0	0	CCGTTAGAGCATTTGTTGAAAATGCTTTCCTTGCTCCATGTGATGACTCT	CABCFGDEEFFEFHGHGGFFGDIGIJFIFHHGHEIFGHBCGHDIFBE9GI	RX:Z:TTA

test/fastq/fastq.tst

@@ -84,6 +84,9 @@ P r2-q.sam $tview -i fastq_aux r2.fa
 P name2.sam $tview -i fastq_name2 name2.fq
 P name2-q.sam $tview -i fastq_name2 name2.fa
 
+# UMI barcodes
+P UMI.sam $tview -i fastq_umi=RX UMI.fq
+
 # --------------------
 # Writing
 
@@ -114,3 +117,6 @@ P r1.fq $tview -f -o fastq_aux -o fastq_rnum r1.sam
 P r2.fq $tview -f -o fastq_aux -o fastq_rnum r2.sam
 P r1.fa $tview -F -o fastq_aux -o fastq_rnum r1.sam
 P r2.fa $tview -F -o fastq_aux -o fastq_rnum r2.sam
+
+# UMI barcodes
+P UMI.fq $tview -f -o fastq_rnum -o fastq_umi UMI.sam