Molecule substructural alignment

There are cases when we want to align molecules which only share a fraction of their structure, like ligands in a congeneric series, where a common scaffold or core is common to all ligands in the series. We can perform this alignment using HTMD function: maximalSubstructureAlignment

In this example, we will use two ligands which bind to the Beta Adrenergic receptor. They share a fraction of their scaffold and bind in the same pocket.

Quick Example

from htmd.ui import *
from htmd.molecule.graphalignment import maximalSubstructureAlignment
from htmd.home import home

path = home(dataDir='test-molecule-graphalignment')

# Load the molecule which will be used as reference
reference_ligand = Molecule(os.path.join(path, 'ref_lig.pdb'))
reference_ligand.view(style='Lines', name='Reference')

# Load other ligand from this congeneric series
ligand_to_align = Molecule(os.path.join(path, 'lig2align.pdb'))
ligand_to_align.view(style='Licorice', name='ToAlign')

# Align
lig_aligned = maximalSubstructureAlignment(reference_ligand,ligand_to_align)

# View the aligned molecule
lig_aligned.view(style='Licorice', name='Aligned')

Detailed explanation

First, we load the molecule which will be used as reference, cocrystallized in PDB entry 2Y02 with resname WHJ

from htmd.ui import *
from htmd.molecule.graphalignment import maximalSubstructureAlignment

reference_crystal = Molecule('2Y02')
reference_crystal.filter('protein or resname WHJ')
reference_ligand = reference_crystal.copy()
# Extracts the ligand from the protein
reference_ligand.filter('resname WHJ')
# There are two ligands with the same resname, so we select one of them
reference_ligand.filter('residue 1')

Now, let’s extract the other ligand from this congeneric series, cocrystallized in PDB entry 2Y03 with resname 5FW

ligand_to_align = Molecule('2Y03')
ligand_to_align.filter('resname 5FW')
# Again, there are two residues with the same name, and
# residue 1 happens to be already aligned with the reference ligand, so we use residue 0
ligand_to_align.filter('residue 0')

Let’s take a look at both molecules as they are now. The reference ligand is displayed with narrower bonds.

reference_ligand.view(style='Lines',name='Reference')
ligand_to_align.view(style='Licorice',name='ToAlign')

Now we align the extracted ligand to the crystal of the molecule we are using as reference. Then, we can see the result. Again, the reference ligand is displayed with narrower bonds.

lig_aligned = maximalSubstructureAlignment(reference_ligand,ligand_to_align)
lig_aligned.view(style='Licorice',name='Aligned')

Finally, we can see how the aligned ligand looks together with the protein.

reference_crystal.view(sel='protein',style='NewCartoon',name='Protein')