单细胞转录组分析做到最后你手里是不是总有一堆“基因集”可能是从文献里找到的疾病特征基因可能是自己差异分析得到的通路基因或者是像MSigDB这样的标准数据库。看着这些基因列表一个最直接的问题就来了在我这个单细胞数据里哪些细胞高表达了这些基因这个问题听起来简单但实际操作起来远不止是“把基因表达量加一加”那么简单。细胞间的表达量差异巨大基因数量不一直接加和会严重偏向高表达基因。更关键的是单细胞数据的稀疏性大量零值让传统的富集分析方法如GSEA常常水土不服。这就是AUCell算法要解决的核心痛点。它不是一个新潮的概念而是单细胞分析中一个极其实用、几乎成为标准流程的评分工具。很多人以为它只是个“打分函数”但它的精妙之处在于它通过计算“曲线下面积”AUC巧妙地绕过了表达量绝对值的干扰专注于基因在单个细胞内的相对排序从而对稀疏、嘈杂的单细胞数据表现出惊人的鲁棒性。如果你正在做单细胞项目并且想知道“我的干细胞特征、免疫细胞特征、代谢通路活性在哪些细胞里更强”那么这篇文章就是为你准备的。我将不仅带你一步步跑通AUCell更会深入解释为什么是AUCell而不是其他方法它的结果到底该怎么解读那些看似简单的参数背后藏着哪些“坑”以及如何将AUCell评分巧妙地融入你的下游分析如聚类、拟时序、细胞通讯为你的故事提供强有力的数据支撑。1. AUCell 要解决的远不止一个“打分”问题在深入代码之前我们必须先理解AUCell背后的设计哲学。这决定了你能否正确使用和解读它的结果。想象一个场景你有一个包含100个基因的“干细胞特征基因集”。细胞A高表达了其中5个关键转录因子表达量都很高细胞B则低表达了另外95个看家基因表达量普遍较低但都有。如果简单求和或取均值细胞B的分数可能会更高但这显然违背了生物学直觉——我们更关心那些高表达核心特征基因的细胞。AUCell的聪明之处在于它完全摒弃了比较细胞间的绝对表达量。相反它只关注单个细胞内部。对于每一个细胞AUCell做这样一件事构建排名列表将这个细胞内所有基因的表达量从高到低排序。定位基因集在这个细胞的私人排名列表中找到你的目标基因集比如那100个干细胞基因都排在什么位置。计算AUC值计算这些基因的排名所构成的“曲线下面积”。你可以把它理解为在这个细胞中目标基因集里的基因是否更多地聚集在表达量最高的头部区域。这个AUC值就是AUCell评分。它的值在0到1之间。分数越接近1说明在这个细胞里你关心的那组基因整体表达排名越靠前分数接近0则说明这些基因在这个细胞里表达排名靠后。这种方法带来了几个关键优势对稀疏性免疫零值多没关系只要目标基因的相对排名高就行。对表达量尺度不敏感不同细胞类型、不同测序深度导致的表达量整体差异不会直接影响评分。专注于“模式”而非“总量”更符合生物学上“特征激活”的概念。所以当你拿到AUCell的评分矩阵时你看到的不是表达量的加和而是每个细胞对于特定基因集的“响应强度”或“特征活性”。这是构建细胞功能状态图谱的关键一步。2. 核心概念AUCell 输入与输出的本质理解输入和输出是正确使用任何工具的前提。AUCell的核心流程可以概括为下图所示graph TD A[输入: 单细胞表达矩阵] -- B[输入: 目标基因集列表] B -- C{AUCell 核心计算}; C -- D[为每个细胞构建基因表达量排名]; D -- E[计算目标基因集在该排名中的 AUC 值]; E -- F[输出: 细胞 x 基因集评分矩阵]; F -- G{下游分析}; G -- H[可视化: 评分分布/UMAP投射]; G -- I[分析: 识别高分细胞群]; G -- J[整合: 驱动聚类/拟时序分析];2.1 输入是什么表达矩阵 (Expression Matrix) 必须是规整的数值矩阵行是基因列是细胞。通常是经过标准预处理质控、归一化、对数转换后的数据。AUCell 接受多种格式如matrix、dgCMatrix(稀疏矩阵) 或SingleCellExperiment/Seurat对象。基因集 (Gene Sets) 一个列表对象。例如在R中是一个list每个元素是一个字符向量包含属于同一个基因集的基因名。基因名必须与表达矩阵的行名基因名匹配。# 示例创建两个基因集 my_gene_sets - list( “Stemness_Signature” c(“SOX2”, “NANOG”, “POU5F1”, “MYC”, “KLF4”), “Inflammation_Response” c(“IL6”, “TNF”, “NFKB1”, “STAT3”, “JUN”) )2.2 输出是什么一个矩阵或者整合到单细胞对象中的一个新维度。矩阵形式 行是细胞列是基因集。每个值就是该细胞对该基因集的AUCell评分0-1。整合到对象中 以Seurat为例评分会被添加到Seurat对象的meta.data中作为细胞的一列新注释方便后续调用和绘图。2.3 AUCell 评分意味着什么需要牢记AUCell评分是细胞内的相对度量只能用于比较同一基因集在不同细胞间的活性差异。分数0.5 是一个随机期望值。意味着目标基因在该细胞中的排名分布和随机抽取差不多。分数 0.5 意味着这些基因在该细胞中整体排名高于随机预期表明该特征可能有活性。分数越高 特征活性越强的可能性越大。重要提醒 不要直接比较不同基因集之间的分数大小因为不同基因集的大小基因数量不同其AUC值的分布基线也不同。基因集越大其AUC值的随机分布越集中。比较“干细胞评分0.8”和“凋亡评分0.6”哪个更高是没有意义的。正确的做法是分别看每个基因集评分的细胞间分布识别高评分细胞群。3. 环境准备在 R 中搭建 AUCell 分析流程我们将使用 R 语言进行演示因为其生态在单细胞分析中最为成熟。确保你已安装以下包# 安装必备包如果尚未安装 if (!requireNamespace(“BiocManager”, quietly TRUE)) install.packages(“BiocManager”) BiocManager::install(“AUCell”) # AUCell 核心算法包 BiocManager::install(“GSEABase”) # 用于处理基因集特别是读取 .gmt 文件 install.packages(“Seurat”) # 单细胞分析主流平台 install.packages(“ggplot2”) # 绘图 install.packages(“dplyr”) # 数据操作 # 加载库 library(AUCell) library(GSEABase) library(Seurat) library(ggplot2) library(dplyr)版本建议 AUCell 版本建议在 1.20.0 以上Seurat 建议使用 4.x 或 5.x 版本。运行packageVersion(“AUCell”)可查看当前版本。数据准备 你需要一个预处理好的单细胞对象。这里假设你有一个名为seurat_obj的 Seurat 对象已经完成了NormalizeData(),FindVariableFeatures(),ScaleData(),RunPCA(),RunUMAP()等标准步骤。4. 核心流程四步走从数据到评分AUCell 的分析流程非常清晰主要分为四步。4.1 第一步提取表达矩阵并构建排名AUCell 需要输入一个数值矩阵。我们从 Seurat 对象中提取经过归一化如LogNormalize的表达矩阵。# 提取表达矩阵。使用 GetAssayData 并指定 slot“data” 获取归一化数据。 expr_matrix - GetAssayData(seurat_obj, slot “data”) # 得到一个稀疏矩阵 dgCMatrix # expr_matrix 是一个基因 x 细胞的矩阵 # 查看矩阵维度 dim(expr_matrix) # [1] 20000 5000 # 例如20000个基因5000个细胞 # 构建基因表达排名。这是AUCell最耗计算的一步但每个排名只需计算一次。 # 它会为每个细胞计算所有基因的表达量排名。 cells_rankings - AUCell_buildRankings( exprMatrix expr_matrix, nCores 4, # 根据你的电脑核心数设置可以加速计算 plotStats TRUE # 绘制排名分布的统计图检查是否正常 )运行AUCell_buildRankings后会生成一个cells_rankings对象。plotStatsTRUE会生成一张图展示每个细胞中基因排名的分布通常应该看到一条平滑的曲线如果出现异常如大量平局排名可能需要检查数据。4.2 第二步准备目标基因集基因集可以来自多种渠道。这里演示三种常见方式方式一手动创建列表# 如前所述手动定义基因集列表 my_signatures - list( “Cell_Cycle_G2M” c(“TOP2A”, “MKI67”, “BIRC5”, “CCNB2”, “UBE2C”), “Hypoxia” c(“VEGFA”, “BNIP3”, “SLC2A1”, “PGK1”, “CA9”), “My_Own_Pathway” c(“GENE1”, “GENE2”, “GENE3”) )方式二从 MSigDB 的 .gmt 文件读取# 假设你从 MSigDB官网下载了文件 c2.cp.kegg.v2023.1.Hs.symbols.gmt gmt_file - “path/to/your/c2.cp.kegg.v2023.1.Hs.symbols.gmt” gene_sets - getGmt(gmt_file) # 转换为 AUCell 需要的列表格式 gene_sets_list - geneIds(gene_sets) # gene_sets_list 就是一个包含多个通路基因集的列表方式三从已有分析结果获取# 例如从你的差异表达分析结果中提取某个细胞簇的上调基因作为特征 cluster_markers - FindMarkers(seurat_obj, ident.1 “Cluster1”, min.pct 0.25) top_genes - rownames(cluster_markers)[1:50] # 取前50个差异基因 my_signatures$“Cluster1_Signature” - top_genes关键检查 必须确保基因集中的基因名与表达矩阵的行名完全匹配大小写、符号一致。# 检查并过滤掉表达矩阵中不存在的基因 for (set_name in names(my_signatures)) { genes_in_matrix - my_signatures[[set_name]] %in% rownames(expr_matrix) cat(set_name, “:”, sum(genes_in_matrix), “/”, length(genes_in_matrix), “genes found\n”) # 可以选择只保留存在的基因 my_signatures[[set_name]] - my_signatures[[set_name]][genes_in_matrix] }4.3 第三步计算 AUC 值这是核心计算步骤使用上一步构建的排名和准备好的基因集。# 计算每个细胞对每个基因集的AUC值 cells_AUC - AUCell_calcAUC( geneSets my_signatures, rankings cells_rankings, nCores 4, # 并行计算 aucMaxRank ceiling(0.05 * nrow(cells_rankings)) # 最重要的参数之一见下文解释 )参数aucMaxRank详解 这是AUCell最关键的参数没有之一。它决定了在计算AUC时只考虑每个细胞中表达量排名前aucMaxRank的基因。为什么需要它在单细胞数据中大量基因表达为0或极低如果考虑所有基因排名AUC值会主要由这些低表达基因的尾部排名决定使得评分区分度下降。限制只关注高表达区域能放大有生物学意义的信号。如何设置通常设置为所有基因数的前5% (0.05 * nrow(rankings))。这是一个经验值你可以尝试2%-10%。可以通过AUCell_exploreThresholds函数来辅助评估。4.4 第四步提取结果并整合计算完成后我们从cells_AUC对象中提取评分矩阵。# 提取评分矩阵 auc_matrix - getAUC(cells_AUC) dim(auc_matrix) # 矩阵维度为基因集数量 x 细胞数量 # 转置一下变成细胞 x 基因集更方便后续操作 auc_matrix - t(auc_matrix)整合到 Seurat 对象 这是推荐的做法便于统一管理。# 将评分矩阵添加到 Seurat 对象的 metadata 中 # 确保细胞顺序一致AUCell默认可能按细胞名排序Seurat对象也是 cell_names - colnames(seurat_obj) auc_matrix - auc_matrix[cell_names, ] # 按Seurat对象细胞顺序重排 for (sig_name in colnames(auc_matrix)) { seurat_obj[[sig_name]] - auc_matrix[, sig_name] } # 检查是否添加成功 head(seurat_objmeta.data)5. 结果可视化与解读让评分“说话”算出分数只是开始如何解读和展示才是关键。5.1 基础可视化在 UMAP/tSNE 上着色最直观的方式是将AUCell评分映射到降维聚类图上。# 绘制单个基因集评分在UMAP上的分布 FeaturePlot(seurat_obj, features “Cell_Cycle_G2M”, # 替换为你的基因集名称 cols c(“lightgrey”, “blue”)) # 颜色从低到高 # 使用更专业的颜色梯度并添加标题 FeaturePlot(seurat_obj, features “Hypoxia”, cols viridis::viridis(10), # 使用viridis色系对色盲友好 order TRUE) # orderTRUE 将高分细胞画在最上层避免被遮盖 ggtitle(“Hypoxia Signature Activity”)5.2 评分分布分析识别“阳性”细胞AUCell评分是连续的但有时我们需要定义一个阈值来区分“高活性”和“低活性”细胞。# 方法1使用AUCell内置函数自动探索阈值 cells_assignment - AUCell_exploreThresholds(cells_AUC, plotHist TRUE, assignCells TRUE) # 这个函数会为每个基因集生成一个直方图并建议一个阈值如“跳出双峰分布的谷底”。 # 结果存储在 cells_assignment 列表中。 # 查看为“Cell_Cycle_G2M”基因集分配的细胞 g2m_cells - cells_assignment$“Cell_Cycle_G2M”$assignment # 这是一个细胞名向量代表被认定为该特征“阳性”的细胞。 # 方法2手动设定阈值例如取前20%的细胞 scores - seurat_obj$“Cell_Cycle_G2M” threshold - quantile(scores, probs 0.8) # 取80%分位数即前20% seurat_obj$“G2M_high” - scores threshold # 在UMAP上查看两类细胞 DimPlot(seurat_obj, group.by “G2M_high”, cols c(“grey”, “red”)) ggtitle(“Top 20% Cells for G2M Signature”)5.3 多基因集比较与热图展示当有多个基因集时可以观察它们在不同细胞群中的活性模式。# 提取所有细胞的评分 auc_scores - FetchData(seurat_obj, vars c(“Cell_Cycle_G2M”, “Hypoxia”, “My_Own_Pathway”)) # 按细胞聚类取平均分观察不同簇的特征活性 avg_scores - seurat_objmeta.data %% group_by(seurat_clusters) %% # 假设你的聚类信息列名为‘seurat_clusters’ summarise(across(c(“Cell_Cycle_G2M”, “Hypoxia”, “My_Own_Pathway”), mean)) # 绘制热图 library(pheatmap) score_mat - as.matrix(avg_scores[, -1]) # 去掉第一列簇名 rownames(score_mat) - avg_scores$seurat_clusters pheatmap(score_mat, cluster_rows TRUE, cluster_cols TRUE, scale “column”, # 按列基因集标准化便于比较不同基因集在各簇的相对高低 color colorRampPalette(c(“navy”, “white”, “firebrick3”))(50), main “Average AUCell Score per Cluster”)这张热图能清晰告诉你哪个细胞簇更富集细胞周期特征哪个簇更处于缺氧状态。6. 进阶应用将 AUCell 评分融入下游分析AUCell评分不仅仅是用来画图的它可以作为强有力的特征输入到下游分析中。6.1 驱动聚类分析传统的聚类基于所有基因的表达变化。你可以尝试加入AUCell评分作为补充特征让聚类结果更倾向于反映功能状态。# 将AUCell评分矩阵作为新的“assay”加入Seurat对象 # 首先确保评分矩阵的细胞顺序与对象一致 auc_matrix - t(getAUC(cells_AUC))[, colnames(seurat_obj)] # 细胞 x 基因集 # 创建一个新的Assay seurat_obj[[“AUC”]] - CreateAssayObject(data t(auc_matrix)) # Assay要求基因 x 细胞 # 设置默认assay为AUC并运行PCA和UMAP DefaultAssay(seurat_obj) - “AUC” seurat_obj - ScaleData(seurat_obj) seurat_obj - RunPCA(seurat_obj, features rownames(seurat_obj[[“AUC”]])) seurat_obj - RunUMAP(seurat_obj, dims 1:10) # 基于功能评分进行聚类 seurat_obj - FindNeighbors(seurat_obj, dims 1:10) seurat_obj - FindClusters(seurat_obj, resolution 0.5) # 可视化基于功能特征的聚类 DimPlot(seurat_obj, reduction “umap”, label TRUE)比较基于基因表达的聚类和基于功能评分的聚类可以发现哪些功能模块定义了新的细胞亚群。6.2 拟时序分析中的功能动力学如果你在使用Monocle3或Slingshot做拟时序分析可以将AUCell评分作为细胞的状态特征观察其沿轨迹的变化。# 假设你已经有了一个cds对象Monocle3 library(monocle3) # 将AUCell评分添加到cds的colData中 cds_col_data - colData(cds) for (sig in colnames(auc_matrix)) { cds_col_data[[sig]] - auc_matrix[rownames(cds_col_data), sig] } colData(cds) - cds_col_data # 绘制评分沿轨迹的平滑曲线 plot_genes_violin(cds, group_cells_by “pseudotime_bin”, # 将拟时序分箱 genes c(“Stemness_Signature”, “Differentiation_Signature”), ncol 2)6.3 细胞-细胞通讯分析中的配体-受体活性在分析细胞通讯如CellChat时你不仅关心配体受体对的表达更关心其下游通路的活性。AUCell可以在这里大显身手。# 假设你有一个从CellChat或其他资源获取的“NF-kB通路靶基因集” nfkb_target_genes - c(“IL6”, “TNF”, “CXCL8”, “ICAM1”, …) # 你的基因列表 # 为每个细胞计算NF-kB通路活性评分 # … (使用AUCell计算步骤同上) … # 然后在分析细胞通讯时你可以问 # “在配体高表达的发送细胞中其NF-kB通路活性是否也更高” # “在受体高表达的接收细胞中NF-kB通路活性是否被激活” # 这可以通过相关性分析或分组比较来实现。7. 常见问题、陷阱与排查指南AUCell用起来简单但想用对、用好必须避开以下几个常见的“坑”。问题现象可能原因排查方式解决方案所有细胞的评分都差不多没有区分度1.aucMaxRank参数设置过大如用了默认值。2. 基因集质量差基因不特异或表达极低。3. 输入的表达矩阵未正确归一化如还是counts。1. 检查aucMaxRank值尝试设置为基因总数的1%-10%。2. 用VlnPlot或FeaturePlot查看基因集中部分基因的表达情况。3. 确认expr_matrix来自slot“data”。1. 调整aucMaxRank并重新计算。2. 重新筛选或寻找更特异的基因集。3. 确保使用归一化后的数据。评分出现大量 NA 或计算错误1. 基因集中的基因在表达矩阵中一个都没找到。2. 基因名不匹配大小写、符号版本。1. 运行第4.2步的“关键检查”代码。2. 使用rownames(expr_matrix)[1:10]和my_signatures[[1]][1:10]对比。1. 修正基因名或使用基因ID转换工具。2. 确保使用一致的基因标识符如Symbol或Ensembl ID。运行AUCell_buildRankings内存不足或极慢细胞数或基因数过多如 10万细胞。监控内存使用。对于超大矩阵排名构建是内存和计算密集型步骤。1. 增加nCores并行。2.对细胞进行降采样先在小样本上测试流程和参数。3. 考虑使用计算集群。不同基因集的评分分布差异巨大这是正常现象。基因集大小直接影响AUC值的理论分布。大基因集AUC值更集中小基因集更分散。使用boxplot或hist查看不同基因集评分的分布。不要直接比较绝对值关注每个基因集内部细胞间的相对高低。使用分位数或Z-score标准化后再比较。我的“阳性”细胞在图上看起来不连续、很散1. 阈值设置不合理太高或太低。2. 该功能特征本身就是跨细胞类型的、或呈梯度变化。3. 降维图UMAP本身不能完美分离所有功能状态。1. 用AUCell_exploreThresholds检查阈值建议的直方图。2. 观察评分在已知细胞类型注释上的分布VlnPlot。1. 调整阈值或使用连续的评分进行分析如相关性分析。2. 结合其他证据如关键标记基因表达综合判断。8. 最佳实践与工程化建议要让AUCell分析真正成为可重复、可信赖的研究的一部分你需要遵循以下实践基因集质量控制是重中之重来源可靠优先使用权威数据库MSigDB, GO, KEGG或高质量文献中验证过的基因集。大小适中避免过小5个基因结果不稳定或过大500个基因信号可能被稀释的基因集。10-200个基因是较理想的区间。特异性确保基因集与你要研究的问题高度相关。使用从同一数据集的差异分析中得到的基因集时要格外小心避免循环论证。参数aucMaxRank需要理性选择不要盲目使用默认值。用AUCell_exploreThresholds函数观察不同aucMaxRank下评分分布的变化。对于关注强表达信号的分析如核心转录因子可以用更小的aucMaxRank如前1%。对于关注广泛、微弱信号的分析如某些代谢通路可以适当调大。在同一项研究中对所有基因集使用相同的aucMaxRank参数以保证可比性。结果解读必须结合生物学背景AUCell高分提示特征活性高但并非绝对证据。必须与已知的细胞类型标记基因表达进行交叉验证。对于新发现的“高评分细胞群”要回到原始表达数据检查是哪些基因驱动了高分。考虑使用多种评分算法如AddModuleScore, UCell, singscore进行交叉验证特别是对于关键结论。代码的可重复性与文档在脚本开头明确记录AUCell包版本、R版本和所有参数尤其是aucMaxRank。保存关键的中间结果如cells_rankings对象和最终的评分矩阵避免重复计算。为每个基因集标注清晰的名称和来源。融入分析流程将AUCell评分计算封装成函数方便在多个数据集或不同预处理条件下重复使用。考虑将评分结果作为标准输出的一部分整合到你的单细胞分析Pipeline中。单细胞数据分析的魅力在于从海量数据中提炼生物学故事。AUCell提供了一把精准的尺子让你能够定量地衡量每个细胞的“功能状态”。它不再让你停留在“这个基因高不高”的层面而是让你能回答“这个细胞群体是否具备某种功能特性”这样的高阶问题。从理解其基于排名的核心思想开始到谨慎地准备基因集和设置参数再到将评分结果可视化并与下游分析深度整合每一步都需要清晰的生物学问题驱动和技术细节把控。当你熟练运用AUCell后你会发现它不仅是流程中的一个步骤更是你探索单细胞世界功能景观的一双眼睛。 
长河君 · 2026/9/4 7:49:58 · 阅读 1.2k 单细胞转录组分析：AUCell算法原理、应用与实战指南
单细胞转录组分析做到最后你手里是不是总有一堆“基因集”可能是从文献里找到的疾病特征基因可能是自己差异分析得到的通路基因或者是像MSigDB这样的标准数据库。看着这些基因列表一个最直接的问题就来了在我这个单细胞数据里哪些细胞高表达了这些基因这个问题听起来简单但实际操作起来远不止是“把基因表达量加一加”那么简单。细胞间的表达量差异巨大基因数量不一直接加和会严重偏向高表达基因。更关键的是单细胞数据的稀疏性大量零值让传统的富集分析方法如GSEA常常水土不服。这就是AUCell算法要解决的核心痛点。它不是一个新潮的概念而是单细胞分析中一个极其实用、几乎成为标准流程的评分工具。很多人以为它只是个“打分函数”但它的精妙之处在于它通过计算“曲线下面积”AUC巧妙地绕过了表达量绝对值的干扰专注于基因在单个细胞内的相对排序从而对稀疏、嘈杂的单细胞数据表现出惊人的鲁棒性。如果你正在做单细胞项目并且想知道“我的干细胞特征、免疫细胞特征、代谢通路活性在哪些细胞里更强”那么这篇文章就是为你准备的。我将不仅带你一步步跑通AUCell更会深入解释为什么是AUCell而不是其他方法它的结果到底该怎么解读那些看似简单的参数背后藏着哪些“坑”以及如何将AUCell评分巧妙地融入你的下游分析如聚类、拟时序、细胞通讯为你的故事提供强有力的数据支撑。1. AUCell 要解决的远不止一个“打分”问题在深入代码之前我们必须先理解AUCell背后的设计哲学。这决定了你能否正确使用和解读它的结果。想象一个场景你有一个包含100个基因的“干细胞特征基因集”。细胞A高表达了其中5个关键转录因子表达量都很高细胞B则低表达了另外95个看家基因表达量普遍较低但都有。如果简单求和或取均值细胞B的分数可能会更高但这显然违背了生物学直觉——我们更关心那些高表达核心特征基因的细胞。AUCell的聪明之处在于它完全摒弃了比较细胞间的绝对表达量。相反它只关注单个细胞内部。对于每一个细胞AUCell做这样一件事构建排名列表将这个细胞内所有基因的表达量从高到低排序。定位基因集在这个细胞的私人排名列表中找到你的目标基因集比如那100个干细胞基因都排在什么位置。计算AUC值计算这些基因的排名所构成的“曲线下面积”。你可以把它理解为在这个细胞中目标基因集里的基因是否更多地聚集在表达量最高的头部区域。这个AUC值就是AUCell评分。它的值在0到1之间。分数越接近1说明在这个细胞里你关心的那组基因整体表达排名越靠前分数接近0则说明这些基因在这个细胞里表达排名靠后。这种方法带来了几个关键优势对稀疏性免疫零值多没关系只要目标基因的相对排名高就行。对表达量尺度不敏感不同细胞类型、不同测序深度导致的表达量整体差异不会直接影响评分。专注于“模式”而非“总量”更符合生物学上“特征激活”的概念。所以当你拿到AUCell的评分矩阵时你看到的不是表达量的加和而是每个细胞对于特定基因集的“响应强度”或“特征活性”。这是构建细胞功能状态图谱的关键一步。2. 核心概念AUCell 输入与输出的本质理解输入和输出是正确使用任何工具的前提。AUCell的核心流程可以概括为下图所示graph TD A[输入: 单细胞表达矩阵] -- B[输入: 目标基因集列表] B -- C{AUCell 核心计算}; C -- D[为每个细胞构建基因表达量排名]; D -- E[计算目标基因集在该排名中的 AUC 值]; E -- F[输出: 细胞 x 基因集评分矩阵]; F -- G{下游分析}; G -- H[可视化: 评分分布/UMAP投射]; G -- I[分析: 识别高分细胞群]; G -- J[整合: 驱动聚类/拟时序分析];2.1 输入是什么表达矩阵 (Expression Matrix) 必须是规整的数值矩阵行是基因列是细胞。通常是经过标准预处理质控、归一化、对数转换后的数据。AUCell 接受多种格式如matrix、dgCMatrix(稀疏矩阵) 或SingleCellExperiment/Seurat对象。基因集 (Gene Sets) 一个列表对象。例如在R中是一个list每个元素是一个字符向量包含属于同一个基因集的基因名。基因名必须与表达矩阵的行名基因名匹配。# 示例创建两个基因集 my_gene_sets - list( “Stemness_Signature” c(“SOX2”, “NANOG”, “POU5F1”, “MYC”, “KLF4”), “Inflammation_Response” c(“IL6”, “TNF”, “NFKB1”, “STAT3”, “JUN”) )2.2 输出是什么一个矩阵或者整合到单细胞对象中的一个新维度。矩阵形式 行是细胞列是基因集。每个值就是该细胞对该基因集的AUCell评分0-1。整合到对象中 以Seurat为例评分会被添加到Seurat对象的meta.data中作为细胞的一列新注释方便后续调用和绘图。2.3 AUCell 评分意味着什么需要牢记AUCell评分是细胞内的相对度量只能用于比较同一基因集在不同细胞间的活性差异。分数0.5 是一个随机期望值。意味着目标基因在该细胞中的排名分布和随机抽取差不多。分数 0.5 意味着这些基因在该细胞中整体排名高于随机预期表明该特征可能有活性。分数越高 特征活性越强的可能性越大。重要提醒 不要直接比较不同基因集之间的分数大小因为不同基因集的大小基因数量不同其AUC值的分布基线也不同。基因集越大其AUC值的随机分布越集中。比较“干细胞评分0.8”和“凋亡评分0.6”哪个更高是没有意义的。正确的做法是分别看每个基因集评分的细胞间分布识别高评分细胞群。3. 环境准备在 R 中搭建 AUCell 分析流程我们将使用 R 语言进行演示因为其生态在单细胞分析中最为成熟。确保你已安装以下包# 安装必备包如果尚未安装 if (!requireNamespace(“BiocManager”, quietly TRUE)) install.packages(“BiocManager”) BiocManager::install(“AUCell”) # AUCell 核心算法包 BiocManager::install(“GSEABase”) # 用于处理基因集特别是读取 .gmt 文件 install.packages(“Seurat”) # 单细胞分析主流平台 install.packages(“ggplot2”) # 绘图 install.packages(“dplyr”) # 数据操作 # 加载库 library(AUCell) library(GSEABase) library(Seurat) library(ggplot2) library(dplyr)版本建议 AUCell 版本建议在 1.20.0 以上Seurat 建议使用 4.x 或 5.x 版本。运行packageVersion(“AUCell”)可查看当前版本。数据准备 你需要一个预处理好的单细胞对象。这里假设你有一个名为seurat_obj的 Seurat 对象已经完成了NormalizeData(),FindVariableFeatures(),ScaleData(),RunPCA(),RunUMAP()等标准步骤。4. 核心流程四步走从数据到评分AUCell 的分析流程非常清晰主要分为四步。4.1 第一步提取表达矩阵并构建排名AUCell 需要输入一个数值矩阵。我们从 Seurat 对象中提取经过归一化如LogNormalize的表达矩阵。# 提取表达矩阵。使用 GetAssayData 并指定 slot“data” 获取归一化数据。 expr_matrix - GetAssayData(seurat_obj, slot “data”) # 得到一个稀疏矩阵 dgCMatrix # expr_matrix 是一个基因 x 细胞的矩阵 # 查看矩阵维度 dim(expr_matrix) # [1] 20000 5000 # 例如20000个基因5000个细胞 # 构建基因表达排名。这是AUCell最耗计算的一步但每个排名只需计算一次。 # 它会为每个细胞计算所有基因的表达量排名。 cells_rankings - AUCell_buildRankings( exprMatrix expr_matrix, nCores 4, # 根据你的电脑核心数设置可以加速计算 plotStats TRUE # 绘制排名分布的统计图检查是否正常 )运行AUCell_buildRankings后会生成一个cells_rankings对象。plotStatsTRUE会生成一张图展示每个细胞中基因排名的分布通常应该看到一条平滑的曲线如果出现异常如大量平局排名可能需要检查数据。4.2 第二步准备目标基因集基因集可以来自多种渠道。这里演示三种常见方式方式一手动创建列表# 如前所述手动定义基因集列表 my_signatures - list( “Cell_Cycle_G2M” c(“TOP2A”, “MKI67”, “BIRC5”, “CCNB2”, “UBE2C”), “Hypoxia” c(“VEGFA”, “BNIP3”, “SLC2A1”, “PGK1”, “CA9”), “My_Own_Pathway” c(“GENE1”, “GENE2”, “GENE3”) )方式二从 MSigDB 的 .gmt 文件读取# 假设你从 MSigDB官网下载了文件 c2.cp.kegg.v2023.1.Hs.symbols.gmt gmt_file - “path/to/your/c2.cp.kegg.v2023.1.Hs.symbols.gmt” gene_sets - getGmt(gmt_file) # 转换为 AUCell 需要的列表格式 gene_sets_list - geneIds(gene_sets) # gene_sets_list 就是一个包含多个通路基因集的列表方式三从已有分析结果获取# 例如从你的差异表达分析结果中提取某个细胞簇的上调基因作为特征 cluster_markers - FindMarkers(seurat_obj, ident.1 “Cluster1”, min.pct 0.25) top_genes - rownames(cluster_markers)[1:50] # 取前50个差异基因 my_signatures$“Cluster1_Signature” - top_genes关键检查 必须确保基因集中的基因名与表达矩阵的行名完全匹配大小写、符号一致。# 检查并过滤掉表达矩阵中不存在的基因 for (set_name in names(my_signatures)) { genes_in_matrix - my_signatures[[set_name]] %in% rownames(expr_matrix) cat(set_name, “:”, sum(genes_in_matrix), “/”, length(genes_in_matrix), “genes found\n”) # 可以选择只保留存在的基因 my_signatures[[set_name]] - my_signatures[[set_name]][genes_in_matrix] }4.3 第三步计算 AUC 值这是核心计算步骤使用上一步构建的排名和准备好的基因集。# 计算每个细胞对每个基因集的AUC值 cells_AUC - AUCell_calcAUC( geneSets my_signatures, rankings cells_rankings, nCores 4, # 并行计算 aucMaxRank ceiling(0.05 * nrow(cells_rankings)) # 最重要的参数之一见下文解释 )参数aucMaxRank详解 这是AUCell最关键的参数没有之一。它决定了在计算AUC时只考虑每个细胞中表达量排名前aucMaxRank的基因。为什么需要它在单细胞数据中大量基因表达为0或极低如果考虑所有基因排名AUC值会主要由这些低表达基因的尾部排名决定使得评分区分度下降。限制只关注高表达区域能放大有生物学意义的信号。如何设置通常设置为所有基因数的前5% (0.05 * nrow(rankings))。这是一个经验值你可以尝试2%-10%。可以通过AUCell_exploreThresholds函数来辅助评估。4.4 第四步提取结果并整合计算完成后我们从cells_AUC对象中提取评分矩阵。# 提取评分矩阵 auc_matrix - getAUC(cells_AUC) dim(auc_matrix) # 矩阵维度为基因集数量 x 细胞数量 # 转置一下变成细胞 x 基因集更方便后续操作 auc_matrix - t(auc_matrix)整合到 Seurat 对象 这是推荐的做法便于统一管理。# 将评分矩阵添加到 Seurat 对象的 metadata 中 # 确保细胞顺序一致AUCell默认可能按细胞名排序Seurat对象也是 cell_names - colnames(seurat_obj) auc_matrix - auc_matrix[cell_names, ] # 按Seurat对象细胞顺序重排 for (sig_name in colnames(auc_matrix)) { seurat_obj[[sig_name]] - auc_matrix[, sig_name] } # 检查是否添加成功 head(seurat_objmeta.data)5. 结果可视化与解读让评分“说话”算出分数只是开始如何解读和展示才是关键。5.1 基础可视化在 UMAP/tSNE 上着色最直观的方式是将AUCell评分映射到降维聚类图上。# 绘制单个基因集评分在UMAP上的分布 FeaturePlot(seurat_obj, features “Cell_Cycle_G2M”, # 替换为你的基因集名称 cols c(“lightgrey”, “blue”)) # 颜色从低到高 # 使用更专业的颜色梯度并添加标题 FeaturePlot(seurat_obj, features “Hypoxia”, cols viridis::viridis(10), # 使用viridis色系对色盲友好 order TRUE) # orderTRUE 将高分细胞画在最上层避免被遮盖 ggtitle(“Hypoxia Signature Activity”)5.2 评分分布分析识别“阳性”细胞AUCell评分是连续的但有时我们需要定义一个阈值来区分“高活性”和“低活性”细胞。# 方法1使用AUCell内置函数自动探索阈值 cells_assignment - AUCell_exploreThresholds(cells_AUC, plotHist TRUE, assignCells TRUE) # 这个函数会为每个基因集生成一个直方图并建议一个阈值如“跳出双峰分布的谷底”。 # 结果存储在 cells_assignment 列表中。 # 查看为“Cell_Cycle_G2M”基因集分配的细胞 g2m_cells - cells_assignment$“Cell_Cycle_G2M”$assignment # 这是一个细胞名向量代表被认定为该特征“阳性”的细胞。 # 方法2手动设定阈值例如取前20%的细胞 scores - seurat_obj$“Cell_Cycle_G2M” threshold - quantile(scores, probs 0.8) # 取80%分位数即前20% seurat_obj$“G2M_high” - scores threshold # 在UMAP上查看两类细胞 DimPlot(seurat_obj, group.by “G2M_high”, cols c(“grey”, “red”)) ggtitle(“Top 20% Cells for G2M Signature”)5.3 多基因集比较与热图展示当有多个基因集时可以观察它们在不同细胞群中的活性模式。# 提取所有细胞的评分 auc_scores - FetchData(seurat_obj, vars c(“Cell_Cycle_G2M”, “Hypoxia”, “My_Own_Pathway”)) # 按细胞聚类取平均分观察不同簇的特征活性 avg_scores - seurat_objmeta.data %% group_by(seurat_clusters) %% # 假设你的聚类信息列名为‘seurat_clusters’ summarise(across(c(“Cell_Cycle_G2M”, “Hypoxia”, “My_Own_Pathway”), mean)) # 绘制热图 library(pheatmap) score_mat - as.matrix(avg_scores[, -1]) # 去掉第一列簇名 rownames(score_mat) - avg_scores$seurat_clusters pheatmap(score_mat, cluster_rows TRUE, cluster_cols TRUE, scale “column”, # 按列基因集标准化便于比较不同基因集在各簇的相对高低 color colorRampPalette(c(“navy”, “white”, “firebrick3”))(50), main “Average AUCell Score per Cluster”)这张热图能清晰告诉你哪个细胞簇更富集细胞周期特征哪个簇更处于缺氧状态。6. 进阶应用将 AUCell 评分融入下游分析AUCell评分不仅仅是用来画图的它可以作为强有力的特征输入到下游分析中。6.1 驱动聚类分析传统的聚类基于所有基因的表达变化。你可以尝试加入AUCell评分作为补充特征让聚类结果更倾向于反映功能状态。# 将AUCell评分矩阵作为新的“assay”加入Seurat对象 # 首先确保评分矩阵的细胞顺序与对象一致 auc_matrix - t(getAUC(cells_AUC))[, colnames(seurat_obj)] # 细胞 x 基因集 # 创建一个新的Assay seurat_obj[[“AUC”]] - CreateAssayObject(data t(auc_matrix)) # Assay要求基因 x 细胞 # 设置默认assay为AUC并运行PCA和UMAP DefaultAssay(seurat_obj) - “AUC” seurat_obj - ScaleData(seurat_obj) seurat_obj - RunPCA(seurat_obj, features rownames(seurat_obj[[“AUC”]])) seurat_obj - RunUMAP(seurat_obj, dims 1:10) # 基于功能评分进行聚类 seurat_obj - FindNeighbors(seurat_obj, dims 1:10) seurat_obj - FindClusters(seurat_obj, resolution 0.5) # 可视化基于功能特征的聚类 DimPlot(seurat_obj, reduction “umap”, label TRUE)比较基于基因表达的聚类和基于功能评分的聚类可以发现哪些功能模块定义了新的细胞亚群。6.2 拟时序分析中的功能动力学如果你在使用Monocle3或Slingshot做拟时序分析可以将AUCell评分作为细胞的状态特征观察其沿轨迹的变化。# 假设你已经有了一个cds对象Monocle3 library(monocle3) # 将AUCell评分添加到cds的colData中 cds_col_data - colData(cds) for (sig in colnames(auc_matrix)) { cds_col_data[[sig]] - auc_matrix[rownames(cds_col_data), sig] } colData(cds) - cds_col_data # 绘制评分沿轨迹的平滑曲线 plot_genes_violin(cds, group_cells_by “pseudotime_bin”, # 将拟时序分箱 genes c(“Stemness_Signature”, “Differentiation_Signature”), ncol 2)6.3 细胞-细胞通讯分析中的配体-受体活性在分析细胞通讯如CellChat时你不仅关心配体受体对的表达更关心其下游通路的活性。AUCell可以在这里大显身手。# 假设你有一个从CellChat或其他资源获取的“NF-kB通路靶基因集” nfkb_target_genes - c(“IL6”, “TNF”, “CXCL8”, “ICAM1”, …) # 你的基因列表 # 为每个细胞计算NF-kB通路活性评分 # … (使用AUCell计算步骤同上) … # 然后在分析细胞通讯时你可以问 # “在配体高表达的发送细胞中其NF-kB通路活性是否也更高” # “在受体高表达的接收细胞中NF-kB通路活性是否被激活” # 这可以通过相关性分析或分组比较来实现。7. 常见问题、陷阱与排查指南AUCell用起来简单但想用对、用好必须避开以下几个常见的“坑”。问题现象可能原因排查方式解决方案所有细胞的评分都差不多没有区分度1.aucMaxRank参数设置过大如用了默认值。2. 基因集质量差基因不特异或表达极低。3. 输入的表达矩阵未正确归一化如还是counts。1. 检查aucMaxRank值尝试设置为基因总数的1%-10%。2. 用VlnPlot或FeaturePlot查看基因集中部分基因的表达情况。3. 确认expr_matrix来自slot“data”。1. 调整aucMaxRank并重新计算。2. 重新筛选或寻找更特异的基因集。3. 确保使用归一化后的数据。评分出现大量 NA 或计算错误1. 基因集中的基因在表达矩阵中一个都没找到。2. 基因名不匹配大小写、符号版本。1. 运行第4.2步的“关键检查”代码。2. 使用rownames(expr_matrix)[1:10]和my_signatures[[1]][1:10]对比。1. 修正基因名或使用基因ID转换工具。2. 确保使用一致的基因标识符如Symbol或Ensembl ID。运行AUCell_buildRankings内存不足或极慢细胞数或基因数过多如 10万细胞。监控内存使用。对于超大矩阵排名构建是内存和计算密集型步骤。1. 增加nCores并行。2.对细胞进行降采样先在小样本上测试流程和参数。3. 考虑使用计算集群。不同基因集的评分分布差异巨大这是正常现象。基因集大小直接影响AUC值的理论分布。大基因集AUC值更集中小基因集更分散。使用boxplot或hist查看不同基因集评分的分布。不要直接比较绝对值关注每个基因集内部细胞间的相对高低。使用分位数或Z-score标准化后再比较。我的“阳性”细胞在图上看起来不连续、很散1. 阈值设置不合理太高或太低。2. 该功能特征本身就是跨细胞类型的、或呈梯度变化。3. 降维图UMAP本身不能完美分离所有功能状态。1. 用AUCell_exploreThresholds检查阈值建议的直方图。2. 观察评分在已知细胞类型注释上的分布VlnPlot。1. 调整阈值或使用连续的评分进行分析如相关性分析。2. 结合其他证据如关键标记基因表达综合判断。8. 最佳实践与工程化建议要让AUCell分析真正成为可重复、可信赖的研究的一部分你需要遵循以下实践基因集质量控制是重中之重来源可靠优先使用权威数据库MSigDB, GO, KEGG或高质量文献中验证过的基因集。大小适中避免过小5个基因结果不稳定或过大500个基因信号可能被稀释的基因集。10-200个基因是较理想的区间。特异性确保基因集与你要研究的问题高度相关。使用从同一数据集的差异分析中得到的基因集时要格外小心避免循环论证。参数aucMaxRank需要理性选择不要盲目使用默认值。用AUCell_exploreThresholds函数观察不同aucMaxRank下评分分布的变化。对于关注强表达信号的分析如核心转录因子可以用更小的aucMaxRank如前1%。对于关注广泛、微弱信号的分析如某些代谢通路可以适当调大。在同一项研究中对所有基因集使用相同的aucMaxRank参数以保证可比性。结果解读必须结合生物学背景AUCell高分提示特征活性高但并非绝对证据。必须与已知的细胞类型标记基因表达进行交叉验证。对于新发现的“高评分细胞群”要回到原始表达数据检查是哪些基因驱动了高分。考虑使用多种评分算法如AddModuleScore, UCell, singscore进行交叉验证特别是对于关键结论。代码的可重复性与文档在脚本开头明确记录AUCell包版本、R版本和所有参数尤其是aucMaxRank。保存关键的中间结果如cells_rankings对象和最终的评分矩阵避免重复计算。为每个基因集标注清晰的名称和来源。融入分析流程将AUCell评分计算封装成函数方便在多个数据集或不同预处理条件下重复使用。考虑将评分结果作为标准输出的一部分整合到你的单细胞分析Pipeline中。单细胞数据分析的魅力在于从海量数据中提炼生物学故事。AUCell提供了一把精准的尺子让你能够定量地衡量每个细胞的“功能状态”。它不再让你停留在“这个基因高不高”的层面而是让你能回答“这个细胞群体是否具备某种功能特性”这样的高阶问题。从理解其基于排名的核心思想开始到谨慎地准备基因集和设置参数再到将评分结果可视化并与下游分析深度整合每一步都需要清晰的生物学问题驱动和技术细节把控。当你熟练运用AUCell后你会发现它不仅是流程中的一个步骤更是你探索单细胞世界功能景观的一双眼睛。