""" ============================================================================= 第 17 讲 · 共享模块：一个"合成蛋白质宇宙"——家族、进化树、数据库的冗余、聚类 ============================================================================= 对应讲义：lectures/lecture17_biology_data.md §0.4 ----------------------------------------------------------------------------- 【这个世界长什么样】 ----------------------------------------------------------------------------- 蛋白质 = 长度 L = 40 的氨基酸序列（20 种氨基酸）。 家族（family）= 一组有共同祖先的序列。每个家族有一个"适应度景观"（能量函数，越低越好）： E(s) = Σ_i h_i(s_i) + Σ_{(i,j)∈接触} J_ij(s_i, s_j) - h_i：位点 i 对每种氨基酸的偏好（有的位点高度保守，例如活性中心；有的位点随便）； - J_ij：三维结构里相互接触的位点对之间的耦合——一个位点突变了，和它接触的位点也要跟着变（"共进化"）。 这就是生物信息学里的 Potts 模型（EVcouplings、GREMLIN、AlphaFold 之前的接触预测都基于它）。 进化：从祖先序列出发，沿一棵进化树分叉；每条分支上做若干次"突变 + 自然选择" （随机改一个位点，按 Metropolis 规则接受：变好总是接受，变差以 exp(−ΔE) 的概率接受）。 于是家族内部有亚家族的结构，序列之间的相似度从 30% 到 100% 不等。 数据库：家族大小是长尾的（少数家族有几千条序列，多数只有几十条）；另外，一些"被研究得多"的分支 被反复测序，留下大量几乎相同的序列——这就是 UniProt 的冗余，UniRef100 / 90 / 50 要解决的问题。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 protein_world.py # 自检：数据库的规模和 UniRef 式的聚类（约 3 秒） 依赖 numpy 和第 1 讲的 console.py。 """ import os import sys import numpy as np _HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.join(_HERE, "..", "lecture01")) from console import print_header, print_observations, row # noqa: E402 AA = "ACDEFGHIKLMNPQRSTVWY" Q = 20 # 氨基酸种类 L = 40 # 序列长度 # ============================================================================= # 1. 家族：能量函数与进化 # ============================================================================= class Family: """ 一个蛋白质家族：能量函数 E(s) = Σ h_i(s_i) + Σ J_ij(s_i, s_j)，一个祖先序列，以及沿进化树产生的成员。 【参数】rng：随机数生成器；n_contacts：接触对的个数；conserve：保守位点的比例； j_sd：耦合的强度（越大，接触位点之间的共进化信号越强） 【思路】 - h：每个位点先抽一个 N(0, 0.5²) 的偏好；保守位点再给一种"正确"的氨基酸 −3 的奖励（能量越低越好）。 - J：随机选 n_contacts 对相距 ≥ 3 的位点作为"接触"，每对有一个 20 × 20 的耦合矩阵 N(0, j_sd²)。 - 祖先：从随机序列出发做 400 步贪心的 Metropolis，得到一个适应度较高的序列。 """ def __init__(self, rng, n_contacts=30, conserve=0.3, j_sd=0.8): self.h = rng.normal(0, 0.5, (L, Q)) cons = rng.random(L) < conserve self.h[np.flatnonzero(cons), rng.integers(0, Q, cons.sum())] -= 3.0 pairs = set() while len(pairs) < n_contacts: i, j = sorted(rng.choice(L, 2, replace=False)) if j - i >= 3: pairs.add((i, j)) self.pairs = np.array(sorted(pairs)) self.J = rng.normal(0, j_sd, (len(self.pairs), Q, Q)) # 每个位点的"邻居"列表：用于快速计算单点突变的能量变化 self.nbr = [[] for _ in range(L)] for k, (i, j) in enumerate(self.pairs): self.nbr[i].append((k, j, 0)) self.nbr[j].append((k, i, 1)) anc = rng.integers(0, Q, (1, L)) self.ancestor = self.evolve(anc, 400, rng, temp=0.3)[0] def energy(self, S): """一批序列 S (n, L) 的能量 (n,)。""" E = self.h[np.arange(L), S].sum(axis=1) for k, (i, j) in enumerate(self.pairs): E += self.J[k, S[:, i], S[:, j]] return E def delta(self, S, pos, new): """ 把每条序列 S[n] 的位点 pos[n] 改成 new[n] 时能量的变化 ΔE (n,)。 【思路】只有 h_pos 和与 pos 接触的 J 会变。为了向量化，对每个可能的位点分组处理。 【边界】new 等于原来的氨基酸时 ΔE = 0。 """ n = len(S) old = S[np.arange(n), pos] dE = self.h[pos, new] - self.h[pos, old] for p in np.unique(pos): m = pos == p for k, other, side in self.nbr[p]: o = S[m, other] if side == 0: dE[m] += self.J[k, new[m], o] - self.J[k, old[m], o] else: dE[m] += self.J[k, o, new[m]] - self.J[k, o, old[m]] return dE def evolve(self, S, steps, rng, temp=1.0): """ 对每条序列做 steps 次"随机突变 + Metropolis 选择"。 【参数】temp：选择的"温度"（越小选择越强） 【返回】新的序列数组（不修改输入） 【数据联系】真实的进化 = 突变 + 选择 + 漂变。天然序列因此"编码"了适应度景观： 保守的位点说明它重要，一起变化的位点对说明它们在三维结构里相互接触。 蛋白质语言模型（ESM）和 AlphaFold 的 MSA 输入，学的都是这些信号。 """ S = S.copy() n = len(S) for _ in range(steps): pos = rng.integers(0, L, n) new = rng.integers(0, Q, n) dE = self.delta(S, pos, new) acc = (dE <= 0) | (rng.random(n) < np.exp(-np.maximum(dE, 0) / temp)) S[np.arange(n)[acc], pos[acc]] = new[acc] return S def members(self, n, rng, depth=6, steps_per_level=10): """ 沿一棵二叉进化树产生 n 条成员序列。 【思路】从祖先出发，每一层每条序列分裂成两条，各自突变 steps_per_level 步；depth 层之后得到 2^depth 个"物种"，再从它们里面（带最后一点突变）抽 n 条。 越晚分开的两条序列越相似——这给出了家族内部的层级结构（亚家族、近缘种、菌株）。 """ S = self.ancestor[None, :] for _ in range(depth): S = self.evolve(np.repeat(S, 2, axis=0), steps_per_level, rng) pick = rng.integers(0, len(S), n) return self.evolve(S[pick], 2, rng), pick def make_database(n_families, rng, min_size=20, zipf=1.2, max_size=1500, redundant_frac=0.3): """ 一个"序列数据库"：n_families 个家族，家族大小服从长尾分布；再给一部分分支加上大量近乎相同的序列。 【参数】zipf：家族大小的长尾程度；redundant_frac：被"反复测序"的家族的比例 【返回】dict：seqs (N, L)、fam (N,) 家族编号、families（Family 对象的列表） 【数据联系】UniProt 里，被研究得多的物种（人、小鼠、大肠杆菌、病原体）和热门的蛋白质有大量几乎相同的条目； 宏基因组让少数大家族变得极其庞大。直接在原始数据库上训练，模型的"注意力"会被这些大家族占据。 """ fams, seqs, fam_id = [], [], [] sizes = np.minimum(max_size, (min_size * (1 + rng.pareto(zipf, n_families)))).astype(int) for f in range(n_families): fam = Family(rng) S, _ = fam.members(sizes[f], rng) if rng.random() < redundant_frac: # 某个分支被反复测序：从 5 条序列出发，复制出很多只差 0–2 个位点的"菌株" base = S[rng.integers(0, len(S), 5)] extra = np.repeat(base, max(sizes[f] // 2, 10), axis=0) mut = rng.random(extra.shape) < 0.02 extra[mut] = rng.integers(0, Q, mut.sum()) S = np.concatenate([S, extra]) fams.append(fam) seqs.append(S) fam_id.append(np.full(len(S), f)) return dict(seqs=np.concatenate(seqs), fam=np.concatenate(fam_id), families=fams) # ============================================================================= # 2. 序列相似度与聚类（UniRef 式） # ============================================================================= def identity(a, B): """序列 a (L,) 和一批序列 B (n, L) 的相同位点比例（等长、不比对）。""" return (B == a).mean(axis=1) def greedy_cluster(S, threshold, rng): """ 贪心聚类（CD-HIT / MMseqs2 的简化版）：随机顺序遍历序列，和已有的某个代表序列的相似度 ≥ threshold 就并入它的簇，否则自己成为一个新的代表。 【参数】S：(n, L)；threshold：相似度阈值（UniRef90 = 0.9，UniRef50 = 0.5） 【返回】(labels (n,), reps (代表序列的下标)) 【思路】UniRef100 合并完全相同的序列；UniRef90 把相似度 ≥ 90% 的序列聚成一簇；UniRef50 把 UniRef90 的代表 再按 50% 聚类。每一级都大幅减少冗余。真实的工具用比对（序列长度不同）和快速的预过滤；这里序列等长，直接逐位比较。 【边界】复杂度是 O(n × 簇数 × L)。只在每个家族内部聚类（不同家族之间的相似度约等于随机的 5%）。 """ order = rng.permutation(len(S)) labels = -np.ones(len(S), int) reps = [] rep_seqs = np.zeros((0, L), dtype=S.dtype) for idx in order: if len(reps): sim = identity(S[idx], rep_seqs) j = int(np.argmax(sim)) if sim[j] >= threshold: labels[idx] = j continue labels[idx] = len(reps) reps.append(idx) rep_seqs = np.vstack([rep_seqs, S[idx]]) return labels, np.array(reps) def cluster_database(db, threshold, rng): """对整个数据库按家族分别聚类，返回全局的簇编号 (N,) 和每个簇的代表序列下标。""" labels = np.zeros(len(db["seqs"]), int) reps = [] offset = 0 for f in np.unique(db["fam"]): idx = np.flatnonzero(db["fam"] == f) lab, rp = greedy_cluster(db["seqs"][idx], threshold, rng) labels[idx] = lab + offset reps.extend(idx[rp]) offset += len(rp) return labels, np.array(reps) # ============================================================================= # 自检 # ============================================================================= if __name__ == "__main__": print_header("protein_world 自检：一个 150 个家族的合成序列数据库，以及 UniRef 式的聚类") rng = np.random.default_rng(0) db = make_database(150, rng) N = len(db["seqs"]) sizes = np.bincount(db["fam"]) uniq = len(np.unique(db["seqs"], axis=0)) c90, _ = cluster_database(db, 0.9, np.random.default_rng(1)) c50, _ = cluster_database(db, 0.5, np.random.default_rng(1)) w = [34, 14] print(row(["统计", "数值"], w)) for name, val in [("序列总数（原始数据库）", N), ("去掉完全相同的序列（UniRef100 式）", uniq), ("按 90% 相似度聚类的簇数（UniRef90 式）", c90.max() + 1), ("按 50% 相似度聚类的簇数（UniRef50 式）", c50.max() + 1), ("最大的家族 / 家族大小的中位数", f"{sizes.max()} / {int(np.median(sizes))}")]: print(row([name, val], w)) same = [identity(db["seqs"][i], db["seqs"][db["fam"] == db["fam"][i]]).mean() for i in range(0, N, max(N // 200, 1))] print(f"\n 同一家族内两条序列的平均相似度：{np.mean(same):.0%}（不同家族之间约 {1 / Q:.0%}）") print_observations([ f"原始数据库 {N} 条序列，按 50% 相似度聚类后只剩 {c50.max() + 1} 个簇：大部分序列是冗余的。" "真实的 UniRef：UniRef100 / 90 / 50 一级比一级小；ESM-2 在 UniRef50 的簇上均匀采样，再从簇里随机取一条 UniRef90 的序列。", f"家族大小是长尾的（最大 {sizes.max()}，中位数 {int(np.median(sizes))}）：直接在原始数据库上训练，大家族会占据大部分训练量。", ])