Software accompanying Feinauer, Christoph, Barthelemy Meynard-Piganeau, and Carlo Lucibello. "Interpretable pairwise distillations for generative protein sequence models." PLOS Computational Biology 18.6 (2022): e1010219
Please cite the paper if you use this code (see below).
- Linux (code has been tested on Ubuntu and Manjaro, not on anything else)
- Julia version 1.6 or higher
- Conda 4.11 or higher
- Git LFS
- Code assumes at various places that a GPU is present
- If you want to reproduce the contact prediction results, hmmer needs to be installed.
Make sure git-lfs is installed by running
$ git lfs install --skip-repoon the shell.
Then, execute
$ git clone git@github.com:christophfeinauer/PairwiseDistillations.git && cd PairwiseDistillationsThis downloads the necessary code files and the data.
Please make sure that the git-lfs objects were pulled correctly by executing in the repository directory
head ./data/YAP1_HUMAN_1_b0.5.a2mThis should show you protein sequences. If not, something has gone wrong.
In the repository directory, run
$ julia --project=. -e 'using Pkg; Pkg.instantiate(;verbose=true)'In the repository directory, run
$ conda env create -f environment.ymlThen, activate it using
$ conda activate PairwiseDistillationsAll subsequent instructions assume that you are in the correct environment.
The EVMutation code comes from the plmc repostory, which is used as a submodule here.
In the repository directory, run
$ git submodule updateto pull the plmc submodule. Then, enter the subdirectory of plmc and build it with
$ pushd evmutation/plmc/ && make all && popdNOTE: The plmc code also allows to compile with parallelization enabled (see readme in the subdirectory). This can be achieved by substituting make all in the command above with make all-openmp. However, the parallel implementation crashed our machine (not just the program, but the complete server) repeatedly. We did not invest any time in debugging it and just used the single-core version.
The complete pipeline is
1. Preprocess Data
2. Train Original Model and Create Samples
3. Extract Independent and Paiwise Models using Samples
4. Evaluate (calculate energies)
5. Plot
There is a bash script called run_all.sh in the repository directory which contains the complete pipeline. However, due to the complex nature of the pipeline it is recommended to run the following commands one by one.
In the repository directoy, run
./data/preprocess_data.shThis creates several new files in the data directory:
- For every
.a2mfile, it creates a set of training and testing sequences (suffix.trainand.test) - For every test sequence, it calculates the mean and minimum Hamming distance (suffix
.train_test_hamming) - For every mutational dataset, it creates an MSA that contains the mutated sequence and the corresponding experimental value in the annotation (suffix
.exp)
The following steps train the original models and create samples from them (except for EVMutation). The later scripts ignore as far as possible cases where original models are missing, so you can also train for example only ArDCA models or just a subset of the VAE models and go on to the later stages.
The ArDCA code comes from the ArDCA repository. It is installed as a registered Julia package (see above).
In the repostitory directory, run
ardca/train_ardca_models.shThis calls the training code found in ardca/ardca.jl on all training datasets and saves the model in models/. The naming of the models follows the scheme {DATASET_FILENAME}.train.ardca.jld. Training is done using threads. The default is to set the number of threads equal to the number of cores in the machine - if that is too much, it can be changed by modifying the option -t $(nproc -all) in the julia invocation in ardca/train_ardca_models.sh.
The code also creates samples from the T, U and M distribution. They are saved in samples/ and follow the naming scheme {MODEL_FILENAME}.samples_{M,T,U}.h5.
Most of the VAE code comes from the repository accompanying the paper "Deciphering protein evolution and fitness landscapes with latent space models". A reduced and modified version is included directly as repository files here.
In the repository directory, run
vae/train_vae_models.shThe code creates the same files as for ArDCA, with the same naming scheme except that the filenames for the models and samples also contain the settings for the latent dimension, the number of hidden units and the weight decays.
Note that creating the samples for the VAE can create a very long time (up to several hours for a single model), depending on the machine/GPU. Most of that time is not spent in creating the sequences, butevaluating their enegies/negative log probabilities due to the fact that we use a very large number of ELBO samples (5000). If you want to quickly reproduce the result you can try reducing this number by modfying the line elbo_samples=5000 in vae/train_vae_models.sh.
In the repository directory, run
evmutation/train_evmutation.shThis trains EVMutation models on all train datasets. The output of EVMutation is a binary file foraatm, which is then transformed into the same format as we use for the extracted pairwise models by the script evmutation/plmc_to_extracted.jl.
After having trained the original models, run in the repository directory
extracted/extract_couplings.shThis extracts independent and pairwise models from the samples and saves them to the models/ folder, adding a extracted_{D} tag, where D is the distribution used (M, T or U). Note that while we also have code that extracts models using samples from the T (Training) distribution, these typically do not perform well and are not reported in the paper.
The extraction code is in extracted/extract_couplings_from_samples.py.
After having trained the original models and extracted pairwise and independent models, run
evaluations/evaluate_all.shin the repository directory. This evaluates the the energies/log probabilities of all models found in the models/ folder on the datasets that end in .test or .exp. These are plain text files, adding a eval_TYPE tag to the filenames, where TYPE indicates whether the values correspond to test or experimental (mutational) sequences.
The naming scheme is cumulative, so for example the evluation file called YAP1_HUMAN_1_b0.5.a2m.train.vae.latentdim_5.weightdecay_0.1.numhiddenunits_40.pt.extracted_PW_M.h5.eval_test would contain the
- Energy values on the test set...
- ... of the pairwise model extracted using samples from the original model distribution
- ... of a VAE model trained with
40hidden units, weight decay set to0.1and a latent dimension of5 - ... trained on the file
YAP1_HUMAN_1_b0.5.a2m.train.
Independent models have IND instead of PW in the name. Note also that original ArDCA and VAE models report the log probability and not the eneriges. Extracted models and EVMutation models on the other hand report the energy. When comparing the evaluations of these model types the sign of one has to be switched.
The script for using extracted models for contact prediction can be found in extracted/predict_contacts.sh. Note that hmmer needs to be installed for this to work.
Plotting scripts can be found in plots/.
This is a list of repositories from which some of the code in this repository has been taken. Please cite appropriately if you use code from this repository.
- plmc (EVMutation): https://github.com/debbiemarkslab/plmc
- ArDCA: https://github.com/pagnani/ArDCA.jl
- VAE: https://github.com/xqding/PEVAE_Paper
Please use these citations if you use this code. They will be updated as soon as the paper is accepted.
Feinauer, Christoph, Barthelemy Meynard-Piganeau, and Carlo Lucibello. "Interpretable pairwise distillations for generative protein sequence models." PLOS Computational Biology 18.6 (2022): e1010219.
@article{feinauer2022interpretable, title={Interpretable pairwise distillations for generative protein sequence models}, author={Feinauer, Christoph and Meynard-Piganeau, Barthelemy and Lucibello, Carlo}, journal={PLOS Computational Biology}, volume={18}, number={6}, pages={e1010219}, year={2022}, publisher={Public Library of Science San Francisco, CA USA} }