perf: build HMMER with AlphaFold 3's jackhmmer speed-up patches

docker/alphafold2.dockerfile

@@ -44,6 +44,19 @@ RUN set -eux; \
     cmake --build build -j"$(nproc)"; \
     cmake --install build
 
+# ---------------------------------------------------------------------------
+# HMMER 3.4 with AlphaFold 3's jackhmmer patches (docker/hmmer/build_hmmer.sh),
+# built like hhblits above.
+# ---------------------------------------------------------------------------
+FROM ubuntu:20.04 AS hmmer-builder
+RUN set -eux; \
+    apt-get update; \
+    DEBIAN_FRONTEND=noninteractive apt-get install -y --no-install-recommends \
+      build-essential ca-certificates curl patch; \
+    rm -rf /var/lib/apt/lists/*
+COPY docker/hmmer/ /tmp/hmmer/
+RUN bash /tmp/hmmer/build_hmmer.sh /opt/hmmer
+
 FROM nvidia/cuda:${CUDA}-cudnn8-runtime-ubuntu20.04 AS base
 ARG CUDA
 ARG MMSEQS_VERSION
@@ -126,6 +139,14 @@ RUN set -eux; \
       *) printf '%s\n' "$out" | head -3; echo "ERROR: hhblits with the realignment fix missing"; exit 1 ;; \
     esac
 
+# Overwrite conda's HMMER the same way. Only the patched build has --seq_limit.
+COPY --from=hmmer-builder /opt/hmmer/bin/ /opt/conda/bin/
+RUN set -eux; \
+    case "$(/opt/conda/bin/jackhmmer -h)" in \
+      *--seq_limit*) echo "patched HMMER in place" ;; \
+      *) echo "ERROR: patched HMMER missing"; exit 1 ;; \
+    esac
+
 #RUN micromamba run -n base python -m pip install --no-cache-dir "openmm==8.1.1"
 RUN python -m pip install --no-cache-dir "setuptools<82" # setuptools>82 breaks pdbfixer at relaxation
 

docker/alphafold3.dockerfile

@@ -60,22 +60,10 @@ COPY --from=ghcr.io/astral-sh/uv:0.9.24 /uv /uvx /bin/
 RUN uv venv ${VIRTUAL_ENV}
 
 # ---------------------------------------------------------------------
-# HMMER (with seq_limit patch)
-# ---------------------------------------------------------------------
-RUN mkdir /hmmer_build /hmmer && \
-    wget http://eddylab.org/software/hmmer/hmmer-3.4.tar.gz -O /hmmer_build/hmmer-3.4.tar.gz && \
-    echo "ca70d94fd0cf271bd7063423aabb116d42de533117343a9b27a65c17ff06fbf3  /hmmer_build/hmmer-3.4.tar.gz" | sha256sum -c - && \
-    tar -xzf /hmmer_build/hmmer-3.4.tar.gz -C /hmmer_build
-
-COPY alphafold3/docker/jackhmmer_seq_limit.patch /hmmer_build/
-
-RUN cd /hmmer_build && \
-    patch -p0 < jackhmmer_seq_limit.patch && \
-    cd /hmmer_build/hmmer-3.4 && \
-    ./configure --prefix=/hmmer && \
-    make -j$(nproc) && \
-    make install && \
-    cd easel && make install && \
+# HMMER with AlphaFold 3's jackhmmer patches (docker/hmmer)
+# ---------------------------------------------------------------------
+COPY docker/hmmer/ /hmmer_build/
+RUN bash /hmmer_build/build_hmmer.sh /hmmer && \
     rm -rf /hmmer_build
 
 # ---------------------------------------------------------------------

docker/hmmer/SHA256SUMS

@@ -0,0 +1,4 @@
+ca70d94fd0cf271bd7063423aabb116d42de533117343a9b27a65c17ff06fbf3  hmmer-3.4.tar.gz
+df9e3ae35ad1659921d96ebfca67a9616a7a467ddde2be18a56f9bd3edb38c41  jackhmmer_seq_limit.patch
+b2e18b9d24b3b0ff0b6dc5a1c24bcd26ad71a099688c321731aca3fde0ef3090  easel_fast_fasta_parsing.patch
+401b6ecb1ef3558b6b9a2b13070b5b278a716b7d7bb71a6182c736da5b0c48ee  jackhmmer_parallel_fasta.patch

docker/hmmer/build_hmmer.sh

@@ -0,0 +1,63 @@
+#!/usr/bin/env bash
+# Build HMMER 3.4 with AlphaFold 3's jackhmmer patches and install it into PREFIX.
+#
+# The patches are docker/*.patch from google-deepmind/alphafold3@03f4f48 (Apache-2.0):
+#   jackhmmer_seq_limit.patch       adds --seq_limit, used by AlphaFold 3's jackhmmer wrapper
+#   easel_fast_fasta_parsing.patch  faster FASTA parsing (EddyRivasLab/easel#83)
+#   jackhmmer_parallel_fasta.patch  every worker reads its own part of the database
+#                                   (EddyRivasLab/hmmer#353)
+set -euo pipefail
+
+prefix="${1:?usage: build_hmmer.sh PREFIX}"
+here="$(cd "$(dirname "$0")" && pwd)"
+work="$(mktemp -d)"
+trap 'rm -rf "$work"' EXIT
+cd "$work"
+
+url=https://eddylab.org/software/hmmer/hmmer-3.4.tar.gz
+if command -v curl >/dev/null; then
+  curl -fsSL "$url" -o hmmer-3.4.tar.gz
+else
+  wget -q "$url" -O hmmer-3.4.tar.gz
+fi
+cp "$here"/*.patch .
+sha256sum -c "$here/SHA256SUMS"
+
+tar -xzf hmmer-3.4.tar.gz
+# In AlphaFold 3's order: both jackhmmer patches edit src/jackhmmer.c.
+for p in jackhmmer_seq_limit easel_fast_fasta_parsing jackhmmer_parallel_fasta; do
+  patch -p0 --forward --fuzz=0 < "$p.patch"
+done
+(
+  cd hmmer-3.4
+  ./configure --prefix="$prefix"
+  make -j"$(nproc)"
+  make install
+  make -C easel install
+)
+
+# --cpu 0 reads the database serially; --cpu 4 takes the patched path where each
+# worker reads its own 512 KB chunks. Copies of the query at the start, middle and
+# end of a ~6 MB database must give the same hits and target count both ways.
+bin="$prefix/bin"
+"$bin/esl-shuffle" -G --amino -N 1 -L 400 --seed 1 | sed '1s/.*/>query/' > query.fa
+"$bin/esl-shuffle" -G --amino -N 10000 -L 300 --seed 2 > a.fa
+"$bin/esl-shuffle" -G --amino -N 10000 -L 300 --seed 3 > b.fa
+{
+  sed '1s/.*/>q_start first copy/' query.fa
+  cat a.fa
+  sed '1s/.*/>q_mid second copy/' query.fa
+  cat b.fa
+  sed '1s/.*/>q_end third copy/' query.fa
+} > db.fa
+for cpu in 0 4; do
+  "$bin/jackhmmer" --cpu "$cpu" -N 2 -o "out.$cpu" --tblout "tbl.$cpu" -A "sto.$cpu" \
+    query.fa db.fa
+  grep -v '^#' "tbl.$cpu" > "hits.$cpu"
+  grep -E '^(Target sequences|Passed .* filter)' "out.$cpu" >> "hits.$cpu"
+done
+grep -q '^Target sequences: *20003 ' hits.4
+test "$(grep -c '^q_' hits.4)" -eq 3
+cmp hits.0 hits.4
+cmp sto.0 sto.4
+echo "HMMER 3.4 with jackhmmer patches installed in $prefix; threaded and serial searches agree"

docker/hmmer/easel_fast_fasta_parsing.patch

@@ -0,0 +1,166 @@
+--- hmmer-3.4/easel/esl_sqio_ascii.c
++++ hmmer-3.4/easel/esl_sqio_ascii.c
+@@ -62,7 +62,7 @@
+ static int  loadmem  (ESL_SQFILE *sqfp);
+ static int  loadbuf  (ESL_SQFILE *sqfp);
+ static int  nextchar (ESL_SQFILE *sqfp, char *ret_c);
+-static int  seebuf   (ESL_SQFILE *sqfp, int64_t maxn, int64_t *opt_nres, int64_t *opt_endpos);
++static int  seebuf   (ESL_SQFILE *sqfp, ESL_SQ *sq, int64_t maxn, int64_t *opt_nres, int64_t *opt_endpos);
+ static void addbuf   (ESL_SQFILE *sqfp, ESL_SQ *sq, int64_t nres);
+ static void skipbuf  (ESL_SQFILE *sqfp, int64_t nskip);
+ static int  read_nres(ESL_SQFILE *sqfp, ESL_SQ *sq, int64_t nskip, int64_t nres, int64_t *opt_actual_nres);
+@@ -786,13 +786,13 @@
+   if ((status = ascii->parse_header(sqfp, sq)) != eslOK) return status; /* EMEM, EOF, EFORMAT */
+
+   do {
+-    if ((status = seebuf(sqfp, -1, &n, &epos)) == eslEFORMAT) return status;
+-    if (esl_sq_GrowTo(sq, sq->n + n) != eslOK) return eslEMEM;
+-    addbuf(sqfp, sq, n);
++    if (sq->n >= sq->salloc - 2 && esl_sq_GrowTo(sq, sq->salloc * 2) != eslOK) return eslEMEM;
++    if ((status = seebuf(sqfp, sq, sq->salloc - 2 - sq->n, &n, &epos)) == eslEFORMAT) return status;
++    ascii->bpos = epos;
+     ascii->L   += n;
+     sq->eoff   = ascii->boff + epos - 1;
+     if (status == eslEOD)     break;
+-  } while ((status = loadbuf(sqfp)) == eslOK);
++  } while (epos < ascii->nc || (status = loadbuf(sqfp)) == eslOK);
+
+   if      (status == eslEOF)
+     {
+@@ -886,7 +886,7 @@
+
+   ascii->L       = 0;
+   do {
+-    status = seebuf(sqfp, -1, &n, &epos);
++    status = seebuf(sqfp, NULL, -1, &n, &epos);
+     ascii->L += n;
+     sq->eoff = ascii->boff + epos - 1;
+     if (status == eslEFORMAT) return status;
+@@ -986,7 +986,7 @@
+   if ((status = ascii->skip_header(sqfp, sq)) != eslOK) return status; /* EOF, EFORMAT */
+
+   do {
+-    if ((status = seebuf(sqfp, -1, &n, &epos)) == eslEFORMAT) return status;
++    if ((status = seebuf(sqfp, NULL, -1, &n, &epos)) == eslEFORMAT) return status;
+     if (esl_sq_GrowTo(sq, sq->n + n) != eslOK) return eslEMEM;
+     addbuf(sqfp, sq, n);
+     ascii->L   += n;
+@@ -2247,9 +2247,9 @@
+  * the need for the <maxn> limit, which ReadWindow() uses.)
+  */
+ static int
+-seebuf(ESL_SQFILE *sqfp, int64_t maxn, int64_t *opt_nres, int64_t *opt_endpos)
++seebuf(ESL_SQFILE *sqfp, ESL_SQ *sq, int64_t maxn, int64_t *opt_nres, int64_t *opt_endpos)
+ {
+-  int     bpos;
++  int64_t bpos;
+   int64_t nres  = 0;
+   int64_t nres2 = 0;/* an optimization for determining lastrpl from nres, without incrementing lastrpl on every char */
+   int     sym;
+@@ -2257,19 +2257,38 @@
+   int     lasteol;
+   int     status  = eslOK;
+
+-  ESL_SQASCII_DATA *ascii = &sqfp->data.ascii;
++  ESL_DSQ dsq_inmap[128];
++  ESL_SQASCII_DATA * restrict ascii = &sqfp->data.ascii;
++  const unsigned char * restrict buf = (const unsigned char *) ascii->buf;
++  const ESL_DSQ * restrict inmap = sqfp->inmap;
++  ESL_DSQ * restrict dst = (sq == NULL) ? NULL : (sq->dsq != NULL ? sq->dsq + sq->n + 1 : (ESL_DSQ *) (sq->seq + sq->n));
++  if (sq != NULL && sq->dsq != NULL) {
++    if (sqfp->do_digital && sqfp->abc == sq->abc) {
++      if (sqfp->inmap['*'] == '*') sqfp->inmap['*'] = sq->abc->inmap['*'];
++    } else {
++      for (int c = 0; c < 128; c++) dsq_inmap[c] = (sqfp->inmap[c] <= 127) ? sq->abc->inmap[sqfp->inmap[c]] : sqfp->inmap[c];
++      inmap = dsq_inmap;
++    }
++  }
+
+   lasteol = ascii->bpos - 1;
+   if (maxn == -1) maxn = ascii->nc; /* makes for a more efficient test. nc is a guaranteed upper bound on nres */
+
+   for (bpos = ascii->bpos; nres < maxn && bpos < ascii->nc; bpos++)
+   {
+-      sym = ascii->buf[bpos];
++      sym = buf[bpos];
+       //printf ("nres: %d, bpos: %d  (%d)\n", nres, bpos, sym);
+       if (!isascii(sym)) ESL_FAIL(eslEFORMAT, ascii->errbuf, "Line %" PRId64 ": non-ASCII character %c in sequence", ascii->linenumber, sym); 
+-      x   = sqfp->inmap[sym];
++      x   = inmap[sym];
+
+-      if      (x <= 127) nres++;
++      if      (x <= 127) {
++        if (dst != NULL) {
++          ESL_DSQ * restrict d = dst + nres;
++          *d = x;
++          for (int64_t end = ESL_MIN(ascii->nc, bpos + maxn - nres) - 1; bpos < end && (sym = buf[bpos + 1], (sym | (x = inmap[sym & 127])) <= 127); ) { bpos++; *++d = x; }
++          nres = (d - dst) + 1;
++        } else nres++;
++      }
+       else if (x == eslDSQ_EOL) 
+       {
+          if (ascii->curbpl != -1) ascii->curbpl += bpos - lasteol;
+@@ -2299,6 +2318,7 @@
+     else if (x != eslDSQ_IGNORED) ESL_FAIL(eslEFORMAT, ascii->errbuf, "inmap corruption?");
+   }
+
++  if (sq         != NULL) sq->n      += nres;
+   if (ascii->curbpl != -1) ascii->curbpl += bpos - lasteol - 1;
+   if (ascii->currpl != -1) ascii->currpl += nres - nres2;
+   if (opt_nres   != NULL) *opt_nres   = nres;
+@@ -2464,11 +2484,11 @@
+   int     status      = eslOK;
+
+   ESL_SQASCII_DATA *ascii = &sqfp->data.ascii;
+-  status = seebuf(sqfp, nskip+nres, &n, &epos);
++  status = seebuf(sqfp, NULL, nskip+nres, &n, &epos);
+   while (status == eslOK && nskip - n > 0) {
+     nskip   -= n;
+     if ((status = loadbuf(sqfp)) == eslEOF) break;
+-    status = seebuf(sqfp, nskip+nres, &n, &epos);
++    status = seebuf(sqfp, NULL, nskip+nres, &n, &epos);
+   }
+
+   if         (status == eslEOF) { 
+@@ -2489,7 +2509,7 @@
+       actual_nres += n;
+       nres        -= n;
+       if ((status = loadbuf(sqfp)) == eslEOF) break;
+-      status = seebuf(sqfp, nres, &n, &epos);
++      status = seebuf(sqfp, NULL, nres, &n, &epos);
+     }
+
+
+@@ -3034,7 +3054,10 @@
+   pos = 0;
+   while (status == eslOK && ! isspace(c))
+   {
+-      sq->name[pos++] = c;
++      char * restrict d = sq->name + pos; const char * restrict s = ascii->buf + ascii->bpos;
++      int64_t k = 0, end = ESL_MIN(ascii->nc - ascii->bpos, sq->nalloc - 1 - pos) - 1;
++      for (d[0] = c; k < end && (unsigned char)s[k + 1] > ' '; ) { k++; d[k] = s[k]; }
++      pos += k + 1; ascii->bpos += k;
+       if (pos == sq->nalloc-1) { ESL_RALLOC(sq->name, tmp, sq->nalloc*2); sq->nalloc*=2; }
+       status = nextchar(sqfp, &c); 
+   }
+@@ -3048,7 +3071,10 @@
+   pos = 0;
+   while (status == eslOK && c != '\n' && c != '\r' && c != 1)
+   {
+-      sq->desc[pos++] = c;
++      char * restrict d = sq->desc + pos; const char * restrict s = ascii->buf + ascii->bpos;
++      int64_t k = 0, end = ESL_MIN(ascii->nc - ascii->bpos, sq->dalloc - 1 - pos) - 1;
++      for (d[0] = c; k < end && (unsigned char)s[k + 1] > '\r'; ) { k++; d[k] = s[k]; }
++      pos += k + 1; ascii->bpos += k;
+       if (pos == sq->dalloc-1) { ESL_RALLOC(sq->desc, tmp, sq->dalloc*2); sq->dalloc*= 2; }
+       status = nextchar(sqfp, &c); 
+   }
+@@ -3337,7 +3363,7 @@
+   if ((status = ascii->parse_header(&sqfp, sq)) != eslOK) return status; /* EOF, EFORMAT */
+
+   do {
+-    if ((status = seebuf(&sqfp, -1, &n, &epos)) == eslEFORMAT) return status;
++    if ((status = seebuf(&sqfp, NULL, -1, &n, &epos)) == eslEFORMAT) return status;
+     if (esl_sq_GrowTo(sq, sq->n + n) != eslOK) return eslEMEM;
+     addbuf(&sqfp, sq, n);
+     ascii->L   += n;