Boosting the Full Potential of PyMOL with Structural Biology Plugins

1The first version of the PyMOL graphical environment for modeling and visualization of molecules, which dates back to 2000, was developed by Warren Lyford DeLano and was at first distributed by DeLano Scientific LLC. After the unexpected death of Warren DeLano in 2009, PyMOL was acquired by Schrödinger Inc. (NY, USA), which continued to support the open-source vision of the original author, but also provided commercial support for maintenance and access to additional features. Currently, PyMOL is a cross-platform, open-source but proprietary software program, maintained, developed and supported by Schrödinger Inc., which also reserves licensing rights. The major changes since Schrödinger Inc. acquisition released with PyMOL version 2 (PyMOL 2), feature a new Graphical User Interface (GUI) that makes use of PyQt5, and replaced Tcl/Tk. PyMOL 2 packages and dependencies management are now on Anaconda, since version 2.3, PyMOL started to be developed and released in Python 3. Such novelties made it possible for PyMOL to be continuously up to date, arriving at the latest available version 2.5, which was released in October 2021. The use of efficient and high-performing graphical libraries, i.e., Open Graphics Library (OpenGL) and the GL Shading Language (GLSL), underlies the advanced images and light renderings of PyMOL, which in turn permits deep control over the setup of final image appearance. Indeed, PyMOL is widely exploited by the scientific community for creating high-quality images and videos that accurately depict the molecular structure being presented.

2One of the main features of PyMOL is the ease for external developers to boost its functionality via scripts and plugins. In this way, PyMOL is programmatically accessible through its editable command-line and Application Programming Interface (API), which make PyMOL features easily exploitable by external software. As a consequence, sequence analysis, molecular docking, molecular dynamics, structure-function relationships analysis, protein structure prediction and virtual screening are all well-known structural bioinformatics approaches for which PyMOL plugins have been developed over the years. In addition to plugins, some useful tools are inherently implemented in PyMOL. For example, PyMOL handles molecular building and editing by providing a builder menu to draw small chemical compounds, peptides and nucleic-acids sequences. The builder menu comes with complete hydrogen and charge fixing tools and with 3D-model cleaning strategies. In the same context, PyMOL supports the possibility to mutate a protein residue or a nucleotide base with the ‘Mutagenesis Wizard’. The ‘Sculpting Wizard’ gives access to a real-time energy minimization function implemented in PyMOL, which is especially useful when dealing with chemical editing.

3PyMOL is an ideal candidate for developing advanced Computational Protein Analysis (CPA) pipelines. Its two-pronged nature of molecular graphics viewer and programmatically accessible interface make it particularly suited to easily develop Python-based software tools for advanced investigations of biological macromolecules. A paper in 2017 outlined PyMOL plugins and features that are useful in computational drug design. Some of these plugins are still actively maintained and are worth citing: CAVER 3.0, a tool for the analysis of transport pathways in dynamic protein structures; PYROSETTA, an interactive interface to the powerful Rosetta molecular modeling suite; Dynamics, which adds Gromacs-based molecular dynamics simulation features to PyMOL, and the APBS Electrostatics plugin, which implements the Adaptive Poisson–Boltzmann Solver is suited for the analysis of solvation and electrostatics in protein structures. The latter, which was has been reviewed in as a standalone plugin, is now available along with PyMOL installation. In addition to those that were previously available, several other PyMOL plugins have been released in recent years. The main aim of this review is to provide an overview on the newly released PyMOL plugins, which are summarized in Table 1, and their scope in CPA. These are introduced by dividing the CPA realm into three categories: (i) Protein Sequences and Structures Analyses (PSSAs); (ii) Protein-Ligand Interactions (PLI); (iii) Protein Dynamics (PD).

4PSSAs are crucial in a variety of biological research areas, particularly in the prediction and modeling of protein structures, which include similarity searches, alignments of sequences and structures, evolutionary and structural comparison, and homology modeling. PyMOL is widely used for PSSAs due to the visual and editing support that is provided for the sequences (i.e., protein, RNA and DNA sequences) of the loaded associated 3D structures.

5PSSAs include a wide range of often-combined methods such as similarity searches, alignments of sequences and structures, evolutionary comparisons and homology modeling. For example, the template-based prediction of protein 3D quaternary structures in complex with ligands and other heteroatoms requires: (i) databases search for homologous protein (s) with known 3D structure (s); (ii) multiple sequence and structure alignments tools, and (iii) algorithms for protein structure prediction and quality assessment. The PyMod plugin provides a GUI-based and easy-to-use environment for PSSAs by implementing, in an integrated manner, a wide range of algorithms, e.g., (PSI-)Blast, MUSCLE, ClustalW, Clustal Omega, WebLogo 3, ESPript 3.0, CAMPO, SCR-Find, PsiPred and MODELLER modules for 3D model building, sequence-structure alignment and quality assessment.

6The latest version of PyMod, i.e., PyMod 3, runs on Windows, MacOS, and Linux on both the open source and incentive PyMOL builds provided by Schrödinger Inc. Alignments, sequences or structures can be imported in PyMod, e.g., protein structures that are loaded in PyMOL. Any imported object is processed by PyMod to retrieve the associated information, e.g., heteroatoms and water molecules, and to display its sequence in the dedicated window. The latter provides the interface for the many PyMod functionalities, and is designed to visually support manual editing of sequences and alignments. The plugin so constituted is widely used for template-based protein structure predictions (Figure 1a–d), alignments, evolutionary conservation analyses, building of phylogenetic trees and database searches. Protein structure analyses now often rely on AlphaFold (AF) accurate predictions; nonetheless, some underlying dynamics features of proteins may be unveiled by comparing experimental and AF-predicted models. Such a hybrid approach was applied in a recent study on the human trans-3-hydroxy-L-proline dehydratase, in which the PyMod support for structural analyses, e.g., structure divergence plots, facilitated the comparisons between experimental and predicted 3D protein structures in different conformational states.

7Protein residue network (PRN) methods model proteins according to the graph theory to obtain a reduction of the complexity of protein structures into more simple descriptors (residues as vertices/nodes and the interactions between them as edges). Analyses such as communication pathways, allostery effects and networks of interactions can all be addressed by making use of PRN-based analyses. The plugin Protein Graph Analyser (pyProGA) offers a GUI for several PRN methods, i.e., centrality calculation, graph partitions evaluation, shortest path identification, ego graph, network differential analysis, binding energy computing, singular value decomposition technique and pair interaction energies analysis. Graphs solely computed from structural data (PDB file format) are distance-based networks (D-PRN). Otherwise, if an output from either the fragment molecular orbital or an Amber calculation (“.prmtop” file) is provided, the computed graphs are based on residue-pair interaction energies (PIE-PRN). The results can be inspected directly in the pyProGA window and/or in PyMOL (Figure 1e). The pyProGA window accommodates plot visualization, while PyMOL is used for color mapping the results on the analyzed proteins (Figure 1f,g). A PRN-based metric, the residue folding degree, has been proposed as a means to measure the propensity of the protein backbone to form secondary structures in Molecular Dynamics (MD) simulations. The residue folding degree analysis relies on conformational state descriptors, namely subgraph centrality, which can be computed from the pyProGA plugin. Additional metrics computable by pyProGA are appropriate for the identification of residues mediating key interactions between proteins, as showed with the de novo synthesized TIM barrel protein (PDB-ID: 5BVL), for which singular value decomposition (SVD) analysis was applied to identify protein-protein interactions between monomers.

8Membrane protein 3D structures are hard to solve, due to their inherent association with lipidic macromolecules that makes the protein solubilization a crucial point. Despite the challenges, many protocols have been developed for solving accurate 3D structures of membrane proteins, but are still cumbersome, expensive and time-consuming. The Small-angle X-ray scattering (SAXS) technique is preferred because of less stringent sample preparation requirements. The output of a SAXS experiment needs to be interpreted according to computational modeling approaches. A PyMOL plugin, MPBuilder, has been developed for assisting the modeling of membrane proteins on SAXS-derived data. Two tabs, for building and refinement (Figure 1j–k), guide the user through the analysis. The required inputs are a file storing SAXS data and the choice of an appropriate protein-lipidic assembly to guide the model building. MPBuilder greatly exploits PyMOL embedded functions, such as translation, rotation, atomic distance computation and atoms selection, as well as protein and lipidic visualization (Figure 1l). The implementation in the PyMOL framework makes it accessible for regular use, even by non-experts. For instance, MPBuilder was used for the interpretation of SAXS data of the phosphoethanolamine enzyme that is active at the lipidic interface, namely lipid A of pathogenic gram-negative bacteria as an antibiotic resistance mechanism.

9Water solvation is well known for its influence in protein conformation, pharmacodynamics, rational drug design and protein structure prediction. Given the relevance of water-mediated interactions, a PyMOL plugin called ProBiS H2O has been developed to infer conserved water binding sites by exploiting information derived from experimentally solved water-containing structures. Compared to other computationally expensive methods, such as molecular dynamics-based methods, it provides an easy and fast way for conserved water site identification. The PDB-ID of the protein of interest is the unique mandatory input of ProBiS H2O, starting from which similar proteins are retrieved, while taking into consideration a user-defined percentage of sequence identity. Afterwards, local superimposition, sampling and clustering of water molecules are carried out and a conservation score is computed for each cluster. In PyMOL, water molecules are visualized and colored according to the computed conservation score. Recently, a new version of ProBiS H2O was released, called ProBiS H2O MD. To account for lesser known proteins, hence proteins lacking enough experimental data for using the standard ProBiS H2O pipeline, this new version of the plugin optimizes the algorithm extending experimental data with MD trajectories. The latter must be a file (i.e., a “.dcd” file) derived from MD simulations of the protein of interest in a water box. At the end of the analysis, from within the plugin window, alternative visualizations of the clusters can be chosen and the amino acids encompassing the chosen binding site can be visualized (Figure 1h,i). An interesting application of ProBiS H2O MD elucidated the influence of local hydration on the PP1-Src kinase ligand binding. Alternatively, interacting water molecules can be inferred by making use of molecular docking engines, as it is done by Waterdock 2.0. This tool implements Autodock Vina and uses optimized parameters to dock water molecules in protein structures (Figure 2a–c). The influence of water-mediated interactions on protein-ligand complexes was recently investigated by making use of Waterdock 2.0 (in comparison with other software intended for the same aim), for estimating the binding free-energy of human acetylcholinesterase in a complex with tacrine at different solvation states.

10The functional and structural characterization of protein folding and conformational variability often employs structure superposition and comparison. PyMOL’s inherent alignment method is sequence-based, i.e., it makes use of a sequence alignment followed by a structural superimposition. However, for an efficient structure superimposition, related proteins with very low sequence identity must rely on structure-based alignment algorithms such as DALI, CE, SSAP and iPBA. The iPBAVizu plugin implements iPBA and simply requires a minimum of two 3D protein structures loaded in the PyMOL workspace to compute the structural alignment. The inferred sequence alignment is then visualized in the plugin window and the 3D structures are superimposed in PyMOL (Figure 2d). Many analyses are now relying on iPBAVizu for efficient protein structure-based superimposition. For instance, iPBAVizu guided the structural comparison of the primary contributors to SARS-CoV-2 antigenicity for the in-silico design of a multi-epitope vaccine candidate.

11Multiple Sequence Alignments (MSAs) have been used for decades to obtain structural, functional and evolutionary information on proteins. MSA-based analyses have applications in structure prediction, protein conformation dynamics, analysis of folding pathways, identification of binding sites, inference of interaction interfaces and prediction of interacting partners. Some methods that fall into the group of Direct Coupling Analyses (DCA) are used to depict the inherent evolutionary information of MSAs of homologous proteins. Co-evolutionary couplings among residue pairs are quantified in DCA with a score, called ‘Direct Information’ (DI), which is computed for each pair of residues. However, to obtain a comprehensive understanding of such analyses, the results must be processed and mapped on 3D structures. The PyMOL plugin DCA-MOL interactively visualizes coevolutionary residue-residue interactions in contact maps and 3D structures. By providing DI-containing files (Figure 2e), an MSA “.fasta” file and the 3D structure of the protein of interest, DCA-MOL maps the contacts onto the 3D structures loaded in PyMOL workspace and shows contacts and distances maps in its dedicated window (Figure 2f). The interactions are captured either in short (<8.0 Å) or long distances (>8.0 Å) and, when mapped in the same chain (intra-chain interactions), can support protein structure prediction analyses. However, since DCA data of concatenated MSAs can be interpreted, DCA-MOL is not limited to intra-chain analyses, but also inter-chain interactions can be analyzed. On the other hand, by providing more 3D structures of the same protein, the mapping of the analyzed coupling residues can be easily switched by DCA-MOL from one structure to another. In the case of multiple structures, provided by either a MD-trajectory or an experimental source, such visualization can help in rationalizing conformational variabilities. These case studies are deeply investigated in the tutorials provided by the developers (https://dca-mol.cent.uw.edu.pl; accessed on 14 October 2022).

12Protein-ligand interactions, prediction and analyses, comprise a set of widely used methods in Structure-Based Drug Discovery (SBDD), i.e., Molecular Docking (MDo) algorithms, re-scoring metrics, post-processing analyses and Virtual Screening (VS). Driven by the need to have in-hand and easy-to-access environments for SBDD, PyMOL has become a platform for assisting PLI analyses.

13The widespread use of protein-ligand docking has led to the development of a wide range of search algorithms, scoring functions, procedures and post-processing analysis tools that create diverse strategies and parameters to be set up by researchers. Although many tools are available for accomplishing each individual task, a lack of integration slows down and hampers the process. A newly developed PyMOL plugin, DockingPie, addresses such limitations by providing an interoperable implementation of the many tools that are needed to carry out each step of a MDo process, from the preparation of input files to the analyses of the results. Currently, the docking engines supported by DockingPie are Smina, Vina, RxDock and ADFR (Figure 3a), as well as their peculiar protocols such as the assignment of side chain flexibility or the setup of pharmacophoric restraints. From within DockingPie, a user can: (i) compute RMSD; (ii) visualize the results by means of predicted affinity-based plots, and (iii) carry out consensus scoring analyses. The input required is any ligand or protein currently loaded in the PyMOL workspace and, for ease of use, DockingPie handles the installation of all external dependencies. How to use the many functionalities of DockingPie is described in different tutorials (available at https://github.com/paiardin/DockingPie/wiki/Tutorials; accessed on 14 October 2022). It has been demonstrated that combining the results of different docking programs with the application of a consensus scoring algorithm can increase the success rates in VS processes. In this context, due to the implementation of different docking engines and consensus scoring metrics, DockingPie is a tool that lets the user collect and re-score the results in a single integrated environment (Figure 3b–d), which is particularly useful in VS campaigns.

14Hot-spots of interactions between proteins and ligands, namely protein regions with potential positive contribution to the ligand binding free energy, are usually inspected at the early stages of SBDD. The FTMap algorithm, which is available as a web server, probes the entire protein surface for detecting hot-spots of interactions, and returns the associated physico-chemical properties. However, the raw nature of FTMap output files prevents the addressing of biologically relevant conclusions. For these reasons, the DRUGpy plugin has been developed to analyze FTMap output data (Figure 3g), and evaluate the level of druggability of the corresponding hot-spots of interaction. DRUGpy quickly locates and defines pockets that FTMap data predict would bind drug-like compounds with high or low affinity (druggable sites), and also displays how protein conformational flexibility affects the target’s druggability. If a protein-ligand complex of an identified hot-spot is available, and is superimposed on the hot-spot, DRUGpy can compute fractional overlap to guide rational ligand design (Figure 3h). DRUGpy was used for the comparative analyses of two homologous enzymes, widely validated as drug targets, trypanothione reductase and glutathione reductase. DRUGpy fractional overlap analysis of the druggable hot spots of both enzymes added significant physico-chemical considerations towards the design of species-specific inhibitors.

15MDo analyses, especially when carried out on a large-scale, return a huge number of binding conformations to be analyzed. This is very common for flexible, ensemble and consensus docking protocols, making necessary the use of a tool for post-processing analyses of docked conformations. Indeed, some MDo engines directly implement a filtering tool for discarding redundant poses. However, the fact that a lot of tools do not consider symmetric poses leads to a misinterpretation of molecule orientations and similarities. To overcome this issue, PoseFilter implements two methods for post-processing analyses of docked complexes, i.e., root mean square deviations (RMSD) and interaction fingerprints, in a way that takes into account symmetric molecules. The correct interpretation of symmetric poses is a hurdle when ligands are docked in symmetrical binding sites, as usually occurs when a ligand binding pocket lies between adjacent monomers of oligomeric molecules. In analyzing these situations, PoseFilter gradually rotates the ligand and carries out an ‘all to all’ comparison, while discarding redundant poses. Either protein-ligand complexes or separated molecules can be provided (Figure 3e), and its use has been shown for (homo-) trimeric, dimeric and tetrameric proteins. Since PyMOL is widely used for inspecting protein-ligand complexes, PoseFilter provides a quick and functional tool for analyzing MDo results (Figure 3f).

16The dynamic behavior of biomolecules is preferably studied with visual support. Indeed, many algorithms for PD have been implemented in different molecular graphics environments. Whereas the length and size of MD trajectories can hamper storing and visualization, PyMOL was adapted to handle the parsing and rendering of big MD-data files. Such novelty has become a step towards the use of PyMOL either for carrying out simulations or for the post-processing of PD-derived data.

17The output of MD simulations is usually in the form of trajectory files that can be further interpreted, visualized and analyzed. Many python modules are available for molecular trajectory file interpretation, yet a free GUI for accessing such features is still missing. Geo-Measures (Figure 4a) was born as a plugin for providing easy access to MDTraj, a python library of trajectory analyses modules. Given a topology and trajectory file (file formats supported by PyMOL are feasible), the following analyses can be carried out: (i) Cα and dihedral angle calculations; (ii) Cα triangle area; (iii) probability density function (PDF) of trajectory frequencies; (iv) root mean square deviation (RMSD) along trajectories; (v) radius of gyration (Rg); (vi) free energy landscape (FEL); (vii) principal component analysis (PCA); (viii) Ramachandran map; (ix) root mean square fluctuation (RMSF); (x) secondary structures definition; (xi) distance calculations, and (xii) trajectory paths visualization. Plots of the results can be visualized and exported as images from a dedicated window of Geo-Measures (Figure 4b–f) The authors provided practical examples of Geo-Measures use on hemoglobin and Ecto-5′-nucleotidase. In a recent work on the mechanistic role of receptors-bound quercetin in hepatoprotection, the bound state stability, the compactness and conformational fluctuations were investigated analyzing the MD trajectories of the complex with GeoMeasures. In this application, Geo-Measures was pivotal to easily compute and plot the RMSD, Rg, RMSF and FEL. Moreover, Geo-Measures was the method of choice for analyzing more flexible and scattered trajectories, due to the support for PCA, which is usually applied in these cases.

18Many environments for MD are available; however, the setup of programs for MD still requires detailed knowledge in the field. Enlighten2, an easy-to-install interface to the Amber force-field, provides direct access to MD simulations. Protocols to set up MD simulations, which have already been tested and used in the previous version of Enlighten, together with the GUI, make MD accessible to all. Enlighten2 offers a GUI to several functionalities for molecule preparation (Figure 4g) and MD simulation set-up (Figure 4h). For instance, Enlighten2 implements Antechamber and propka to, respectively, compute ligands parametrization and to modify ligands and protein residues according to user-defined pH values. Furthermore, Enlighten2 supports the possibility to set up how to handle solvation and co-factors. When using Enlighten2, each step that modifies the protein’s properties is promptly visualized in PyMOL. Following molecule preparation, a dedicated object is created in PyMOL to store and visualize the solvent cap (Figure 4i). Moreover, the results are loaded in PyMOL as multiple separated objects. The interoperability of PyMOL and Enlighten2 is also emphasized in the tutorials provided on the Enlighten2 websites (https://enlighten2.github.io/tutorial1; accessed on 14 October 2022), where it is explained how to analyze mutant models of two enzyme-ligand complexes created with the PyMOL ‘Mutagenesis Wizard’. Providing easy access to MD simulations, also for non-expert users, it is now often a method of choice to support experimental data, as has been proposed for the engineering of glycan-binding proteins (GBP).

19The analysis of PD data, either derived from MD simulations or coarse-grained Elastic Network Models (ENMs), is often a complex task; hence dimensionality reduction algorithms have taken hold. Even though many algorithms for deciphering biological relevant conclusions from PD have been developed, their availability as standalone command-line tools, or integration in commercial software, limits a free and easy application. In this context, a recently released tool, MODE-TASK offers an array of implements for analyzing PD-derived data. MODE-TASK can be exploited to carry out the following analyses: (i) PCA (i.e., standard, kernel or incremental PCA); (ii) Normal Mode Analysis (NMA); (iii) Multidimensional Scaling (MDS); (iv) t-Distributed Stochastics, and (v) Neighbor Embedding (t-SNE). Furthermore, MODE-TASK incorporates a new algorithm, namely a coarse graining technique, which was validated for being less computationally expensive than those previously available. PyMOL has been chosen for hosting the GUI of MODE-TASK as a plugin, known as pyMODE-TASK (Figure 4j,k). Such integration makes the analysis of proteic large-scale motions easily accessible to PyMOL users. The pyMODE-TASK site provides introductory tutorials (https://pymode-task.readthedocs.io/en/latest/pyMODE-TASK_usage.html; accessed on 14 October 2022) on the use of this plugin, and its application in NMA and PCA of proteins (Figure 4l–n). PCA when applied to MD is pivotal to spot differences in the analysis of the same biological system in different states, as described for the destabilizing mutation of renin-angiotensinogen system, and the effect of phosphorylation on toll/interleukin-1 receptor domain-containing adapter protein (TIRAP)-mediated signaling. In both works, pyMODE-TASK was exploited to compute and plot the principal components of the system.

20It is clear that an interoperable use of software and tools can greatly facilitate the development of pipelines to approach CPA. Several algorithms with applications in this field have been developed accompanied by a command-line interface only, limiting their use to experts in computational biology. Moreover, given the computational capacity and precision, alongside existing ones, an increasing number of in-silico approaches will be available, placing the scientific community in front of a vast choice of alternatives. In this context, PyMOL is positioned as an integrative environment for bridging molecular visualization and CPA. Especially, PyMOL eligibility as a versatile platform, must be acknowledged in relation to accessing its framework through plugins. In addition, its most outstanding features, namely excellent image rendering, programmatic accessibility and inherent chemical editing tools, have contributed to making it an integral part of structural bioinformatics. This review highlights how the use of plugins that exploit and, at the same time, boost the use of PyMOL, can help create an environment to seamlessly carry out CPA.