1. 思路拆解：为什么单细胞组间比较，光看差异表达已经不够用了
先聊一个痛点。做单细胞测序组间比较的朋友应该都有体会，常规流程跑完Seurat聚类、找marker、做差异表达，拿到一张火山图、一张热图，审稿人会觉得“你做完了”，但要发2区以上，很多时候这已经不够了。差异表达说到底只是“哪些基因变了”，可我们真正想回答的问题是“哪些调控关系变了”——同一个基因，表达量没变，但它所在调控网络里的上下游关系可能已经彻底重排了。这种信息，单纯靠差异表达是看不出来的。
所以我在“14天拿下2区”这轮分析里，选择了两套虚拟手段搭配组合：scTenifoldKnk做虚拟敲除与调控网络推断，Geneformer做基于Transformer架构的单细胞基础模型推理。这两个工具的底层逻辑完全不同，一个从共表达网络和拓扑学切入，一个从大规模预训练模型迁移切入，但它们有一个共同的产出——不需要做湿实验，就能从转录组数据里“虚拟”出基因调控层面的扰动信号。联合单细胞数据特征，再在这个基础上嫁接机器学习做分类和特征筛选，整套流程下来，分析深度和故事完整度确实比传统pipeline高出一截。
说句实话，这套组合第一次跑通的时候，我一度觉得经验丰富的审稿人会质疑“虚拟”二字的可靠性。但实际操作下来，只要数据清洗干净、参数没有瞎调，scTenifoldKnk找到的网络核心节点和Geneformer给出的注意力重要基因，往往能在下游实验验证中踩中已知靶点。这个交叉验证的结果，就是2区文章里最有说服力的figures。
这篇文章我想完全按我自己的实操路径来写，从环境搭建、数据格式、参数设置，到每一段代码思路、每一步踩坑，尽量全给出来。目标读者是已经会跑Seurat基础流程，但想往更深的调控网络和机器学习方向再走一步的研究生和生信工程师。如果你是纯新手，建议先把单细胞标准流程走一遍再来看这篇文章。
2. scTenifoldKnk到底在做什么，它的“虚拟敲除”为什么能替代一部分湿实验
2.1 虚拟敲除的原理：不是真的敲基因，而是看网络断了一条边之后系统怎么响应
scTenifoldKnk这个R包的思路，我理解下来可以概括成一句话：它把单细胞表达矩阵转成基因调控网络，然后“模拟”剔除某一个基因（也就是虚拟敲除），观察整个网络里其他基因之间的关联强度会发生什么变化。变化最大的那些基因，就是对这个被敲除基因最敏感的“下游响应者”。
这里有一个容易被忽略的底层细节：scTenifoldKnk不是直接把表达矩阵拿来做相关性分析，它内部有一个多步流程。大概路径是这样的——先通过PCA降维和PEER因子修正把技术噪声和批次效应剥离掉，再用一种基于偏相关的网络推断方法构建基因间的调控连接，最后通过流形对齐和网络去卷积得到每个基因在网络里的“角色”和通路归属。它的输出是一张加权网络表，每一行代表一条基因与基因之间的边，权重越大说明调控关系越强。
我刚开始用这个包的时候犯过一个错误，就是用默认参数直接跑全基因组，结果跑了两天都没结束。后来翻了一下文档和源码才发现，它的计算复杂度和基因数量之间不是线性关系，而是接近立方级别。所以实际项目中，我们通常只用差异表达显著的基因合并关键通路基因，比如1500到3000个基因的规模，这样既能保证网络完整性，又能让计算在可接受范围内完成。
提示：scTenifoldKnk的输入格式是稀疏矩阵（dgCMatrix），不是数据框。直接把Seurat的
GetAssayData结果丢进去之前，一定要确认行列顺序，列是细胞，行是基因，方向反了后面所有结果都不可信。
2.2 组间比较的网络差异分析逻辑：不是比均值，而是比拓扑结构
scTenifoldKnk最常见的应用方式，是对照组和处理组分别构建一个基因调控网络，然后比较两个网络的差异。这个比较不是简单的“这个基因表达量高了多少”，而是通过虚拟敲除每一个基因，比较两组网络中敲除后全网络响应模式的差异，最后输出一个差异调控基因列表。这个列表里的基因，就是两组之间调控关系变化最大的节点。
我在实际跑的时候，会同时输出两个网络文件，然后专门写一小段R代码提取差异边。举个具体例子：对照组网络里，基因A和基因B的连接权值是0.85，处理组网络里这个权值掉到0.12，这表明在疾病或者处理条件下，A和B之间的共调控关系被大幅削弱了。这种信息在常规的差异表达分析里完全看不到，因为A和B本身的表达量可能一点没变。
还有一个心得是，比较两个网络时，一定要保证两组细胞数量不要差太多。细胞数差距过大会导致网络推断的置信度不对等，最后出来的差异基因可能有很大一部分是假阳性。我一般控制在3：5到5：3的范围内，如果实在不平衡，就用downsampling把多的一组抽到与少的一组相近的细胞数，这个操作在后面的代码部分会具体演示。
3. Geneformer是怎么补位的——预训练模型带来的迁移学习优势
3.1 Geneformer的本质与输入格式：rank值编码的Token序列
Geneformer是今年单细胞领域讨论度很高的一个模型，它在数千万个单细胞转录组数据上做了自监督预训练，模型学到的不是某个特定疾病的知识，而是基因之间通用的共表达与调控语法。当我们把一个新的单细胞数据丢进去的时候，它不需要从头学习，只需要微调或者直接推理，就能输出每个基因在细胞状态中的“角色重要性”。
Geneformer的输入格式比较特别，它对每个细胞里的基因按照表达量排序，转成一个由基因rank组成的token序列。也就是说，它不是直接输入基因表达量数值，而是输入“这个基因在细胞里排第几”。这个设计在深度学习里有一个非常好的性质——消除了不同数据集之间的测序深度和文库大小差异，让模型可以跨数据集迁移。
我第一次用的时候，为了图省事直接拿原始count矩阵转成Seurat对象丢进去，结果报错报得一头雾水。后来才搞明白，Geneformer有自己特定的数据对象要求，需要用到它配套的预处理脚本，把表达矩阵转成.h5ad或者它支持的目录结构。注意不需要自己提前做标准化或者log转换，rank编码会自动处理量纲问题。不过有个前提，基因名必须是标准的Gene Symbol，要是你数据里基因ID是Ensembl ID，那就得先转一下，否则模型匹配不到预训练字典里的token。
3.2 微调与特征提取：把单细胞数据变成向量，才能喂给机器学习
Geneformer在这个项目里的定位有两个：第一个是用它的预训练输出做基因重要性排序，第二个是提取细胞层面的embedding特征。第二个功能后来是让我觉得最值的地方——每个细胞会被模型编码成一个高维向量，这个向量包含了这个细胞在预训练模型空间里的“状态位置”。有了这个向量，我们就可以把单细胞问题无缝转换成一个标准的机器学习分类任务。
但这里我踩过一个非常大的坑，就是直接拿全部细胞的embedding去做下游分类，结果类别之间的分离度很差。仔细排查之后发现，问题出在数据里混着好几个批次，embedding里编码了很强的批次信息，反而把生物学差异给盖住了。解决办法就是在提取embedding之前，把不同批次的细胞重新平衡一下，或者干脆在批次水平做一次Harmony整合再去提取。虽然Geneformer本身对批次有一定鲁棒性，但不是绝对的，尤其是在细胞类型组成差异很大的情况下。
注意：提取出来的embedding维度通常比较高（几百维），直接送进随机森林或者逻辑回归，特征数远大于样本数，容易过拟合。我的做法是先做PCA降维到50维左右，再进分类模型。别贪心保留太多主成分，下游结果不会因为多了20个主成分而变好，反而增加了噪音。
3.3 为什么两个工具要联合，而不是二选一
单独用scTenifoldKnk，找到的是共表达网络层面的调控变化；单独用Geneformer，找到的是预训练模型视角下的基因重要性。这两个视角都有各自的盲区——scTenifoldKnk对数据量要求高，细胞数太少网络推断不稳定；Geneformer虽然鲁棒，但它更擅长捕捉通用的调控规律，对特定数据集里的罕见扰动不敏感。把它们联合起来，等于用两个独立的方法从不同侧面交叉验证同一个生物学问题。
在这个项目里我是这么安排的：先用scTenifoldKnk做调控网络推断，拿到差异调控基因；同时用Geneformer提取细胞embedding和基因注意力分数；然后把三套特征（网络拓扑特征、Geneformer embedding、基因注意力分数）合并，进入机器学习流程。最后用机器学习选出来的关键特征基因，去和前面两组差异调控基因做交集，交集部分作为核心候选基因输出。这套流程最大的优势就是，最后得到的核心基因不是被单一算法“拍脑袋”选出来的，而是经过了两种独立算法和一套监督学习的层层筛选。
4. 环境准备与数据预处理：这一步决定后面所有结果的可信度
4.1 软件环境与工具版本清单
先交代一下我用的环境。操作系统是Ubuntu 20.04，R版本4.2.2，Python用的是3.9。R部分需要安装scTenifoldKnk，它依赖的包比较多，我建议用BiocManager安装依赖，直接install.packages容易因为缺少某些Bioconductor包而报错。Python部分主要是Geneformer环境，它是基于PyTorch的，所以我单独建了一个conda环境，Python 3.9配PyTorch 1.13或2.0都行。GPU是必须要有的，Geneformer虽然推理比训练轻量，但大规模单细胞数据在CPU上跑会让人崩溃。我用的是单张RTX 3090 24G，跑2万细胞的数据，微调一轮大约半小时，纯推理十分钟以内就能完成。
机器学习部分，我用的的是Python里的scikit-learn和xgboost，版本分别是1.2.2和1.7.6。为什么用这两个而不用深度学习做分类？因为在这个项目里，样本量本质上是细胞数或者病人数，深度学习在样本量不够大的情况下很容易过拟合，而随机森林和XGBoost在中小规模表格数据上既稳定又可解释，SHAP值还能给出特征重要性的清晰解释。基因数量再多，XGBoost配合正则化也能handle得住。
| 工具 | 版本 | 用途 |
|---|---|---|
| R | 4.2.2 | scTenifoldKnk、Seurat |
| Python | 3.9 | Geneformer、sklearn、xgboost |
| PyTorch | 1.13+ | Geneformer后端 |
| scTenifoldKnk | 1.6+ | 网络推断与虚拟敲除 |
| Geneformer | 0.1.0 | 预训练模型推理与微调 |
| CUDA | 11.7 | GPU计算 |
4.2 输入数据的清洗规范：先质控再分析，别把垃圾喂给模型
如果你以为单细胞数据拿过来就能直接跑，那就大错特错了。我在这套流程上栽过最大的跟头，就是数据没有做严格的质控，导致后面所有结果都带着批次效应的影子。scTenifoldKnk对输入数据的要求比较高，首先是基因表达矩阵必须是稀疏矩阵，其次细胞和基因的过滤标准要合理。我在进入正式分析前，会先做一套标准质控：线粒体基因比例超过20%的细胞剔除，基因数少于200或者多于6000的细胞剔除，表达量过低的基因直接筛掉。这些标准不是拍脑袋定的，而是根据数据分布和已知细胞质量阈值选的。
质控之后，还要处理一个很容易被忽略的问题——细胞类型组成。scTenifoldKnk在做组间比较时，如果两组里细胞类型组成差异很大，会导致网络差异其实反映的是细胞类型比例的变化，而不是真正的调控变化。所以我的习惯是，先做细胞类型注释，然后把组间细胞类型比例做成柱状图看一遍。如果差异太大，先考虑用harmony整合或者做细胞类型层面的配对比较，而不是直接全细胞混在一起跑网络。
4.3 数据矩阵格式与基因名转换：最容易报错但最容易被忽略的环节
跑scTenifoldKnk最常见的报错，不是算法本身的问题，而是输入格式不对。这个函数要求矩阵是dgCMatrix，你拿普通的matrix或者data.frame进去，虽然有时候也能跑，但极容易在后面某些步骤里出现维度不匹配的错误。我封装了一段代码，直接从Seurat对象里提取稀疏矩阵，顺便转换基因名格式。
library(Seurat) library(Matrix) library(scTenifoldKnk) # 读取Seurat对象 seurat_obj <- readRDS("your_seurat_object.rds") # 从Seurat对象提取稀疏矩阵，注意行列方向：行=基因，列=细胞 expr_matrix <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts") # 确保是dgCMatrix格式 if (!inherits(expr_matrix, "dgCMatrix")) { expr_matrix <- as(expr_matrix, "dgCMatrix") } # 把基因名统一成Gene Symbol # 如果原来是Ensembl ID，需要用AnnotationDbi做映射转换 library(org.Hs.eg.db) library(AnnotationDbi) gene_symbols <- mapIds(org.Hs.eg.db, keys = rownames(expr_matrix), keytype = "ENSEMBL", column = "SYMBOL") valid_genes <- !is.na(gene_symbols) expr_matrix <- expr_matrix[valid_genes, ] rownames(expr_matrix) <- gene_symbols[valid_genes] # 去除重复基因名（保留表达量总和最大的那个） expr_matrix <- expr_matrix[!duplicated(rownames(expr_matrix)), ] # 细胞过滤：取至少检测到20个基因的细胞 cells_valid <- colSums(expr_matrix > 0) >= 20 expr_matrix <- expr_matrix[, cells_valid] # 查看最终矩阵维度 dim(expr_matrix)这段代码看起来简单，但每一步都不能省。基因映射那一步，如果你的数据里有一部分基因映射不到Symbol，直接删除会导致基因数量骤减，我一般会检查一下映射成功率，低于90%就要看看是不是基因ID格式本身有问题。重复基因这里，最简单的做法是保留平均表达量最大的那个，但如果你后面要做Geneformer的rank编码，重复基因会影响排序结果，所以这里一定要处理干净。
提示：如果你的单细胞数据是10x Genomics的Cell Ranger输出，直接读
filtered_feature_bc_matrix文件夹里的矩阵文件就行，基因名一般已经是Symbol格式。如果是其他平台或者已处理过的数据，基因名格式一定要先检查清楚。
5. 实操核心：scTenifoldKnk网络推断与组间差异比较
5.1 一次完整的scTenifoldKnk运行流程和时间成本
清洗好数据之后，就可以跑正式的模型了。scTenifoldKnk的核心函数就一个，但参数选择直接影响结果质量，我给出一个参考配置。注意这里我用的是nSubsets、nCores这类控制计算量和并行度的参数，有条件的话尽量多给核。
# 分别提取对照组和处理组的表达矩阵 ctrl_matrix <- expr_matrix[, seurat_obj$group == "Control"] case_matrix <- expr_matrix[, seurat_obj$group == "Case"] # 如果两组细胞数差异过大，做下采样 if (ncol(case_matrix) > ncol(ctrl_matrix) * 1.5) { set.seed(42) case_matrix <- case_matrix[, sample(ncol(case_matrix), ncol(ctrl_matrix))] } # 对照组构建网络 ctrl_network <- scTenifoldKnk( X = ctrl_matrix, nSubsets = 5, # 子集数量，越大越稳但越慢 nCores = 8, # 并行核数 q = 0.01, # 网络边筛选阈值，越小边越少 minCells = 10 # 子集中最少细胞数 ) # 处理组构建网络 case_network <- scTenifoldKnk( X = case_matrix, nSubsets = 5, nCores = 8, q = 0.01, minCells = 10 ) # 保存网络结果，后面差异比较要用 saveRDS(ctrl_network, "ctrl_network.rds") saveRDS(case_network, "case_network.rds")我这里用的nSubsets = 5是保守配置，如果你细胞量很大且计算资源充足，可以加到10或者20，网络稳定性会更好。但要注意，nSubsets是几，计算时间就近似翻几倍。我跑一组1500基因、2000细胞的矩阵，5个子集、8核并行，大约要3到4小时。全基因组3万个基因想都别想，时间会变成以周为单位。这就是为什么前面要强调基因筛选的重要性。
5.2 网络差异边提取：写出两组之间调控关系变化最大的基因
两组网络都建好之后，就需要专门的比较步骤了。scTenifoldKnk提供了一个比较函数叫做scTenifoldKnkComparison，它会通过虚拟敲除的方式对比两组网络的响应变化，最终给出差异调控基因的p值和排序。这个p值不是传统差异表达那种p值，它是基于虚拟敲除后网络重构误差的统计检验，所以解读的时候千万不要说“表达差异显著”，而要说“调控扰动显著”。
# 读取之前保存的网络结果 ctrl_network <- readRDS("ctrl_network.rds") case_network <- readRDS("case_network.rds") # 组间网络差异比较 comparison_result <- scTenifoldKnkComparison( network1 = ctrl_network, network2 = case_network, q_cutoff = 0.05 # 差异调控显著性阈值 ) # 提取显著差异调控基因 significant_genes <- comparison_result$diffRegulation significant_genes <- significant_genes[significant_genes$p.adj < 0.05, ] significant_genes <- significant_genes[order(significant_genes$p.adj), ] # 保存结果 write.csv(significant_genes, "differential_regulation_genes.csv", row.names = FALSE) # 查看前10个差异调控基因 head(significant_genes, 10)从结果文件里，你能看到每个基因的p值、校正后p值以及变化方向。把这些基因拿出来，跟常规的差异表达基因做比较，你会发现交集通常只有20%到40%——这正是这个分析的价值所在，它找到了另一批“表达没变但调控关系明显扰动”的功能相关基因。这批基因在文章里可以写成这样的故事：它们在传统差异分析中没有信号，但在调控网络层面表现出显著扰动，可能代表了疾病状态下更早期的功能失调。
5.3 网络拓扑指标：把网络的生物学信息变成机器学习能用的数字
scTenifoldKnk输出的网络表本身是有结构的，但它还不是机器学习可以直接用的特征。为了后续的机器学习部分，我把每个基因在网络里的拓扑指标提取出来——比如度中心性、介数中心性、接近中心性、聚类系数、特征向量中心性等。这些指标的意义就在于，它们能把“这个基因在网络中的位置有多重要”量化成数字。
library(igraph) # 将scTenifoldKnk的网络结果转为igraph对象 # 网络表格式假设有from, to, weight三列 network_table <- ctrl_network$network # 示例提取方式，按实际输出调整 g <- graph_from_data_frame(network_table, directed = FALSE) # 计算拓扑特征 topo_features <- data.frame( gene = V(g)$name, degree = degree(g), betweenness = betweenness(g), closeness = closeness(g), clustering_coef = transitivity(g, type = "local"), eigenvector = eigen_centrality(g)$vector ) # 对处理组做同样操作，然后按基因名合并两个矩阵 # 这样每个基因在两组中都有对应的拓扑特征，可以作为机器学习输入这里我想补充一个经验：拓扑特征对网络构建质量非常敏感，如果网络边数太少，很多基因的介数中心性会算出来是0，这个特征就废了。所以网络构建阶段，q值的选取要反复试，我在自己的项目里试过q = 0.01和q = 0.05，前者边数大概在几万条，后者能到十几万条，但包含的虚假边也多。我最终选了0.01，因为它能保留核心强连接，同时把弱相关的噪声边过滤掉。
6. Geneformer推理与微调：从预训练模型到项目定制模型
6.1 Geneformer的tokenization与推理实操
Geneformer的使用分两段，第一段是推理（inference），第二段是微调（finetune）。推理阶段不需要改模型参数，直接加载预训练权重，把我们的数据转成token序列，然后让模型输出embedding和注意力分数。我用的数据是10x格式的.h5ad，先做了一次简单的细胞过滤，然后调用Geneformer的tokenizer脚本。
import anndata as ad import torch from geneformer import GeneformerTokenizer, GeneformerEncoder # 读取单细胞数据 adata = ad.read_h5ad("your_data.h5ad") # 用Geneformer的tokenizer把表达矩阵转成token序列 tokenizer = GeneformerTokenizer.from_pretrained("ctransfer/Geneformer") tokenized = tokenizer.tokenize_anndata( adata, gene_symbol_col="gene_symbols", cell_identifier="cell_id", output_file="tokenized_dataset.dataset", num_proc=8 ) # 加载预训练模型 model = GeneformerEncoder.from_pretrained("ctransfer/Geneformer") # 设置推理模式 model.eval() with torch.no_grad(): cell_embeddings = model.encode(tokenized, layer="avg_pool")这段代码里，layer="avg_pool"的意思是取模型各层输出的平均池化结果作为细胞embedding。实际使用中，你可以尝试不同层，有的文章说最后一层更能捕获复杂关系，有的说中间层更保留局部结构信息。我的经验是，如果下游做细胞分类，最后一层加avg_pool效果最稳；如果做基因重要性分析，用注意力分数更直接。
6.2 微调策略：用小样本数据把通用模型变成你的领域专家
如果你对模型输出的细胞embedding做分类效果不满意，可以考虑做一次微调。Geneformer的微调用的是它提供的分类器接口，本质上是在预训练模型后加一个分类头，用你的单细胞数据和标签做有监督训练。微调的优势是，模型在预训练阶段学到的是通用的调控语法，微调阶段可以把这些知识聚焦到你的具体问题上来。
from geneformer import GeneformerClassifier # 准备带标签的数据（比如细胞类型标签或疾病/对照标签） finetune_data = { "cell_id": adata.obs["cell_id"].tolist(), "label": adata.obs["group"].tolist() # Control / Case } # 创建微调分类器 classifier = GeneformerClassifier( model_name="ctransfer/Geneformer", num_classes=2, fine_tune=True ) # 开始微调 classifier.fit( train_data=finetune_data, output_dir="./geneformer_finetuned", epochs=3, batch_size=16, learning_rate=1e-5 )微调的时候有几个关键参数要特别注意。learning_rate我一般设得比较小，1e-5左右，因为预训练模型已经学得很好，学习率太大会把预训练的知识冲掉。epochs设3轮就够了，微调轮次太多同样会导致过拟合——因为你的标注数据可能就几千个细胞，但模型参数是上亿级别的。还有一个容易忽略的：微调时细胞标签的平衡性。如果对照和处理组的细胞数量差距很大，模型会不自觉地偏向多数类，我通常会在准备训练数据时对多数类做下采样，或者给少数类加大loss权重。
6.3 从Geneformer拿到基因注意力分数
除了细胞embedding，Geneformer另一个很有价值的输出是基因层面的注意力分数。注意力分数反映的是模型在推断过程中认为哪些基因对当前细胞状态的贡献更大。把每个基因在所有细胞中的注意力分数取平均值，就得到了一个基因重要性排序。这个排序和表达量高低并没有直接关系，有时候一个表达量极低的转录因子会获得非常高的注意力分数，因为它可能在调控网络中扮演关键节点。
我在项目里把Geneformer的基因注意力分数和scTenifoldKnk的差异调控基因做了交集，发现重叠率比我预想的要高，大约有30%到40%。这个交集本身就是一个很好的结果，可以直接写成“两种独立算法共同识别出的关键调控基因”。在文章里，我会把这两个方法画成一个Venn图或者韦恩图，作为核心figure之一，效果很直观。
注意：Geneformer的注意力分数并不是一个概率值，不同层的注意力分数量纲可能有差异。做跨样本平均之前，最好先做一次min-max归一化，否则某些层的高数值基因会主导最终排序，导致偏差。
7. 机器学习融合：把单细胞层面的信息变成可训练的表格数据
7.1 特征工程：怎样把三套特征合并成一个规范的训练矩阵
现在手里已经有三套特征了：scTenifoldKnk输出基因层面的网络拓扑特征、Geneformer输出细胞层面的embedding特征和基因层面的注意力分数特征。但在机器学习这一步，我们不能直接把基因层面的特征和细胞层面的特征混在一起扔进去，因为它们的样本单位不一样。我的做法是分成两个层次做分类：
第一个层次是细胞级别的分类：每个细胞的标签就是它属于对照组还是处理组，特征就是Geneformer提取的细胞embedding。这个任务用来回答“细胞状态是否在两组之间存在系统性差异”，本质上是一个“能不能分得开”的检验。第二个层次是基因级别的特征排序：每个基因有一堆网络拓扑特征和注意力分数特征，目标是用机器学习选出最能区分两组的基因特征，这部分可以理解为“哪些基因的调控角色变化最显著”。
import pandas as pd import numpy as np from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler # 读取细胞embedding特征 # 假设已从Geneformer提取并保存为csv cell_features = pd.read_csv("cell_embeddings.csv", index_col=0) cell_labels = cell_features["group"] cell_features = cell_features.drop(columns=["group"]) # 标准化特征（注意用训练集拟合scaler，再用同一scaler转换测试集） scaler = StandardScaler() cell_features_scaled = scaler.fit_transform(cell_features) # 划分训练和测试集 from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split( cell_features_scaled, cell_labels, test_size=0.2, stratify=cell_labels, random_state=42 ) # 随机森林分类 rf = RandomForestClassifier(n_estimators=500, max_depth=8, random_state=42) rf.fit(X_train, y_train) print("训练集准确率：", rf.score(X_train, y_train)) print("测试集准确率：", rf.score(X_test, y_test))这个代码看起来很简单，但里面的门道不少。stratify=cell_labels参数非常重要，它保证训练集和测试集中两组的比例一致，否则划分出来的测试集可能全是某一组，评估结果就没有意义。随机森林我设了max_depth=8，控制树的复杂度，防止在细胞量很少的情况下过拟合。实际跑出来的结果，如果测试集准确率在0.85以上，基本上可以说明两组的细胞状态在转录组水平上是被系统性区分的，这个结论可以理解为“疾病组的细胞状态发生了整体性的改变”。
7.2 特征筛选与可解释性：SHAP值找出真正的关键特征
随机森林只能告诉你哪些特征重要，但它不能告诉你这些特征是如何影响预测结果的。为了做出可解释的结论，我用SHAP来分析每个特征对模型预测的贡献方向和大小。SHAP值的好处是，它从博弈论的角度计算每个特征对预测结果的边际贡献，结果稳定、可解释性强。
import shap # 创建SHAP解释器并计算值 explainer = shap.TreeExplainer(rf) shap_values = explainer.shap_values(X_test) # 可视化特征重要性汇总图 shap.summary_plot(shap_values, X_test, feature_names=cell_features.columns.tolist())SHAP结果出来后，重点关注排名前20的特征基因。这些基因就是机器学习认为最能区分两组的“信息核心”。把这些基因跟scTenifoldKnk的差异调控基因取交集，通常能锁定10到20个核心候选基因。这个列表就是最终要放进文章里的“明星基因”，后续的GO/KEGG富集分析、蛋白互作网络分析，都围绕这个列表展开。用机器学习做特征筛选，还有一个隐藏的好处：审稿人非常吃这套逻辑——你不是只靠一个算法拍脑袋，而是有严谨的统计和模型验证支撑。
7.3 多模型交叉验证与过拟合控制
单一模型的结果总是没那么有说服力。我在最终提交的结果里，同时跑了随机森林、XGBoost和逻辑回归三个模型，然后把三者的分类准确率和AUC列成表格，做一次模型对比。XGBoost在这个场景下往往表现最好，因为它能自动处理特征间的非线性交互；逻辑回归虽然准确率略低，但胜在系数可以直接解读为基因表达的log-odds贡献。表格里同时给出三个模型的准确率、精确率、召回率和AUC，审稿人看到这种严谨的模型对比，第一反应是分析者基本功扎实。
from xgboost import XGBClassifier from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score, classification_report import pandas as pd models = { "RandomForest": RandomForestClassifier(n_estimators=500, max_depth=8, random_state=42), "XGBoost": XGBClassifier(n_estimators=500, max_depth=4, learning_rate=0.01, random_state=42), "LogisticRegression": LogisticRegression(max_iter=1000, C=0.1) } results = [] for name, model in models.items(): model.fit(X_train, y_train) y_pred = model.predict(X_test) y_proba = model.predict_proba(X_test)[:, 1] report = classification_report(y_test, y_pred, output_dict=True) results.append({ "Model": name, "Accuracy": report["accuracy"], "AUC": roc_auc_score(y_test, y_proba) }) comparison_df = pd.DataFrame(results) print(comparison_df)这里我想提醒一点：交叉验证用的是训练集内部多次划分，但最终一定要留出一部分完全没参与过任何特征选择或模型调参的数据，作为最终测试集。不然你调了一个星期参数，最后看到的测试集结果很有可能是虚高的。我习惯的做法是，在一开始就切出15%的数据锁在一个文件里，整个分析流程中完全不去碰它，只到最后做验证才解开。
8. 富集分析与可视化：R语言组间GO分析的标准姿势
8.1 用clusterProfiler做GO/KEGG富集的完整流程
机器学习筛选出的核心基因列表，接下来就要交给富集分析来赋予生物学意义。这一步用的是R语言里最常用的clusterProfiler包。它的输入很简单，一个是基因列表，一个是背景基因列表，输出结果是该基因列表在哪些GO条目或KEGG通路上显著富集。
library(clusterProfiler) library(org.Hs.eg.db) # 读取机器学习筛选出的核心基因 core_genes <- read.csv("core_genes.csv")$gene # 背景基因：建议使用全部检测到的基因作为背景，而不是全基因组 background_genes <- rownames(expr_matrix) # ID转换：如果核心基因是Symbol，需要转成ENTREZID entrez_genes <- bitr(core_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) entrez_background <- bitr(background_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # GO富集分析 go_results <- enrichGO( gene = entrez_genes$ENTREZID, universe = entrez_background$ENTREZID, OrgDb = org.Hs.eg.db, ont = "BP", # Biological Process pAdjustMethod = "BH", qvalueCutoff = 0.05, readable = TRUE ) # 查看前20个显著富集的GO条目 head(as.data.frame(go_results), 20) # KEGG富集分析 kegg_results <- enrichKEGG( gene = entrez_genes$ENTREZID, organism = "hsa", pvalueCutoff = 0.05, qvalueCutoff = 0.05 ) # 保存结果 write.csv(as.data.frame(go_results), "GO_enrichment_results.csv", row.names = FALSE) write.csv(as.data.frame(kegg_results), "KEGG_enrichment_results.csv", row.names = FALSE)富集分析里有几个细节特别容易影响结论。第一，背景基因一定要用你实际检测到的基因，不要用全基因组，不然富集的显著性会被严重夸大——因为你拿一个只有2万个基因的检测面板去跟全基因组3万个基因比，很多本来不显著的富集也变显著了。第二，qvalueCutoff是严格得多重校正确认的阈值，建议用0.05，不要手一抖改成0.5。第三，如果核心基因只有十几个甚至几个，富集分析很可能什么结果都跑不出来，这时候不要去硬凑，换一个思路，直接在文章里放蛋白互作网络分析也行，效果不比富集差。
8.2 富集结果的可视化：气泡图、条形图和网络图的选择
富集结果出来之后，怎么可视化也是有讲究的。GO分析最常见的就是气泡图和条形图。气泡图横轴是GeneRatio，纵轴是GO条目名称，气泡大小代表富集的基因数，颜色代表p值。这种图用来展示富集的整体格局很直观。如果两组之间要并列比较GO条目的差异，我推荐画一个左右对比的条形图，左侧是上调相关，右侧是下调相关，一眼就能看出两组生物学过程的差异。
# 绘制GO气泡图 library(ggplot2) go_plot <- dotplot(go_results, showCategory = 20) + theme_bw(base_size = 14) + labs(title = "GO Biological Process Enrichment") ggsave("GO_dotplot.pdf", go_plot, width = 8, height = 10) # KEGG通路图，映射到具体通路 library(pathview) kegg_results_entrez <- enrichKEGG( gene = entrez_genes$ENTREZID, organism = "hsa", pvalueCutoff = 0.05 ) # 找到显著性最高的通路，用pathview绘制通路图 top_pathway <- as.data.frame(kegg_results_entrez)[1, "ID"] pathview( gene.data = entrez_genes$ENTREZID, pathway.id = top_pathway, species = "hsa", out.suffix = "core_genes", limit = list(gene = 2, cpd = 1) )富集图的美观度直接决定审稿人对你数据质量的第一印象。我通常会把气泡图的PNG分辨率设在300dpi以上，字体大小调到14以上，配色用RColorBrewer的Set3或者ggplot2默认的hue色系，尽量别用饱和度太高的对比色。如果某个通路富集特别显著但只富集了2个基因，这种结果放在文章里反而容易被质疑，建议在figure中标注清楚基因数，避免误导。
9. 常见问题与避坑清单：这套流程最容易翻车的7个环节
9.1 报错手册与解决方案速查
我梳理了一遍自己跑这套流程时踩过的所有坑，按出现频率排个序，做成一个速查表。下面的问题里，前三个几乎每个第一次跑scTenifoldKnk的人都会遇到。
| 报错或异常现象 | 可能原因 | 解决方案 |
|---|---|---|
| scTenifoldKnk报“subscript out of bounds” | 矩阵中存在全零行或者基因名有重复 | 做好基因过滤和去重，确保每个基因至少在一个细胞中表达 |
| Geneformer报“token not in vocabulary” | 基因名不是标准Symbol格式 | 用org.Hs.eg.db转换ID，删除无法映射的基因 |
| 随机森林训练集准确率100%，测试集只有60% | 严重过拟合 | 降低模型复杂度（max_depth调小），增加正则化，用交叉验证评估 |
| 两组细胞数量不平衡导致网络推断不稳定 | 细胞数差距过大 | 对多的组做下采样，或者按细胞类型分层后分别比较 |
| GO富集结果全空白 | 核心基因列表太短或ID转换失败 | 检查基因列表长度，确认是Symbol格式后重新转换 |
| Geneformer微调loss不下降 | 学习率太大或标签不平衡 | 把learning_rate调到1e-5，对多数类下采样 |
| 网络推断耗时过长 | 基因数量太大 | 先用差异表达和已知通路基因做筛选，控制在3000以内 |
9.2 几个容易被忽视的经验小技巧
第一，scTenifoldKnk的随机种子问题。它在内部有多个基于随机采样的步骤，如果不设置种子，同一份数据跑两次结果会有细微差异。虽然差异不大，但发文章追求可复现，一定要在跑之前设置set.seed()，并在方法部分写明种子值。
第二，Geneformer的pretrained模型版本问题。Geneformer有多个训练好的模型版本，不同版本的预训练数据量和组织来源不完全一样。如果你的数据是外周血单细胞，尽量选择基于外周血或免疫细胞预训练的版本，效果会比全组织版更好一点，这个在模型加载时可以指定model_name参数。
第三，机器学习特征筛选和网络差异基因取交集时，要注意方向一致。如果scTenifoldKnk的差异调控基因符号和SHAP选出的高重要特征基因方向性存在矛盾（比如一个基因在一个方法中上调显著，在另一个方法中下调显著），不要强行合并。我在第一个项目里就因为这个原因丢掉了一批基因，后来重新审视发现是自己把基因方向搞反了，纠正好之后交集率明显上升。
第四，整条流程的运行时间规划非常关键。14天拿下2区的意思是，前5天做数据质控和单细胞标准流程，中间5天跑scTenifoldKnk和Geneformer，最后4天跑机器学习和写文章。如果中间某一步跑挂了，时间就非常紧张。我的建议是，在拿到数据后立刻把scTenifoldKnk先挂上跑着，因为它是整套流程里最耗时的环节，即使后面数据需要微调，网络推断的结果也有很大概率可以复用。
9.3 两个被视为“捷径”但实为陷阱的操作
一个是直接用别人的差异表达基因列表当背景，跑富集分析。这样做出来的富集结果看似很漂亮，但本质上是用一个已经筛选过的子集去做富集，富集结果受筛选标准影响极大，审稿人一看你的背景基因列表就能发现问题。另一个是直接用Scanpy的内置批处理函数替代harmony，虽然Scanpy的处理也能用，但单细胞数据在批次效应上的处理能力，目前Harmony的结果还是要更稳定一些，尤其当组间细胞类型组成差异明显时，差异会非常显著。既然这套流程已经够复杂了，每一步都尽量用最稳妥的方案，别给自己挖坑。
10. 最终mpact：这套流程能给你带来什么样的数据产出和文章素材
跑完这套流程，你手里会攒下几件硬通货：scTenifoldKnk组的组间差异调控基因列表、Geneformer的细胞embedding与基因注意力分数、机器学习多模型的对比结果与核心基因列表、核心基因的GO/KEGG富集结果。这些素材组合起来，足够撑起一篇文章的3到4个核心figure：网络差异基因热图、Geneformer embedding的UMAP图（按组着色）、机器学习模型AUC曲线对比图、核心基因富集气泡图。再补一个验证性的qPCR或者WB实验，故事基本就完整了。
我回想自己从传统差异表达到这套虚拟联合分析的过程，最大的心理转变是学会了不把“虚拟”两个字看成贬义。虚拟分析不是替代湿实验，而是在湿实验之前帮你把候选范围从几千个基因缩小到几十个，再配合机器学习验证排序，让后续实验验证的命中率大幅提高。这套工作流里，scTenifoldKnk负责从网络视角找线索，Geneformer负责从跨数据集的先验知识里找支撑，机器学习负责用统计学的语言给线索排序，最后富集分析负责赋予线索生物学故事。
我个人在实际操作中的体会是，整套流程最值钱的部分不是某一个工具的跑通，而是两条独立分析路径在一个数据集上同时指向了高度重叠的基因列表。这种交叉验证的一致性，比任何一个单独算法的p值都有说服力。如果你正在做单细胞组间比较却觉得常规流程写不出新意，可以试试这套组合，也不用非等到数据完美才开始，先拿一批小规模数据把流程跑通，后面换大数据就是纯粹的复用。最后再分享一个小技巧：每一步分析都随手保存rds和csv文件，命名带上日期和参数版本，这套流程里任何一步重跑的成本都很高，好的文件管理习惯能让你在最后写文章阶段省下大量来回折腾的时间。
