Module 13: Topological Analysis#

Welcome back, Apprentice Master. In Module 12: Building, Repairing and Auditing Systems, you explored system quality auditing and structural repair with msm.build. Now we turn to non-spatial features: Topological Analysis with msm.topology.

In computational biophysics, a molecular system contains properties that are entirely independent of 3D spatial coordinates—such as covalent graph connectivity, connected molecular components, covalent path distances, and primary amino acid sequence ordering. While spatial coordinates change continuously during a molecular dynamics trajectory, the underlying chemical topology and residue sequence remain invariant.

The msm.topology function module provides specialized tools to inspect and analyze these non-spatial properties:

  • get_covalent_blocks(): Partitions connected covalent components and analyzes sub-blocks under hypothetical bond cuts.

  • get_covalent_paths(): Measures shortest topological bond paths between specified atom pairs.

  • get_sequence_alignment(): Aligns residue sequences across different molecular systems.

  • get_sequence_identity(): Computes percentage sequence identity between matching chains or systems.

1. Covalent Blocks and Paths#

Let’s begin by importing MolSysMT and loading our T4 Lysozyme demonstration system.

import molsysmt as msm
from molsysmt import systems

# Load T4 Lysozyme as native MolSys
lysozyme = msm.convert(systems['T4 lysozyme L99A']['181l.bcif.gz'], to_form='molsysmt.MolSys')

A covalent block is a set of atoms connected continuously through covalent bonds. When called without parameters, msm.topology.get_covalent_blocks() partitions the system into independent connected components:

# Partition system into connected covalent blocks
blocks = msm.topology.get_covalent_blocks(lysozyme)
n_components = msm.get(lysozyme, n_components=True)

print(f"Covalent blocks count: {len(blocks)}")
print(f"System components count: {n_components}")

Covalent blocks count: 141
System components count: 141

You can also use remove_bonds to simulate hypothetical bond cuts. For example, cutting an N-CA backbone bond splits a connected protein component into two sub-blocks:

# Select the N and CA atom indices of group index 20
n_atom = msm.select(lysozyme, selection='group_index == 20 and atom_name == "N"')[0]
ca_atom = msm.select(lysozyme, selection='group_index == 20 and atom_name == "CA"')[0]

# Partition blocks after hypothetically cutting the N-CA bond
blocks_after_cut = msm.topology.get_covalent_blocks(lysozyme, remove_bonds=[[n_atom, ca_atom]])

print(f"Blocks count before bond cut: {len(blocks)}")
print(f"Blocks count after N-CA cut : {len(blocks_after_cut)}")

Blocks count before bond cut: 141
Blocks count after N-CA cut : 142

To measure topological distances (the minimum number of covalent bonds separating two atoms), use msm.topology.get_covalent_paths():

# Calculate topological covalent path between atom 0 and atom 1
path_result = msm.topology.get_covalent_paths(lysozyme, path=[[0, 1]])
print(f"Covalent path array between atom 0 and 1:\n{path_result}")

Covalent path array between atom 0 and 1:
[[0]
 [1]]

Hint

msm.topology.get_covalent_blocks(): Returns sets of atoms remaining connected when specified bonds are cut. See API doc: molsysmt.topology.get_covalent_blocks().

2. Sequence Alignment and Identity#

Besides chemical bond networks, msm.topology includes sequence comparison utilities. You can align residue sequences and calculate sequence identity percentages between different molecular systems.

Let’s compare T4 Lysozyme (181L) against a related T4 Lysozyme variant structure (PDB ID 5X33):

# Load reference T4 Lysozyme variant
t4_variant = 'pdb_id:5X33'

# Compute sequence alignment
alignment = msm.topology.get_sequence_alignment('pdb_id:181L', reference_molecular_system=t4_variant)
print(f"Aligned T4 Lysozyme (181L) seq: {alignment[0][:60]}...")
print(f"Aligned T4 Variant  (5X33) seq: {alignment[1][:60]}...")

# Calculate sequence identity percentage
seq_id = msm.topology.get_sequence_identity('pdb_id:181L', reference_molecular_system=t4_variant)
print(f"Sequence Identity: {seq_id[0]:.2f}%")

Aligned T4 Lysozyme (181L) seq: ------------------------------------------------------------...
Aligned T4 Variant  (5X33) seq: SNTFIPLLAMILLSVSMVVGLPGNTFVVWSILKRMRKRSVTALMVLNLALADLAVLLTAP...

Sequence Identity: 52.65%

🏆 Challenge 13: The Sequence & Topology Master#

  1. Load T4 Lysozyme (pdb_id:181L) and the variant structure (pdb_id:5X33).

  2. Partition T4 Lysozyme into covalent blocks using msm.topology.get_covalent_blocks().

  3. Compute the sequence alignment between both structures using msm.topology.get_sequence_alignment().

  4. Verify that their sequence identity percentage is around 52.65% using msm.topology.get_sequence_identity().

Topological analysis bridges chemistry, graph theory, and bioinformatics. In Module 14: Comparing Systems, we will explore 3D structural alignment and RMSD calculations.