PyAutoFEP

PyAutoFEP: an automated FEP workflow for GROMACS integrating enhanced sampling methods

PyAutoFEP is a tool to automate Free Energy Perturbations (FEP) calculations to estimate Relative Free Energies of Binding (RFEB) of small molecules to macromolecular targets. It automates the generation of perturbation maps, the building of dual-topologies for ligand pairs, the setup of MD systems, and the analysis. Distinctively, PyAutoFEP supports multiple force fields, integrates enhanced sampling methods, and allows flexible λ windows schemes. Furthermore, it aims to be as flexible as possible, giving the user great control over all steps. Nevertheless, reasonable defaults and automation are provided, so that PyAutoFEP can be used by non experts. PyAutoFEP is written in Python3 and uses GROMACS.

Announcements

Commit SHA1 TBD (20.06.2023)
This commit sets a new default value: output_collecttype=python. Previously,the default collect script was a binary created using PyInstaller. This was done because back in 2020 not every cluster had a modern enough Python3 version, so this could save users some time in setting up an env just for a simple script. Nowadays Python3 is everywhere and keeping the old binary as default makes no sense. Note that output_collecttype=bin is still allowed, but I expect to drop support in the upcoming updates.

Commit f20b3c0 (11.01.2023)
This commit changes the poses reading. Before it, prepare_dual_topology.py would read poses for every molecule in the input molecules, although only the poses for which the A states are needed would ever be used. Now, only the poses needed will be read. This speeds up execution, specially when using pose_loader=superimpose.

Commit c6e427c (06.11.2022)
This commit fixes the restraint code. Restraints can be selected by passing mdp_substitution to prepare_dual_topology.py, eg:

mdp_substitution = define: -DPOSRES_PROTEIN_CA

The following restraints are automatically generated and can bu used in this manner.

Molecule	`make_ndx` selection restraint	Affected atoms	Define
Ligand	LIG	All atoms	POSRES_LIG
Protein	Protein-H	Protein's heavy atoms	POSRES_PROTEIN
Protein	C-alpha	Protein's Cα	POSRES_PROTEIN_CA
Protein	Backbone	Protein's backbond	POSRES_PROTEIN_BB
Water	Water	Water molecules	POSRES_WATER

Commit f8324a3 (22.09.2022)
This commit contains a massive backport of upcoming charge-changing perturbations to the topology reading code. This new code broadens the topology support, allowing for the following.

Topologies containing A and B states. Note, however, that B states will not be used to prepare the dual topology! Dual topology is built from distinct molecules. In case the input topology has both A and B states, prepare_dual_topology.py will merely check that states A and B match (eg, as generated by SwissParam), and use data from state A.
New code for concomitant perturbing charge and VdW parameter of atoms. Not really used yet, but paves way for charge-changing perturbations and perturbing atoms in rings.
General topology coverage. Online dihedrals and more terms are supported.

Commit 190bef6 (06.06.2022)
This commit brings three new major features:

OpenBabel 3.X is now supported, along with backward compatibility with OpenBabel 2.4. Because OpenBabel 3 series is the one actively maintained upstream, all users are recommended to upgrade.
Reading of poses from Vina (and derivatives) .pdbqt files. prepare_dual_topology.py now accepts pose_loader=vina which reads .pdbqt files. Poses in a ligand file can be selected setting poses_advanced_options cluster_docking_data (see manual for more info).
Automated parameterization of input ligands. prepare_dual_topology.py can now run external small molecule parameterization tools automatically (using parameterize and related options, see manual for further info). Currently, the only parameterization software supported is AcPYPE, but I hope others can be added soon.

Furthermore, this update fixes some bugs and improves 3D MCS code, which is now more robust. Two defaults changed: in prepare_dual_topology.py, gmx_bin_run=gmx_mpi is now the default, as a MPI-enabled gmx executable is required and GROMACS compiles a gmx_mpi when using MPI. Also, generate_perturbation_map.py now will not save full_graph.svg and result_plot.svg unless the verbosity is at least INFO.

As always, please, fill issues should you experience any problems.

Commit fe41f7d (19.01.2022)
This commit introduces a bunch of new features and code changes. Even though I tested newly implemented and rewritten code, things may break. Please, fill issues should you experience any problem. Main changes:

Modifications to the atom matching functions in PyAutoFEP were done to implement support for user-supplied atom maps.
User-supplied atom maps are also available for the superimpose pose loader
Selecting 3D MCS for merge_topologies.merge_topologies is now supported (3D MCS in generate_perturbation_map.py and superimpose loader coming soon)
3D MCS code rewritten for better support for ligand pairs with multiple atom matches

Requirements

Common GNU programs: Bash, awk, tar
GROMACS 2016 or newer
GNU parallel on the run node, used during the rerun step (See manual for details).
Python 3.6+
rdkit 2019.03+
networkx 2.X (1.X versions are not supported)
alchemlyb 0.6.0 & pymbar 3.0.5 OR alchemlyb 0.3.0 & pymbar 3.0.3 (Because of choderalab/pymbar#419)
openbabel 2.4 or 3.X (sparsely used, mainly to load receptor files in prepare_perturbation_map.py.)
matplotlib (required only in analyze_results.py, optional in generate_perturbation_map.py)
numpy

Optional requirements. The following are not required to run basic calculations in PyAutoFEP, but are needed for specific functions.

biopython (allows sequence alignment when reading initial pose data)
mdanalysis (allows use of atom selection language in some contexts)
pytest (required to run Python tests)
packaging (used to compare package versions, falling back to distutils)

Install

To install PyAutoFEP, please, clone this repository using

git clone https://github.com/luancarvalhomartins/PyAutoFEP.git

Required dependencies can be installed using Anaconda, except for pymbar and alchemlyb, which must be installed using pip.

# Create a conda environment and activate it.
conda create -n PyAutoFEP
conda activate PyAutoFEP

# Install stuff
conda install -c conda-forge rdkit
conda install -c conda-forge openbabel
conda install matplotlib networkx pip
# Use pip to install pymbar and alchemlyb
pip install pymbar alchemlyb==0.6.0

Documentation

Manual

PyAutoFEP manual describes in detail its functions and options.

Tutorials

Farnesoid X receptor tutorial - A tutorial using the Farnesoid X receptor and a series of rigid binders is available.

(We plan to add more tutorials, covering specific aspects of PyAutoFEP in the near future.)

Issues & pull requests

Issues and pull requests are welcome. When filling a GitHub issue, please include as much details as possible. Inputs and verbose outputs are also useful (if available/relevant). And thanks for reporting bugs!

Roadmap

Aside from bug squashing, I am currently working on charged perturbations and in a GUI.

Charged pertubations are being implemented using the alchemical co-ion method, which seems to be both most general and simpler to code. Briefly, and random water molecule will be perturbed to an ion of the opposite charge as the ligand perturbation. At first, only charge differences of +1 and -1 will be supported (this should cover most of the use cases, anyway). Code for regular, non-charged perturbations will not be affected.
A PyQt5 GUI for PyAutoFEP is being written. So far, a perturbation map editor was mostly implemented and a ligand table is beign worked on. The GUI will be frontend to the scripts, so that no function will depend on using the first. The GUI development is low priority right now, so this is not making into the tree anytime soon.

Further goals

Covalent perturbations
Automated cycle-closure histeresis (and likely other analysis as well)
Support for peptides as ligands
More tutorials
A website

Citation

If PyAutoFEP is used in scientific publications, please cite:

LC Martins, EA Cino, RS Ferreira. PyAutoFEP: An Automated Free Energy Perturbation Workflow for GROMACS Integrating Enhanced Sampling Methods. Journal of Chemical Theory and Computation. 2021 17 (7), 4262-4273. LINK

Legal notice

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses

lmmpf/PyAutoFEP