E1 Phase 0: HFSP rebuild, pair deduplication, and a uniform embedding cohort - #5
Open
tsenoner wants to merge 6 commits into
Open
E1 Phase 0: HFSP rebuild, pair deduplication, and a uniform embedding cohort#5tsenoner wants to merge 6 commits into
tsenoner wants to merge 6 commits into
Conversation
Mahlich et al. 2018 (Bioinformatics 34:i304-i312), whose HFSP score this pipeline reproduces, set --alignment-mode 3, --num-iterations 3, --e-profile 1e-10 and -e 1e-3 but never --max-seqs -- i.e. MMseqs2's default of 300. The sequence table on disk is at exactly 300; the script's stated default of 1000 did not describe the data it produced. Foldseek stays at 1000, matching its own table. HFSP is defined only against the sequence search, so binding the structural search to 300 would buy no methodological fidelity while invalidating every alntmscore probe. Both caps bite and both belong in Methods: 234,571 of 542,238 sequence queries (43.3%) sit exactly at 300, contributing 68.7% of all rows, versus 11.8% of queries / 22.9% of rows for Foldseek at 1000. The prefilter ranks by score, so the dropped hits are the remote ones -- the pair set is enriched for high-identity pairs inside large families. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
…l pair
--num-iterations 3 makes the search a profile search, which is directional:
a pair may be reported as (A,B), as (B,A), or as both, with different
alignments. Nothing collapsed them, so 61.5% of unordered sequence pairs
(74.9% of structural ones) carried double weight while the rest carried
single weight.
That weighting is not random. Duplicated pairs are measurably the MORE
similar ones -- mean fident 0.5158 vs 0.4142, mean alntmscore 0.8304 vs
0.6319, median E-value 1.34e-59 vs 8.74e-35. The raw table therefore
inflates exactly the sequence-similarity signal the reviewers questioned.
It also feeds conflicting targets to the symmetric probe arms
(linear_distance, euclidean) for 30.9% of sequence and 44.9% of structural
pairs, while the concatenating arms never see the conflict -- biasing the
capacity ladder.
The two sides get deliberately different rules:
sequence -- keep the single lowest-E-value alignment. HFSP is non-linear
in (PIDE, L), so averaging would fabricate an alignment that
never existed and desync fident/nident/mismatch.
structure -- mean alntmscore. It is a lone scalar, and symmetric by
definition: verified on pairs whose proteins differ in
length, 41.7% carry identical scores both ways versus 12.6%
matching a target-normalised model, so Foldseek #312 does
not apply here. Max would re-bias the duplicated subset up.
--no-dedupe reproduces the legacy directional table and writes it to
merged_protein_similarity_nodedup.parquet, so it can never overwrite the
canonical one and the HFSP delta stays separately attributable.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
The 15 embedding arms do not cover the same proteins, and src/shared/datasets.py:99-105 drops a pair when EITHER protein is missing, so coverage loss is quadratic. esm2_3b was trained and scored on 558,947 test pairs where ten other arms got 872,572, yet is published at rank #10. A ranking is only a ranking if every row was scored on the same data, so this figure is the honest disclosure of the input to that claim. Measured: only 422,972 of 542,238 proteins (78.00%) are in every arm, in five tiers -- 542,238 complete; 542,237 (one outlier, Q9H195, the longest protein in the cohort at 13,477 aa, which ProtT5 could not embed and prottucker/random_1024 inherit); 540,881 (default --max_seq_len 2000); 526,871 (ESM-1b's 1022-token cap, which CLEAN inherits by construction); 435,298 (esm2_3b, an interrupted run, being completed). Once esm2_3b lands the intersection rises to 526,871 and stops there -- ESM-1b's cap then binds, which is why topping arms up can never yield a uniform test set and intersecting is the only construction that can. freeze/embedding_key_coverage.json (2.2 KB) is committed and is the default source, so the figure redraws offline: the .h5 files live on LRZ/Zenodo and enumerating 542k HDF5 keys costs ~3 min per file over GPFS. Same philosophy as the pair index -- the repo pins the manifest, the image is a reproducible build product. Integrity is asserted on write: the patterns must reconstruct counts, universe and intersection exactly. Set maths is separated from rendering and pinned by tests; --h5-dir rebuilds from the HDF5 files, caching a <stem>.keys.txt sidecar stamped with (size, mtime) so a changed file invalidates its own cache. Exposed as `plm figures coverage-upset`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
esm2_3b.h5 holds 435,298 of the cohort's 542,238 proteins. The gap is an
interrupted run, not a filter: present and missing length distributions
are indistinguishable (median 295 vs 296, mean 353.3 vs 354.2, p95 860 vs
864, >1024 aa 2.8% in both) and the cohort's longest protein was embedded.
So an append-mode resume is the correct fix -- _partition_cohort filters
already-present proteins before any forward pass.
Three traps handled:
--max_seq_len 6000: the DEFAULT IS 2000 and silently skips longer
sequences, and the file already holds 1,092 proteins above it. The
missing set maxes at 5,412, so 6000 covers all of them.
no --bf16: the existing 435,298 vectors are float32; mixing precision
inside one arm is the confound the reviewers objected to.
a key-count assertion instead of an exit code: the append writer counts
a skip as a success (embedding_generation.py:965), so a misconfigured
resume writes nothing, logs success and exits 0.
It seeds a working copy so the md5-verified Zenodo deposit is never
appended to. Partition measured with --test-only on 2026-08-05: MIG starts
same-day where H100 was 2026-08-18 and A100-80x4 2026-08-19; shorter
--time moved none of them, so queue position rather than backfill binds.
MIG needs BOTH --qos=mig and the typed --gres=gpu:3g.20gb:1.
Also gitignores docs/manuscript_changes.md, the internal ledger of changes
that must reach the manuscript text (origin is public; same rule as D-8).
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
The 15 arms do not cover the same proteins, and _load_and_filter_data drops a pair when EITHER protein is missing from that arm's HDF5, so the loss is quadratic in coverage. Measured on the published grid, esm2_3b was scored on 558,947 test pairs where ten other arms got 872,572, yet is reported at rank #10. A ranking whose rows were scored on different data is not a ranking. Exclusion rather than completion, because one gap cannot be closed: ESM-1b's learned positional embeddings cap at 1022 tokens, so ~2.8% of the cohort can never be embedded by it and CLEAN inherits exactly that set by construction. Topping arms up can therefore never yield a uniform test set; intersecting is the only construction that can. The single longest protein (Q9H195, 13,477 aa, which ProtT5 could not embed and prottucker and random_1024 inherit) is the smallest instance of the same rule. A load-time filter over a committed id list, NOT deletion from the .h5 files: those are the md5-verified Zenodo deposit, deletion is irreversible, and each file would stop matching its published checksum. The freeze stores the EXCLUDED ids because that half is ~34x smaller (15,367 vs 526,871). No freeze committed yet => empty exclusion => this is a no-op, so adopting the cohort is an explicit act rather than an accident. The list itself is generated once esm2_3b finishes: doing it now would bake in 119,266 exclusions instead of 15,367 and discard 103,899 proteins for nothing. exclusion_summary separates "removed" from "excluded but already absent here", so the filter cannot look like it did more work on the deficient arms than it did. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
One pass over each arm does three jobs at once, because each is a full
rewrite anyway:
CONSOLIDATE 15 arms live in data/zenodo_17469268/ while
run_experiments.py:49 globs <cohort>/embeddings/*.h5 -- it would find 2
files, print "Found 2 embedding files" and report success (blocker B2).
The directory is also named after the VERSION DOI, while the concept DOI
17469267 is the one to cite, and a new version is imminent.
CUT to <=2000 aa the cohort becomes 540,881 of 542,238. Not an invented
threshold: it is the pipeline's own --max_seq_len default, and
esm3_open/esmc_300m/esmc_600m already sit at exactly 540,881.
REPACK `del f[key]` unlinks without reclaiming space, so deleting in
place would ship a new deposit at the old size with dead blocks inside.
Measured result (LRZ, 2026-08-06), now committed as
freeze/embedding_key_coverage_cohort2k.json: 13 arms hold the full 540,881
and the membership structure collapses to TWO patterns -- 526,871 present
everywhere, 14,010 missing only clean/esm1b. Before the cut it was seven
patterns and a 422,972 intersection (78.00%); after, 97.41%.
clean/esm1b stay at 526,871 by design. Cutting to ESM-1b's 1022-token cap
would cost every arm a further 14,010 proteins to accommodate one model's
positional limit, so the discrepancy is documented instead. Because a pair
needs both proteins, 97.41% protein coverage is ~94.9% of pairs for those
two arms -- disclose wherever either appears in a ranking.
esm2_3b is sourced from the completed working copy, not the stale deposit.
Its 51 CUDA-OOM failures are ALL >2000 aa (min 2,414, max 5,412), so the
cut removes exactly the proteins that failed: esm2_3b lands at 540,881.
Also applies the measured HDF5 fix to load_h5_keysets: these are
superblock-v0 symbol-table groups, ~0.52 metadata reads per key at 356 B,
and on GPFS an identical 356 B read costs 3.0 us within an 8.4 MB span vs
175 us across 4.62 GB. The core driver reads the file in one sequential
pass instead -- measured 66.8x on esm1b.h5 (161.7 s -> 2.4 s), key sets
byte-identical. Size-guarded at 4.0e9: a 3.53 GB file peaked at 3.70 GB RSS
under the login node's 4 GiB per-user cgroup while a 4.62 GB one was
SIGKILLed, and it must not run concurrently.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Phase 0 of the resubmission. Rebuilds the pair table on a corrected HFSP score, removes a
systematic weighting bias, and makes all embedding arms cover the same proteins.
What changed
HFSP + search parameters.
--max-seqsis pinned to 300, the value that actually producedthe pair table on disk and the one Mahlich et al. 2018 used (their Methods set
--alignment-mode 3,--num-iterations 3,--e-profile 1e-10,-e 1e-3and never--max-seqs). The script had claimed 1000, which did not describe the data. Foldseek stays at1000, matching its own table.
Pair deduplication. The profile search is directional, so a pair may be reported as (A,B),
as (B,A), or both. Nothing collapsed them, so 61.5% of unordered sequence pairs carried double
weight. That weighting is not random — duplicated pairs are measurably the more similar ones
(mean fident 0.5158 vs 0.4142; median E-value 1.34e-59 vs 8.74e-35), so the raw table inflated
exactly the sequence-similarity signal under review. It also fed conflicting targets to the
symmetric probe arms for 30.9% of sequence and 44.9% of structural pairs while the
concatenating arms never saw the conflict, biasing the capacity ladder.
Sequence side keeps the lowest-E-value alignment whole (HFSP is non-linear in (PIDE, L), so
averaging would fabricate an alignment that never existed); structure side takes the mean
alntmscore, which is symmetric by definition — verified on pairs whose proteins differ inlength, 41.7% carry identical scores both ways versus 12.6% matching a target-normalised model.
Uniform cohort. The 15 arms did not cover the same proteins, and a pair is dropped when
either protein is missing, so the loss is quadratic. The cohort is now proteins of ≤2000
residues (540,881 of 542,238) — the pipeline's own
--max_seq_lendefault, which three armswere already built to. Membership collapses from seven patterns to two: 526,871 present
everywhere, 14,010 missing only clean/esm1b. Intersection 78.00% → 97.41%.
clean/esm1b stay at 526,871 by design; cutting to ESM-1b's 1022-token cap would cost every arm
a further 14,010 proteins to accommodate one model. The discrepancy is documented instead.
Verification
Every predicted number was stated before the run and matched exactly: MMseqs side
68,976,426 retained without dedup and 40,916,356 with; Foldseek 185,683,143 / 98,893,186;
totals 187,450,209 / 99,803,629.
--no-dedupereproduces the legacy table under a distinctfilename so the HFSP delta stays separately attributable.
Suite: 1180 passed, 2 skipped (was 1160). New tests cover the dedup rules, the cohort filter,
and the coverage set maths.
Also here
plm figures coverage-upsetplus a committed coverage freeze, so the figure redraws offline(enumerating 542k HDF5 keys costs ~3 min/file over GPFS).
load_h5_keysetsuses thecoredriver — measured 66.8x on esm1b.h5 (161.7s → 2.4s),size-guarded at 4.0e9 against the login node's 4 GiB per-user cgroup.
🤖 Generated with Claude Code
https://claude.ai/code/session_011cdSkVB9XdhjZP4smvX44Z