Make the pipeline's outputs reproducible

[mrubash1]

Running the pipeline twice on identical inputs produces different maps. There are three independent causes, all pre-existing on main. This PR fixes each.

This came out of designing the validation for a possible foldseek version bump: a before/after output diff turns out to have no power, because the outputs already differ between two runs of the same version.

1. PCA was unseeded and used a randomized solver

dim_reduction.py constructed PCA(n_components=n_components, **kwargs) — no random_state, and svd_solver left at its default of 'auto'. main() passes random_state to calculate_TSNE and calculate_UMAP but never to calculate_PCA.

'auto' selects the randomized solver once the matrix exceeds 500 rows or columns, and that solver draws its projection matrix from the global numpy random state:

N= 500 -> svd_solver = full
N= 501 -> svd_solver = randomized

Both configured plotting modes (pca_tsne, pca_umap) run a 30-component PCA and feed the result into t-SNE/UMAP, so seeding those two did not make anything reproducible — their input was changing. Measured on a 600-protein matrix:

two unseeded PCA runs identical? False
max abs diff across runs      : 2.52
std of last PC                : 0.35

The trailing components are several times their own standard deviation of pure run-to-run noise, and all 30 are handed downstream.

Fixed by threading the existing random_state through and requesting svd_solver="full", which is exact and deterministic. Verified afterwards that two runs are bit-identical and that the result no longer depends on column order (which matters for cause 2).

2. The similarity matrix's column order was hash-randomized

reading_data returned targets as a set, and pivot_foldseek_results writes the header and every row by iterating it. Python salts string hashing per process, so the column order changed on every run, and nothing sets PYTHONHASHSEED.

Rows were already sorted() while columns were not, so all_by_all_tmscore_pivoted.tsv was a matrix whose row i and column i were different proteins, with the 1.0 diagonal scattered through it.

Harmless on its own for Euclidean row distances, but not harmless in combination with cause 1, because the randomized solver is sensitive to column order. Fixed by sorting, which also makes the column order match the row order.

3. Which proteins reached the map was decided by hash randomization

aggregate_hits accumulates accessions into a set and writes them in iteration order. download_pdbs then truncates that file with accessions[:max_structures] (default 5000).

So whenever there were more hits than max_structures — the situation the parameter exists for — the proteins that survived were an arbitrary sample chosen by hash order. Same accessions, three processes:

seed=1 -> first 3 kept: ['A0A286Q506', 'P60713', 'D7RIF5']
seed=2 -> first 3 kept: ['D7RIF5', 'P60713', 'A0A286Q506']
seed=3 -> first 3 kept: ['A0A850ZFV5', 'Q6QAQ1', 'P60713']

Fixed by sorting.

This makes the choice reproducible, not principled. The retained hits are now the first max_structures in accession order, which is stable but still not biologically meaningful. Ranking by significance (best e-value across queries for BLAST hits, TM-score for Foldseek hits) would be the better rule, but it changes which proteins appear in a map rather than only making the existing selection stable, so it belongs in its own PR with its own discussion.

Verification

• PCA: two runs bit-identical, and invariant to a random column permutation (max diff < 1e-9).
• Ordering: both fixes produce identical output across PYTHONHASHSEED 1/2/3.
• Unit tests pass; cluster-mode DAG resolves unchanged (16 rules).
• ruff check, ruff format --check, snakefmt --check clean under Python 3.9.

What was not verified

I have not measured the end-to-end effect on a real map, i.e. how much a published figure changes between two runs today. The PCA measurement above is on a synthetic 600×600 matrix. The direction is not in doubt, but the magnitude on real data is unquantified.

Rebased onto main after #103 merged. The conflict was confined to dim_reduction.py; #103's structure is kept and the seed added on top. calculate_PCA on main still takes no random_state, so the defect this PR fixes is live and unaffected by #103.

test_pipeline_in_cluster_mode and the unit tests pass on the rebased branch, as does lint under CI's pinned versions. (The earlier note here, that make test failed on this base for the envs/analysis.yml drift, no longer applies: #103 pinned matplotlib and setuptools, and main is green.)

Implications

• Any map of more than ~500 proteins cannot currently be regenerated exactly from its inputs, so figures cannot be reproduced from archived inputs alone.
• Any before/after comparison — a foldseek bump, a BLAST database change, a dependency update — is confounded until this lands. That includes the nr vs refseq_protein question in Fix search-mode reliability for BLAST, UniProt, AFDB, and small maps #103.
• After this, the right check is a self-diff: run the pipeline twice at the same versions and confirm the outputs are identical. That is the precondition for any version-bump validation being meaningful.

[ahmedhosny]

Looks good. Three independent sources of run-to-run map variation, each with a matching fix:

1. PCA now gets the existing random_state and svd_solver="full", so maps with >500 proteins no longer draw a randomized projection from the global numpy RNG. calculate_PCA is only called from main() in this file, and the new argument is threaded the same way as t-SNE/UMAP.
2. Similarity-matrix columns are sorted so they match the already-sorted rows instead of depending on PYTHONHASHSEED.
3. Aggregated accessions are sorted before download_pdbs truncates to max_structures, so which proteins survive is stable.

Agree that (3) makes the truncation reproducible rather than biologically ranked; accession order is the right change for this PR, and ranking by e-value/TM-score belongs in a follow-up.