DeepEST

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:

1. the fine-tuned structure module [1], and
2. 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

---

Install

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.

Quick start

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_B31

This 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.

Setup

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.py

data/ now holds every input, and nothing else needs configuring as long as you run from the repository root.

User specific paths

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:

1. an environment variable, DEEPEST_DATA or DEEPEST_BASELINE_OUT
2. user.yaml at the repository root
3. a repository-relative default, <repo>/data or <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 effect

Relative paths in user.yaml resolve against the repository root rather than your working directory. user.yaml is gitignored.

Reproduce results

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.csv

The 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_B31

The 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.

Repository map

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 scripts

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

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

Config fields

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.

Third-party code

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.

Citation and licence

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.

Data sources

Expression: PATHOgenex [2]. Structures: AlphaFold DB [3]. Annotations: UniProt [4], accessed 12 July 2023, and the GO ontology released 7 October 2022.

References

[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).