TCGA+单细胞+空间转录组：生信高分论文全流程拆解

这套“TCGA数据 + 单细胞手动注释 + 双虚拟（敲除/过表达）+ 分子亚型 + 细胞轨迹 + 空间转录组 + 模型综合分析”的组合，是目前生信论文里非常典型的高分套路。它解决的实际问题只有一个：从公共数据出发，把一个基因、一群细胞或一条通路的故事讲完整。比如你想看 STMN1 在宫颈癌里的表达和细胞来源，常规 TCGA 只能告诉你表达量高低，补上单细胞注释就能知道它主要在哪种细胞里表达，再补上轨迹和空间信息，就能解释这种表达差异是怎么形成、在哪形成。

这篇文章面向三类人：准备做生信数据挖掘的医学研究生、想复现高分文章分析逻辑的科研人员，以及要把公共数据做成预测模型的算法工程师。最值得关注的点是：这套流程真正难的不是某个软件跑不跑得通，而是每一步都在做取舍——注释准不准、分组合不合理、模型有没有过拟合。下面按我实际落地这套分析时的顺序拆开讲。

1. 这类综合分析最值钱的不在“跑通”，而在每一步的取舍

1.1 先确定这套流程解决什么科研问题

标题里的每个词分别对应一个分析目标：

• TCGA 数据：提供大样本的转录组、突变、临床随访信息，用来做表达差异、预后关联和模型训练。
• 单细胞（手动注释）：解决肿瘤异质性问题，明确目标基因在哪些细胞类型里表达，而不是只看组织块里的平均值。
• 双虚拟（敲除/过表达）：不依赖实验，直接在计算层面模拟目标基因下调或上调后，细胞状态和通路会发生什么变化。
• 亚型：把样本或细胞划分成不同的功能亚群，用来解释为什么同样一个基因在不同病人里的影响不一样。
• 细胞轨迹：描述细胞从一种状态到另一种状态的连续变化，回答“谁从哪来、要往哪去”。
• 空间转录组：把基因表达放回组织位置里去看，回答“高表达区域旁边是什么细胞，有没有空间上的相互作用”。
• 模型综合分析：把前面的特征收敛成一个可判断的东西，通常是风险评分、分型标签或预测概率，方便临床解读。

这七个部分不是堆砌，而是一条证据链。TCGA 给出相关性，单细胞给出细胞层级，双虚拟给出功能方向，亚型和轨迹给出状态变化，空间数据给出位置证据，模型把相关性转成可用的输出。

1.2 分析链路里的一个关键判断点

整条链路里最容易翻车的位置，不是最后建模，而是中间的数据衔接。TCGA 是 bulk 转录组，单细胞是单细胞分辨率，空间转录组是带坐标的矩阵。这三类数据虽然都叫表达矩阵，但维度、稀疏度、分布都不一样。

不少人直接拿 Seurat 的标准化结果去做 limma 差异分析，或者拿 monocle 的拟时序结果去和 TCGA 的生存数据合并，最后发现分组对不上、细胞名混乱。原因是这三类数据有不同的对象粒度：TCGA 的样本是人，单细胞的样本是细胞，空间数据的样本是 spot 或 bin。分析之前要把对象粒度统一，不能混着算。

我一般会在项目目录里单独建一个 metadata 文件夹，里面放三张表：样本信息表、细胞信息表、spot 信息表。每次合并结果前先检查 ID 是否有重复、是否有缺失，再往下走。这一步看起来简单，但能省掉后面大量报错排查时间。

2. 数据准备：先把 TCGA 和单细胞输入整理干净

2.1 TCGA 数据怎么下载、怎么转成矩阵

TCGA 数据现在推荐从 GDC 官网或官方 API 下载，也可以用第三方 R 包，比如 TCGAbiolinks 或 cBioPortal 的导出功能。第一次接触的人最容易卡在两个地方：

一是下载下来的文件不是直接可用的矩阵。GDC 默认下载的是每个样本一个表达文件，需要合并成矩阵。使用 TCGAbiolinks::GDCquery 和 GDCdownload 下载后，可以用 GDCprepare 读取，得到的是 SummarizedExperiment 对象，再通过 assay() 提取表达矩阵。

二是样本名和样本类型对不上。TCGA 的样本 barcode 第 14、15 位代表组织来源码，“01”是肿瘤，“11”是癌旁。做差异分析前要按 barcode 拆分，不能把所有样本混成一个组。

一个通用处理流程：

r复制# 示例：R 环境下整理 TCGA 表达矩阵
library(TCGAbiolinks)
query <- GDCquery(
  project = "TCGA-CESC",          # 以宫颈癌为例
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "STAR - Counts"
)
GDCdownload(query)
GDCprepare(query)

下载后建议做四件事：过滤低表达基因、把 Ensembl ID 换成 Symbol、保留蛋白编码基因、检查缺失值。常见环境中基因数会在 2 万到 3 万之间，过滤后通常剩 1.5 万到 2 万。

2.2 单细胞数据格式与 QC 前检查

单细胞数据一般下载的是 10x Genomics 的格式：barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz。如果拿到的是 h5 文件，也可以用 Read10X_h5() 读取。

QC 前先看三样东西：

• cellranger 输出的 metrics_summary.csv，里面有测序饱和度、细胞数、中位基因数。
• 是否需要过滤双细胞。常见做法是用 DoubletFinder 或 scrublet。
• 样本来源是否混杂了多个病人。如果有多个样本，要保留 orig.ident 信息，后面合并时要用。

数据格式检查通常这样开始：

r复制library(Seurat)
counts <- Read10X(data.dir = "path/to/filtered_feature_bc_matrix")
obj <- CreateSeuratObject(counts = counts, project = "cervical", min.cells = 3, min.features = 200)

2.3 环境准备：硬件、软件和依赖

这套分析不是单机一键完成。TCGA 的大矩阵在普通笔记本上能跑，但单细胞数据如果超过 10 万细胞，Seurat 的默认流程会明显变慢，更不用说空间转录组和轨迹分析。

一个比较稳妥的配置参考：

没有 GPU 也能做，只是部分步骤慢。Seurat 的 PCA/UMAP 直接用 CPU 跑，10 万细胞的 UMAP 可能要十分钟以上。空间转录组如果用深度学习方法做反卷积或配准，GPU 才有明显优势。

3. 单细胞质控与手动注释，决定后面所有结论的可信度

3.1 质控参数不是死值，得看你自己的数据

单细胞分析最常见的质控指标是三个：每个细胞的基因数（nFeature）、UMI 总数（nCount）和线粒体基因比例（percent.mt）。很多教程直接写“过滤 nFeature < 200 和 percent.mt > 20%”，我一般不建议直接照抄。因为不同组织类型里这三个指标的合理范围差异很大。

比如肿瘤样本里的上皮细胞往往 UMI 很高，而 T 细胞、中性粒细胞的基因数天然偏低。如果你直接按一个固定阈值过滤，可能把真实的免疫细胞全部丢掉。

更稳的做法是分两步：

第一步，先过滤明显异常值：基因数为 0、线粒体比例超过 50% 的细胞，这类通常是破损细胞或空液滴。
第二步，用 VlnPlot 看每个样本的分布，再根据分布折点确定具体阈值。

质控代码通常长这样：

r复制obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 20)

3.2 聚类之后先读 marker，再动手注释

质控完就是标准化、PCA、UMAP、聚类。这些步骤在 Seurat 里是标准流程，跑起来很快，真正的重点是注释。自动注释可以先用 SingleR 或 CellTypist 给一个粗略结果，但不能直接拿这个当最终答案，尤其是肿瘤样本，免疫细胞、基质细胞、恶性上皮细胞经常被误判。

手动注释的基本逻辑是：先看每个 cluster 的高表达 marker，再结合已知生物学背景判断它是什么细胞。以宫颈癌为例：

• 上皮细胞：KRT5、KRT17、KRT8、EPCAM
• T 细胞：CD3D、CD3E、CD2
• B 细胞：MS4A1、CD79A
• 巨噬细胞：LYZ、CD68、AIF1
• 成纤维细胞：COL1A1、DCN、LUM
• 内皮细胞：PECAM1、VWF、CDH5

单看一个 marker 不够，要同时看几个。比如 CD3D 阳性且 CD14 阴性的群更可能是 T 细胞；如果一群细胞同时高表达 KRT 和 CD68，那更可能是双细胞或者技术污染，要回到原始矩阵检查。

3.3 手动注释的质量判断标准

手动注释做到什么程度才算合格？我总结三个可检查的标准：

• 每个主要细胞类型至少有 3 个已知 marker 支持，不是只看一两个基因。
• 注释后的细胞类型分布符合组织逻辑，比如宫颈癌样本里上皮细胞和成纤维细胞通常占比最高，内皮和免疫细胞比例因人而异。
• 注释结果在 UMAP 上有可解释性，同一类的细胞应该相对集中，而不是散成一团。

如果某个 cluster 注释不了，不要硬起名字。可以暂时标记为“Unknown”或“Low quality”，在后续分析里排除，或者补充更权威的参考数据重新映射。

做完手动注释后，常见的一个动作是“看 stmn1 在宫颈癌的表达情况”，这时你已经可以回答得很具体：STMN1 在宫颈癌组织里主要表达于哪类细胞，在恶性上皮细胞里的比例是多少，而在成纤维细胞里几乎不表达。这个结论只有做单细胞注释之后才能得出。

4. 双虚拟分析（敲除/过表达）怎么落地，结果如何判断

4.1 “双虚拟”到底是什么

“双虚拟”在生信分析里指的是不进行真实实验，而是在计算环境中模拟基因功能的变化。最常见的是两种：

• 虚拟敲除：模拟目标基因表达量降低或完全缺失后，细胞内的转录调控、通路活性发生什么变化。
• 虚拟过表达：模拟目标基因表达量升高后，哪些下游通路被激活。

因为用的是真实测序数据，这种模拟并不是凭空捏造，而是通过比较目标基因高表达组和低表达组的系统差异，间接推断功能。它不是证明因果，而是给出强关联和可验证的假设。

4.2 常用实现路径：表达分层、回归调控、工具预测

实现虚拟敲除和过表达，我常用三种方式，各有适用场景。

第一种是表达分层，最简单。把细胞或样本按照目标基因表达量排序，取前 30% 作为高表达组模拟过表达，后 30% 作为低表达组模拟敲除，然后做差异表达、GO/KEGG 富集、GSEA。这个方案不需要额外工具，结果是可解释的，但存在过拟合风险，分组阈值不同可能改变结论。

第二种是用调控网络/转录因子活性推断。通过 SCENIC、pySCENIC 或 PROGENy 推断目标基因下游的转录因子和通路，然后比较敲除和过表达后的调控子活性差异。这种方式更接近生物学机制，但工具安装和运行时间都比较长，Python 环境里 pySCENIC 的依赖相对复杂。

第三种是通过已有的数据库或工具做扰动预测，比如基于基因调控网络做 in silico perturbation。这类工具对数据质量要求高，不太适合初学者直接作为主分析，可以当作补充验证。

4.3 敲除与过表达结果的差异怎么解释

如果 STMN1 虚拟敲除后，细胞周期相关通路显著下调，而上皮间质转化通路在过表达组里显著上调，那么这个故事就非常顺：STMN1 可能与细胞增殖和转移表型有关。但要注意，这种关联可能来自批次效应或细胞组成差异，尤其当目标基因在不同细胞类型里表达差异很大时。

我一般会加一步去混淆分析。把免疫细胞和上皮细胞分开做双虚拟分析，或者用线性回归把细胞类型比例作为协变量校正。否则得到的结果可能是“细胞类型差异”，而不是“基因功能效应”。

4.4 验证：in silico 结果不能替代实验，但能筛选

双虚拟分析的定位是假设筛选，不是实验替代。普通 TCGA 数据只能说明 STMN1 高表达和较差预后相关，虚拟敲除可以进一步指出可能的机制通路，但最终结论需要湿实验验证。

在实际复现高分文章的时，这里通常会画一个机制图，图里一半是计算分析得到的通路变化，一半是文献报道的已知功能。我做的时候会先在本地跑通通路富集，再对比 KEGG、Hallmark 里收集到的参考条目，看看有没有互相印证。如果计算结果显示一个完全没有文献支撑的奇怪通路，先检查输入数据和分组代码，不要急着当新发现。

5. 分子亚型与细胞轨迹：只看聚类不够，还要看关系

5.1 亚型聚类怎么和前面的单细胞注释衔接

分子亚型通常是把 TCGA 样本按表达谱聚类，得到几个亚群。常见方法有 ConsensusClusterPlus、NMF、k-means。这一步容易犯的错是：只拿全部基因做聚类，没有过滤高度变异的基因或没有去除批次效应。

更好的做法是结合单细胞结果来选基因。比如在单细胞注释里发现目标基因主要表达于某个上皮亚群，就用这个亚群的 marker 基因去对 TCGA 样本做分型，这样亚型结果和细胞类型的对应关系更清晰，后面的机制解释更有依据。

分型完成后要做三件事：

• 确定最优簇数：看累积 CDF 曲线下降趋势和 consensus 矩阵是否清晰。
• 检查亚型间目标基因表达是否有差异。
• 比较亚型间的生存曲线。

5.2 细胞轨迹分析选 Monocle 还是 CytoTRACE

细胞轨迹分析的作用是把离散的细胞群排序成一个连续过程。做之前先想清楚：你要描述的是哪一种连续变化？是上皮细胞恶性转化过程、T 细胞耗竭过程，还是成纤维细胞向肌成纤维细胞转变？目标不同，选的方法就不同。

主要选择可以这样判断：

如果只是想看目标基因在细胞状态变化中的趋势，CytoTRACE 是最快的，几分钟就能出结果。如果想还原一个分支事件，比如某群细胞分化成两个方向，Monocle 3 更合适。Monocle 2 很多人习惯用，但遇到大数据时很容易内存爆掉，处理 10 万细胞要特别小心。

5.3 亚型、轨迹和预后的联动判断

轨迹分析的输出是一根轴，代表细胞状态的变化，每个细胞或 spot 在这根轴上有一个位置。这个位置和 TCGA 的预后信息怎么连接？常见做法是计算每个样本“轨迹高分细胞比例”或“分化末端 signature 分数”，然后和生存数据合并。

这种联动是整篇分析里最需要解释清楚的地方。如果 STMN1 在轨迹末端表达升高，同时高表达 STMN1 的亚型预后更差，那就可以说明：STMN1 高表达更可能来自终末状态的恶性细胞，而不是早期分化细胞。这个逻辑链在审稿人眼里比较完整。

但注意：轨迹分析是计算推断，不是命运图谱实验。论文里写作要谨慎，不能直接写“证明细胞分化到末端后才开始表达 STMN1”，可以写“提示可能存在这种过渡关系”。

6. 空间转录组：能加就加，但要先想清楚补什么信息

6.1 空间转录组补的是“位置”信息

空间转录组和单细胞相比，最大的增量是位置。同样是 STMN1 高表达，在单细胞里只知道是上皮细胞表达，在空间里能知道这个高表达区域紧挨着成纤维细胞区域，旁边还有一圈免疫细胞，那这就引出了一个更具体的假设：高表达 STMN1 的恶性细胞可能和微环境里的成纤维细胞有空间互作，从而影响免疫浸润。

这个假设可以继续用 CellChat 做配受体分析验证。但空间数据本身噪声大、每个 spot 可能是多个细胞的混合物，表达矩阵稀疏程度比单细胞还高。

6.2 空间转录组与单细胞整合

空间转录组分析最常用的方案还是 Seurat 的 Spatial 流程。大致的整合路径是：

1. 读取空间数据，检查 spot 数、基因数、线粒体比例。
2. 用标准流程做 PCA 和 UMAP。
3. 找空间高变基因（SVG），常用的工具是 FindSpatiallyVariableFeatures 或 Moran's I。
4. 将单细胞数据作为参考，用 Deconvo 方法对空间 spot 反卷积，得到每个 spot 的细胞类型比例。
5. 可视化目标基因的表达位置，并做邻近分析。

反卷积是我觉得最值得花时间调的部分。同一个空间数据集，用不同反卷积工具得到的结果差异可能很大。我常用的策略是同时跑 Cell2location 和 SPOTlight，如果两个结果一致，结论才敢放出来；如果差异大，优先检查参考单细胞的注释质量。

6.3 什么时候可以不做空间整合

不是每个课题都需要空间转录组。如果你的公共数据和目标基因已经足够支撑核心结论，加上空间数据反而会被审稿人追问“空间结果和前面的单细胞结论有什么增量”。这种情况下，不加更安全。

判断标准是这样：如果单细胞轨迹结果显示目标基因在某个亚群中高表达，而空间数据能显示这个亚群在组织里的位置，并且和免疫排斥或基质重塑相关，那就有明显增量。如果空间结果只显示“高表达区域集中在一个区域”，没有和细胞互作、配受体、邻近关系等分析配套，那就是一张漂亮的图，但证据级别有限。

7. 模型综合分析：把前面所有信息收敛到临床问题

7.1 特征怎么选：从差异表达到高维特征

前面几步得到的信息很多：TCGA 差异基因、单细胞 marker、双虚拟分析中的通路、亚型特征、轨迹 signature。直接把这些全部扔进模型，维度爆炸，还会导致严重过拟合。

我一般先把特征分成三类：

• 基因层面的特征，比如目标基因及共表达模块。
• 通路层面的特征，比如 GSVA/Hallmark 通路分数。
• 细胞层面的特征，比如免疫细胞浸润比例、某种细胞亚群的比例。

建模时优先用通路分数和细胞比例作为特征。原因很简单：基因层面的噪声很大，反复横跳的基因太多；通路分数更稳定，细胞比例更接近生物学状态。

7.2 建模：LASSO、随机森林、Cox 怎么搭

如果是生存数据，经典组合是：先用单因素 Cox 筛选候选特征，再用 LASSO-Cox 压缩变量，最后用多因素 Cox 建立风险评分。如果是分类任务，比如区分两个亚型，可以用随机森林加逻辑回归。

一个比较常见的流程：

r复制# 示例：LASSO-Cox 变量筛选（R）
library(glmnet)
x <- as.matrix(train_expr)
y <- survival::Surv(time, status)
fit <- cv.glmnet(x, y, family = "cox", alpha = 1)
coef_min <- coef(fit, s = "lambda.min")

参数选择上，alpha = 1 对应 LASSO，适合大量变量中挑少数关键变量；alpha = 0 是 Ridge；alpha 在 0 到 1 之间是 Elastic Net。LASSO 更受欢迎，因为它会直接把很多系数压成 0，给一个清晰的风险基因列表。

用了 LASSO 并不代表万事大吉。判断标准要看三个数：训练集 C-index 或 AUC、验证集的对应指标、校准曲线。我在实际审查模型时会优先怀疑以下情况：训练集 AUC 0.95，验证集只有 0.6，这个模型大概率是过度拟合了，特征数量也要检查。

7.3 验证：包括外部验证和风险分层

模型的验证分为内部和外部。内部验证包括时间依赖 AUC、C-index、校准曲线，常见用 timeROC 或 survcomp。外部验证最好找 GEO 或 ICGC 数据库里的另一个独立队列，如果不能做外部验证，至少要留出一部分数据做验证集。

除了数值化指标，还要把模型转成可以解释的临床分层。患者可以按风险评分中位数或最佳截断值分成高低风险组，展示两条生存曲线的分离程度。这个图是审稿人最关心的部分之一。

7.4 模型“综合分析”的常见误区

最大误区是把所有特征堆到一起，得到 AUC 很高就认为模型好。我见过很多分析，把 TCGA 训练集里揉进去了几十个免疫浸润比例、几十条通路分数，模型指标好看，但实际临床中这些变量很难收集。一个好的模型应该满足三个条件：输入变量可测、系数相对稳定、在外部数据中不会明显衰退。

第二个误区是不区分“相关性”和“预测力”。目标基因 STMN1 和相关通路可以做机制故事，但它们加入模型以后不一定提高预测能力。建模前可以做一次变量重要度排序，不想放进最终模型的变量，可以作为讨论内容。

第三个误区是使用全数据集反复调参。这样会导致隐式过拟合，等到拿测试集验证时，结果已经不等于模型真实表现。正确做法是先用训练集完成特征筛选、参数调整、模型选择，测试集只跑一次。

8. 全流程实操清单与排查顺序

8.1 实操顺序：先跑通再调参

如果从零开始复现这套分析，我的建议是不要按文章模块的完成顺序走。先把每个模块的最小流程跑通，再统一调到高质量图表的参数。

具体顺序是：

1. TCGA 差异分析，确认目标基因在肿瘤和癌旁之间有差异。
2. 单细胞质控和注释，得到细胞类型。
3. 目标基因在单细胞里的表达分布。
4. 双虚拟分析，得到通路和富集结果。
5. 分子亚型和轨迹，补充状态变化和分群信息。
6. 空间数据，作为位置证据。
7. 建模，把所有表达结果收敛成模型。

先跑通指的是：单细胞很粗糙的过滤、注释只到主要细胞大类、轨迹不考虑分支。这时候的目标是管线不报错、输出文件齐全。跑通后再回头调优才有意义。

8.2 各类报错的排查优先级

这套分析链条长，报错种类也多。我一般按这个顺序排查：

1. 先看路径问题。R 或 Python 报“cannot open file”或“No such file”，大多数时候是路径带中文、相对路径工作目录不对、文件扩展名不对。
2. 再看格式问题。表达矩阵的行列有没有反、基因 ID 是 Symbol 还是 Ensembl ID、是否有重复基因名。
3. 再看依赖版本。Seurat 升级到 5.0 以后，很多旧代码的 pData、idents 调用方式都有变化；Monocle 2 和 Monocle 3 的输入要求也不一样。
4. 最后看资源占用。内存不足在空间转录组和轨迹分析里非常常见，尤其是 10 万细胞以上用 Monocle 2，很容易卡死。

8.3 如果想复现这类 Q1 区 9.7 分文章，核心要盯住什么

“9.7 分”不是一个可以由分析流程直接换来的结果，它背后通常是完整的证据链设计加实验验证。但这套计算框架确实值得复现，因为它把几乎所有主流生信分析模块串成了闭环。

复现时最该盯住的不是某个包用得多熟练，而是下面四点：

• 每个模块的输出是否能在下一个模块里被使用。TCGA 的特征基因能不能映射到单细胞注释里，单细胞亚群能不能映射到空间反卷积里。
• 每个分组变量是否有明确的定义。细胞类型、风险高低、亚型类别，这些定义不一致，整个故事就会塌。
• 每个结论是否有至少两种方法或两个数据集交叉验证。
• 图表是否直接服务于那条核心故事线，而不是为了展示分析量大而堆图。

这套分析的完整落地时间，根据不同电脑配置和熟练程度，从一周到一个月都可能。第一次做，不需要追求每张图都达到发表水平，先把通路的输入输出和判断标准搞清楚。后面再优化时，主要是调注释、调参数、调整模型特征，而不是推翻重跑。