Python计算rmsd完全指南(附Python代码)

一、什么是rmsd？

RMSD 全称为 root-mean-square deviation，是一种用于测量同一分子不同构象之间结构相似度的指标。在药物设计、蛋白结构分析等领域得到广泛应用，通过计算不同构象之间每个原子在三维空间中的坐标差异，以此衡量两个结构的相似性。

二、rmsd的公式和计算方法

RMSD 的计算公式如下：

RMSD = sqrt( 1/n  * Σ (RMSD i)^2)

其中，n代表原子数，RMSDi代表第i个原子的RMSD。

计算步骤如下：

1、读取两个分子的三维坐标信息。

2、对应原子的坐标分别取出，计算它们之间的欧氏距离。

3、将不同原子间欧氏距离的平方和求平均数。

4、对求得的平均数开方即得到RMSD。

下面是一个计算rmsd的示例代码：

import numpy as np

def calc_rmsd(coords1, coords2):
    """计算坐标数据的RMSD值"""
    num_atoms = len(coords1)
    diff = coords1 - coords2
    rmsd = np.sqrt(np.sum(diff**2) / num_atoms)
    return rmsd

# 比较两个PDB文件中第一个model的坐标，并输出RMSD值
from Bio.PDB import Superimposer, PDBParser

pdb_parser = PDBParser(QUIET=True)
structure1 = pdb_parser.get_structure("X", "protein1.pdb")
structure2 = pdb_parser.get_structure("Y", "protein2.pdb")

model1 = structure1[0]
model2 = structure2[0]

coords1 = np.array([atom.get_coord() for atom in model1.get_atoms()])
coords2 = np.array([atom.get_coord() for atom in model2.get_atoms()])

print(calc_rmsd(coords1, coords2))

三、如何优化rmsd计算的速度

计算RMSD需要进行大量的计算，因此速度是一个重要的考虑因素。下面介绍几种优化RMSD计算速度的方法：

1、使用 numpy

numpy可以帮助我们更快速地处理坐标信息，实现向量化计算。将坐标转换为numpy数组后，使用向量化计算即可大大提高计算速度。如前面的示例代码所示，np.sum和np.sqrt函数可以直接对数组进行运算，速度会更快。

2、使用 Bio.PDB.Superimposer

使用Bio.PDB.Superimposer功能，我们可以在同一坐标系中对两个分子进行结构比较，然后计算RMSD指标。这个函数是C语言实现的，所以它非常快。

import numpy as np
from Bio.PDB import Superimposer, PDBParser

pdb_parser = PDBParser(QUIET=True)
structure1 = pdb_parser.get_structure("X", "protein1.pdb")
structure2 = pdb_parser.get_structure("Y", "protein2.pdb")

model1 = structure1[0]
model2 = structure2[0]

coords1 = np.array([atom.get_coord() for atom in model1.get_atoms()])
coords2 = np.array([atom.get_coord() for atom in model2.get_atoms()])

# 实例化 Superimposer 对象
superimposer = Superimposer()

# 设置参考目标，即目标坐标为 coords1
superimposer.set_atoms(model1.get_atoms(), model2.get_atoms())

# 获取旋转矩阵和平移向量
rot, tran = superimposer.rotran

# 用旋转矩阵和平移向量，计算出最小平方根距离
diff = coords1 - coords2
coords2_new = np.dot(diff, rot) + tran
rmsd = np.sqrt(np.sum((coords1 - coords2_new)**2) / len(coords1))

print(rmsd)

三、使用OpenMP进行多线程并行

多线程并行是另一个提高计算速度的方法。OpenMP是一种并行编程接口，能够使得多个线程同时执行同一个代码段。使用OpenMP时，只需在循环前添加两行指令，即可实现多线程并行计算。但是要注意，在Python中使用OpenMP需要安装numba模块，并启用自动并行化功能。

import numpy as np
from numba import jit, prange

@jit(nopython=True, parallel=True)
def calc_rmsd(coords1, coords2):
    """计算坐标数据的RMSD值"""
    num_atoms = len(coords1)
    diff = coords1 - coords2
    rmsd = np.sqrt(np.sum(diff**2) / num_atoms)
    return rmsd

# 比较两个PDB文件中第一个model的坐标，并输出RMSD值
from Bio.PDB import PDBParser

pdb_parser = PDBParser(QUIET=True)
structure1 = pdb_parser.get_structure("X", "protein1.pdb")
structure2 = pdb_parser.get_structure("Y", "protein2.pdb")

model1 = structure1[0]
model2 = structure2[0]

coords1 = np.array([atom.get_coord() for atom in model1.get_atoms()])
coords2 = np.array([atom.get_coord() for atom in model2.get_atoms()])

print(calc_rmsd(coords1, coords2))

四、如何可视化rmsd计算结果

我们可以通过绘制坐标点的三维空间图，来观察分子间结构的变化。本示例中，我们使用matplotlib库来绘制三维图形。下面是一个简单的示例代码：

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

from Bio.PDB import PDBParser

pdb_parser = PDBParser(QUIET=True)
structure1 = pdb_parser.get_structure("X", "protein1.pdb")
structure2 = pdb_parser.get_structure("Y", "protein2.pdb")

model1 = structure1[0]
model2 = structure2[0]

coords1 = np.array([atom.get_coord() for atom in model1.get_atoms()])
coords2 = np.array([atom.get_coord() for atom in model2.get_atoms()])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x1, y1, z1 = np.split(coords1, 3, axis=-1)
x2, y2, z2 = np.split(coords2, 3, axis=-1)

ax.scatter(x1, y1, z1, c='r', marker='o')
ax.scatter(x2, y2, z2, c='b', marker='^')

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

五、如何使用rmsd进行分子对齐

RMSD不仅可以用于衡量两个结构的相似性，还可以用于分子对齐。下面是一个简单的实现示例，使用Bio.PDB.Superimposer对两个分子进行对齐，并输出坐标信息和RMSD值：

from Bio.PDB import PDBParser, Superimposer

pdb_parser = PDBParser(QUIET=True)
structure1 = pdb_parser.get_structure("X", "protein1.pdb")
structure2 = pdb_parser.get_structure("Y", "protein2.pdb")

model1 = structure1[0]
model2 = structure2[0]

coords1 = [atom.get_coord() for atom in model1.get_atoms()]
coords2 = [atom.get_coord() for atom in model2.get_atoms()]

super_imposer = Superimposer()
super_imposer.set_atoms(model1.get_atoms(), model2.get_atoms())
super_imposer.apply(structure2[0].get_atoms())

rmsd = super_imposer.rms
print(f"RMSD: {rmsd}")

# Output aligned coordinates
print("Aligned Coordinates:")
for a1, a2 in zip(model1.get_atoms(), model2.get_atoms()):
    print(f"{a1}  {a2.get_coord()}")