<jupyter_start><jupyter_text>QuickstartThis is a very short guide on how to use ProLIF to generate an interaction fingerprint for a ligand-protein complex from an MD simulation.Let's start by importing MDAnalysis and ProLIF to read our input files:<jupyter_code>import MDAnalysis as mda
import prolif as plf
# load trajectory
u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)
# create selections for the ligand and protein
lig = u.atoms.select_atoms("resname LIG")
prot = u.atoms.select_atoms("protein")
lig, prot<jupyter_output><empty_output><jupyter_text>MDAnalysis should automatically recognize the file type that you're using from its extension. Click [here](https://userguide.mdanalysis.org/stable/quickstart.html) to learn more about loading files with MDAnalysis, and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.Next, lets make sure that our ligand was correctly read by MDAnalysis.⚠ This next step is crucial if you're loading a structure from a file that doesn't explicitely contain bond orders and formal charges, such as a PDB file or most MD trajectory files. MDAnalysis will infer those from the atoms connectivity, which **requires all atoms including hydrogens to be present in the input file**.ProLIF molecules are built on top of RDKit and are compatible with its drawing code. Let's have a quick look at our ligand:<jupyter_code>from rdkit import Chem
from rdkit.Chem import Draw
# create a molecule from the MDAnalysis selection
lmol = plf.Molecule.from_mda(lig)
# cleanup before drawing
mol = Chem.RemoveHs(lmol)
mol.RemoveAllConformers()
Draw.MolToImage(mol, size=(400,200))<jupyter_output><empty_output><jupyter_text>We can do the same for the residues in the protein (I'll only show the first 20 to keep the notebook short):<jupyter_code>pmol = plf.Molecule.from_mda(prot)
frags = []
# to show all residues, simply use `for res in pmol:`
for i in range(20):
    res = pmol[i]
    mol = Chem.RemoveHs(res)
    mol.RemoveAllConformers()
    frags.append(mol)
Draw.MolsToGridImage(frags,
                     legends=[str(res.resid) for res in pmol], 
                     subImgSize=(200, 140),
                     molsPerRow=4,
                     maxMols=prot.n_residues)<jupyter_output><empty_output><jupyter_text>If only one of the two molecules was not read correctly, it might be because your input file doesn't contain bond information (typically the case for PDB files) so MDAnalysis has to guess it. This is usually done automatically, but when it's done on two different atom selections it will fail, so one has to explicitely tell MDAnalysis to do bond detection for the whole system.  To fix this, simply add `guess_bonds=True` when creating the Universe: `u = mda.Universe("file.pdb", guess_bonds=True)`.In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of MDAnalysis molecules to RDKit, but workarounds exist:* Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:  `u = mda.Universe("file.pdb", guess_bonds=True, vdwradii={"H": 1.05, "O": 1.48})`* If you have separate PDB files for your ligand and protein, you can also directly use RDKit:```pythonmol = Chem.MolFromPDBFile("prot.pdb", removeHs=False)prot = plf.Molecule(mol)mol = Chem.MolFromPDBFile("lig.pdb", removeHs=False)lig = plf.Molecule(mol)fp = plf.Fingerprint()fp.run_from_iterable([lig], prot)``` Everything looks good, we can now compute a fingerprint:<jupyter_code># use default interactions
fp = plf.Fingerprint()
# run on a slice of frames from begining to end with a step of 10
fp.run(u.trajectory[::10], lig, prot)<jupyter_output><empty_output><jupyter_text>The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. Alternatively, you can pass a list of residues like so:```pythonfp.run(..., residues=["TYR38.A", "ASP129.A"])```Or simply use `fp.run(..., residues="all")` to use all residues in the `prot` selection.The best way to access our results is to export our interaction fingerprints to a Pandas DataFrame. By default, the resulting DataFrame only keeps track of residues and interaction types that were seen in at least one of the frames in your trajectory. If needed, you can access the full results with `fp.to_dataframe(drop_empty=False)`.<jupyter_code>df = fp.to_dataframe()
# show only the 10 first frames
df.head(10)<jupyter_output><empty_output><jupyter_text>If the generated DataFrame is empty and you are using a PDB file, please manually call `u.atoms.guess_bonds()` right after creating the Universe with MDAnalysis, and rerun the code up to this point. This problem typically happens when the PDB file is missing all or some of the `CONECT` records to indicate bonds, which leads to structures with atoms but no bonds. By design, those structures cannot match any of the patterns that are used to define interactions, hence no interaction is detected, resulting in an empty dataframe.<jupyter_code># drop the ligand residue column since there's only a single ligand residue
df = df.droplevel("ligand", axis=1)
df.head(5)
# show all pi-stacking interactions
df.xs("PiStacking", level="interaction", axis=1).head(5)
# show all interactions with a specific protein residue
df.xs("ASP129.A", level="protein", axis=1).head(5)
# or more simply
df["ASP129.A"].head(5)<jupyter_output><empty_output><jupyter_text>Here's a simple example to plot the interactions over time<jupyter_code>import seaborn as sns
import pandas as pd

# reorganize data
data = df.reset_index()
data = pd.melt(data, id_vars=["Frame"], var_name=["residue","interaction"])
data = data[data["value"] != False]
data.reset_index(inplace=True, drop=True)

# plot
sns.set_theme(font_scale=.8, style="white", context="talk")
g = sns.catplot(
    data=data, x="interaction", y="Frame", hue="interaction", col="residue",
    hue_order=["Hydrophobic", "HBDonor", "HBAcceptor", "PiStacking", "CationPi", "Cationic"],
    height=3, aspect=0.2, jitter=0, sharex=False, marker="_", s=8, linewidth=3.5,
)
g.set_titles("{col_name}")
g.set(xticks=[], ylim=(-.5, data.Frame.max()+1))
g.set_xticklabels([])
g.set_xlabels("")
g.fig.subplots_adjust(wspace=0)
g.add_legend()
g.despine(bottom=True)
for ax in g.axes.flat:
    ax.invert_yaxis()
    ax.set_title(ax.get_title(), pad=15, rotation=60, ha="center", va="baseline")
# calculate the occurence of each interaction on the trajectory
occ = df.mean()
# restrict to the frequent ones
occ.loc[occ > 0.3]
# regroup all interactions together and do the same
g = (df.groupby(level=["protein"], axis=1)
       .sum()
       .astype(bool)
       .mean())
g.loc[g > 0.3]<jupyter_output><empty_output><jupyter_text>You can also compute a Tanimoto similarity between each frame:<jupyter_code>from rdkit import DataStructs
bvs = fp.to_bitvectors()
tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bvs[0], bvs)
tanimoto_sims<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10082597550555397, 'fraction_numerical': 0.010253489034463116, 'mean_word_length': 5.2761394101876675, 'pattern_counts': {'":': 2, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '566012', 'license': 'Apache-2.0', 'n_tokens_mistral': 2304, 'n_tokens_neox': 2136, 'n_words': 792, 'path': 'docs/notebooks/quickstart.ipynb', 'repo_name': 'ale94mleon/ProLIF', 'size': '7022'}