单细胞分析十八般武艺16：AUCell

      单细胞测序技术的发展日新月异，新的分析工具也层出不穷。每个工具都有它的优势与不足，在没有权威工具和流程的单细胞生信江湖里，多掌握几种分析方法和工具，探索数据时常常会有意想不到的惊喜。

往期专题

单细胞转录组基础分析专题

单细胞转录组高级分析专题

看过我以前帖子的朋友应该对基于GSVA包的基因集评分方法有所了解，也可能在做SCENIC分析时接触了AUCell，他们之间有什么区别呢？最大的不同在于GSVA算法是多年前为表达谱芯片和bulkRNA测序数据开发的，并没有考虑到10X Genomics 技术平台（或其他高通量单细胞技术平台）产生的稀疏性数据，AUCell则是针对单细胞数据开发的基因集评分方法。

AUCell简介

安装AUCell

AUCell收录在Bioconductor平台，安装起来比较简单。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("AUCell")
# 查看官方教程
browseVignettes("AUCell")

基本原理

AUCell为每个细胞基于表达值对基因进行排序（rank），表达值相同的基因随机赋予顺序，然后用rank值代替表达值进行后续分析。它用恢复曲线（recovery curve）下的面积代表富集分数。曲线的横坐标用rank值（如下图的1-500），每个rank值对应的基因集合记作topRankSet，纵坐标对应topRankSet在评分基因集中的基因数量。例如，rank100对应的曲线横坐标是100，其topRankSet包含此细胞排名前100的基因，假设topRankSet有15个基因在评分基因集中，那么rank100对应的曲线纵坐标就是15。

计算AUC值的横坐标终点对应AUCell_calcAUC函数中的aucMaxRank参数，默认值为数据集中所有基因的前5%。这可以更快地完成大数据集上运算，并减少排名底部噪音的影响（例如很多基因的表达值为0，但是也有rank值）。对于基因数量（例如10x数据中的nfeature_count）比较高的数据集，或表达值比较高的数据集，最好提高该阈值。

分析流程

AUCell的分析流程包含以下三个步骤：

1. 将基因表达矩阵转换为rank矩阵；

2. 计算恢复曲线下面积；

3. 为每个基因集设置活性阈值。

注意事项

AUC的值并不是绝对的，它取决于细胞类型（如细胞尺寸、转录本数量）、测序技术（如检测灵敏度）和基因集本身。通常不能获得一个只在感兴趣细胞中完全“打开”，而在其他细胞中“关闭”的理想基因集。对于基因数量较小的基因集，很可能得到AUC=0的细胞。较大的基因集（100-2k基因）可能更稳定，结果更准确。作者在官方教程中强调：基因集活性阈值的自动选择并非详尽无遗，因此他们强烈建议用户检查AUC直方图，并在需要时手动选择阈值。

数据实测

测试数据

使用SeuratData包提供的pbmc3k数据集，细胞类型已经注释好了。

library(AUCell)
library(Seurat)
library(SeuratData)
library(tidyverse)
library(patchwork)
rm(list = ls())

### 整理演示数据
# 演示数据可能需要下一行代码提前安装
#InstallData("pbmc3k.SeuratData")  
pbmc <- LoadData("pbmc3k.SeuratData")
pbmc <- UpdateSeuratObject(pbmc)
pbmc <- NormalizeData(pbmc)
pbmc <- SCTransform(pbmc)
pbmc <- RunPCA(pbmc)
ElbowPlot(pbmc, ndims = 50)
pc.num <- 1:15
pbmc <- pbmc %>% RunTSNE(dims=pc.num) %>% RunUMAP(dims=pc.num)
pbmc <- pbmc %>% FindNeighbors(dims=pc.num) %>% FindClusters()
pbmc <- pbmc[,!is.na(pbmc$seurat_annotations)]
p <- DimPlot(pbmc, group.by = "seurat_annotations", label = T)
ggsave("1601_celltype.png", p, width = 8, height = 6)

### 提取细胞类型signature
Idents(pbmc) <- "seurat_annotations"
celltype.markers <- FindAllMarkers(pbmc, assay = "RNA", only.pos = T)
signatures <- lapply(levels(pbmc), function(x){
  rownames(subset(celltype.markers, cluster==x&p_val_adj<0.05))
})

AUCell评分

### AUCell评分并分配活性阈值
# 提取矩阵
exprMatrix <- pbmc@assays$RNA@data
# 表达矩阵转rank矩阵
cells_rankings <- AUCell_buildRankings(exprMatrix)
#    min      1%      5%     10%     50%    100% 
# 212.00  341.00  458.85  554.00  819.00 2413.00
#以上数据意味着排名2413以后的基因没有价值

# 参数aucMaxRank的设置最好参考上面的结果，此处先不调整默认参数
cells_AUC <- AUCell_calcAUC(signatures, cells_rankings, nCores = 10, 
                            aucMaxRank=nrow(cells_rankings)*0.05)
# 给每个基因集分配活性阈值
set.seed(123)
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, assign=TRUE)

自动分配的活性阈值，效果确实值得商榷，如下面两图：

AUCell判定所有细胞对应Naive CD4 T的signature都是激活的

CD14+ Mono出现了明显的双峰，阈值判定结果相对可靠一些

结果展示

先直接用AUC评分画个热图，看看是否与细胞类型匹配。

library(pheatmap)
aucMat <- getAUC(cells_AUC)
pheatmap(aucMat, show_colnames = F, scale = "row")

再来看看使用活性阈值将auc矩阵转为二进制矩阵后的效果。

### AUC矩阵二进制转换
# 自动选择的阈值有些不合适，人工调整一下
T.adj <- c("Naive CD4 T" = 0.52,
           "Memory CD4 T" = 0.20,
           "CD14+ Mono" = 0.235,
           "B" = 0.13,
           "CD8 T" = 0.165,
           "FCGR3A+ Mono" = 0.11,
           "NK" = 0.095,
           "DC" = 0.035,
           "Platelet" = 0.04)
# auc矩阵转为只有0和1的二进制矩阵，1代表激活          
aucBin <- aucMat
for(i in rownames(aucMat)){
  aucBin[i,] <- ifelse(aucMat[i,]>T.adj[i], 1, 0)
}
pheatmap(aucBin, show_colnames = F)

性能对比

这次测试选择已知细胞类型的topMarkers作为评分基因集，主要是为了直观地对比几种常见评分算法。效果孰优孰劣，大家自行判断

AddModuleScore

Seurat包中的AddModuleScore函数，第一步是计算表达矩阵中所有基因的平均表达值，然后根据平均表达值把所有基因分成24个（默认值）背景基因集。假设评分基因集为ScoreSet，分割的背景基因集为BinSets。第二步是为ScoreSet中的每个基因找到对应的BinSet，从其中随机挑选100个（默认值）基因作为背景基因。如果ScoreSet有10个基因，那么算法会为它随机挑选最多1000个背景基因（选完之后会去重），那么最后的评分是用这10个基因的均值减去那1000个基因的均值。

sco <- AddModuleScore(pbmc, features = signatures, name = "M", assay = "RNA")
mat <- t([e-mail][,c(paste0("M",1:9))])
rownames(mat) <- levels(pbmc)
pheatmap(mat, show_colnames = F, scale = "row")

ssGSEA

首先用表达值对每个样本（细胞）的基因进行排名，并用排名顺序代替表达值用于后续分析。然后用样本基因在评分基因集内和基因集外的经验累积分布函数之间的差值生成富集分数。通俗地讲就是按基因排名顺序，逐个检查基因是否在评分基因集中；如果在基因集内富集分数加分，反之扣分，并且排名越靠前加分权重越大。最后，根据基因表达值的极差对富集分数进行标准化。

library(GSVA)
expr <- as.matrix(pbmc@assays$RNA@data)
mat1 <- gsva(expr, signatures, method="ssgsea", parallel.sz=18L)
pheatmap(mat1, show_colnames = F, scale = "row")

GSVA

GSVA算是ssGSEA升级版，主要变化是在排序之前会用核函数对表达矩阵进行归一化，不同的数据类型归一化算法不一样，可以通过kcdf参数指定。

mat2 <- gsva(expr, signatures, method="gsva", kcdf="Gaussian", parallel.sz=18L)
pheatmap(mat2, show_colnames = F, scale = "row")

交流探讨

如果您阅读此文有所疑惑，或有不同见解，亦或其他问题，欢迎添加我的微信探讨。我工作的公司2016年开始单细胞测序服务，至今已完成近万例样本的单细胞测序，服务质量经过了10X Genomics公司的权威认证。欢迎大家联系Kinesin洽谈单细胞测序及数据分析业务！