WASP: Protein Functional Annotation using AlphaFold structures

Welcome to the official repository for the paper WASP: A pipeline for functional annotation based on AlphaFold structural models!

WASP, Whole-proteome Annotation via Structural-homology Prediction, is a python-based software designed for comprehensive organism annotation at the whole-proteome level based on structural homology.

WASP is a user-friendly command-line tool that only requires the NCBI taxonomy ID of the organism of interest as an input. Using the computational speed of Foldseek [1], WASP generates a graphical representation of reciprocal hits between the organism protein query and the AlphaFold database [2, 3], enabling downstream robust functional enrichment and statistical testing. WASP annotates uncharacterised proteins using multiple functional descriptors, including GO terms, Pfam domains, PANTHER family classification and CATH superfamilies, Rhea IDs and EC numbers. Additionally, WASP provides a module to map native proteins to orphan reactions in genome-scale models based on structural homology.

Table of Contents

1. Installation
   • External requirements
   • Quickstart
2. Usage
   • Whole-proteome annotation
   • GEM gap-filling module
3. Additional utils
   • Download custom proteins
   • Predict custom strucrures
4. References

1. Installation

To ensure full reproducibility and easy setup, WASP packages all external dependencies (Foldseek, gsutil) via Conda.

✅ Recommended: Install via Conda

Ensure you have Conda or Miniconda installed. This will automatically set up Python 3.10, Foldseek 9 or later, gsutil, and the WASP pipeline in an isolated environment.

# Clone the repository
git clone https://github.com/gioodm/wasp-proteins-annotation.git
cd wasp-proteins-annotation

# Create the environment and install all dependencies
conda env create -f environment.yml

# Activate the environment
conda activate wasp-env

The manuscript results were obtained using Python 3.10.14 and Foldseek 8-ef4e960.

2. Usage

2.1 Whole-proteome annotation

Usage

First, identify the NCBI taxonomy ID of your organism of interest. You can find the ID at NCBI Taxonomy. You can use AlphaFold DB to check how many structures are linked to that ID by searching the same ID on the AFDB search bar.

Once you have identified the taxid, run the pipeline as follows (e.g., for organism S. cerevisiae S288c, taxid: 559292):

wasp-run -t 559292

This requires gsutil installed.
Note: On the first run, the AlphaFold DB clustered at 50% will need to be downloaded using Foldseek. This can take some time depending on your machine's performance and internet speed. Ensure you have sufficient storage space to host the database. After the initial setup, WASP annotation will take up to 3 hours for iteration for a proteome of approximately 6000 proteins.

To annotate proteins without predicted structures, provide a FASTA file instead. WASP downloads the ProstT5 weights required by Foldseek on the first FASTA run and performs a structure-based sequence search. This mode does not require gsutil:

wasp-run --fasta proteins.fasta

Use exactly one input option: --taxid, --fasta, or --structures. The output directory is named after the input filename without its extension when using a FASTA file or structure archive.

Additional parameters can be customised, including:

• -e evalue_threshold: set the evalue threshold (default: 10e-10)
• -b bitscore threshold: set the bitscore threshold (default: 50)
• -n max_neighbours: set the max number of neighbours (default: 10)
• -s step: set step to add to max neighbours (n) in additional iterations (default: 10)
• -i iterations: set number of iterations to perform (default: 3)
• -m max sequence identity: set maximum sequence identity (0.0 to 1.0) to exclude sequence homologs (default: 1.0)

Usage examples:

wasp-run -e 1e-50 -b 200 -n 5 -i 5 -t 559292
wasp-run -s 5 -t 559292
wasp-run -f proteins.fasta -e 1e-50 -b 200

To use a custom structural dataset (e.g., a newly sequenced genome or a set of proteins from different organisms), create a tar archive containing protein structures in .cif.gz or .pdb.gz format. Pass its path with --structures; no gsutil download is performed:

wasp-run -p proteins.tar
# equivalent to: wasp-run --structures proteins.tar

2.2 GEM gap-filling module

Pre-processing

Some pre-processing steps are required to obtain a standardized input file for the WASP pipeline. Potential modifications to the Python scripts might be needed depending on the GEM format.

• find-orphans: identifies orphan reactions in the GEM (accepted extensions: .xml, .sbml, .json, and .mat) using the Python3 cobrapy module. Annotation in different formats (accepted: BiGG, Rhea, EC number, KEGG, PubMed, MetaNetX) present in the model is retrieved - input example at: example_files/gap_filling/iYLI649.xml, output example at: example_files/gap_filling/iYLI649_orphans.txt.

• run-rxn2code: each reaction annotation is mapped to the corresponding Rhea reaction ID and/or EC number when available. The output file includes: 1. reaction id, 2. reaction extended name, 3. Rhea IDs, 4. EC numbers. If MetaNetX codes are present, the reac_xref.tsv file (retrieved from MetaNetX) is needed - input example at: example_files/gap_filling/iYLI649_orphans.txt, output example at: example_files/gap_filling/iYLI649_gaps.txt.

The final gaps_file.txt must contain all 4 columns; if the reaction extended name is not present in the model, an empty column should be present. If no Rhea/EC IDs are identified, empty columns should be present instead.

Usage

wasp-gem [-h] [-e evalue_threshold] [-b bitscore_threshold] [-t tmscore] taxid gaps_file.txt

An example of the final output can be found at example_files/gap_filling/iYLI649_hits.txt

3 Additional utils

3.1 Downloading custom protein structures from AFDB

The afdb-download utility allows you to fetch protein structures from the AlphaFold Database (AFDB) using either a FASTA file containing protein sequences or a list of UniProt IDs.

Usage:

afdb-download [-h] --input_file INPUT_FILE --output_dir OUTPUT_DIR [--num_cores NUM_CORES]

3.2 Predicting structures with ColabFold

A Python script (predict_structures.py) is provided to run ColabFold predictions on a list of UniProt IDs.

Prerequisites:

• ColabFold must be installed and accessible in your PATH.
• Ensure your system meets hardware requirements (GPU recommended for large predictions).

Usage:

python3 predict_structures.py -i <uniprot_ids.txt> -o <output_dir> [-c <cores>]  

Best practices & recommendations

For large-scale predictions, combine proteins into a single FASTA file to reduce overhead, but split into smaller batches (50-100 proteins) to avoid memory issues. Use parallel processing (-c 8 for 8 cores) and submit as a batch job (e.g., SLURM) for heavy workloads.

Pre-download FASTA sequences to avoid failures. Test with one protein first to verify setup. If predictions fail, reduce --num-recycle (default: 3) or check GPU memory.

For custom workflows, modify the script (e.g., add --amber or --templates). See the ColabFold docs for advanced options.

4 References

[1] van Kempen M, Kim SS, Tumescheit C, Mirdita M, Lee J, Gilchrist CLM, et al. Fast and accurate protein structure search with Foldseek. Nature Biotechnology. 2023. doi: https://doi.org/10.1038/s41587-023-01773-0

[2] Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate protein structure prediction with AlphaFold. Nature. 2021;596(7873):583-9. doi: https://doi.org/10.1038/s41586-021-03819-2

[3] Varadi M, Bertoni D, Magana P, Paramval U, Pidruchna I, Radhakrishnan M, et al. AlphaFold Protein Structure Database in 2024: providing structure coverage for over 214 million protein sequences. Nucleic Acids Research. 2023;52(D1):D368–D375. doi: http://dx.doi.org/10.1093/nar/gkad1011