1. 从零开始：为什么我们要给CD8T细胞做“功能评分”？
如果你刚刚接触单细胞测序数据分析，看到“功能基因集评分”这个词可能会有点懵。这到底是在干什么？简单来说，这就像给一群学生（CD8T细胞）做一次“综合能力测评”。我们手头有一份单细胞测序数据，里面记录了成千上万个CD8T细胞各自表达了哪些基因。但基因成千上万，我们不可能一个个去看。这时候，前辈们通过大量研究，已经总结出了一些“考题集”——也就是功能基因集。比如，一套关于“细胞毒性”的考题，里面包含了穿孔素、颗粒酶等基因；另一套关于“耗竭”的考题，则包含了PD-1、TIM-3等基因。
给细胞做功能评分，就是拿着这些“考题集”，去计算每个细胞在这些“综合科目”上的得分。得分高的细胞，说明它在这个功能上比较活跃。比如，一个细胞在“细胞毒性”基因集上得分很高，那它很可能就是一个战斗力强悍的“杀手T细胞”。这么做的巨大好处是，我们能从更高维度、更功能化的视角来理解细胞群体。不再局限于看单个基因的表达高低，而是看细胞整体的功能状态。这对于在肿瘤免疫、自身免疫病等领域，深入理解CD8T细胞亚群的特性和作用，简直是神器。
我刚开始做分析的时候，总喜欢盯着差异表达基因列表看，看得头晕眼花还得不出整体结论。后来用了基因集评分方法，一下子豁然开朗，能清晰地看到不同病人或不同治疗组之间，细胞功能状态的整体差异。今天，我就手把手带你复现一篇顶刊里的经典分析图，并把这个流程掰开揉碎了讲清楚。你会发现，从数据准备、评分计算到可视化出图，每一步都有需要注意的细节和“坑”，我会把我踩过的那些都告诉你。
2. 实战准备：搭建你的分析环境与理解数据
工欲善其事，必先利其器。在开始敲代码之前，我们需要把环境和数据都准备好。这个过程看似简单，但却是后续一切分析的基础，这里出问题，后面全白搭。
2.1 核心R包安装与加载
我们的整个分析流程将完全在R语言环境中完成，主要依赖Seurat这个单细胞分析的金牌工具包。我强烈建议你使用RStudio来操作，体验会好很多。首先，我们需要安装并加载必要的R包。如果你还没安装，可以把下面这行代码的注释去掉运行。
# 安装核心包（如果尚未安装）
# install.packages(c("Seurat", "tidyverse", "pheatmap", "scales"))
安装好后，在每个分析脚本的开头，我们都需要加载它们：
# 加载必要的库
library(Seurat) # 单细胞分析核心
library(tidyverse) # 数据操作和绘图，包含ggplot2、dplyr等
library(pheatmap) # 绘制漂亮的热图
library(scales) # 提供数据缩放函数，如rescale
这里有个小经验：tidyverse是一个“包集合”，加载它之后，你就可以用非常优雅和一致的语法（比如管道操作符 %>%）来处理数据了，效率能提升不少。pheatmap是我个人最喜欢的热图包，因为它默认出来的图就很好看，自定义选项也丰富。
2.2 理解并加载你的单细胞数据
接下来是加载数据。原始文章里用的是load(“CD8_Obj.Rdata”)，这里隐藏了一个关键前提：你已经有了一个处理好的、只包含CD8T细胞的Seurat对象。什么叫“处理好的”？通常意味着你已经完成了：
- 数据质控：过滤了低质量细胞和过多线粒体基因的细胞。
- 标准化与归一化。
- 降维聚类：使用PCA、t-SNE或UMAP，并且已经完成了细胞聚类。
- 细胞类型注释：并且从中提取出了CD8T细胞亚群，保存为独立的Seurat对象。
这一步是最大的分水岭。 很多新手会卡在这里，因为他们手上的数据可能还是原始的基因表达矩阵。你需要确保你的CD8_Obj已经完成了上述基础分析。你可以通过运行 CD8_Obj 来查看对象的基本信息，确认细胞数量和基因数量是否符合预期。
# 假设你的CD8T细胞Seurat对象已经保存为Rdata文件
load(“CD8_Obj.Rdata”)
# 查看对象摘要
CD8_Obj
加载成功后，我们先做一点准备工作：确保细胞聚类信息seurat_clusters被正确地设置为因子（factor），并设为当前对象的标识（Ident）。这步很重要，因为后续我们按聚类来汇总分数时，需要它是因子类型。
# 将seurat_clusters设置为因子类型
CD8_Obj$seurat_clusters <- as.factor(CD8_Obj$seurat_clusters)
# 设置当前细胞标识为聚类ID
Idents(CD8_Obj) <- CD8_Obj$seurat_clusters
3. 核心步骤一：读取与整理功能基因集
有了细胞数据，我们还需要“考题集”——功能基因集。原始文章的作者很贴心，在补充材料里提供了他们使用的19个功能基因集。我已经帮你整理成了一个CSV文件CD8-markers.csv。这个文件是“宽格式”的，每一列代表一个功能，每一行是该功能下的基因，列名就是功能名称。
3.1 读取CSV文件并处理格式
直接读取这个文件时要注意，有时CSV第一行可能是备注，我们需要跳过。用read_csv函数（来自tidyverse）可以轻松处理。
# 设置输出图片的目录，养成好习惯
figurePath <- "./figures"
dir.create(figurePath, showWarnings = FALSE, recursive = TRUE)
# 读取标记基因CSV文件，跳过可能存在的标题行
CD8_markers <- read_csv("CD8-markers.csv", skip = 1)
现在CD8_markers是一个数据框（DataFrame）。但Seurat的评分函数AddModuleScore需要的是一个列表（list），列表的每个元素是一个字符向量（即基因名），元素的名字就是功能名称。所以我们需要进行格式转换。这里有个坑：基因集名称里有时会有“/”这样的字符，这在后续做图或命名时可能会被误解为路径，所以我们最好把它替换掉，比如换成“:”。
# 初始化一个空列表，用于存放转换后的基因集
marker.list <- list()
# 循环处理CSV的每一列（即每一个功能基因集）
for(col_name in colnames(CD8_markers)) {
# 从该列中提取所有非空且非NA的基因名
genes <- CD8_markers[[col_name]]
genes <- genes[!is.na(genes) & genes != ""]
# 清理功能名称：将“/”替换为“:”
clean_name <- gsub("/", ":", col_name)
# 如果这个基因集里有基因，就把它添加到列表里
if(length(genes) > 0) {
marker.list[[clean_name]] <- genes
}
}
# 查看一下我们转换好的基因集列表的前几个
head(marker.list, 3)
运行完这段代码，marker.list就应该是一个包含19个元素的命名列表了。你可以检查一下每个元素里基因的数量，确保没有空的基因集被错误地导入。
4. 核心步骤二：使用AddModuleScore进行基因集评分
这是整个流程的技术核心。我们将使用Seurat包中的AddModuleScore函数来计算每个细胞在每个功能基因集上的“得分”。这个函数的原理是什么？我打个比方：它并不是简单地把基因集中所有基因的表达值加起来。那样的话，高表达基因多的细胞永远得分高，不公平。
AddModuleScore的做法更聪明：它会为每个功能基因集，从整个表达矩阵中随机抽取一组数量相同、表达水平分布类似的“控制基因”（由参数ctrl指定，比如5组）。然后用功能基因集的平均表达值，减去这些控制基因集的平均表达值，得到最终的分数。这样就消除了基因表达量基线的影响，分数更能反映该功能基因集相对于细胞自身背景的“特异性活跃程度”。
# 进行基因集评分，ctrl参数控制用于背景对比的随机基因集数量
CD8_Obj <- AddModuleScore(CD8_Obj,
features = marker.list,
ctrl = 5, # 通常5-100之间，需要根据数据量调整
name = "FunctionScore")
运行完后，新的评分结果被添加到了CD8_Obj@meta.data这个细胞的元数据表中。但是，默认的列名不太友好，是FunctionScore1， FunctionScore2……我们需要把它们换成我们熟悉的基因集名称。
# 将默认的评分列名重命名为我们基因集的实际名称
for(i in 1:length(marker.list)){
old_name <- paste0("FunctionScore", i)
new_name <- names(marker.list)[i]
# 在元数据中找到旧列名，替换为新列名
colnames(CD8_Obj@meta.data)[colnames(CD8_Obj@meta.data) == old_name] <- new_name
}
现在，你可以查看CD8_Obj@meta.data的列名，应该能看到“Naive”、“Cytotoxicity”这些功能名称已经作为新列出现了。每一行是一个细胞，每一列是这个细胞在该功能上的得分。
5. 核心步骤三：整理评分矩阵，为热图做准备
我们现在有了每个细胞的分数，但热图通常展示的是细胞群体（比如每个聚类）的平均功能状态。所以我们需要把数据“聚合”一下：计算每个细胞聚类（cluster）在每个功能基因集上的平均得分。此外，原文的热图将功能分成了四大类：分化状态、免疫功能、代谢和凋亡。我们需要按照这个逻辑来整理和排序我们的数据行。
5.1 定义功能分类与排序
首先，我们按照原文的定义，把19个功能归到四个大类里，并创建一个最终的行名顺序向量。
# 定义四个大类的具体功能（名称必须与marker.list中的完全一致）
Differentiation <- c("Naive", "Activation:Effector function", "Exhaustion")
Function <- c("TCR Signaling", "Cytotoxicity", "Cytokine:Cytokine receptor", "Chemokine:Chemokine receptor", "Senescence", "Anergy", "NFKB Signaling", "Stress response", "MAPK Signaling", "Adhesion", "IFN Response")
Metabolism <- c("Oxidative phosphorylation", "Glycolysis", "Fatty acid metabolism")
Apoptosis <- c("Pro-apoptosis", "Anti-apoptosis")
# 合并成最终的行顺序向量
MarkerNameVector <- c(Differentiation, Function, Metabolism, Apoptosis)
5.2 计算聚类平均得分矩阵
接下来，我们创建一个空矩阵，行是功能，列是细胞聚类，然后通过双层循环，计算每个功能在每个聚类里的平均得分。
# 获取唯一的聚类ID
cluster_ids <- levels(CD8_Obj$seurat_clusters) # 注意这里用levels，能保证因子顺序
# 创建一个空矩阵
FunctionScoreMatrix <- matrix(0,
ncol = length(cluster_ids),
nrow = length(MarkerNameVector))
# 设置行名和列名
rownames(FunctionScoreMatrix) <- MarkerNameVector
colnames(FunctionScoreMatrix) <- paste0("CD8_", cluster_ids) # 给列名加个前缀
# 双层循环，计算均值
for(ci in 1:ncol(FunctionScoreMatrix)) { # 循环每个聚类
for(ri in 1:nrow(FunctionScoreMatrix)) { # 循环每个功能
# 从元数据中提取所有细胞在当前功能上的得分
score_vector <- CD8_Obj@meta.data[[ MarkerNameVector[ri] ]]
# 计算当前聚类（ci）内细胞的平均得分
cluster_mean_score <- mean(score_vector[CD8_Obj$seurat_clusters == cluster_ids[ci]])
# 填入矩阵
FunctionScoreMatrix[ri, ci] <- cluster_mean_score
}
}
这里有个非常重要的步骤：缩放（Scaling）。 直接计算出来的平均得分，不同功能之间可能数值范围差异很大。为了在热图上公平地展示，我们需要对每一行（每个功能）的数据进行缩放，通常缩放到[-1, 1]或[0, 1]的区间。这里我们使用scales包的rescale函数，缩放到[-1, 1]。
# 对矩阵的每一行进行缩放（-1到1）
FunctionScoreMatrix <- t(apply(FunctionScoreMatrix, 1, rescale, to=c(-1, 1)))
缩放后，数值1代表该功能在该聚类中表达最强（相对于其他聚类），-1代表最弱。这样，热图的颜色对比就非常清晰了。
5.3 调整聚类顺序
有时候，我们不想按照聚类ID的数字顺序（0,1,2…）来展示，可能想根据其他知识（比如发育轨迹、细胞状态）来手动排序。这时可以预先定义一个顺序向量orderC。
# 假设我们想调整聚类顺序，例如想把CD8_5和CD8_2放到前面
orderC <- c(5, 2, 0, 1, 3, 4, 6, 7, 8, 9, 10) # 请根据你的实际聚类ID修改
orderC_with_prefix <- paste0("CD8_", orderC)
# 按照新的列顺序重排矩阵
FunctionScoreMatrix <- FunctionScoreMatrix[, orderC_with_prefix]
6. 核心步骤四：绘制发表级功能热图
数据矩阵准备好了，最后一步就是把它变成一张直观、美观、可以放进论文里的热图。我们使用pheatmap包，它能提供丰富的注释和美化选项。
6.1 准备热图注释信息
原文的热图左边有颜色条，表示功能属于哪个大类（分化、功能、代谢、凋亡）。我们需要为这个注释创建一个数据框。
# 创建行注释信息（功能大类）
signatureType_row <- data.frame(
Signature.type = c(
rep("Differentiation", length(Differentiation)),
rep("Function", length(Function)),
rep("Metabolism", length(Metabolism)),
rep("Apoptosis", length(Apoptosis))
)
)
# 行注释的行名必须与评分矩阵的行名一致
rownames(signatureType_row) <- MarkerNameVector
6.2 自定义颜色与断点
默认的热图颜色可能对比不够鲜明。我们可以自定义一个从蓝到白到红的渐变色，并精细控制颜色变化的断点。这能让正值（红色）和负值（蓝色）的对比非常清晰。
# 定义颜色断点：从-1到0分成10段，从0到1分成10段
my.breaks <- c(seq(-1, 0, length.out=11), seq(0.1, 1, length.out=10))
# 定义对应的颜色：-1到0用蓝色到白色，0到1用白色到红色
my.colors <- c(
colorRampPalette(colors = c("#6DCCFD", "white"))(10), # 蓝-白渐变
colorRampPalette(colors = c("white", "#FD9AA0"))(10) # 白-红渐变
)
6.3 绘制并保存热图
万事俱备，调用pheatmap函数。这里参数比较多，我解释几个关键的：
cluster_rows/cols = F：不对行和列进行聚类排序，保持我们预设的顺序。gaps_row = c(3, 14, 17)：在行方向（功能）的第3行后、第14行后、第17行后画分隔线，正好对应四个大类之间的分界。border_color = “NA”：去掉每个单元格的边框，让图看起来更干净。
pheatmap(FunctionScoreMatrix,
show_colnames = T,
show_rownames = T,
annotation_row = signatureType_row, # 添加行注释
gaps_row = c(3, 14, 17), # 在大类之间加分隔线
cluster_rows = F,
cluster_cols = F,
breaks = my.breaks, # 自定义断点
color = my.colors, # 自定义颜色
border_color = "NA", # 无边框
fontsize = 8, # 字体大小
width = 5, # 输出PDF的宽度（英寸）
height = 3.8, # 输出PDF的高度（英寸）
filename = file.path(figurePath, "CD8_FunctionScore_heatmap.pdf")) # 保存
运行完这行代码，你会在之前创建的./figures文件夹里找到一个PDF格式的热图。这张图已经非常接近发表水平了。你可以把它导入到Adobe Illustrator或Inkscape这类矢量图软件中，微调字体、图例位置，或者添加最终的图表标题，就能直接用于你的文章或报告中。
7. 避坑指南与进阶思考
走完整个流程，你应该已经成功复现了分析。但在实际项目中，我遇到过不少问题，这里分享几个最常见的“坑”和解决思路。
第一个坑：评分结果全是负值或正值不明显。 这很可能是因为ctrl参数设置不合理。ctrl太小（比如1），背景噪音控制不足；太大则可能过度校正。我的经验是，对于万级别的细胞数，从5、10、20、50开始尝试，观察评分分布的变化。也可以尝试Seurat的另一函数AddModuleScore的替代方案，如AUCell或UCell包，它们用了不同的算法，有时效果更好。
第二个坑：热图看起来“糊成一团”，没有明显模式。 这可能有两个原因。一是功能基因集本身不适合你的数据，比如用的是小鼠的基因集来分析人的数据。二是细胞聚类分得不好，或者聚类分辨率不合适，导致每个聚类内部细胞异质性仍然很大。这时需要回头检查你的聚类分析，尝试不同的分辨率重新聚类。
第三个坑：如何选择或自定义基因集？ 顶刊的基因集是很好的起点，但未必完全适用于你的研究场景。你可以从MSigDB、GO、KEGG等权威数据库下载通路基因集，也可以从与你疾病相关的关键文献中收集。自己构建时，切记基因集不能太大（一般建议15-500个基因），太大容易失去特异性；也不能太小，太小容易受噪声影响。
进阶思考：评分之后做什么？ 功能评分不是一个终点，而是一个强大的新起点。你可以：
- 比较组间差异：将病人分为应答组与非应答组，比较两组间CD8T细胞在“耗竭”评分上是否有显著差异。
- 关联临床指标：计算每个样本的“细胞毒性”平均分，与患者的生存时间或肿瘤负荷做相关性分析。
- 轨迹分析：将功能评分作为细胞的一个属性，投射到拟时序轨迹上，观察功能如何随时间动态变化。
我自己的一个项目中，就是通过功能评分发现，对免疫治疗无效的患者，其肿瘤内CD8T细胞不仅“耗竭”评分高，关键的是“细胞毒性”与“耗竭”评分呈现强烈的共高表达，这是一种功能失调的状态，而不仅仅是单纯的衰竭。这个发现仅看差异基因是很难直观得出的。所以，熟练掌握这套流程，真正理解其背后的生物信息学逻辑，能为你挖掘单细胞数据打开一扇新的大门。
258
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
