From 5b5a34eb47f0cbd82a2195c2a96796af80b93a5a Mon Sep 17 00:00:00 2001 From: Dima Molodenskiy Date: Mon, 19 Oct 2026 15:36:45 +0200 Subject: [PATCH] perf: build HMMER with AlphaFold 3's jackhmmer speed-up patches AlphaFold 3 (google-deepmind/alphafold3@03f4f48) now builds HMMER 3.4 with two more patches on top of --seq_limit: - easel_fast_fasta_parsing.patch: Easel copies residues while it scans the read buffer instead of in a second pass (EddyRivasLab/easel#83). Reported 2-2.5x faster jackhmmer. - jackhmmer_parallel_fasta.patch: on a rewindable FASTA every worker opens the database itself and claims 512 KB chunks, so jackhmmer is no longer limited by one reader thread (EddyRivasLab/hmmer#353). Reported 98 s -> 38 s at 8 CPUs and 94 s -> 18 s at 24 CPUs on UniRef90. Both upstream issues are open with no maintainer response yet. The three patches are vendored in docker/hmmer with their SHA256 sums (Apache-2.0, from that commit; the seq_limit patch is byte-identical to the one in the alphafold3 submodule). docker/hmmer/build_hmmer.sh downloads HMMER 3.4, checks every checksum, applies the patches with --fuzz=0 in AlphaFold 3's order, installs into a prefix, and then checks that a serial (--cpu 0) and a threaded (--cpu 4) jackhmmer give the same hits, alignment and target count on a generated ~6 MB database. - AF3 image: builds HMMER with the script instead of applying only the seq_limit patch from the submodule. - AF2 image: builds it in an ubuntu:20.04 stage and overwrites conda's HMMER binaries, as is already done for hhblits; the build fails if the installed jackhmmer lacks --seq_limit. - docs: how to install the same build into a conda environment. Co-Authored-By: Claude Opus 5.5 --- docker/alphafold2.dockerfile | 21 ++ docker/alphafold3.dockerfile | 20 +- docker/hmmer/SHA256SUMS | 4 + docker/hmmer/build_hmmer.sh | 63 ++++++ docker/hmmer/easel_fast_fasta_parsing.patch | 166 +++++++++++++++ docker/hmmer/jackhmmer_parallel_fasta.patch | 221 ++++++++++++++++++++ docker/hmmer/jackhmmer_seq_limit.patch | 32 +++ docs/backend_installation.md | 8 + 8 files changed, 519 insertions(+), 16 deletions(-) create mode 100644 docker/hmmer/SHA256SUMS create mode 100755 docker/hmmer/build_hmmer.sh create mode 100644 docker/hmmer/easel_fast_fasta_parsing.patch create mode 100644 docker/hmmer/jackhmmer_parallel_fasta.patch create mode 100644 docker/hmmer/jackhmmer_seq_limit.patch diff --git a/docker/alphafold2.dockerfile b/docker/alphafold2.dockerfile index 5757142b..cfa8748b 100644 --- a/docker/alphafold2.dockerfile +++ b/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 diff --git a/docker/alphafold3.dockerfile b/docker/alphafold3.dockerfile index 8392697d..cd8a3725 100644 --- a/docker/alphafold3.dockerfile +++ b/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 # --------------------------------------------------------------------- diff --git a/docker/hmmer/SHA256SUMS b/docker/hmmer/SHA256SUMS new file mode 100644 index 00000000..d3a01e8b --- /dev/null +++ b/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 diff --git a/docker/hmmer/build_hmmer.sh b/docker/hmmer/build_hmmer.sh new file mode 100755 index 00000000..fff95f27 --- /dev/null +++ b/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" diff --git a/docker/hmmer/easel_fast_fasta_parsing.patch b/docker/hmmer/easel_fast_fasta_parsing.patch new file mode 100644 index 00000000..0f64418e --- /dev/null +++ b/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 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; diff --git a/docker/hmmer/jackhmmer_parallel_fasta.patch b/docker/hmmer/jackhmmer_parallel_fasta.patch new file mode 100644 index 00000000..562595c7 --- /dev/null +++ b/docker/hmmer/jackhmmer_parallel_fasta.patch @@ -0,0 +1,221 @@ +--- hmmer-3.4/src/jackhmmer.c ++++ hmmer-3.4/src/jackhmmer.c +@@ -2,6 +2,7 @@ + */ + #include + ++#include + #include + #include + #include +@@ -34,6 +35,10 @@ + typedef struct { + #ifdef HMMER_THREADS + ESL_WORK_QUEUE *queue; ++ ESL_SQFILE *dbfp; ++ off_t *next_offset; ++ off_t file_size; ++ int sstatus; + #endif + P7_BG *bg; + P7_PIPELINE *pli; +@@ -549,6 +554,7 @@ + info[i].bg = p7_bg_Clone(bg); + #ifdef HMMER_THREADS + info[i].queue = queue; ++ info[i].dbfp = NULL; + #endif + } + +@@ -765,7 +771,12 @@ + /* Cleanup - prepare for successful exit + */ + for (i = 0; i < infocnt; ++i) +- p7_bg_Destroy(info[i].bg); ++ { ++ p7_bg_Destroy(info[i].bg); ++#ifdef HMMER_THREADS ++ esl_sqfile_Close(info[i].dbfp); ++#endif ++ } + + #ifdef HMMER_THREADS + if (ncpus > 0) +@@ -1633,14 +1644,60 @@ + } + + #ifdef HMMER_THREADS ++#define FASTA_CHUNK_SIZE (512 * 1024) ++ + static int + thread_loop(ESL_THREADS *obj, ESL_WORK_QUEUE *queue, ESL_SQFILE *dbfp) + { +- int status = eslOK; +- int sstatus = eslOK; ++ int status = eslOK; ++ int sstatus = eslOK; + int eofCount = 0; + ESL_SQ_BLOCK *block; + void *newBlock; ++ int nworkers = esl_threads_GetWorkerCount(obj); ++ int w; ++ ++ if (dbfp->format == eslSQFILE_FASTA && esl_sqfile_IsRewindable(dbfp)) ++ { ++ off_t next_offset = 0; ++ if (fseeko(dbfp->data.ascii.fp, 0, SEEK_END) != 0) p7_Fail("fseeko() failed"); ++ off_t file_size = ftello(dbfp->data.ascii.fp); ++ if (file_size < 0) p7_Fail("ftello() failed"); ++ ++ for (w = 0; w < nworkers; ++w) ++ { ++ WORKER_INFO *info = (WORKER_INFO *) esl_threads_GetData(obj, w); ++ if (info->dbfp == NULL && ++ esl_sqfile_OpenDigital(info->om->abc, dbfp->filename, dbfp->format, p7_SEQDBENV, &info->dbfp) != eslOK) ++ p7_Fail("Failed to open target sequence database %s for reading\n", dbfp->filename); ++ info->next_offset = &next_offset; ++ info->file_size = file_size; ++ info->sstatus = eslEOF; ++ } ++ ++ esl_threads_WaitForStart(obj); ++ esl_threads_WaitForFinish(obj); ++ ++ off_t err_offset = file_size; ++ sstatus = eslEOF; ++ for (w = 0; w < nworkers; ++w) ++ { ++ WORKER_INFO *info = (WORKER_INFO *) esl_threads_GetData(obj, w); ++ if (info->sstatus != eslEOF && info->file_size < err_offset) ++ { ++ sstatus = info->sstatus; ++ err_offset = info->file_size; ++ } ++ } ++ if (sstatus == eslEFORMAT) ++ { ++ ESL_SQ *tmpsq = esl_sq_CreateDigital(dbfp->abc); ++ esl_sqfile_Position(dbfp, 0); ++ while ((sstatus = esl_sqio_Read(dbfp, tmpsq)) == eslOK) esl_sq_Reuse(tmpsq); ++ esl_sq_Destroy(tmpsq); ++ } ++ return sstatus; ++ } + + esl_workqueue_Reset(queue); + esl_threads_WaitForStart(obj); +@@ -1699,6 +1756,114 @@ + + info = (WORKER_INFO *) esl_threads_GetData(obj, workeridx); + ++ if (info->dbfp != NULL) ++ { ++ ESL_SQ *dbsq = esl_sq_CreateDigital(info->om->abc); ++ ESL_SQASCII_DATA *ascii = &info->dbfp->data.ascii; ++ off_t file_size = info->file_size; ++ off_t cur_roff = -1; ++ off_t x0; ++ ++ while ((x0 = __atomic_fetch_add(info->next_offset, FASTA_CHUNK_SIZE, __ATOMIC_RELAXED)) < file_size) ++ { ++ off_t x1 = ESL_MIN(x0 + FASTA_CHUNK_SIZE, file_size); ++ status = eslOK; ++ ++ if (x0 <= cur_roff) ++ { ++ if (cur_roff >= x1) continue; ++ } ++ else if (x0 == 0) ++ { ++ if (ascii->boff != 0 || ascii->bpos != 0 || ascii->nc == 0) ++ status = esl_sqfile_Position(info->dbfp, 0); ++ char prev = 0; ++ while (status == eslOK && ascii->bpos < ascii->nc && isspace((unsigned char) ascii->buf[ascii->bpos])) ++ { ++ prev = ascii->buf[ascii->bpos++]; ++ if (ascii->bpos == ascii->nc) ++ status = esl_sqfile_Position(info->dbfp, ascii->boff + ascii->nc); ++ } ++ if (status == eslOK && ascii->boff + ascii->bpos >= x1 && ++ ascii->buf[ascii->bpos] == '>' && (prev == '\n' || prev == '\r')) ++ { ++ cur_roff = ascii->boff + ascii->bpos; ++ continue; ++ } ++ } ++ else ++ { ++ off_t pos = x0 - 1; ++ while (pos < x1 && (status = esl_sqfile_Position(info->dbfp, pos)) == eslOK) ++ { ++ off_t limit = ESL_MIN((off_t) ascii->nc, x1 - ascii->boff); ++ const char *p = ascii->buf + 1; ++ const char *end = ascii->buf + limit; ++ while (p < end && (p = memchr(p, '>', end - p)) != NULL) ++ { ++ if (p[-1] == '\n' || p[-1] == '\r') break; ++ ++p; ++ } ++ if (p != NULL && p < end) { ascii->bpos = (int) (p - ascii->buf); break; } ++ if (ascii->boff + limit >= x1 || ascii->nc <= 1) { pos = x1; break; } ++ pos = ascii->boff + ascii->nc - 1; ++ } ++ if (status == eslOK && pos >= x1) { cur_roff = -1; continue; } ++ } ++ if (status != eslOK) ++ { ++ if (status != eslEOF) ++ { ++ info->sstatus = status; ++ info->file_size = x0; ++ __atomic_store_n(info->next_offset, file_size, __ATOMIC_RELAXED); ++ } ++ break; ++ } ++ ++ cur_roff = -1; ++ while ((status = esl_sqio_Read(info->dbfp, dbsq)) == eslOK) ++ { ++ off_t next_roff = dbsq->eoff + 1; ++ ++ p7_pli_NewSeq(info->pli, dbsq); ++ p7_bg_SetLength(info->bg, dbsq->n); ++ p7_oprofile_ReconfigLength(info->om, dbsq->n); ++ ++ p7_Pipeline(info->pli, info->om, info->bg, dbsq, NULL, info->th); ++ ++ esl_sq_Reuse(dbsq); ++ p7_pipeline_Reuse(info->pli); ++ ++ if (next_roff >= x1) ++ { ++ if (next_roff >= file_size) break; ++ if (ascii->bpos == 0) ++ { ++ if ((status = esl_sqfile_Position(info->dbfp, next_roff - 1)) != eslOK) break; ++ ascii->bpos = 1; ++ } ++ if (ascii->buf[ascii->bpos - 1] == '\n' || ascii->buf[ascii->bpos - 1] == '\r') ++ { ++ cur_roff = next_roff; ++ break; ++ } ++ } ++ } ++ if (status != eslOK && status != eslEOF) ++ { ++ info->sstatus = status; ++ info->file_size = x0; ++ __atomic_store_n(info->next_offset, file_size, __ATOMIC_RELAXED); ++ break; ++ } ++ } ++ ++ esl_sq_Destroy(dbsq); ++ esl_threads_Finished(obj, workeridx); ++ return; ++ } ++ + status = esl_workqueue_WorkerUpdate(info->queue, NULL, &newBlock); + if (status != eslOK) p7_Fail("Work queue worker failed"); + diff --git a/docker/hmmer/jackhmmer_seq_limit.patch b/docker/hmmer/jackhmmer_seq_limit.patch new file mode 100644 index 00000000..9249a1f6 --- /dev/null +++ b/docker/hmmer/jackhmmer_seq_limit.patch @@ -0,0 +1,32 @@ +--- hmmer-3.4/src/jackhmmer.c ++++ hmmer-3.4/src/jackhmmer.c +@@ -73,6 +73,7 @@ static ESL_OPTIONS options[] = { + { "--noali", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "don't output alignments, so output is smaller", 2 }, + { "--notextw", eslARG_NONE, NULL, NULL, NULL, NULL, NULL, "--textw", "unlimit ASCII text output line width", 2 }, + { "--textw", eslARG_INT, "120", NULL, "n>=120", NULL, NULL, "--notextw", "set max width of ASCII text output lines", 2 }, ++ { "--seq_limit", eslARG_INT, NULL, NULL, NULL, NULL, NULL, "--seq_limit", "if set, truncate all hits after this value is reached", 2 }, + /* Control of scoring system */ + { "--popen", eslARG_REAL, "0.02", NULL, "0<=x<0.5",NULL, NULL, NULL, "gap open probability", 3 }, + { "--pextend", eslARG_REAL, "0.4", NULL, "0<=x<1", NULL, NULL, NULL, "gap extend probability", 3 }, +@@ -298,6 +299,7 @@ output_header(FILE *ofp, ESL_GETOPTS *go + if (esl_opt_IsUsed(go, "--noali") && fprintf(ofp, "# show alignments in output: no\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (esl_opt_IsUsed(go, "--notextw") && fprintf(ofp, "# max ASCII text line length: unlimited\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (esl_opt_IsUsed(go, "--textw") && fprintf(ofp, "# max ASCII text line length: %d\n", esl_opt_GetInteger(go, "--textw")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); ++ if (esl_opt_IsUsed(go, "--seq_limit") && fprintf(ofp, "# set max sequence hits to return: %d\n", esl_opt_GetInteger(go, "--seq_limit")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (esl_opt_IsUsed(go, "--popen") && fprintf(ofp, "# gap open probability: %f\n", esl_opt_GetReal (go, "--popen")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (esl_opt_IsUsed(go, "--pextend") && fprintf(ofp, "# gap extend probability: %f\n", esl_opt_GetReal (go, "--pextend")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (esl_opt_IsUsed(go, "--mx") && fprintf(ofp, "# subst score matrix (built-in): %s\n", esl_opt_GetString (go, "--mx")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); +@@ -674,6 +676,13 @@ serial_master(ESL_GETOPTS *go, struct cf + /* Print the results. */ + p7_tophits_SortBySortkey(info->th); + p7_tophits_Threshold(info->th, info->pli); ++ /* Limit the number of hits if specified. */ ++ if (esl_opt_IsOn(go, "--seq_limit")) ++ { ++ int seq_limit = esl_opt_GetInteger(go, "--seq_limit"); ++ info->th->N = ESL_MIN(info->th->N, seq_limit); ++ } ++ + p7_tophits_CompareRanking(info->th, kh, &nnew_targets); + p7_tophits_Targets(ofp, info->th, info->pli, textw); if (fprintf(ofp, "\n\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + p7_tophits_Domains(ofp, info->th, info->pli, textw); if (fprintf(ofp, "\n\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); diff --git a/docs/backend_installation.md b/docs/backend_installation.md index 4336a1db..94ede473 100644 --- a/docs/backend_installation.md +++ b/docs/backend_installation.md @@ -109,6 +109,14 @@ Notes: python -m pip install -e . --no-deps ``` +## Faster jackhmmer (optional) + +The Docker images build HMMER 3.4 with AlphaFold 3's jackhmmer patches: faster FASTA parsing, and threads that each read their own part of the database. To use the same build in a conda environment, overwrite its HMMER: + +```bash +bash docker/hmmer/build_hmmer.sh "$CONDA_PREFIX" +``` + ## Cluster validation On EMBL, the AF2 and AF3 functional tests already have default database roots baked into the test files, so the only environment variable you need for the standard cluster runs is: