1. 从“形状”入手:为什么基因和癌症数据需要TDA?

大家好,我是老张,在AI和生物信息交叉领域摸爬滚打了十来年。今天想和大家聊聊一个听起来有点“高冷”,但实际威力巨大的工具——拓扑数据分析,也就是TDA。如果你正在处理基因表达数据、癌症样本分类,或者任何高维、复杂的生物数据,感觉用传统方法总是隔靴搔痒,抓不住那些微妙的关键信息,那这篇文章可能就是为你准备的。

我们先从一个简单的比喻开始。想象一下,你面前有一大团乱麻(这就是你的高维数据,比如几万个基因在数百个样本中的表达量)。传统的分析方法,比如主成分分析或者各种聚类算法,就像是用手去捋这团麻线,试图把它捋成几股清晰可见的粗绳。它们能帮你看到最明显的几个大类,比如把癌症样本和正常样本大致分开。这当然很有用,但问题在于,癌症不是铁板一块。比如同样是肺癌,驱动基因突变不同、对药物反应不同的亚型,可能就像乱麻中几根颜色、质地都略有不同,但又紧紧缠绕在一起的细丝。用手捋(传统降维/聚类)很容易就把这些细微但至关重要的差别给忽略或抹平了。

TDA的思路完全不同。它不急于把数据“压扁”或强行分成几大块,而是像一个精密的扫描仪,去探究这团数据本身固有的“形状”和“结构”。这个“形状”,在数学上对应的就是拓扑性质:比如数据中是否存在“空洞”、“隧道”、“分支点”。这些结构特征,往往对应着数据中隐藏的、连续变化的规律或者多个微小但稳定的亚群。在基因和癌症研究中,一个“空洞”可能意味着某种特定的基因表达模式在样本空间中天然缺失,这或许指向一种罕见的、但生存预后截然不同的癌症亚型;一条“隧道”可能连接着两种看似不同的状态,揭示了癌细胞从一种亚型向另一种亚型演变的潜在路径。

我刚开始接触TDA时,也被它的数学背景吓到过,但实际用起来你会发现,它的核心思想非常直观。它不做“硬切割”,而是通过一种叫Mapper的算法,像用一张带重叠网格的渔网去“打捞”你的数据,在每个网格局部进行精细聚类,再根据重叠部分把局部聚类连接起来,最终形成一张展现数据整体拓扑结构的网络图。这张图,就是理解你数据“形状”的钥匙。它能帮你发现那些样本量很少、但生物学意义重大的稀有亚型,也能揭示不同亚型之间非线性的过渡关系,这些都是传统方法很难甚至无法做到的。

2. 实战拆解:Mapper算法如何“编织”数据的拓扑网络?

上一节我们打了比方,说TDA像用一张特殊的网去打捞数据。这一节,我们就来亲手拆解一下这张“网”是怎么织出来的,也就是Mapper算法的具体步骤。别担心,我会尽量避开复杂的公式,用操作流程和代码示例带你走一遍。

### 2.1 第一步:选择一个“透镜”去看数据

Mapper的第一步,是选择一个“滤镜函数”或者叫“透镜”。这步的目的是把高维数据先映射到一个更低维、更容易处理的空间,为我们后续的“织网”提供一个视角。这个透镜不是随便选的,它决定了你将从哪个角度去观察数据的形状。

常见的透镜有:

  • 主成分:使用PCA的第一或前几个主成分。这是最常用的方法之一,能抓住数据最大方差的方向。
  • 数据密度:计算每个数据点周围的密度。这对于发现密度差异大的簇特别有效。
  • 自定义函数:你可以用任何有意义的指标。比如在癌症研究中,用一个与患者生存期高度相关的基因表达评分作为透镜,你就能直接看到数据在“预后”这个维度上的拓扑结构。

在Python的Giotto-tda库中,这一步非常灵活。你甚至可以用UMAP或t-SNE的嵌入结果作为透镜,虽然它们本身是非线性降维,但作为Mapper的输入透镜,能带来更丰富的视角。

# 示例:使用PCA第一主成分作为透镜
from sklearn.decomposition import PCA
import numpy as np

# 假设 X 是你的基因表达矩阵,形状为 (n_samples, n_genes)
lens = PCA(n_components=1).fit_transform(X)  # 输出形状 (n_samples, 1)
# lens 现在就是每个样本的一维投影值

### 2.2 第二步:覆盖与划分——设置网格的“分辨率”和“重叠度”

有了透镜投影后的一维(或二维)数据,我们就要在上面铺“网格”了。这里有两个超参数至关重要,直接决定了你最终拓扑图的精细程度:

  1. resolution(分辨率/间隔数):把透镜值的范围分成多少个区间。比如分成10个区间。这个数越大,网格越密,捕捉的细节越多,但图也可能更复杂。
  2. overlap(重叠比例):相邻区间之间重叠的比例。比如设为0.2(20%)。这个参数是Mapper的灵魂! 正是重叠部分的存在,使得数据点可以同时属于两个相邻的区间,从而在后续步骤中将不同的局部簇连接起来,形成复杂的拓扑结构。没有重叠,你得到的就只是一堆孤立的簇,而不是一个网络。

在Giotto-tda中,我们可以这样设置:

from gtda.mapper import CubicalCover

cover = CubicalCover(
    n_intervals=10,    # 分辨率:分成10个区间
    overlap_frac=0.3   # 重叠度:30%的重叠
)

### 2.3 第三步:局部聚类——在每个网格里“精耕细作”

现在,对于上一步划分出的每一个区间(比如第i个区间),我们找到所有透镜值落在这个区间内的原始高维数据点(这叫原像 f^{-1}(U_i)),然后在这个数据子集内部进行聚类。你可以使用任何你熟悉的聚类算法,比如DBSCAN、K-Means、层次聚类等。

这一步是关键:全局上我们只用了一个简单的透镜,但在每个局部区间内,我们却是在原始高维空间中进行聚类。这意味着,我们既利用了透镜提供的全局引导视角,又没有丢失高维数据本身的丰富信息。DBSCAN在这里特别受欢迎,因为它不需要预先指定簇的个数,能基于密度自然发现形状各异的簇。

from sklearn.cluster import DBSCAN

# 在每个区间内进行聚类
clusterer = DBSCAN(eps=0.5, min_samples=3)
# eps和min_samples需要根据你的数据特性进行调整

### 2.4 第四步:构建拓扑网络——连接就是发现

经过上一步,每个区间内部都产生了若干个簇。我们把每个簇变成一个拓扑图的“节点”。然后,我们检查任意两个节点:如果它们包含至少一个共同的数据点(这得益于第二步设置的重叠区间!),我们就在这两个节点之间连一条“边”。

最终,所有节点和边就构成了一张拓扑网络图。这张图就是你的数据“形状”的可视化呈现。一个节点代表一群在局部高维空间中相似的数据点(样本),节点的大小通常代表该簇中样本的多少。连接两个节点的边,则代表这两群样本之间存在“过渡”或“相似”的样本,可能意味着生物学上的连续状态或中间表型。

3. 癌症研究中的“火眼金睛”:TDA如何发现隐藏的亚型?

理论说了这么多,TDA在真实的基因和癌症研究中到底能干什么?我结合自己参与过和文献中看到的一些案例,给大家讲讲它如何扮演“火眼金睛”的角色。

### 3.1 案例一:在乳腺癌数据中揪出罕见但凶险的亚群

几年前,我们团队分析一个公开的乳腺癌基因表达数据集,里面包含上千个样本,已有经典的PAM50分型(Luminal A, Luminal B, HER2-enriched, Basal-like等)。我们用常规的聚类和t-SNE可视化,结果确实能把几大亚型分开,但图里总有一些点散落在主要集群之间或边缘,说不清道不明。

后来我们尝试了TDA。我们选择一组与细胞周期和免疫应答相关的基因签名作为透镜,用Mapper生成了拓扑网络。结果非常震撼:在主要对应于Basal-like亚型的大节点旁边,我们发现了两个非常小、但连接结构独特的节点。这两个节点里的样本,在PAM50分型里也被标记为Basal-like,但在我们的拓扑图中,它们通过很长的“路径”才与主节点相连,且彼此之间有直接连接。

深入分析这两个小节点的样本,我们发现它们的突变谱很特别:一个富集了特定的DNA损伤修复基因突变,另一个则表现出强烈的上皮-间质转化特征。更重要的是,临床数据回溯显示,这两个小节点患者的无进展生存期显著短于主Basal-like群体。TDA成功地从传统分类认为的“同一类”里,剥离出了预后更差的稀有亚型。这对于后续开展精准治疗,比如为这些患者优先考虑PARP抑制剂或特定靶向药,提供了关键线索。

### 3.2 案例二:描绘癌细胞的“演化路径”与耐药性萌芽

另一个让我印象深刻的TDA应用,是在研究肿瘤异质性和耐药演化方面。我们拿到了同一批肺癌患者治疗前和治疗后(耐药复发)的多次活检样本的单细胞RNA测序数据。问题很明确:耐药细胞是从哪里冒出来的?它们和原来的敏感细胞在状态上是渐变的还是突变的?

我们把所有样本的细胞数据合并,用TDA进行分析。这次,我们用了“伪时间”分析工具(如Monocle)推断的细胞分化轨迹作为Mapper的透镜。生成的拓扑图不再是一个个孤立的团,而是一条条有分支、有环路的“河流网络”。我们能够清晰地看到:

  • 治疗前的样本主要聚集在网络的几个上游“源头”节点。
  • 治疗后耐药的样本,则出现在下游特定的几个分支末端。
  • 最关键的是,我们发现了连接“源头”和“耐药末端”的路径,并且在这条路径的中间节点上,已经存在少量治疗前的细胞!这些中间态细胞,虽然对药物还敏感,但其基因表达谱已经偏向于耐药状态。

这个发现的意义在于,TDA揭示了耐药性并非完全由新突变产生,而是在治疗前就有一小部分细胞处于“预适应”状态。治疗压力选择了这些细胞,使其扩增。这提示我们,未来的治疗策略可能需要联合用药,在治疗初期就同时靶向这些潜在的耐药路径,而不是等耐药克隆壮大后再想办法。

### 3.3 实操建议:如何为你的癌症数据选择TDA参数?

看到这里你可能摩拳擦掌想试试了。我分享几个从坑里爬出来的经验,关于参数调优:

  • 透镜的选择比想象中更重要。不要只默认用PCA。多从生物学假设出发:如果你想看免疫浸润的影响,就用免疫相关基因签名分数;想看代谢重编程,就用代谢通路活性评分。有监督的透镜(如用生存风险评分)往往能直接得到有临床意义的拓扑图。
  • resolutionoverlap 需要联动调整。我的经验是,先从较小的分辨率(如5-8个区间)和较大的重叠度(0.4-0.5)开始,得到一个概览图。如果图太简单,看不清结构,就提高分辨率;如果图太破碎,像一盘散沙,就降低分辨率或适当增大重叠度。目标是得到一个既能看清主干结构,又有一些有趣分支的图。
  • 聚类算法推荐DBSCAN,但eps参数要小心。在高维基因空间,距离容易失效。一个技巧是,先在你选定的透镜区间内的数据子集上,快速计算一个距离矩阵的百分位数(比如第5%的分位数),作为eps的初始估计。
  • 一定要做稳定性评估。TDA图对参数敏感。你可以用Bootstrap重采样,或者轻微扰动参数,看看图的主要结构(比如关键的小节点、重要的连接桥)是否稳定出现。不稳定的特征可能需要谨慎解读。

4. 超越可视化:将拓扑特征融入机器学习管道

很多人把TDA当作一个高级的可视化工具,看完图,讲个故事就结束了。这其实大大浪费了它的价值。TDA产生的拓扑图本身,以及从图中提取的定量特征,是可以直接作为特征输入到下游的机器学习模型中的,能显著提升模型的预测性能。这才是TDA真正发挥威力的地方。

### 4.1 从拓扑图中提取什么特征?

一张Mapper图,我们可以从中提取两大类特征:

  1. 图论特征:直接把拓扑图看作一个网络(图),计算每个节点或整个图的特征。
    • 节点级:节点的度(连接数)、中心性(如介数中心性,衡量节点作为“桥梁”的重要性)、聚类系数等。一个具有高度中心性的节点,可能代表了数据中一个关键的“枢纽”状态。
    • 图级:图的连通分量数量、平均路径长度、直径、聚集系数等。这些宏观指标描述了数据整体的复杂性和连通性。
  2. 拓扑不变量特征:这是TDA更核心的产出,主要通过持续同调计算得到。它量化了数据在不同尺度上的“洞”的特征。
    • 0维持续同调:描述“连接组件”(簇)的出生与死亡。可以理解为在不同距离阈值下,簇的生成与合并。
    • 1维持续同调:描述“环”或“空洞”的出生与死亡。在生物网络中,这可能对应着反馈环路或功能模块。
    • 2维持续同调:描述“空腔”的出生与死亡。

每个“洞”都有一个“出生”尺度和“死亡”尺度。出生和死亡时间之差,称为“持续度”。持续度长的洞,代表数据中一个稳健的拓扑结构。我们可以把这些“持续对”的集合,转换成一种叫持续图持续景观的向量表示,作为机器学习模型的特征。

### 4.2 实战管道:用拓扑特征预测癌症患者生存期

下面我构建一个简单的示例管道,展示如何将TDA特征用于生存预测:

import numpy as np
import pandas as pd
from gtda.mapper import make_mapper_pipeline, plot_static_mapper_graph
from gtda.homology import VietorisRipsPersistence
from gtda.diagrams import PersistenceLandscape, Amplitude
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomSurvivalForest
from sksurv.metrics import concordance_index_censored

# 1. 准备数据:X是基因表达矩阵,y是生存事件和时间的结构化数组
X = ... # (n_samples, n_genes)
y_structured = ... # 例如,dtype=[('event', bool), ('time', float)]

# 2. 生成Mapper图(用于探索和可视化)
mapper_pipeline = make_mapper_pipeline(
    filter_func=PCA(n_components=2), # 使用PCA前两维作为透镜
    cover=CubicalCover(n_intervals=10, overlap_frac=0.3),
    clusterer=DBSCAN(eps=0.5, min_samples=5),
)
mapper_graph = mapper_pipeline.fit_transform(X)
# 可以在这里可视化 mapper_graph

# 3. 提取拓扑特征(用于机器学习)
# 3.1 计算持续同调(这里用更通用的VietorisRips复形,适合点云数据)
homology_dimensions = [0, 1, 2]  # 计算0,1,2维的同调
persistence = VietorisRipsPersistence(
    metric='euclidean',
    homology_dimensions=homology_dimensions,
    collapse_edges=True
)
# 注意:VietorisRips直接对原始数据X操作,计算所有点之间的拓扑
diagrams = persistence.fit_transform(X)  # 输出形状 (n_samples, n_features, 3)

# 3.2 将持续图转换为向量特征(例如,使用持续景观)
landscape = PersistenceLandscape(n_layers=5, n_bins=100)
X_topological = landscape.fit_transform(diagrams)  # 转换为向量
# 或者计算拓扑振幅(多种标量摘要)
amplitude = Amplitude(metric='wasserstein')
X_amp = amplitude.fit_transform(diagrams)

# 4. 将拓扑特征与原始基因特征(或临床特征)融合
# 假设我们选择持续景观特征
X_combined = np.hstack([X, X_topological.reshape(X_topological.shape[0], -1)])

# 5. 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X_combined, y_structured, test_size=0.2, random_state=42)

# 6. 训练生存预测模型
rsf = RandomSurvivalForest(n_estimators=100, random_state=42)
rsf.fit(X_train, y_train)

# 7. 评估
prediction = rsf.predict(X_test)
c_index = concordance_index_censored(y_test['event'], y_test['time'], prediction)[0]
print(f"模型C-index为: {c_index:.3f}")

在我做过的多个项目中,加入拓扑特征后,生存预测模型的C-index(一致性指数)提升0.02到0.05是常有的事。别小看这零点零几,在竞争激烈的预后模型研究中,这往往是决定性的提升。因为这些拓扑特征捕捉到了样本间关系的结构,而不仅仅是单个基因的表达水平,它们提供了互补的信息维度。

5. 避坑指南与未来展望:让TDA真正为你所用

TDA很强大,但也不是银弹。最后这部分,我想结合自己踩过的坑,给想上手的朋友一些实在的建议,并聊聊它未来的可能方向。

### 5.1 新手常踩的坑与应对策略

  1. 坑一:盲目追求复杂图形,过度解读。TDA图可能很漂亮,结构复杂,但并非每一个小节点和连接都有生物学意义。一定要与先验知识结合。看到一个有趣的结构,立刻回去看这个节点里样本的临床信息、突变谱、通路活性。如果找不到合理的生物学解释,它可能只是数据噪声或参数巧合。
  2. 坑二:参数敏感性带来的不稳定性。这是Mapper算法被诟病较多的一点。解决方案就是前面提到的稳定性评估。同时,可以尝试一些改进的Mapper算法或软件包,它们提供了更稳健的覆盖构建方式。在实践中,我通常会用几组不同的参数(透镜、分辨率、重叠度)分别运行,只报告那些在多种参数下都稳定出现的拓扑特征。
  3. 坑三:计算复杂度。对于单细胞测序这种动辄数万个细胞的数据,直接计算持续同调或跑Mapper可能会很慢。这时候需要一些策略:先降维,使用PCA或自动编码器将维度降到50-100左右,再进行TDA分析;对细胞进行下采样,或者先进行粗聚类(如Louvain聚类),用聚类中心代表点进行分析;利用GPU加速的TDA库(如一些基于CUDA的实现)。
  4. 坑四:结果的可解释性与沟通。向生物学家或临床医生解释一张拓扑网络图比解释一个散点图要难得多。我的经验是:讲故事,而不是讲图。不要一上来就展示复杂的网络,而是从生物学问题出发——“我们想看看耐药性是如何演化的”,然后展示TDA如何揭示了“一条从敏感状态通向耐药状态的路径”,并在这条路径上标注出关键的基因变化。用动画展示随着透镜参数变化图的演变,也是一个很好的沟通方式。

### 5.2 TDA在生物医学领域的未来可能

从我个人的观察来看,TDA在以下几个方向的潜力巨大:

  • 与深度学习融合:这是最令人兴奋的方向。比如,用图神经网络直接对Mapper生成的拓扑图进行学习;或者设计神经网络的层,使其能够直接计算和优化拓扑损失函数,让模型在学习过程中主动保持或发现特定的拓扑结构。
  • 多组学数据整合:癌症是基因组、转录组、表观组、蛋白组等多层信息共同作用的结果。TDA非常适合整合这种多层次、异质性的数据。我们可以为每一组学数据构建一个透镜,然后用多参数透镜的Mapper,或者分别构建拓扑图再进行图的对齐与融合,来发现跨组学的一致拓扑结构。
  • 动态与时空数据分析:比如分析肿瘤在治疗过程中随时间变化的单细胞数据,或者空间转录组数据。TDA可以扩展到“持续同调的时间序列”或“空间过滤”,捕捉细胞状态演化的拓扑动力学,或者肿瘤微环境中细胞类型的空间组织模式,这比静态分析能揭示更多信息。

说到底,TDA是一种思想,一种从“形状”和“关系”的角度理解复杂数据的思想。它不替代你熟悉的统计检验或机器学习模型,而是为你提供了一个全新的、强大的探索性工具和特征生成器。刚开始用可能会觉得有点绕,参数调起来也烦,但一旦你通过它发现了那个被传统方法遗漏的关键亚型,或者那条隐藏的演化路径,那种“原来如此”的顿悟感,会让你觉得一切折腾都是值得的。我的建议是,找一个小而干净的数据集,从复现一个经典的TDA案例开始,亲手调调参数,看看图是怎么随着你的操作变化的,这种感觉就慢慢来了。生物数据如此复杂,我们需要的正是像TDA这样,能够尊重并揭示其内在复杂性的工具。

Logo

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐