📚 课程目录第3章距离矩阵与接触图

距离矩阵与接触图：从三维坐标到接触图的计算

嗯，这一章咱们来聊聊蛋白质结构里最基础、也最核心的计算——从三维坐标到接触图。说白了，就是把那些原子在空间里的位置信息，转换成一张能直接看到残基之间亲密关系的图。

我刚开始做这个的时候，总觉得坐标文件读进来直接算距离就行了呗。结果第一次算出来，画出来的图跟论文上完全对不上——有的地方该有的接触没有，不该有的却冒出一大片。后来才明白，这里头有不少坑。

距离矩阵是什么？

先讲清楚距离矩阵这个概念。假设你的蛋白质有N个残基，每个残基选一个代表原子（通常是Cα或者Cβ），那么这N个点两两之间的距离，就构成了一个N×N的对称矩阵。说白了，就是一张“谁和谁离得远、谁和谁离得近”的表。

矩阵的(i, j)元素 = 残基i和残基j代表原子之间的欧氏距离
因为是距离，所以对角线全是0（自己和自己距离为0）
矩阵是对称的：d(i,j) = d(j,i)

我习惯把距离矩阵当成蛋白质结构的“指纹”。你看，不管这个蛋白质怎么旋转平移，距离矩阵里每个元素的值是不变的！这就是为什么很多结构比对算法都基于距离矩阵——它天然具有旋转平移不变性。

实际小技巧：千万别用所有原子去算距离矩阵。通常只取Cα（主链碳原子）就够了，这样矩阵大小是残基数×残基数，非常简洁。用全原子的话，几千几万个原子算下来，矩阵能把你内存撑爆——我曾经用全原子算一个2000残基的蛋白，MATLAB直接崩了。

从距离矩阵到接触图

距离矩阵是连续的数值，接触图是二值化的。我们需要一个阈值——比如8Å（埃）。如果两个残基之间的距离小于这个阈值，我们就把它们标记为“接触”（1），否则就是“不接触”（0）。

你看，就这么简单。但这里的门道可不少。

参数	典型值	说明
阈值（threshold）	8 Å	最常用，也被称为Cα接触距离
残基对范围	\|i - j\| >= 4	避免相邻残基的天然近邻干扰
代表原子	Cα 或 Cβ	Cα用于主链，Cβ用于侧链
矩阵对称性	自动对称	接触关系无方向性

为什么要排除相邻残基？你想想看，肽链里相邻的两个残基，距离本来就不到4Å，算上它们，接触图里对角线附近会有一圈全是1，根本看不出有价值的远距离接触信息。所以标准的做法是忽略序列上相差少于4个残基的对。

我曾经掉过的坑：有些工具默认不排除相邻残基，画出来的接触图对角线附近全是深色块，乍一看好像接触很多，实际上全是“伪接触”。后来我写代码时都会手动加一个条件：if abs(i - j) >= 4，才纳入计算。

计算流程：一步一步来

这里我把从PDB文件到接触图的完整流程拆解一下：

读取PDB文件，提取所有残基的Cα原子坐标
构建N×3的坐标矩阵（N个残基，每个有x,y,z）
计算两两欧氏距离，得到距离矩阵
应用阈值和序列间隔过滤，生成二值接触矩阵
可视化——用热图展示接触图

核心代码大概长这样（Python + NumPy实现）：

import numpy as np
from Bio.PDB import PDBParser

# 读取PDB，提取Cα坐标
parser = PDBParser(QUIET=True)
structure = parser.get_structure("protein", "1crn.pdb")
coords = []
for model in structure:
    for chain in model:
        for residue in chain:
            if residue.get_id()[0] != ' ':
                continue
            try:
                ca = residue['CA']
                coords.append(ca.get_vector().get_array())
            except KeyError:
                continue
coords = np.array(coords)  # shape: (N, 3)

# 计算距离矩阵
n = len(coords)
dist_matrix = np.zeros((n, n))
for i in range(n):
    diff = coords - coords[i]   # 利用广播
    dist_matrix[i] = np.sqrt(np.sum(diff**2, axis=1))

# 生成接触图（阈值8Å，排除序列距离<4）
contact_map = (dist_matrix <= 8.0).astype(int)
for i in range(n):
    for j in range(n):
        if abs(i - j) < 4:
            contact_map[i, j] = 0

这段代码我用了广播技巧，比三层循环快了两个数量级。实际项目里，一个300残基的蛋白，用这方法几十毫秒就能出结果。

关键提醒：上面的代码里，我用了 abs(i - j) < 4 来排除相邻残基。但注意，这里的“序列距离”是以残基编号为准的。如果你的PDB里残基编号不连续（比如1,2,3,5,6），那abs(i-j)就不等于真实序列差。我遇到过一个坑：PDB里有些残基编号跳号，结果我的过滤条件把不该排除的排除了。后来我改用残基索引序列而不是绝对编号来算间隔，才避免了这个bug。

一个完整的计算框架图

我把整个从坐标到接触图的流程画了个图，方便你理解各个模块的关系：

这张图里，我特意用虚线标出了另一条路径——其实你完全可以直接对每对残基算距离，不用先构建坐标矩阵。但用NumPy向量化先把坐标矩阵建好，后面算距离的时候能利用广播加速，处理几百个残基的蛋白非常快。

一些实用细节

接触图计算里，有几个细节我觉得值得专门拿出来讲：

阈值选择：8Å是主流，但不同应用可能有变化。比如一些远程接触预测会用12Å、15Å。我见过有的论文用“每对残基的最小距离”而不是Cα距离，结果完全不同。选什么取决于你关注的是骨架接触还是侧链接触。
接触图的稀疏散度：一般蛋白质的接触图里，约5%~10%的条目是1（接触），其余都是0。这个稀疏性后续在机器学习建模时很重要，直接拿全矩阵喂神经网络会有大量冗余。
对称性与三角：接触图是对称矩阵，所以我们通常只用上半三角或下半三角存储。用numpy的triu或tril就能提取，节省一半空间。

我个人的编码习惯：算完接触图后，立刻用 np.count_nonzero(contact_map) / (n*n/2) 算一下接触密度。如果密度大于20%，很可能是阈值设得太大了，或者没正确排除相邻残基。这个检查能帮你快速发现bug。

再从接触图回到距离矩阵

你可能会问：既然接触图只是距离矩阵二值化的结果，那为啥不直接拿距离矩阵做分析，还要多此一举转成二值的？

嗯，我刚开始也这么想。后来在实际项目中发现了几个关键原因：

二值化的接触图更容易可视化。 距离矩阵是连续灰度，人眼看不出细微差别。接触图只有黑白两色，哪里接触、哪里分离一目了然。
很多算法（比如接触预测、拓扑分析）天然需要0/1输入。 比如神经网络结构预测中，接触图是标准的监督信号，连续距离反而不好用。
接触图对实验误差更鲁棒。 距离矩阵里一个原子位置偏差0.5Å，所有距离都会变，但用8Å阈值一过滤，只要没跨过阈值线，结果就不变。我处理过同源建模的坐标，不同模型主链偏差1-2Å，距离矩阵差异很大，但接触图的核心拓扑基本一致。

你看，这就是为什么接触图在结构生物学里这么受欢迎——它抓住了结构的关键特征，同时扔掉那些不重要的细节。