转录组分析：读懂细胞城市的工作日报

可视化演示

演示RNA-seq数据比对与标准化、自编码器降维聚类发现疾病亚型，以及差异表达与通路富集三个核心环节。

转录组基础：先看懂城市的工作日报

请把每个细胞想象成一座分工精细的城市。两万多个基因就是两万多个岗位：有制造胰岛素的“食品厂工人”，有负责 DNA 修复的“电路检修员”，有传递信号的“外卖小哥”。每个岗位的职责都写成一份岗位说明书，锁在市中心档案室（细胞核）的保险柜里——这份说明书就是 DNA，它从不轻易离开档案室。当城市需要某个岗位开工时，档案室会复印一份说明书交给车间执行，这份复印本就叫 mRNA。此刻城市里所有正在流通的工单加在一起，就是这座城市的“工作日报”——转录组。

RNA-seq 就是清点这份日报的系统：先把城市里所有 mRNA 抽出来，反转录成稳定的 cDNA，上机测序得到一条条短读段，再把读段“贴回”参考基因组，数一数每个基因名下有多少条读段——读段越多，说明这个岗位今天越忙碌。但直接解读原始计数有两个陷阱：长的基因更容易被读到，测序量大的实验所有计数都会整体膨胀。就像统计工单时，印刷字体大的部门天然多占纸张，而日报总页数变多时所有部门看起来都在忙。

于是需要归一化，最常用的指标是 TPM（Transcripts Per Million，每百万转录本数），它同时校正基因长度和测序深度两个偏差：

$$\displaystyle \text{TPM}_i = \frac{r_i / l_i}{\sum_{j=1}^{n} r_j / l_j} \times 10^6$$

1. 计算读段密度：把基因 i 的读段数 r_i 除以基因长度 l_i，得到每碱基上的读段数，消除“长基因天然被多读到”的长度偏差。

$$\displaystyle \frac{r_i}{l_i}$$

2. 归一化为占比：用该基因的密度除以所有基因密度之和，表达它占城市全部工单数的比例，消除测序深度带来的整体放大效应。

$$\displaystyle \frac{r_i / l_i}{\sum_{j=1}^{n} r_j / l_j}$$

3. 放大到百万量级：乘上 10^6，得到“每百万条转录本中该基因占多少条”，数值更直观、便于跨样本比较。

$$\displaystyle \text{TPM}_i = f_i \times 10^6$$

有了 TPM 这把公平的尺子，跨城市、跨实验比较才成为可能。接下来，我们要把健康城市和生病城市的日报摆在一起，逐岗核对谁在疯狂加班、谁已彻底停工——这正是差异表达分析的工作。

差异表达分析：找出加班与停工的岗位

现在手头有两类城市：健康的和生病的。我们的目标是翻遍两万多个岗位的工作记录，找出哪些岗位在生病城市里明显多开工（上调表达）或彻底停工（下调表达）。但城市每天发工单的数量有随机波动，就像昨天全城外卖 3 万单、今天 3.2 万单，不能仅凭数字变化就断定外卖行业变繁荣了。这需要一场统计学审问：眼前这个差别，到底是真实差异，还是随机波动？

经典做法是对每个基因做两样本 t 检验：把两种城市里该基因的表达量看成两组观测值，比较它们的均值差是否显著大于组内波动。公式如下：

$$\displaystyle t_g = \frac{\bar{x}_{g,\text{case}} - \bar{x}_{g,\text{control}}}{\sqrt{\dfrac{s_{g,\text{case}}^2}{n_{\text{case}}} + \dfrac{s_{g,\text{control}}^2}{n_{\text{control}}}}}$$

1. 分子：两组平均表达之差：拿出患病城市与健康城市中该基因的平均表达量相减，反映岗位繁忙度的变化方向与幅度。

$$\displaystyle \bar{x}_{\text{case}} - \bar{x}_{\text{control}}$$

2. 分母：合并标准误：把两个城市的组内波动按各自样本量折算后相加再开方，衡量均值差的“测量噪声”有多大。

$$\displaystyle \sqrt{\frac{s_{\text{case}}^2}{n_{\text{case}}} + \frac{s_{\text{control}}^2}{n_{\text{control}}}}$$

3. 比值即检验统计量：均值差除以标准误得到 t 值，|t| 越大，p 值越小，越说明差异不太可能由随机波动造成。

$$\displaystyle |t_g| = \left|\frac{\bar{x}_{g,\text{case}} - \bar{x}_{g,\text{control}}}{\text{SE}}\right|$$

4. BH 校正控制假阳性：对全部基因的 p 值按秩调整得到 FDR（错误发现率），常与倍数变化阈值组合成双条件筛选。

$$\displaystyle q_{(i)} = \min_{k \ge i}\left\{ \frac{m \cdot p_{(k)}}{k} \right\}$$

经过这一轮“审问”，我们得到一张火山图：横轴是 log2 倍数变化（变化幅度），纵轴是 -log10 p 值（变化可信度），右上角和左上角的点就是值得关注的差异基因。但差异基因往往成百上千个，它们背后有没有隐藏的结构性规律？下一节，我们学习如何给城市画像、分班。

聚类分析：给城市的工作风格分班

假设你是教育局局长，要对 60 座城市的运行手册做调研，但每份手册有两万多个指标——人脑根本看不过来。聪明的做法是先抓“城市风格”的主干：有的城市重工业，有的重服务业，有的重科研。聚类分析就是这个“分班”过程：让工作模式相似的城市自然聚成一个班，全程不需要事先贴标签，属于无监督学习。

转录组数据动辄上万维，直接聚类会受噪声和冗余干扰，所以第一步通常用 PCA 或 UMAP 把高维表达谱压缩到可视化的低维空间，相当于只考几门“代表科目”来概括学生水平。然后在低维空间做 K 均值聚类，它的目标非常直白：找到 K 个聚类中心，让每个样本到所属中心距离的平方和最小：

$$\displaystyle \min_{\{\boldsymbol{\mu}_k\}} \sum_{k=1}^{K} \sum_{i \in C_k} \left\|\mathbf{x}_i - \boldsymbol{\mu}_k\right\|^2$$

1. 随机初始化 K 个中心：先随机放 K 个“班主任”站位，比如 K=3，代表我们想分三个班。

$$\displaystyle \boldsymbol{\mu}_1^{(0)},\ \boldsymbol{\mu}_2^{(0)},\ \boldsymbol{\mu}_3^{(0)}$$

2. 样本归属最近中心：每个样本（城市）划给距离最近的班主任，形成初始分班。

$$\displaystyle i \in C_k \iff \|\mathbf{x}_i - \boldsymbol{\mu}_k\|^2 \le \|\mathbf{x}_i - \boldsymbol{\mu}_{k'}\|^2,\ \forall k' \ne k$$

3. 用类内均值更新中心：每班所有学生的平均位置成为新班主任站位，班级“重心”随成员变化而移动。

$$\displaystyle \boldsymbol{\mu}_k = \frac{1}{|C_k|} \sum_{i \in C_k} \mathbf{x}_i$$

4. 迭代至收敛：重复“归类—更新中心”两步，直到班主任位置不再明显移动，输出稳定分班结果。

聚类让我们看到城市之间存在天然的“风格部落”。但把两万维数据硬压到两三个维度时，线性降维难免丢失非线性细节——有没有更聪明的压缩方式？下一节的自编码器，交给神经网络来学“怎么写摘要”。

自编码器：AI 笔下的城市摘要师

聚类之前我们用 PCA 做机械压缩，但 PCA 是线性变换，而城市的工作风格可能是高度非线性的：岗位之间的协作、代偿、级联都藏在工单数字里。这时可以请出一位“摘要师”——自编码器（Autoencoder）。它像个实习生，读完整份两万行的日报，只允许写 20 个字的摘要；但必须凭摘要尽可能还原出原日报。如果还原得像，说明那 20 个字真的抓住了精髓。

自编码器由编码器 f 和解码器 g 组成：编码器把高维表达向量 x_i 压缩成低维潜变量 z_i，解码器把 z_i 还原成重建向量 \hat{x}_i。训练目标就是最小化重建误差——原文与“摘要重建文”之间的差距：

$$\displaystyle \mathcal{L} = \frac{1}{n} \sum_{i=1}^{n} \left\|\mathbf{x}_i - g\big(f(\mathbf{x}_i)\big)\right\|_2^2$$

1. 编码器压缩：高维表达向量经过多层非线性网络，压成远小于原维度的潜变量 z_i。

$$\displaystyle \mathbf{z}_i = f(\mathbf{x}_i)$$

2. 解码器重建：潜变量再经过对称结构的网络，还原成与原维度相同的重建向量。

$$\displaystyle \hat{\mathbf{x}}_i = g(\mathbf{z}_i)$$

3. 最小化重建误差：用均方误差衡量原文与重建文的差距，梯度回传同时更新编码器和解码器参数。

$$\displaystyle \min_{f,g} \frac{1}{n} \sum_{i=1}^{n} \left\|\mathbf{x}_i - \hat{\mathbf{x}}_i\right\|_2^2$$

4. 潜变量即特征：训练完成后，把 z_i 取出来作为样本的低维特征，可用于聚类、疾病分型，甚至批次效应校正。

摘要师提炼出的低维特征，往往比原始两万维数据更干净、更抗噪。但岗位之间并非各自独立，它们构成了一张协作网——要真正读懂这张网，需要图神经网络登场。

图神经网络：透视岗位协作牵制网

城市里的岗位从来不是孤岛：外卖需求暴增，餐饮店、包装厂、交通调度都会跟着忙起来；某条供应链断裂，上下游部门连锁反应。要把这种“联动关系”写进模型，光靠向量不够，还需要一张图——节点是基因，边是它们之间的调控或共表达关系，这就是基因调控网络。

图神经网络（GNN）的核心操作叫消息传递：每一层，每个节点不仅要更新自己的特征，还要“听取”邻居节点的特征，把自己的状态与邻居状态融合。层数越多，节点能感知的“朋友圈”范围越广。一个典型的消息传递公式如下：

$$\displaystyle \mathbf{h}_v^{(l+1)} = \sigma\left(\mathbf{W}_{\text{self}} \mathbf{h}_v^{(l)} + \sum_{u \in \mathcal{N}(v)} \mathbf{W}_{\text{nei}} \mathbf{h}_u^{(l)}\right)$$

1. 自身特征变换：节点自己上一层的表示乘上权重矩阵，保留“历史存量”信息。

$$\displaystyle \mathbf{W}_{\text{self}} \mathbf{h}_v^{(l)}$$

2. 聚合邻居信息：把每个邻居 u 的特征乘权重后求和，相当于把一圈朋友圈的意见汇总起来。

$$\displaystyle \sum_{u \in \mathcal{N}(v)} \mathbf{W}_{\text{nei}} \mathbf{h}_u^{(l)}$$

3. 融合并激活：自身信息与邻居信息相加，过非线性激活函数 σ（如 ReLU），得到节点的新表示，进入下一层。

$$\displaystyle \sigma\left(\mathbf{W}_{\text{self}} \mathbf{h}_v^{(l)} + \sum_{u \in \mathcal{N}(v)} \mathbf{W}_{\text{nei}} \mathbf{h}_u^{(l)}\right)$$

有了 GNN，我们不仅能看单个基因的异常，还能看“异常如何通过网络扩散”——这比单独看基因更贴近疾病真实机制。不过故事还没完：成百上千个差异基因扎堆出现，会不会它们其实属于同一条流水线？下一节我们找一找“集体加班的部门”。

通路富集：定位集体加班的部门

假设城市日报显示：消防局、医院、交警队、应急管理局同时大量开工。单独看，四个部门各自在加班；但把名单放在一起，立刻明白——城市可能发生了重大事故。生物学里也有类似的“部门”：通路（Pathway）。KEGG、GO 等数据库把参与同一功能的基因编成基因集，就像把员工编入同一个部门目录。

通路富集分析要回答：我们筛出的差异基因，是显著集中出现在某个通路里，还是纯属随机凑巧？这可以用超几何检验定量：想象一个装满 N 个基因的大口袋，其中 M 个属于某通路；随机摸出 n 个差异基因，恰好有 k 个来自该通路的概率服从超几何分布。富集的 p 值就是“摸到 ≥k 个通路基因”的尾概率：

$$\displaystyle p = 1 - \sum_{i=0}^{k-1} \frac{\binom{M}{i} \binom{N-M}{n-i}}{\binom{N}{n}}$$

1. 定义四个关键参数：N 为全部基因数，M 为该通路基因数，n 为差异基因总数，k 为差异基因中属于该通路的个数。

$$\displaystyle N,\ M,\ n,\ k$$

2. 单点概率：任意抽 n 个基因中恰有 i 个通路基因的超几何概率，即“恰好重叠 i 个”的组合数占比。

$$\displaystyle \frac{\binom{M}{i} \binom{N-M}{n-i}}{\binom{N}{n}}$$

3. 累计尾概率：从总概率 1 中减去“重叠少于 k 个”的所有概率，剩下的就是“≥k 个重叠”的富集 p 值，越小越富集。

$$\displaystyle p = \Pr(\text{overlap} \ge k)$$

富集分析把“一堆基因”升华成“一条机制”，给疾病故事提供了叙事主线。现在所有线索都齐了——差异基因、低维特征、调控网络、富集通路——该用它们做点真正有价值的事：给疾病“分型”，走向精准医疗。

疾病分型应用：从一刀切到精准分班

同一个病名的患者，预后可能天差地别：有人对药物敏感，有人完全无效。就像同一本教材教出来的学生，学习风格各不相同。过去医学常按组织形态“一刀切”分型，如今我们希望用转录组数据驱动，把患者精确分班，让每个班匹配最合适的治疗方案，这就是疾病分型（如乳腺癌的 PAM50 亚型、肺癌的免疫亚型）。

一套典型流程是：先做差异表达与富集分析，筛选一组标志基因签名；再用自编码器提取低维特征，配合聚类或分类器把患者划入亚型；最后用生存分析与临床指标评估分型是否有预后价值。给患者打分常用的“签名分数”公式：

$$\displaystyle S_i = \sum_{j=1}^{m} w_j \cdot z_{ij}$$

1. 挑选标志基因集合：从差异基因与富集通路里挑出 m 个最具代表性的基因，构成疾病标签签名。

$$\displaystyle G_{\text{sig}} = \{g_1, g_2, \dots, g_m\}$$

2. 标准化消除量纲：每个基因在不同患者间做 z-score，消除绝对表达量差异，让权重可以跨基因比较。

$$\displaystyle z_{ij} = \frac{x_{ij} - \bar{x}_j}{\sigma_j}$$

3. 加权求和得分数：把患者 i 在各标志基因上的标准化表达乘权重 w_j 后累加，得到分型分数。

$$\displaystyle S_i = w_1 z_{i1} + w_2 z_{i2} + \dots + w_m z_{im}$$

4. 划分并验证亚型：设定阈值或结合聚类结果把患者划入不同亚型，再用生存曲线、AUC 等指标验证分型效果。

$$\displaystyle \text{if } S_i \le \tau \Rightarrow \text{亚型 A},\ \text{否则为亚型 B}$$

从读日报、找差异、聚类画像、神经网络压缩、网络透视、通路定位，到分型落地——最初那座混沌的细胞城市，已经在 AI 的帮助下变得井井有条。这就是计算生物学与精准医学结合的美妙之处：数据驱动，但每一步都映照着真实的生物学机制。

AI 应用场景

Python 演示

import math import numpy as np  def simulate_rnaseq(n_genes=400, n_ctrl=30, n_case=30, seed=42):    """模拟两个组的 RNA-seq 计数矩阵, 并埋入 25 个真实差异基因."""    rng = np.random.default_rng(seed)    ctrl = rng.poisson(lam=50, size=(n_genes, n_ctrl)).astype(float)    case = rng.poisson(lam=50, size=(n_genes, n_case)).astype(float)    truth = rng.choice(n_genes, size=25, replace=False)    case[truth, :] = np.round(case[truth, :] * 2.5)    return ctrl, case, truth  def diff_express(ctrl, case):    """对每个基因做两样本 t 检验, 返回 p 值与 log2 倍数变化."""    n = ctrl.shape[0]    pvals = np.ones(n)    log2fc = np.zeros(n)    for g in range(n):        a, b = ctrl[g] + 1.0, case[g] + 1.0        na, nb = len(a), len(b)        va, vb = a.var(ddof=1), b.var(ddof=1)        t = abs(a.mean() - b.mean()) / np.sqrt(va / na + vb / nb)        pvals[g] = 2 * (1 - 0.5 * (1 + math.erf(t / math.sqrt(2))))        log2fc[g] = np.log2(b.mean() / a.mean())    return pvals, log2fc  def bh_correct(pvals, fdr=0.05):    """BH 方法做多重假设校正, 返回是否显著的真值数组."""    order = np.argsort(pvals)    q = pvals[order] * len(pvals) / np.arange(1, len(pvals) + 1)    q = np.minimum.accumulate(q[::-1])[::-1]    qvals = np.empty_like(pvals)    qvals[order] = np.minimum(q, 1.0)    return qvals < fdr  def kmeans(X, k=2, iters=20, seed=0):    """对样本表达谱做 K 均值聚类, 返回每个样本的簇标签."""    rng = np.random.default_rng(seed)    centers = X[rng.choice(len(X), k, replace=False)]    labels = np.zeros(len(X), dtype=int)    for _ in range(iters):        dists = np.linalg.norm(X[:, None, :] - centers[None, :, :], axis=2)        labels = dists.argmin(axis=1)        for j in range(k):            if np.any(labels == j):                centers[j] = X[labels == j].mean(axis=0)    return labels  def main():    ctrl, case, truth = simulate_rnaseq()    pvals, log2fc = diff_express(ctrl, case)    de_idx = np.where(bh_correct(pvals) & (np.abs(log2fc) >= 1.0))[0]    hits = int(np.isin(de_idx, truth).sum())    print("DE genes:", len(de_idx), "| true hits:", hits)    sig = np.hstack([ctrl[de_idx], case[de_idx]]).T    labels = kmeans(sig, k=2)    acc = max((labels[:30] == 0).mean(), (labels[:30] == 1).mean())    print("Sample clustering accuracy: %.1f%%" % (acc * 100))  if __name__ == "__main__":    main()

1. case[truth, :] = np.round(case[truth, :] * 2.5)把 25 个“真差异基因”在患病组放大 2.5 倍，埋入已知答案供后续验证，就像事先知道哪几个岗位真的在加班。
2. a, b = ctrl[g] + 1.0, case[g] + 1.0给计数加 1 再做对数处理，避免零计数导致除零或无穷大。
3. t = abs(a.mean() - b.mean()) / np.sqrt(va / na + vb / nb)两样本 t 检验统计量：均值差除以合并标准误，是“差异表达”核心审问公式的代码实现。
4. pvals[g] = 2 * (1 - 0.5 * (1 + math.erf(t / math.sqrt(2))))用正态近似计算双侧 p 值，等价于查 t 分布表。
5. q = pvals[order] * len(pvals) / np.arange(1, len(pvals) + 1)把 p 值排序后按秩放大，是 BH 多重假设校正的关键一步。
6. np.minimum.accumulate(q[::-1])[::-1]逆向累积最小值，保证校正后的 q 值单调不减，符合 FDR 的数学定义。
7. dists = np.linalg.norm(X[:, None, :] - centers[None, :, :], axis=2)利用广播机制一次性算出每个样本到所有聚类中心的欧氏距离，对应 K 均值“就近分班”步骤。
8. de_idx = np.where(bh_correct(pvals) & (np.abs(log2fc) >= 1.0))[0]FDR 与倍数变化双阈值筛选差异基因，呼应正文“既要可信度、又要幅度”的筛选逻辑。

常见误区

面试题