Torsion Sampling

[jameseastwood]

Type
System sampling check
Question
Is phase space well sampled?
How
Minima sampled ~60%
Tool
MDA + rdkit
Note
Ideally what you should be looking at here is a combination of the solvent, vacuum, and complex differences and the variation in dihedral sampling across lambda states. I.e. is there a major change happening in my transformation that isn't being sampled across all states -- @IAlibay

[RiesBen]

@jameseastwood I somehow di dnot manage to link this to the PR, that closes this task. ... But there is now a new "task" in the project linked to the pr. shall we close this one then? :)

New Task: #20

[jameseastwood]

No problem, I linked them. Carry on with the PR, and it will close this issue when you merge the PR.

[jthorton]

I recently needed something similar when inspecting the results of an RBFE calculation, see the script below which borrows some of the work from Ben in PR #20. The plots are not as nice as Bens work and there is an issue with building the molecule for state B (see the carbon atom with a negative charge in the stateB example) but this should be a good start for a general workflow.

Examples

StateA

State B

"""Script to analyze ligand torsions in a trajectory for the end states."""

import base64
import io
import MDAnalysis as mda
from MDAnalysis.analysis.dihedrals import Dihedral
from openfe_analysis import FEReader
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
import pandas as pd
import matplotlib.pyplot as plt


# functions from bens PR to openfe-analysis
def get_rotatable_bonds(mol: Chem.Mol) -> list[Chem.Bond]:
    """Function to find all rotatable bonds in a molecule taking symmetry into account.
    Parameters
    ----------
    mol : Chem.Mol
        The rdkit.Chem.Mol object
    Returns
    -------
    List[Chem.Bond]
         List of rdkit.Chem.Bond that were found in a molecule taking symmetry into account.
    """
    RotatableBondSmarts = Chem.MolFromSmarts("[!$(*#*)&!D1]-&!@[!$(*#*)&!D1]")
    find_rbs = lambda x, y=RotatableBondSmarts: x.GetSubstructMatches(y, uniquify=1)
    rbs = find_rbs(mol)
    bonds = [mol.GetBondBetweenAtoms(*inds) for inds in rbs]

    return bonds

def get_torsion_atoms_idx(bond: Chem.Bond) -> list[int]:
    """Function that finds the atomic ids that specify a torsion around the bond of interest.
    Parameters
    ----------
    bond : Chem.Bond
        The bond of interest around which the torsion should be specified.
    Returns
    -------
    List[int]
        List of atomic ids that specify a torsion around the bond of interest.
    """
    bond_atoms = [bond.GetBeginAtom(), bond.GetEndAtom()]
    additional_atom1 = list(filter(lambda x: x.GetIdx() != bond_atoms[1].GetIdx(), bond_atoms[0].GetNeighbors()))[0]
    additional_atom2 = list(filter(lambda x: x.GetIdx() != bond_atoms[0].GetIdx(), bond_atoms[1].GetNeighbors()))[0]
    torsion_atoms = [additional_atom1] + list(bond_atoms) + [additional_atom2]
    torsion_atom_ids = [a.GetIdx() for a in torsion_atoms]

    return torsion_atom_ids

def figure_to_img(figure) -> str:
    """Convert a matplotlib figure to an embeddable HTML image tag.

    Args:
        figure: The figure to convert.

    Returns:
        The HTML image tag.
    """

    with io.BytesIO() as stream:
        figure.savefig(stream, format="svg")
        data = base64.b64encode(stream.getvalue()).decode()

    return f'<img src="data:image/svg+xml;base64,{data}"></img>'

def mol_to_html(mol: str) -> str:
    raw_data = base64.b64encode(mol.encode()).decode()
    return f"<img src='data:image/svg+xml;base64,{raw_data}'/>"

def draw_molecule(mol: Chem.Mol, highlight_atoms: list[int]) -> str:
    """Draw a molecule with highlighted atoms and return as HTML image tag.

    Args:
        mol: The molecule to draw.
        highlight_atoms: The atoms to highlight.

    Returns:
        The HTML image tag.
    """
    AllChem.Compute2DCoords(mol)
    drawer = Draw.rdMolDraw2D.MolDraw2DSVG(300, 300)
    Draw.rdMolDraw2D.PrepareAndDrawMolecule(drawer, mol, highlightAtoms=highlight_atoms)
    drawer.FinishDrawing()
    svg = drawer.GetDrawingText()

    return mol_to_html(svg)


def main():

    for state_id, label, bfactor in zip([0, -1], ["stateA", "stateB"], [0.25, 0.75]):
        print(f"Processing state {state_id} ({label})")
        u = mda.Universe("hybrid_system.pdb", "simulation.nc", format=FEReader, state_id=state_id)

        # get only the atoms in using the bfactor hack
        state_atoms = np.array([atom.ix for atom in u.atoms if atom.bfactor in (bfactor, 0.5)])
        state = u.atoms[state_atoms]
        state_lig = state.select_atoms("resname UNK")
        # guess the bonds for rdkit
        state_lig.guess_bonds()
        # convert to rdkit
        state_rd = state_lig.atoms.convert_to("RDKIT")
        # get the rotatable bonds
        rotatable_bonds = get_rotatable_bonds(state_rd)
        # now we have the torsion index create an atom group for each index and use MDAnalysis to get the diehdral values 
        atom_groups = [u.atoms[get_torsion_atoms_idx(bond)] for bond in rotatable_bonds]
        dihedrals = Dihedral(atom_groups).run()
        # for each dihedral angle we want to make a matplotlib histogram and rdkit image with the dihedral highlighted on the molecule
        # and save the data as a html file using pandas
        all_data = []
        for i, bond in enumerate(rotatable_bonds):
            angle_data = dihedrals.angles[:, i]
            torsion_atoms = get_torsion_atoms_idx(bond)
            fig, axis = plt.subplots(1, 1, figsize=(4.0, 4.0))
            axis.hist(angle_data, bins=36, range=(-180, 180), density=True)
            axis.set_xlabel("Dihedral angle (degrees)")
            axis.set_ylabel("Density")
            axis.set_title(f"Torsion around bond {torsion_atoms[1]}-{torsion_atoms[2]}")
            fig.tight_layout()
            fig_html = figure_to_img(fig)
            plt.close(fig)
            mol_html = draw_molecule(state_rd, torsion_atoms)

            all_data.append({
                "fig": fig_html,
                "mol": mol_html,
            })

        # Save all data to a HTML file
        df = pd.DataFrame(all_data)
        df.to_html(f"torsion_analysis_{label}.html", escape=False, )

if __name__ == "__main__":
    main()