Home

Awesome

<!-- omit in toc -->

Vaxformer: Antigenicity-controlled Transformer for Vaccine Design Against SARS-CoV-2

This repository contains pre-trained models, corpora, indices, and code for pre-training, finetuning, retrieving and evaluating for "Vaxformer: Antigenicity-controlled Transformer for Vaccine Design Against SARS-CoV-2"

<img src="assets/vaxformer.png">

The SARS-CoV-2 pandemic has emphasised the importance of developing a universal vaccine that can protect against current and future variants of the virus. However, research in the field of protein design is limited despite the advancements in AI and our knowledge of the human immune system. In this context, the present study proposes a novel conditional protein Language Model architecture, called Vaxformer, which is designed to produce antigenicity-controlled SARS-CoV-2 spike proteins. We evaluate the generated protein sequences of the Vaxformer model using DDGun protein stability measure, netMHCpan antigenicity score, and a structure fidelity score using the root mean square deviation between the protein folded with AlphaFold and a Covid reference protein to gauge its viability for vaccine development. Our results show that Vaxformer outperforms the existing state-of-the-art Conditional Variational Autoencoder model to generate antigenicity-controlled SARS-CoV-2 spike proteins. These findings suggest promising opportunities for Conditional Transformer models to expand our understanding of vaccine design and their role in mitigating global health challenges.

Authors (equal contribution):

<!-- omit in toc -->

Table of Contents

πŸ› οΈ Setup

Python packages

This codebase requires the following dependencies:

- biopython
- matplotlib
- seaborn
- numpy
- pandas
- pydantic
- python-dotenv
- PyYAML
- tqdm
- wandb
- jupyterlab

We opted in to using conda as our package manager. The following will install all necessary dependencies for a GPU training:

ENV_NAME=vaxformer
conda create -n ${ENV_NAME} python=3.8 -y
conda activate ${ENV_NAME}
conda install pytorch torchvision torchaudio pytorch-cuda=11.7 -c pytorch -c nvidia -y
pip install -r requirements.txt

netMHCpan

Note

netMHCpan only runs on Linux or Darwin machine

Follow this step to setup netMHCpan:

  1. Download https://services.healthtech.dtu.dk/services/NetMHCpan-4.1/
  2. Follow their installation instruction outlined in the netMHCpan-4.1.readme file

DDGun

We included DDGun as a submodule to this repo. To install, you need to:

# Clone the DDGun submodule
git submodule update --init --recursive

# Install DDGun
cd submodules/ddgun
python setup.py

Alphafold 2

  1. Install PyMol
  2. Install alphafold2 or use HPC cluster with alphafold installed (recommended)

Dataset

Vaxformer is trained with a dataset of spike Covid proteins from GI-SAID. You have to have the appropriate GI-SAID credentials to download the dataset. To obtain comparable data splits, we inquired to the author of the previous publication (ImmuneConstrainedVAE).

⌨️ Codebase Structure

.
β”œβ”€β”€ configs                                       # Config files
β”‚Β Β  β”œβ”€β”€ test/                                     # Config files for sampling runs
β”‚Β Β  └── train/                                    # Config files for training runs
β”œβ”€β”€ datasets/                                     # Datasets of sequences and immunogenicity scores
β”œβ”€β”€ scripts                                       # Scripts to start runs
β”‚Β Β  β”œβ”€β”€ netmhcpan/
β”‚   β”‚Β Β  β”œβ”€β”€ netmhcpan_allele_scores_one_file.sh   # Script to run netMHCpan scoring for multiple peptide files
β”‚   β”‚Β Β  β”œβ”€β”€ netmhcpan_allele_scores.sh            # Script to run netMHCpan scoring for one peptide file
β”‚   β”‚Β Β  β”œβ”€β”€ compute_hits_from_peptides.py         # Script to compute netMHCpan hits and scores from peptides
β”‚   β”‚Β Β  └── generate_peptides_from_sequences.py   # Script to generate peptides from sequences
β”‚Β Β  β”œβ”€β”€ evaluation/
β”‚   β”‚Β Β  β”œβ”€β”€ alphafold.sh                          # Script to run the AlphaFold evaluation
β”‚   β”‚Β Β  β”œβ”€β”€ ddGun.sh                              # Script to run the DDGun evaluation
β”‚   β”‚Β Β  β”œβ”€β”€ create_data_for_ddgun.ipynb           # Notebook to create the necessary data for DDGun run
β”‚   β”‚Β Β  β”œβ”€β”€ dd_eval.ipynb                         # Notebook to read the DDGun evaluation and find the top 10 proteins
β”‚   β”‚Β Β  └── evaluation.ipynb                      # Notebook for PCA evaluation
β”‚Β Β  β”œβ”€β”€ baseline/
β”‚   β”‚Β Β  └── naive_bayes.ipynb                     # Notebook for 
β”‚Β Β  β”œβ”€β”€ slurm/                                    # Slurm scripts for training and sampling runs
β”‚Β Β  β”œβ”€β”€ sample.py                                 # Script to run sampling with a model of choice
β”‚Β Β  └── train.py                                  # Script to run training with a model configuration of choice
β”œβ”€β”€ src
β”‚   β”œβ”€β”€ dataset
β”‚   β”‚Β Β  β”œβ”€β”€ __init__.py
β”‚   β”‚Β Β  β”œβ”€β”€ dataset.py                            # Dataset class to prepare and iterate through the dataset
β”‚   β”‚Β Β  └── tokenizer.py                          # Tokenizer class to preprocess the input sequence
β”‚   β”œβ”€β”€ models 
β”‚   β”‚Β Β  β”œβ”€β”€ __init__.py
β”‚   β”‚Β Β  β”œβ”€β”€ baseline.py                           # Naive Bayes baseline model
β”‚   β”‚Β Β  β”œβ”€β”€ lstm.py                               # Conditional LSTM model
β”‚   β”‚Β Β  β”œβ”€β”€ vae.py                                # Conditional VAE model
β”‚   β”‚Β Β  └── vaxformer.py                          # The proposed Vaxformer model
β”‚   β”œβ”€β”€ utils
β”‚   β”‚   β”œβ”€β”€ __init__.py
β”‚   β”‚   β”œβ”€β”€ common_utils.py                       # Utility functions to prepare trainings
β”‚   β”‚   └── model_utils.py                        # Utility functions for modelling purposes
β”‚   β”œβ”€β”€ __init__.py
β”‚   β”œβ”€β”€ configs.py                                # Pydantic configs validator
β”‚   β”œβ”€β”€ constants.py                              # Constants for the training
β”‚   └── trainer.py                                # Trainer class to handle training operations
β”œβ”€β”€ submodules
β”‚   └── ddgun/                                    # The DDGun package that we use for one of the evaluations
β”œβ”€β”€ requirements.txt                              # Necessary Python Packages
β”œβ”€β”€ README.md                                     # You are here

πŸ€– Training

Prepare the dataset

Once you obtained the dataset, you need to first run netMHCpan to all sequences to compute the immunogenicity scores. First, you need to generate 9-mer peptides of each sequence of each dataset split:

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/TRAIN_DATASET.txt

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/VALID_DATASET.txt

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/TEST_DATASET.txt

These runs will create .pep files containing 9-mer peptides for each sequences which then can be passed to the netMHCpan:

cd scripts/netmhcpan/
bash netmhcpan_allele_scores_one_file.sh PATH/TO/TRAIN_DATASET_PEPTIDES_FILE.pep
bash netmhcpan_allele_scores_one_file.sh PATH/TO/VALID_DATASET_PEPTIDES_FILE.pep
bash netmhcpan_allele_scores_one_file.sh PATH/TO/TEST_DATASET_PEPTIDES_FILE.pep

These runs will create .pep.out files which contains the immunogenicity score for each peptides. Finally, we need to reconcile the peptides into sequences and calculate the hits scores of each sequence:

python scripts/netmhcpan/compute_hits_from_peptides.py \
--sequences_filepath=PATH/TO/TRAIN_DATASET.txt \
--peptides_dir=PATH/TO/TRAIN_PEPTIDES_DIR/ \
--peptide_file_prefix=TRAIN_PEPTIDE_FILE_PREFIX

python scripts/netmhcpan/compute_hits_from_peptides.py \
--sequences_filepath=PATH/TO/VALID_DATASET.txt \
--peptides_dir=PATH/TO/VALID_PEPTIDES_DIR/ \
--peptide_file_prefix=VALID_PEPTIDE_FILE_PREFIX

python scripts/netmhcpan/compute_hits_from_peptides.py \
--sequences_filepath=PATH/TO/TEST_DATASET.txt \
--peptides_dir=PATH/TO/TEST_PEPTIDES_DIR/ \
--peptide_file_prefix=TEST_PEPTIDE_FILE_PREFIX

Each of the peptide files that are generated from the previous command would have a prefix. For instance, from the previous commands you obtained 12 files (vaxformer_large_A0101.pep.out, ... , vaxformer_large_B5801.pep.out), then the peptide file prefix is vaxformer_large.

Practically, we can obtain quantiles from the distribution of hits of the training data split. For the sake of simplicity and reproducibility (and sanity check), you can check the Q1 and Q3 that we used to calculate the immunogenicity scores in the src/constants.py which are denoted as IMMUNOGENICITY_Q1 and IMMUNOGENICITY_Q3 respectively.

These quantiles are used to compute the immunogenicity scores of each sequence. Sequences whose scores are lower than the Q1 are considered to have low immunogenicity scores (0), in between Q1 and Q3 are intermediate immunogenicity scores (1), and above Q3 are high immunogenicity scores (3).

Training the model

Once the sequences and immunogenicity scores datasets are obtained, we can run a training process.

python scripts/train.py \
--config_filepath=PATH/TO/TRAIN_CONFIG_FILE

Selections of train config files can be found in the configs/train/ folder.

Both LSTM and Vaxformer can be trained with -small, -base, or -large setting. They differ in terms of the number of hidden layers and their sizes.

πŸ“ Generating sequences

Before any evaluation steps, we need to first generate the sequences with a pretrained model of choice

python scripts/sample.py \
--config_filepath=PATH/TO/TEST_CONFIG_FILE
--num_sequences=2000

Examples of test config files can be found in the configs/test/ folder.

πŸ“ˆ Evaluation

DDGun

  1. Use create_data_for_ddgun.ipynb provided in scripts folder to generate input for ddgun
  2. Install DDGun (https://github.com/biofold/ddgun)
  3. Use ddGun.sh provided in scripts folder
  4. dd_eval.ipynb provides functions necessary to read the output of the DDGun

netMHCpan

To evaluate the immunogenicity level of the generated sequences, we need to first generate 9-mer peptides for each generated sequences

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/FULL_GENERATED_LOW_IMMUNO_SEQUENCES.txt

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/FULL_GENERATED_INTERMEDIATE_IMMUNO_SEQUENCES.txt

python scripts/netmhcpan/generate_peptides_from_sequences.py \
--sequences_filepath=PATH/TO/FULL_GENERATED_HIGH_IMMUNO_SEQUENCES.txt

This run will create 3 .pep files (low, intermediate, and high). Next, we need to run the netMHCpan scoring script:

cd scripts/netmhcpan/
bash netmhcpan_allele_scores.sh PATH/TO/PEP_FILES_DIRECTORY/

This run will create a combined.pep file (in the same directory) which contains the concatenation of low, intermediate, and high peptides (in that order). The combined.pep file will then be passed to the netMHCpan computation. This will results in 12 MHC files that are scored. Notice that this run is different from the train-valid-test data preparation as this is designed to concatenate generated sequences across different immunogenicity scores mainly to understand the distribution of each immunogenicity score.

Finally, we need to reconcile the peptides into sequences and calculate the score quantiles:

python scripts/netmhcpan/compute_hits_from_peptides.py \
--sequences_filepath=PATH/TO/GENERATED_SEQUENCES.fasta \
--peptides_dir=PATH/TO/PEPTIDES_DIR/ \
--peptide_file_prefix=combined

Similar to the train-valid-test data preparation process, each peptide file that is generated would have a prefix. In this case, it would be combined.

AlphaFold 2

  1. Fold the proteins using alphafold.sh provided in scripts folder
  2. Open PyMol, open the proteins and run:
align wuhan_folded, protein_to_evaluate

where wuhan_alphafold denotes the reference protein and protein_to_evaluate denotes the generated protein.

βš–οΈ Results

PCA

<img src="assets/pca_square.png" width=500>

Perplexity

ModelTrainVal.
LSTM-small1.0191.021
LSTM-base1.0161.018
LSTM-large1.0161.019
Vaxformer-small1.0341.128
Vaxformer-base1.0141.043
Vaxformer-large1.0131.014

DDGun & RMSD

<img src="assets/folded_proteins.png" width=900>
ModelRMSD (\AA)$\Delta \Delta G$ (KCAL/MOL)
Random Mutation0.32 $\pm$ 0.23-2.51 $\pm$ 0.31
Naive Bayes$0.59 \pm 0.23$-0.5 $\pm$ 0.30
VAE0.48 $\pm$ 0.28-5.17 $\pm$ 0.51
LSTM-base0.51 $\pm$ 0.17-4.95 $\pm$ 1.09
Vaxformer-large0.67 $\pm$ 0.31-5.45 $\pm$ 0.72

netMHCpan

<table> <tr> <th rowspan="2">Model</th> <th colspan="3">T-statistics</th> <th colspan="3">U-statistics</th> </tr> <tr> <th>Medium-Low</th> <th>High-Medium</th> <th>High-Low</th> <th>Medium-Low</th> <th>High-Medium</th> <th>High-Low</th> </tr> <tr> <td>VAE</td> <td>15.22</td> <td>27.68</td> <td>39.17</td> <td>3.99 Γ— 10<sup>5</sup></td> <td>4.73 Γ— 10<sup>5</sup></td> <td>5.23 Γ— 10<sup>5</sup></td> </tr> <tr> <td>LSTM-base</td> <td>14.23</td> <td>16.83</td> <td>29.75</td> <td>4.45 Γ— 10<sup>5</sup></td> <td>4.37 Γ— 10<sup>5</sup></td> <td>5.16 Γ— 10<sup>5</sup></td> </tr> <tr> <td><strong>Vaxformer-large</strong></td> <td><strong>26.30</strong></td> <td><strong>38.02</strong></td> <td><strong>56.62</strong></td> <td><strong>4.75 Γ— 10<sup>5</sup></strong></td> <td><strong>5.21 Γ— 10<sup>5</sup></strong></td> <td><strong>5.56 Γ— 10<sup>5</sup></strong></td> </tr> </table> <img src="assets/netMHCpan.png" width=900>