前言
今日份看文献《SCENIC: single-cell regulatory network inference and clustering》
发现SCENIC的最后一部分，计算基因集在每一个细胞中的活跃程度调用的是AUCell进行计算的
AUCell输入基因集的一个关键子集是否在每个细胞的表达基因中活跃
AUCell
该软件的地址为：传送门
AUCell根据单细胞的表达矩阵进行计算的，首先我们随机生成了单细胞表达矩阵
# 生成表达矩阵
exprMatrix <- matrix(rnorm(15*45),nrow=45)
colnames(exprMatrix) = c(paste("cell",rep(1:15),sep = '_'))
row.names(exprMatrix) = c(paste("gene",rep(1:45),sep = '_'))
# 设定基因集
genes_1 <- c("gene_1", "gene_2", "gene_3","gene_4", "gene_5", "gene_6")
genes_2 = c("gene_13", "gene_14", "gene_15","gene_10", "gene_11", "gene_12")
geneSets <- list(geneSet1=genes_1,geneSet2=genes_2)
# 读为data.table的对象
exprMat <- data.table::data.table(exprMat, keep.rownames=TRUE)
# 按照基因名进行排序
data.table::setkey(exprMat, "rn") # (reorders rows)
colsNam <- colnames(exprMat)[-1] # 1=row names
# 对每一个细胞的基因表达量分别排序，并且将表达量换作序数
exprMat[, (colsNam):=lapply(-.SD, data.table::frank, ties.method="random", na.last="keep"),
.SDcols=colsNam]
# 重命名
rn <- exprMat$rn
exprMat <- as.matrix(exprMat[,-1])
rownames(exprMat) <- rn
names(dimnames(exprMat)) <- c("genes", "cells")
其中：
exprMat
- rn 列代表的是基因
- 每一个cell 中的数值代表某基因在某细胞中表达量排序的序号（表达量从高到低排序）
计算AUC
aucMaxRank=3
rankings = data.frame(exprMat)
row.names(rankings) = rankings[,1]
rankings = rankings[,-1]
# 定义计算AUC的函数
.AUC.geneSet <- .AUC.geneSet_norm
#### 1. Calculate the AUC for each gene set
gSetName <- NULL
aucMatrix <- sapply(names(geneSets), function(gSetName)
.AUC.geneSet(geneSet=geneSets[[gSetName]], rankings=rankings,
aucMaxRank=aucMaxRank, gSetName=gSetName))
aucMatrix <- t(aucMatrix)
上述结果反应了之前定义geneSet在每一个细胞中的活性，值越高，活性越高
计算AUC的函数
我们以其中一个基因集为例：
# 过滤基因集
geneSet <- unique(geneSets$geneSet2)
nGenes <- length(geneSet)
geneSet <- geneSet[which(geneSet %in% rownames(rankings))]
missing <- nGenes-length(geneSet)
# 获取基因集中的基因在每个cell中表达量排序的序号
gSetRanks <- rankings[which(rownames(rankings) %in% geneSet),,drop=FALSE]
# 定义AUC的阈值
aucThreshold <- round(aucMaxRank)
# 计算最大的AUC值
## x_th 为1:nrow(gSetRanks)的序数
x_th <- 1:nrow(gSetRanks)
## 对序号进行排序，序号小的代表基因表达量高
x_th <- sort(x_th[x_th<aucThreshold])
## y_th为长度为x_th的序数
y_th <- seq_along(x_th)
maxAUC <- sum(diff(c(x_th, aucThreshold)) * y_th)
# Apply by columns (i.e. to each ranking)
## 按列计算每一个cell的 AUC
auc <- apply(gSetRanks, 2, .auc, aucThreshold, maxAUC)
# 计算AUC的函数：
.auc <- function(oneRanking, aucThreshold, maxAUC)
{
x <- unlist(oneRanking)
# 选择满足阈值的基因序号
x <- sort(x[x<aucThreshold])
# y = 1,2,3 ...；作为纵坐标用于计算面积
y <- seq_along(x)
sum(diff(c(x, aucThreshold)) * y)/maxAUC
}
如何理解AUC值表征基因集在cell中的活性
不难发现，AUCell基于基因在cell中表达量的排序（序号）进行AUC的计算。序号小的代表该基因表达量高；反之序号大的代表该基因表达量低
假设我们给的基因集里面的基因个数都是5个基因，，AUC阈值是7，我们举如下几个例子：
[1]. cell 1中该基因集的序号为：4，5，6，7（紫色）
[2]. cell 2中该基因集的序号为：1，2，3，6，7 （绿色）
[3]. cell 1中该基因集的序号为：1，5，7 （黄色）
红色为最大的AUC值
显然绿色AUC值最大，黄色次之，紫色最小；由它们的序数我们也可以得知：紫色的序数靠后，说明基因表达量在该细胞中不算特别高；绿色的序数有三个都属于top3，所以AUC最大；黄色则次之
由此可见AUC的大小反应某基因集的基因在细胞中的表达量大小，AUC越大，该基因集的基因在该细胞中的表达量越大，故越活越
