单细胞轨迹分析利器CytoTRACE：从基因计数预测细胞命运

1. 项目概述：从单细胞数据中“听”见细胞命运的轨迹

如果你手头有一批单细胞转录组数据，看着UMAP图上密密麻麻的细胞群，除了知道它们属于不同的细胞类型，会不会好奇：这些细胞之间有没有“辈分”关系？哪个细胞更“年轻”，哪个更“老”？它们沿着什么路径在分化？这就是单细胞轨迹分析，或者说拟时分析要回答的核心问题。它不满足于静态的细胞分类，而是要重建细胞状态演变的动态过程，就像给细胞拍一部“成长纪录片”。

今天要聊的CytoTRACE，就是这部“纪录片”的一位独特导演。它不像Monocle、PAGA那样依赖复杂的图论或机器学习模型来构建轨迹，而是另辟蹊径，从一个非常直观的生物学假设出发：分化程度越高的细胞，其转录组越特化，表达的基因总数往往会减少。因此，它用每个细胞表达的基因数量（减去一些噪音）作为一个简单的指标，来推断细胞的发育潜能或分化状态。数值越高，代表细胞越“幼稚”、分化潜能越大；数值越低，则代表细胞越“成熟”或特化。这个方法在2018年由Gulati等人发表在《Cell》上，因其原理简单、计算快速、无需先验知识而备受关注，特别适合在分析的早期阶段快速评估细胞的分化层次。

简单来说，CytoTRACE帮你做两件事：第一，给你数据里的每个细胞打一个从0到1的“分化潜力”分数；第二，基于这个分数，对基因进行排序，找出那些可能驱动早期命运决定的“潜力股”基因。对于刚拿到单细胞数据，想快速看看里面有没有分化轨迹、哪个群可能是起点的研究者来说，它是一个非常高效的“侦察兵”。

2. 核心原理与算法逻辑拆解：为什么数基因能预测命运？

CytoTRACE的核心思想，源于发育生物学中的一个经典观察：多能干细胞或祖细胞通常具有更活跃、更广泛的转录活动，而终末分化的细胞则倾向于高表达少数特定功能基因，整体转录的广度下降。CytoTRACE将这一观察量化成了一个可计算的指标。

2.1 算法核心四步走

它的计算流程可以概括为四个关键步骤，理解了这几步，你就能明白结果是怎么来的，以及该如何解读。

第一步：基因表达矩阵的过滤与准备输入是一个标准的单细胞RNA-seq基因表达矩阵（细胞×基因）。首先，算法会进行基础过滤，通常只保留在至少一定比例细胞中表达的基因，以去除大量零表达带来的噪音。这里的一个关键点是，CytoTRACE关注的是基因是否被“检测到”，而不是表达量高低，因此它通常使用原始计数或二值化（表达记为1，不表达记为0）后的数据。

第二步：计算每个细胞的基因表达数量这是最直观的一步。对于每个细胞，计算其表达量大于零的基因数量。我们暂时称这个原始值为“原始基因数”。但是，直接使用这个数会有明显偏差，因为测序深度（每个细胞捕获的mRNA总量）深的细胞，天然会检测到更多基因，这并不一定代表其分化潜能高，可能只是技术因素。

第三步：校正测序深度偏差为了剔除技术噪音，CytoTRACE采用了一种基于排序的局部回归校正方法。它会绘制每个细胞的“原始基因数”与“总UMI数”（或总读数）的散点图。然后，计算一条拟合曲线（通常是Loess回归），这条曲线代表了在给定测序深度下，“预期”的基因数。最后，用细胞的“原始基因数”减去该深度下的“预期基因数”，得到残差。这个残差就是初步校正后的分数，它反映了基因数相对于其测序深度的“超额”部分，这更可能来源于生物学本质。

第四步：平滑与最终CytoTRACE分数计算单细胞数据本身存在技术噪音和生物学异质性。为了得到更稳健的趋势，CytoTRACE引入了细胞间的相似性进行平滑。它首先基于基因表达谱计算细胞间的相似性矩阵（如皮尔逊相关系数），然后对于每个细胞，其最终的CytoTRACE分数是其自身校正后分数与其相似邻居细胞分数的加权平均。这一步相当于让每个细胞的分数“参考”了其所在局部细胞群体的信息，使发育轨迹在降维图上更加连续和平滑。最终分数会被归一化到0到1之间，1代表预测分化潜能最高（最幼稚），0代表最低（最成熟）。

2.2 与主流轨迹分析方法的本质区别

理解CytoTRACE的独特之处，能帮助你在正确场景下选择它。

• 无需预设起点或终点：像Monocle、Slingshot这类方法，通常需要用户指定一个根节点（起始细胞群）或提供细胞的时间序列信息。CytoTRACE完全无监督，它自己从数据中推断“起点”（分数高的细胞群）。
• 基于全局特征而非局部连接：PAGA、Diffusion Map等方法侧重于分析细胞在高维空间中的邻近关系和过渡概率来构建轨迹图。CytoTRACE不直接构建“路径图”，而是先给每个细胞一个全局性的“潜能”标量值。你可以后续将这个值作为颜色映射到UMAP/t-SNE图上观察梯度，或者作为伪时间值输入其他工具进行分支分析。
• 计算速度极快：由于核心是计数和回归校正，避开了复杂的图优化或概率建模，CytoTRACE的计算通常在几分钟内完成，对于大型数据集（数万细胞）非常友好。
• 假设的局限性：它的核心假设（基因数多=更幼稚）在大多数发育、分化场景中成立，但在一些特殊情况下可能失效。例如，在细胞周期活跃的群体中，处于S/G2/M期的细胞由于DNA复制，整体转录本数量会增加，可能导致CytoTRACE分数被高估。又或者，某些终末分化细胞（如高度代谢活跃的肝细胞）可能仍然表达大量基因。因此，它通常需要与其他生物学标记结合验证。

注意：CytoTRACE分数是一个连续的分化状态指标，而不是离散的“伪时间”。伪时间通常将细胞排列在一条或多条具有方向性的路径上，而CytoTRACE分数更接近于一个细胞内在的“干性”或“分化度”标尺。你可以把它看作伪时间分析一个极好的起点或补充。

3. 实战演练：使用R语言完整复现CytoTRACE分析

纸上得来终觉浅，我们直接上手，用一个真实的单细胞数据集（比如一个造血干细胞分化的数据集）来跑通整个流程。这里以R语言环境为例，因为原版CytoTRACE就是一个R包。

3.1 环境准备与数据加载

首先，确保你的R环境已经就绪。CytoTRACE包可以从GitHub安装。

# 安装必要的包 if (!require("devtools")) install.packages("devtools") devtools::install_github("digitalcytometry/cytotrace") # 加载包 library(CytoTRACE) library(Seurat) # 假设我们使用Seurat对象，这是目前最流行的单细胞分析框架 library(ggplot2) library(patchwork) # 用于拼图 # 加载示例数据集，这里我们使用内置数据集或从公开数据库下载 # 例如，我们可以模拟一个数据加载过程，实际中请替换为你的数据 # 假设我们已经有了一个经过标准预处理（QC、标准化、降维、聚类）的Seurat对象 `seurat_obj` # 本示例假设数据已准备好

如果你的数据是Seurat对象，需要从中提取表达矩阵。CytoTRACE要求输入一个矩阵，行是基因，列是细胞。

# 从Seurat对象中提取标准化后的数据（例如log1p(CPM)后的数据），或者原始计数。 # 作者推荐使用去除了批次效应的标准化数据，但并非必须。 # 我们这里使用log标准化后的数据作为输入，这也是CytoTRACE文档中常用的。 expr_matrix <- as.matrix(seurat_obj@assays$RNA@data) # 获取log标准化数据 # 或者使用原始计数，可能对于基因计数更敏感 # expr_matrix <- as.matrix(seurat_obj@assays$RNA@counts) # 确保矩阵是数值矩阵，并且行名是基因名，列名是细胞ID dim(expr_matrix)

3.2 运行CytoTRACE核心分析

运行主函数非常简单。CytoTRACE()函数是核心。

# 运行CytoTRACE分析 # 注意：对于大型数据集（>5000细胞），可以使用`ncores`参数进行并行计算加速 cyt_result <- CytoTRACE(expr_matrix, ncores = 4) # 使用4个CPU核心 # 查看结果结构 names(cyt_result) # 通常会包含： # - `CytoTRACE`: 每个细胞的最终CytoTRACE分数（数值向量） # - `CytoTRACErank`: 细胞按分数从高到低的排名 # - `exprMatrix`: 过滤后的表达矩阵 # - `gcs`: 每个细胞的基因计数特征（Gene Count Signature） # - `filteredCells`: 过滤掉的细胞（如果有） # - `filteredGenes`: 过滤掉的基因

运行完成后，最重要的结果就是cyt_result$CytoTRACE，它是一个以细胞ID为名字、CytoTRACE分数为值的向量。

3.3 结果可视化与解读

将计算出的分数整合回你的Seurat对象，并可视化。

# 将CytoTRACE分数添加到Seurat对象的元数据中 seurat_obj$cytotrace_score <- cyt_result$CytoTRACE[colnames(seurat_obj)] # 可视化：在UMAP图上用颜色深浅表示CytoTRACE分数 p1 <- DimPlot(seurat_obj, reduction = "umap", group.by = "celltype", label = TRUE) + ggtitle("Cell Type") p2 <- FeaturePlot(seurat_obj, features = "cytotrace_score", reduction = "umap") + scale_colour_gradientn(colours = c("blue", "green", "yellow", "red")) + ggtitle("CytoTRACE Score (Red: High/Immature, Blue: Low/Mature)") p1 + p2 # 并排查看细胞类型注释和CytoTRACE分数分布

解读UMAP图：

• 寻找颜色梯度：关注从红色（高分）到蓝色（低分）的连续过渡区域。这很可能就是一条分化轨迹。
• 验证起点：查看你已知的干细胞或祖细胞群（例如，造血系统中的HSC（造血干细胞）群）是否被赋予了最高的CytoTRACE分数（显示为红色/黄色）。这是验证分析是否合理的第一步。
• 观察分支：如果细胞类型形成多个簇，观察每个簇内部是否有一个从高到低的分数梯度。这可能意味着每个簇代表一条独立的分化路径。

除了整体分数，CytoTRACE还能预测与高分化潜能相关的基因。

# 获取预测的“干细胞性”相关基因（即表达与CytoTRACE分数正相关最强的基因） # 使用`plotCytoTRACE`函数或直接提取结果 # 我们可以计算基因分数（gene score） gene_scores <- cyt_result$gcs # 基因计数特征，可以近似看作基因的重要性指标 top_genes <- names(sort(gene_scores, decreasing = TRUE))[1:20] print("Top 20 genes predicted to be associated with high differentiation potential:") print(top_genes) # 你可以用这些基因做后续验证，例如查看它们是否在已知的干性基因集中

3.4 高级应用：结合其他轨迹推断工具

CytoTRACE分数可以作为一个优秀的“伪时间”起点，输入到其他更擅长构建复杂分支轨迹的工具中。例如，我们可以用Slingshot。

# 假设我们已经有了在低维空间（如PCA）的嵌入 library(slingshot) # 获取细胞在低维空间的坐标（例如PCA的前50个主成分） pca_coords <- Embeddings(seurat_obj, reduction = "pca")[, 1:50] # 使用CytoTRACE分数作为“起始点”的指引 # 我们假设CytoTRACE分数最高的细胞簇是起点 start_cluster <- seurat_obj$seurat_clusters[which.max(seurat_obj$cytotrace_score)] # 运行Slingshot，指定起始簇 sce <- as.SingleCellExperiment(seurat_obj) # 转换为Slingshot需要的对象 sce <- slingshot(sce, clusterLabels = 'seurat_clusters', reducedDim = 'PCA', start.clus = start_cluster) # 提取Slingshot推断的伪时间 pseudotime <- slingPseudotime(sce) # 可以将多条曲线的伪时间整合或选择第一条曲线进行分析

这种组合策略既利用了CytoTRACE无监督推断起点的优势，又借助了Slingshot在构建复杂分支轨迹方面的能力。

4. 参数调优、注意事项与避坑指南

在实际操作中，直接运行默认参数可能不会总是得到理想结果。下面是一些关键的调参点和常见陷阱。

4.1 关键参数解析

• expr_matrix输入：这是最重要的选择。官方推荐使用标准化但未缩放的数据（如log1p(CPM)）。使用原始计数可能会因测序深度差异过大而引入强烈噪音。绝对避免使用经过ScaleData缩放后的数据（均值为0，方差为1），这会彻底破坏基因计数的生物学意义。
• enableFast参数：默认是TRUE，它会使用一种快速近似算法来计算细胞相似性。对于绝大多数数据集，这已经足够且能极大提升速度。只有在结果非常不连续、且数据量不大时，可以尝试设为FALSE使用精确计算。
• ncores：并行计算核心数。对于超过1万个细胞的数据集，建议设置与CPU核心数相近的值以节省时间。
• 基因过滤：CytoTRACE内部会过滤低表达基因。如果你事先已经进行了严格的基因过滤，可以关注cyt_result$filteredGenes看看是否过滤掉了你关心的基因。通常无需干预。

4.2 常见问题与解决方案

问题1：CytoTRACE分数在我的UMAP图上没有显示出清晰的梯度，而是斑驳的斑点。

• 可能原因1：数据批次效应强烈。不同批次间细胞的技术差异可能掩盖了生物学上的发育梯度。
• 解决方案：在运行CytoTRACE之前，务必使用Harmony、BBKNN或Seurat的IntegrateData等功能进行批次校正。对校正后的整合数据再运行CytoTRACE。
• 可能原因2：细胞类型过于离散，发育轨迹不连续。如果你的数据包含多个完全独立的细胞谱系（如同时有神经元、免疫细胞和上皮细胞），它们之间可能没有连续的过渡状态。
• 解决方案：尝试分群单独分析。先根据广谱的细胞类型注释（如major.celltype）将数据子集化，对每个可能包含连续分化的子集（例如所有的免疫细胞）单独运行CytoTRACE。

问题2：已知的干细胞群（如HSCs）没有得到最高分，反而是某个分化中的细胞群分数最高。

• 可能原因1：细胞周期影响。该高分群可能处于活跃的细胞周期（S/G2/M期），转录活动整体增强。
• 解决方案：计算并回归掉细胞周期评分的影响。在Seurat中，可以使用CellCycleScoring函数，并在运行CytoTRACE时，考虑使用回归了细胞周期效应后的表达矩阵（ScaleData时指定vars.to.regress = c(“S.Score”, “G2M.Score”)，然后取@scale.data？注意：这里需谨慎，最好使用SCTransform的vars.to.regress或直接对原始计数进行细胞周期回归后再标准化）。
• 可能原因2：该数据集不适用于CytoTRACE的核心假设。有些细胞类型的发育可能不伴随基因表达数量的显著下降。
• 解决方案：用已知的干性标记基因（如POU5F1(OCT4),NANOG,SOX2用于多能干细胞；CD34,KIT用于造血干细胞）进行双重验证。如果CytoTRACE高分群也高表达这些标记，那结果可能是合理的；否则，应考虑使用其他轨迹推断方法（如基于扩散图的DPT或基于RNA速度的方法）作为主要工具，将CytoTRACE作为参考。

问题3：运行速度非常慢，甚至内存不足。

• 可能原因：细胞数过多（>5万）或基因数过多。
• 解决方案：
  1. 使用enableFast = TRUE（默认）。
  2. 增加ncores参数充分利用多核。
  3. 在运行前，对表达矩阵进行更激进的基因过滤，例如只保留在所有细胞中表达比例前5000-8000个高变基因。这通常不会丢失关键信息，因为与发育相关的基因大多是高表达的。
  4. 如果可能，先进行细胞亚群的下采样分析，找到感兴趣的区域后再在全数据集上验证。

4.3 结果可信度验证策略

不要盲目相信任何一个计算工具的输出。对于CytoTRACE的结果，建议从以下几个维度交叉验证：

1. 已知标记物验证：检查CytoTRACE预测的高潜能细胞群是否高表达该领域公认的干/祖细胞标记基因。同时，检查低分群是否高表达终末分化标记基因。
2. 与RNA速度结果对比：RNA速度能提供细胞状态变化的向量方向。将CytoTRACE分数作为伪时间，与RNA速度流场叠加在同一个UMAP图上，观察速度向量的方向是否大体上从高分区域指向低分区域。如果方向一致，则结果可信度大增。
3. 轨迹一致性检验：使用其他至少1-2种轨迹推断方法（如Diffusion Map伪时间、Monocle3、PAGA轨迹）独立分析。比较不同方法推断出的“起点”和细胞排序是否具有一致性。如果多种方法指向相似的结论，那么这个轨迹就更可能是真实的生物学信号，而非算法假象。
4. 功能富集分析：对CytoTRACE预测的Top 100个与高潜能正相关的基因进行GO或KEGG富集分析。看看这些基因是否显著富集在“干细胞多能性维持”、“细胞周期”、“DNA复制”等相关通路。这从功能层面提供了支持。

5. 扩展应用场景与代码复现案例

CytoTRACE的应用远不止于经典的发育生物学。结合最新的热点，我们可以探索一些更前沿或实用的分析场景。

5.1 场景一：在肿瘤微环境中鉴定干细胞样癌细胞

肿瘤异质性是癌症研究的核心。肿瘤内部也存在类似干细胞的群体，即癌症干细胞，它们被认为是耐药和复发的根源。CytoTRACE可以用来从肿瘤单细胞数据中识别这些潜在的高恶性、去分化细胞亚群。

分析思路：

1. 从肿瘤单细胞数据中，分离出恶性细胞（通常通过拷贝数变异推断或已知标记）。
2. 仅对恶性细胞子集运行CytoTRACE。
3. 将细胞按CytoTRACE分数从高到低排序。
4. 分析高分细胞群（前10%-20%）特有的基因表达特征，并与已知的癌症干细胞标记（如CD44,ALDH1A1,CD133等）进行比对。
5. 进一步，可以检查这些高分细胞是否富集在特定的空间位置（如有空间转录组数据），或与特定的免疫微环境特征相关。

# 假设 `seurat_tumor` 是只包含肿瘤细胞的Seurat对象 # 运行CytoTRACE cyt_tumor <- CytoTRACE(as.matrix(seurat_tumor@assays$RNA@data), ncores=4) seurat_tumor$cytotrace_score <- cyt_tumor$CytoTRACE[colnames(seurat_tumor)] # 定义“高干细胞潜能”的阈值，例如前20% score_threshold <- quantile(seurat_tumor$cytotrace_score, probs = 0.8) seurat_tumor$csc_potential <- ifelse(seurat_tumor$cytotrace_score >= score_threshold, "High", "Low") # 可视化 DimPlot(seurat_tumor, group.by = "csc_potential", cols = c("grey", "red")) + ggtitle("Putative Cancer Stem-like Cells (CytoTRACE Top 20%)") # 差异表达分析，寻找高潜能细胞的标志物 Idents(seurat_tumor) <- seurat_tumor$csc_potential markers <- FindMarkers(seurat_tumor, ident.1 = "High", ident.2 = "Low", min.pct = 0.25) head(markers, 10)

5.2 场景二：整合分析细胞分化与基因调控网络

CytoTRACE分数可以作为一个关键的连续表型，用于关联分析基因表达动态和调控网络变化。例如，我们可以进行“伪时间”差异表达分析和趋势聚类。

library(tradeSeq) # 用于沿着轨迹的差异表达分析 # 假设我们已经有了Slingshot推断的伪时间曲线 `slingPseudotime(sce)[,1]`，或者直接使用CytoTRACE分数作为伪时间 pseudotime_vector <- seurat_obj$cytotrace_score # 使用CytoTRACE分数 # 为了使用tradeSeq，我们需要一个细胞×伪时间的矩阵，这里简单处理 # 更严谨的做法是用Slingshot的曲线 # 此处演示思路：我们可以用CytoTRACE分数对基因表达做相关性分析来替代复杂的轨迹模型 # 方法：计算每个基因的表达与CytoTRACE分数的相关性 cor_results <- apply(GetAssayData(seurat_obj, slot = "data")[VariableFeatures(seurat_obj)[1:2000], ], 1, function(gene_exp) cor(gene_exp, pseudotime_vector, method = "spearman")) # 得到与分化潜能最正相关和最负相关的基因 top_pos_genes <- names(sort(cor_results, decreasing = TRUE))[1:50] top_neg_genes <- names(sort(cor_results, decreasing = FALSE))[1:50] # 对这些基因进行热图可视化，按CytoTRACE分数排序细胞 DoHeatmap(subset(seurat_obj, downsample = 500), features = c(top_pos_genes[1:10], top_neg_genes[1:10]), group.by = "celltype", cells = order(seurat_obj$cytotrace_score)) + scale_fill_gradient2(low = "blue", high = "red", mid = "white")

这张热图可以直观展示随着分化（分数降低），哪些基因模块被逐渐抑制（蓝色），哪些被逐渐激活（红色）。

5.3 场景三：评估重编程或去分化过程的效率

在细胞重编程（如成纤维细胞诱导为iPS细胞）实验中，CytoTRACE可以用来定量评估重编程效率。通过对比重编程不同时间点的细胞样本，观察整体CytoTRACE分数的提升情况，可以量化去分化的进程。

分析流程：

1. 整合重编程第0天（体细胞）、第7天、第14天、第21天和完全重编程的iPS细胞的单细胞数据。
2. 运行CytoTRACE（注意需先校正批次效应）。
3. 比较各时间点细胞群体的平均CytoTRACE分数或分数分布。
4. 可以观察到分数分布从低分（体细胞）向高分（iPS细胞）移动的过程，并且中间时间点会出现双峰分布，代表部分重编程和完全重编程细胞的混合状态。

这个应用提供了一个超越简单标记基因表达的、全局性的量化指标来监测复杂的细胞身份转换过程。

通过以上从原理、实战到高级应用的拆解，相信你已经对CytoTRACE这个工具有了立体的认识。它不是一个万能的轨迹推断神器，但其简洁的生物学逻辑和快速的运算能力，使其成为单细胞数据分析武器库中一把非常独特的“快刀”。在项目初期，用它来快速扫描数据中的分化潜能梯度，往往能带来意想不到的发现，为后续更精细的分析指明方向。记住，任何计算工具的结果都需要结合坚实的生物学知识进行审慎的解读和验证。