Ctrl+K

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 moleculekit.tools.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

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')

2024-06-12 13:08:39,766 - moleculekit.molecule - INFO - Removed 370 atoms. 4694 atoms remaining in the molecule.
2024-06-12 13:08:39,797 - moleculekit.molecule - INFO - Removed 4640 atoms. 54 atoms remaining in the molecule.
2024-06-12 13:08:39,798 - moleculekit.molecule - INFO - Removed 27 atoms. 27 atoms remaining in the molecule.

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26], dtype=int32)

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')

2024-06-12 13:08:42,724 - moleculekit.molecule - INFO - Removed 4832 atoms. 30 atoms remaining in the molecule.
2024-06-12 13:08:42,726 - moleculekit.molecule - INFO - Removed 15 atoms. 15 atoms remaining in the molecule.

array([15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
      dtype=int32)

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')
On this page