This is the repository with the implementation of DeepEST, a multimodal method to perform protein function prediction for bacterial species. Given the characteristics of their genome, the functional characterization of bacteria cannot rely solely on protein sequences. In fact, it requires the use of data sources capable of capturing different dimensions of the protein functionality, i.e., protein structures and expression-location patterns.
The architecture of DeepEST, visualized in the following figure, is comprised of two main modules:
- the fine-tuned structure module [1], and
- the expression-location module.
Subsequently, the two modalities are integrated through learnable linear combination.
Contents: Install · Quick start · Setup · Reproduce results · Repository map · Preprocessing · Citation · References
git clone https://github.com/BorgwardtLab/DeepEST
cd DeepEST
conda env create -f environment.yml -p ./envs/deepest
conda activate ./envs/deepest
pip install -e .This recreates the environment with the pinned packages used to reproduce the results in the paper: Python 3.10.20, torch 2.0.0+cu118, pytorch-lightning 1.9.4, scikit-learn 1.0.2, numpy 1.21.5 and pandas 1.4.4, plus matplotlib, scipy, openpyxl, captum, biopython and joblib for the analysis code.
environment.yml pins the CUDA build of torch, so the environment is about 10 GB on disk and takes
roughly fifteen minutes to build. A CUDA GPU is used when available. The code can be run on CPU as
well and is sufficient for the example below.
The example data for Borrelia burgdorferi B31, the smallest of the 25 species at 423 proteins, is in the repository. External data (go.obo file and DeepFRI model) need to be fetched with the code below. Run everything from the repository root:
scripts/fetch_reference_data.sh # GO 2022-10-07 and the DeepFRI model, into data/This downloads about 900 MB (the DeepFRI release is one 872 MB tarball, of which only the 142 MB
BP model is kept) and leaves roughly 180 MB in data/. Both downloads are checksum-verified, and
re-running the script skips whatever is already present and verified.
Then you can train a DeepEST model on the data fold 0:
deepest train \
--split 0 \
--splitdir examples/borrelia_b31/splits/ \
--config examples/borrelia_b31/config.yaml \
--expr-loc examples/borrelia_b31/expr_loc.csv \
--structures examples/borrelia_b31/structures.pkl \
--label examples/borrelia_b31/labels.csv \
--conversion_dict examples/borrelia_b31/conversion_dict.txt \
--outdir predictions/ \
--species Borrelia_burgdorferi_B31This will take about 30 seconds of training on a CPU, about a minute end to end. It should report a gene-centric F-max of 0.6730 and a micro term-centric
AUPRC of 0.7326. Precision might vary with hardware and thread count.
examples/borrelia_b31/expected_metrics.json holds the reference values, the tolerance and the machine
they were measured on. Note that this is only a quick test for one species and one fold. They are not paper
results, which average over 25 species and 5 folds.
The published data record is at 10.5281/zenodo.18695434.
| File | Bytes | Contents |
|---|---|---|
features.zip |
17,769,556 | expression-location features per species |
matrix_labels.zip |
2,779,968 | GO label matrices |
all_proteins.pkl |
410,039,073 | DeepFRI structure embeddings, 66,163 proteins |
splits_foldseek.zip |
1,191,071 | structure-aware cross-validation splits |
protein_sequences.zip |
74,443,617 | source sequences, only to regenerate embeddings |
structure.zip |
3,207,332,351 | AlphaFold structures, only to regenerate embeddings |
A training run needs nine inputs:
| Input | From |
|---|---|
all_proteins.pkl |
Zenodo |
matrix_labels/label_matrix_<species>_thr10_new.csv |
Zenodo, matrix_labels.zip |
splits_foldseek/<species>/ |
Zenodo, splits_foldseek.zip |
expr-loc/<species>_exprloc.csv |
derived from Zenodo features.zip, step 4 |
conversion_dicts/conv_<species>.txt |
derived from Zenodo features.zip, step 4 |
train_model/configs/foldseek/config_combined_<species>.yaml |
this repository |
deepfri_terms.pkl |
this repository, at assets/, copied into place in step 3 |
go.obo |
downloaded in step 3 |
deepfri_model.hdf5 |
downloaded in step 3 |
Run these from the repository root, with <zenodo> the directory holding the downloaded record:
# 1. unpack
mkdir -p data
unzip -q <zenodo>/features.zip -d data -x '__MACOSX/*'
unzip -q <zenodo>/matrix_labels.zip -d data -x '__MACOSX/*'
unzip -q <zenodo>/splits_foldseek.zip -d data -x '__MACOSX/*'
cp <zenodo>/all_proteins.pkl data/
# 2. verify
md5sum data/all_proteins.pkl # 1a3faa70097538bd61816ddd2664357a
# 3. download go.obo and the DeepFRI model into the same directory
scripts/fetch_reference_data.sh
# 4. derive expr-loc/ and conversion_dicts/ for all 25 species
python experiments/prepare_all_species.py
python scripts/verify_zenodo_derivations.pydata/ now holds every input, and nothing else needs configuring as long as you run from the
repository root.
Only needed if data/ is somewhere else. Two directories vary between machines: the Zenodo download
and where the baseline run writes. Each resolves in this order, first hit wins:
- an environment variable,
DEEPEST_DATAorDEEPEST_BASELINE_OUT user.yamlat the repository root- a repository-relative default,
<repo>/dataor<repo>/baseline_run
cp user.example.yaml user.yaml # set `data:` to wherever you unpacked it
python deepest_paths.py # prints what is currently in effectRelative paths in user.yaml resolve against the repository root rather than your working directory.
user.yaml is gitignored.
One species, one fold:
deepest train --species Borrelia_burgdorferi_B31 --data data--split defaults to 0, --outdir to predictions/, and --data-root to the current directory,
which is where the config finds data/go.obo, data/deepfri_model.hdf5 and data/deepfri_terms.pkl.
Any path can be given explicitly, and an explicit flag wins:
deepest train --species Borrelia_burgdorferi_B31 --data data \
--split 3 --label /elsewhere/labels.csvThe long form naming all six paths:
deepest train \
--split 0 \
--splitdir data/splits_foldseek/Borrelia_burgdorferi_B31 \
--config train_model/configs/foldseek/config_combined_Borrelia_burgdorferi_B31.yaml \
--expr-loc data/expr-loc/Borrelia_burgdorferi_B31_exprloc.csv \
--structures data/all_proteins.pkl \
--label data/matrix_labels/label_matrix_Borrelia_burgdorferi_B31_thr10_new.csv \
--conversion_dict data/conversion_dicts/conv_Borrelia_burgdorferi_B31.txt \
--data-root . \
--outdir predictions/ \
--species Borrelia_burgdorferi_B31The full table is 125 runs, 25 species by splits 0 to 4. Each species has its own config, because the
number of stress conditions differs between them (expression.input_dim). The batch scripts that ran
these sweeps were written for one cluster and are not part of this release.
| Path | Contents |
|---|---|
src/deepest/ |
the model and training code, which is what pip install installs |
examples/borrelia_b31/ |
the minimal example and its expected output |
train_model/ |
per-species configs and split indices. The .py files re-export deepest under the old module names |
preprocessing/ |
scripts that build the model inputs from raw structures, genomes and RNA-seq |
baselines/ |
Foldseek, DIAMOND/BLAST and DeepFRI comparisons. Re-running them needs foldseek, diamond and blastp on your PATH, plus protein_sequences.zip and structure.zip from Zenodo. The published summary is committed at baselines/baselines_summary.csv, which compare_to_paper.py uses when no fresh run is present |
analysis/ |
the scoring scripts behind the manuscript's analyses |
experiments/ |
the long-running jobs that produce those caches: training ablations, input staging, interpretability sweeps |
analysis/data/ |
cached results of those jobs, so a score can be re-read without re-running them |
scripts/ |
prediction comparison and the Zenodo derivation check |
tests/ |
import, path, scheduler, metric and end-to-end checks |
annotations_hypothetical_proteins/ |
DeepEST predictions for about 7,000 unannotated proteins |
trained_models/ |
where training writes checkpoints, created at the repository root on the first run. About 61 MB per run. Not committed, and runs accumulate rather than overwrite |
assets/ |
phylogeny, contig-simulation inputs, the PATRIC protein-family table |
deepest_paths.py |
resolves the two machine-specific directories. Run it to see what is in effect |
user.example.yaml |
template for user.yaml, the one file a new user may need to edit |
analysis/ holds the scripts behind the manuscript's analyses. Each score_*.py answers one
question against the trained models and writes the answer into analysis/data/:
| Script | Question |
|---|---|
score_noloc_all25.py, score_modality_ablation.py |
what each modality contributes, across all 25 species |
score_cond_ablation.py, score_cond_ablation_genes.py |
what removing a stress condition costs, per species and per gene |
score_contig_fulltest.py, score_contig_perfold.py, score_contig_sl1344.py |
robustness to fragmented draft assemblies |
depth_analysis.py |
the location ablations stratified by GO depth |
analysis/deepest_interpret.py is the shared integrated-gradients library the condition-attribution
scripts build on. The long-running jobs that produce their inputs live in experiments/.
preprocessing/ builds the model inputs from raw structures, genomes and RNA-seq.
NOTE This is not needed to reproduce the paper.
The Zenodo record contains the preprocessing pipeline's outputs. These scripts are only needed for adapting the pipeline to another species.
| Directory | Builds |
|---|---|
scripts/ |
the per-species GO label matrices, the cross-validation splits, and the genomic-location features |
structure_clustering/ |
the Foldseek structure-aware splits behind the paper's main results, clustered at LDDT 0.7. Needs a Foldseek binary, which is not vendored here |
gene_function/ |
the 3x512 DeepFRI-GCN embedding per protein that becomes all_proteins.pkl, plus the AlphaFold and UniProt downloaders |
Each species has its own config under train_model/configs/foldseek/. Three data paths in it
are read at runtime and resolve against --data-root:
| Field | File |
|---|---|
model_params.go_fn |
data/go.obo |
structure.model_file |
data/deepfri_model.hdf5 |
structure.deepfri_terms |
data/deepfri_terms.pkl |
Every other path in the config, and the whole sequence: block, is legacy and unused.
src/deepest/lr_scheduler.py is vendored from lightning-flash under Apache-2.0. It is one class,
copied byte-faithfully with its original attribution. See NOTICE. It replaces a dependency
on a package old enough to be an install risk. tests/test_lr_scheduler.py checks the vendored copy
against a reference sequence captured from upstream and committed as
tests/lr_schedule_reference.json. The step-by-step comparison against lightning-flash itself runs
only if that package happens to be installed, which environment.yml deliberately does not do.
If you use DeepEST in your research, please cite:
@article{muzio2024bacterial,
title={Bacterial protein function prediction via multimodal deep learning},
author={Muzio, Giulia and Adamer, Michael and Fernandez, Leyden and Borgwardt, Karsten and Avican, Kemal},
journal={bioRxiv},
pages={2024--10},
year={2024},
publisher={Cold Spring Harbor Laboratory}
}See LICENSE.md for licensing.
Expression: PATHOgenex [2]. Structures: AlphaFold DB [3]. Annotations: UniProt [4], accessed 12 July 2023, and the GO ontology released 7 October 2022.
[1] Gligorijević, V., et al. Structure-based protein function prediction using graph convolutional networks. Nature Communications 12, 3168 (2021).
[2] Avican, K., et al. RNA atlas of human bacterial pathogens uncovers stress dynamics linked to infection. Nature Communications 12, 3282 (2021).
[3] Varadi, M., et al. AlphaFold Protein Structure Database: massively expanding the structural coverage of protein-sequence space with high-accuracy models. Nucleic Acids Research 50(D1), D439-D444 (2022).
[4] The UniProt Consortium. UniProt: the Universal Protein Knowledgebase in 2023. Nucleic Acids Research 51(D1), D523-D531 (2023).