thoipapy 3.2.0

Step 1: install DVC

If you created the conda environment from environment.yml you already have it, and can skip to step 2. Otherwise, pick whichever you prefer:

# with conda or mamba (recommended, matches environment.yml)
conda install -c conda-forge dvc dvc-http

# or with pip
pip install "dvc[http]"

dvc-http matters: this project’s public data store is served over HTTPS, and plain dvc cannot read an HTTPS remote without it. If you see URL 'https://...' is supported but requires these missing dependencies, this is what is missing.

Check it worked:

dvc --version

Step 2: fetch the data

From the root of your clone:

dvc pull -j 4

Expect roughly 3800 files, 128 MB, and a minute or two on a normal connection.

Use -j 4. The -j flag sets how many files DVC downloads at once. Its default is four times your CPU count, which is enough parallel requests that the data store starts refusing them, and the pull fails with 429 Too Many Requests. -j 4 avoids this. If you still see 429 errors, lower it further with -j 2.

To confirm everything arrived:

dvc status

Data and pipelines are up to date. means the contents of data/ match the checksums recorded in git, so you have exactly the files the authors had.

Fetching only part of the data

The full set is rarely needed. To fetch one directory, name its .dvc file:

dvc pull -j 4 test/regression_data.dvc  # run the test suite                    28 KB
dvc pull -j 4 data/features.dvc         # retrain the model                      18 MB
dvc pull -j 4 data/results.dvc          # trained models and validation output   25 MB
dvc pull -j 4 data/homologues.dvc       # re-extract features from alignments    60 MB

The remaining directories are small: data/Proteins (420 KB), data/Input_data (272 KB) and data/ETRA_data (112 KB). data/Predictions (26 MB) holds the published predictions and is needed only to compare against them.

Common problems

429 Too Many Requests

Too many parallel downloads. Use dvc pull -j 4, or -j 2 on a fast connection.

URL 'https://...' is supported but requires these missing dependencies

dvc-http is not installed. See step 1.

No DVC remote is specified

You are not in the repository root, or .dvc/config is missing. Run dvc remote list; it should print a remote named public.

dvc: command not found

DVC installed into a different environment than the one you have active. Run conda activate thoipapy first.

The pull is slow or stalls

The data store is rate-limited. This is expected and -j 4 handles it. DVC resumes where it stopped, so re-running dvc pull -j 4 after an interruption is safe and will not re-download what you already have.

How old the alignments are, and what that means

The DVC-tracked homologues were downloaded from NCBI nr in May 2020. The modernisation updated the code, not the data. Every alignment-derived feature in data/features/ – conservation, rate4site, the PSSM columns, coevolution, LIPS, and the polarity features that read from the PSSM – descends from that download, and so does the shipped model.

The provenance is recorded per protein and can be checked:

tar -xzOf data/homologues/xml/crystal/1orqC4.surr20.BLAST.xml.tar.gz --wildcards '*details*'
# acc  1orqC4
# download_date  20200518
# database  ncbi_nr

Regenerating the alignments against current databases is a reasonable thing to want, and the pipeline supports it. It has not been shown to produce a better model. Two rebuilds bear on this, one from a September 2026 UniRef90 search and one from the ColabFold MSA server, and neither improved on the 2020 archive. Scores are averaged over the 34 CD-HIT clusters rather than the 40 proteins, because several of the proteins are homologues of each other:

Every interval contains zero, and ColabFold minus UniRef90 is -0.005 [-0.026, +0.016]. A failure to reject is not proof of equivalence, so the honest statement is what the intervals exclude: any ColabFold advantage over UniRef90 above +0.016 AUC, or over nr above +0.012, is ruled out at 95%, while an improvement below about 0.03 AUC would be too small for 40 proteins to detect.

Depth was then tested directly, with everything else held constant: the same server, the same search, the same query and the same folds, differing only in whether the MMseqs2 diversity filter runs. A 2.6x increase in total alignment depth moved ROC AUC by +0.0001 [-0.016, +0.016]. That comparison carries no confound of time or algorithm, and it is the clearest evidence that depth by itself is not the constraint. It does not show that no better alignment exists – a differently built or differently filtered one has not been tried.

The features that disagree most between alignment sources are the coevolution scores, and those are among the features the model relies on most. DImax correlates at Spearman 0.45 between the nr and UniRef90 feature sets and 0.26 between nr and ColabFold, against 1.00 for every sequence-derived feature. A signal that can be rewritten that far without moving accuracy is where the next accuracy work belongs – starting with MI_highest_face, which a review of this comparison found to be a constant heptad mask rather than a coevolution feature at all. The full analysis, including the limits of the comparison, is in docs/homologue_source_comparison.md.

Two consequences worth keeping in mind:

• Do not assume a refresh is an upgrade. If you regenerate, re-run the validation and compare, rather than adopting the new model because its alignments are deeper.

• Do not mix vintages. A feature table with some proteins from 2020 nr and others from a current search is not a coherent dataset. Regenerate all of them or none.

scripts/compare_blast_database_depth.py performs the comparison above between any two data directories, and is how those numbers were produced. The full write-up, including a negative result on protein language model embeddings, is in thoipapy/docs/embeddings_and_alignment_depth.rst.

Where the data comes from

dvc pull fetches from a publicly readable Cloudflare R2 bucket, declared in .dvc/config. It is read-only and anonymous by design: anyone can fetch the data, nobody can overwrite it. Maintainers push with a separate credentialed remote configured in .dvc/config.local, which is not tracked by git.

The data originates from the dataset published with the paper, archived at the Open Science Framework: https://osf.io/txjev/. The OSF archive is the citable record of what the paper used and does not change. The DVC store is the working copy the pipeline reads, and it tracks the code. Where the two differ, OSF is authoritative for the published results.

The protein sets

data/sets/ holds the three protein sets published with the paper, matching data/protein_lists/ in the OSF repository:

set05 is the union of set07 and set08. The copies here carry a few extra annotation columns beyond the OSF versions; the shared columns are identical.

These were Excel workbooks and are now CSV. They are single flat tables, so the workbook format bought nothing and made them undiffable in git and awkward to read without Excel. The same applies to data/protein_names.csv and the ETRA scanning-mutagenesis data in data/ETRA_data/Average_with_interface/. Pipeline outputs under data/results/ remain Excel where a file genuinely holds several tables; the ones that did not are CSV.

Why this changed

THOIPA was built in 2017-2020, when a BLAST query to NCBI returned in a few minutes. The sequence databases have grown enormously since: nr now holds over a billion sequences, and NCBI queues automated clients accordingly. A probe of a three-residue query in August 2026 returned RTOE=11165 – NCBI’s own estimate of 3.1 hours before results would be ready. A prediction that took minutes when the paper was published can now take hours, or be refused outright if the same client queries repeatedly.

That is not a fault in THOIPA and it is not something THOIPA can fix from the client side. NCBI also does not archive old nr releases, so the May 2020 snapshot behind the published numbers cannot be recovered by anyone. nr is the historical reference, not an option.

Retrieving homologues from the ColabFold MSA server

The server runs MMseqs2 with three profile iterations against UniRef30 and an environmental database assembled from BFD, MGnify, MetaEuk and SMAG, then expands each matched cluster back to its members. A search takes 10 to 70 seconds per protein against hours of NCBI queue, needs no local database, and returns an archive containing msa.sh – the exact command line and database versions behind the alignment. The BLAST path recorded no equivalent.

export THOIPA_COLABFOLD_CONTACT_EMAIL=you@example.com

Then switch homologue_source to colabfold in the settings CSV and enable the run_retrieve_homologues_from_colabfold and run_parse_colabfold_a3m_into_csv stages. They replace the blastp search and the XML parse; everything downstream is unchanged, because the homologue CSV they write carries the same columns.

The standalone predictor uses the same setting: run_THOIPA_prediction reads homologue_source from thoipapy/setting/standalone_run_settings.csv and takes the same ColabFold path, so a single prediction no longer has to wait on an NCBI queue either.

A deployment selects its source with an environment variable instead of by editing a file inside the installed package:

export THOIPA_HOMOLOGUE_SOURCE=colabfold
export THOIPA_COLABFOLD_CONTACT_EMAIL=you@example.com

THOIPA_HOMOLOGUE_SOURCE overrides the settings file when set. THOIPA_COLABFOLD_CONTACT_EMAIL is required on that path and a prediction refuses to start without it, because the MSA server asks automated clients to identify themselves. The shipped default remains ncbi, so an upgrade does not change what an existing installation returns.

Licensing and attribution

Nothing here restricts commercial use. ColabFold and MMseqs2 are MIT licensed, and the databases distributed at colabfold.mmseqs.com state that “all files are available under a Creative Commons Attribution 4.0 International (CC BY 4.0) License”. UniProt, from which UniRef30 derives, is also CC BY 4.0. All of these permit commercial use provided the sources are attributed, which means citing ColabFold, MMseqs2 and the databases in anything published from THOIPA output.

The public server publishes no terms of service, and the constraint on it is a fair-use expectation rather than a licence: the compute is donated by KOBIC and the Söding Lab, the operators describe it as a limited shared resource, and they reserve the right to limit access case by case. Treat it as a courtesy that can be withdrawn rather than a service with a guarantee, and keep the self-hosting route below available for anything that outgrows it.

api.colabfold.com is a free academic service that processes a few thousand alignments a day and asks for serial queries from a single IP. A 40-protein set is three orders of magnitude inside that, and this webserver’s traffic has never come near it, so the public server is the default rather than a stopgap. THOIPA submits one protein at a time and never searches in parallel. A deployment with steady traffic should host its own MMseqs2 server and point THOIPA_COLABFOLD_HOST at it; nothing else changes. Self-hosting is not a way to save the disk, though: the two ColabFold databases are 103 GB and 118 GB as downloads, and ColabFold’s own setup guidance puts them nearer 1 TB once built with indexes, against 96 GB for the local UniRef90 BLAST database. Self-hosting buys independence from a shared service, and costs roughly an order of magnitude more disk than the thing it replaces.

Fallback: a local UniRef90 BLAST database

Use this if the MSA server is unavailable, or if you need the pipeline to run with no third-party service in the path at all. It costs about 96 GB of disk and a rebuild every UniProt release, and it is indistinguishable from ColabFold on accuracy (-0.005 ROC AUC, see above), so the disk buys availability rather than quality.

Either way this is a breaking change: predictions will not match the published results. A different database yields a different homologue set, a different alignment, and therefore different conservation, coevolution and polarity features. The scores move.

Setting up UniRef90

UniRef90 clusters UniProt at 90% identity: 121 million clusters against nr’s billion-plus sequences, with little practical loss for this pipeline, which collapses exact-duplicate TMD sequences and runs CD-HIT anyway.

Requirements: about 160 GB of free disk and BLAST+. The script keeps the 32 GB download, the decompressed FASTA and the built database side by side, and deletes none of them, so peak usage is all three at once.

mamba install -c bioconda blast

# Downloads UniRef90 and runs makeblastdb. Idempotent: an existing download is not refetched.
# Expect several hours for the download and 1-2 hours to build.
python scripts/build_local_blast_db.py

export THOIPA_LOCAL_BLAST_DB=/mnt/shared/blastdb/uniref90

Or by hand:

wget https://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz
gunzip uniref90.fasta.gz
makeblastdb -in uniref90.fasta -dbtype prot -out uniref90 -title uniref90 -parse_seqids

With THOIPA_LOCAL_BLAST_DB set, predictions search that database. Unset, they fall back to querying NCBI, which requires THOIPA_NCBI_CONTACT_EMAIL because NCBI expects automated clients to identify themselves. That path is kept for reproducing old runs and is not recommended for new ones: see the queueing times above.

A smaller database for tests

Set DATABASE_TO_BUILD = "swissprot" in the script for a 575k-sequence, ~340 MB database that builds in ten seconds and answers in under one. Use it for tests, CI and development.

Do not use Swiss-Prot for production. It yields roughly a quarter of nr’s alignment depth: 26 hits against 76 unique TMD homologues for 1xioA4, 18 against 101 for 4hksA1, and only 10 survived filtering for glycophorin A. Conservation, rate4site and coevolution rank 4th, 2nd and 6th by mean decrease in impurity, and all derive from that alignment.

A note on database choice

A database filtered to membrane proteins does not help. THOIPA needs homologues of the query protein itself – glycophorin A’s useful hits are glycophorins in other species – so what matters is taxonomic depth within the query’s own family, not enrichment for membrane proteins. Filtering on a transmembrane annotation would also silently drop genuine homologues that merely lack the annotation, and alignment depth is precisely what must not be lost.

NCBI does not host UniRef90, so the remote BLAST path cannot use it. NCBI’s BLAST service offers nr, refseq_protein, swissprot, pdb, env_nr, tsa_nr, landmark and pataa; of those only nr has comparable depth, and it is the one that queues for hours. This is why the homologue search moved to the ColabFold MSA server rather than to another NCBI database.