在读张泽民老师发表的新冠文献的时候看到计算免疫细胞的cytokine score或inflammatory score使用了Seurat包的AddMouduleScore函数，在计算细胞周期的CellCycleScoring函数的原代码中也使用了这个函数。
此功能可用于计算任何基因列表的监督模块评分，非常实用！！！
1. 查看一下这个函数的用法
library(Seurat)
?AddModuleScore
2. 使用
输入数据：Seurat对象和一个gene list。
- 1 准备矩阵
library(tidyverse)
library(Matrix)
library(cowplot)
pbmc <- readRDS("pbmc.rds")
- 2 准备想要计算的gene list
inflammatory_gene <- readxl::read_xlsx("inflammatory_gene.xlsx")
View(inflammatory_gene)
#转换成list
gene <- as.list(inflammatory_gene)
- 3 计算
pbmc <- AddModuleScore(
object = pbmc,
features = gene,
ctrl = 100, #默认值是100
name = 'CD_Features'
)
计算结果保存在pbmc@meta.data[["CD_Features1"]]
得到的score是在每个细胞中算出来的我们感兴趣的基因的表达均值。
其在使用时需要先计算基因集中所有基因的平均值，再根据平均值把表达矩阵切割成若干个bin，然后从切割后的每一个bin随机抽取对照基因（基因集外的基因，默认100个）作为背景值。最后所有的目标基因算一个平均值，所有的背景基因算一个平均值，两者相减就是该gene set的score值。因此，在整合不同样本的情况下，即使使用相同基因集为相同细胞打分，也会产生不同的富集评分。从本质上看它和Zscore的算法很类似，Zscore又称Z值，原是一个统计学概念，表示的是个体测量值X以标准差σ为单位，偏离总体均数μ的距离，即：Z score=（X-μ）/σ。牵扯到统计学的概念不免有些难以理解，简单说它就是处理过的平均值。
colnames(pbmc@meta.data)
# [1] "orig.ident" "nCount_RNA" "nFeature_RNA"
# [4] "percent.mt" "RNA_snn_res.0.5" "seurat_clusters"
# [7] "cell_type" "CD_Features1"
colnames(pbmc@meta.data)[8] <- 'Inflammatory_Score'
VlnPlot(pbmc,features = 'Inflammatory_Score')
绘制箱线图
data<- FetchData(pbmc,vars = c("cell_type","Inflammatory_Score"))
p <- ggplot(data,aes(cell_type,Inflammatory_Score))
p+geom_boxplot()+theme_bw()+RotatedAxis()
绘制分数的分布图
library(ggplot2)
mydata<- FetchData(pbmc,vars = c("UMAP_1","UMAP_2","Inflammatory_Score"))
a <- ggplot(mydata,aes(x = UMAP_1,y =UMAP_2,colour = Inflammatory_Score))+geom_point(size = 1)+scale_color_gradientn(values = seq(0,1,0.2),colours = c('blue','cyan','green','yellow','orange','red'))
a+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
打分是一个连续变化的值，我们也可以直接将其转换为一定的分类变量：
#按照均值分
pbmc[["stage"]] <- ifelse(pbmc@meta.data[,"Inflammatory_Score"] > mean(pbmc@meta.data[,"Inflammatory_Score"]),"High","Low")
rownames(pbmc@meta.data)
#按照75%分
pbmc[["stage"]] <-pbmc@meta.data[,"Inflammatory_Score"] <- ifelse(pbmc@meta.data[,"Inflammatory_Score"] > quantile(pbmc@meta.data[,"Inflammatory_Score"],0.75),"High","Low")
#按照具体数值分
pbmc[["stage"]] <-pbmc@meta.data[,"Inflammatory_Score"] <- ifelse(pbmc@meta.data[,"Inflammatory_Score"] > 0.5,"High","Low")
Idents(pbmc) <- 'Inflammatory_Score'
DimPlot(pbmc,reduction = "umap",label = TRUE,pt.size = 1.5)
3. 主要参数
| 参数 | 意义 |
|---|---|
| object | Seurat object |
| features | A list of vectors of features for expression programs; each entry should be a vector of feature names |
| pool | List of features to check expression levels against, defaults to rownames(x = object) |
| nbin | Number of bins of aggregate expression levels for all analyzed features |
| ctrl | Number of control features selected from the same bin per analyzed feature |
| k | Use feature clusters returned from DoKMeans |
| assay | Name of assay to use |
| name | Name for the expression programs; will append a number to the end for each entry in features (eg. if features has three programs, the results will be stored as name1, name2, name3, respectively) |
| seed | Set a random seed. If NULL, seed is not set. |
| search | Search for symbol synonyms for features in features that don't match features in object? Searches the HGNC's gene names database; see UpdateSymbolList for more details |
| ... | Extra parameters passed to UpdateSymbolList
|
4. 函数代码
AddModuleScore
function (object, features, pool = NULL, nbin = 24, ctrl = 100,
k = FALSE, assay = NULL, name = "Cluster", seed = 1, search = FALSE,
...)
{
if (!is.null(x = seed)) {
set.seed(seed = seed)
}
assay.old <- DefaultAssay(object = object)
assay <- assay %||% assay.old
DefaultAssay(object = object) <- assay
assay.data <- GetAssayData(object = object)
features.old <- features
if (k) {
.NotYetUsed(arg = "k")
features <- list()
for (i in as.numeric(x = names(x = table(object@kmeans.obj[[1]]$cluster)))) {
features[[i]] <- names(x = which(x = object@kmeans.obj[[1]]$cluster ==
i))
}
cluster.length <- length(x = features)
}
else {
if (is.null(x = features)) {
stop("Missing input feature list")
}
features <- lapply(X = features, FUN = function(x) {
missing.features <- setdiff(x = x, y = rownames(x = object))
if (length(x = missing.features) > 0) {
warning("The following features are not present in the object: ",
paste(missing.features, collapse = ", "), ifelse(test = search,
yes = ", attempting to find updated synonyms",
no = ", not searching for symbol synonyms"),
call. = FALSE, immediate. = TRUE)
if (search) {
tryCatch(expr = {
updated.features <- UpdateSymbolList(symbols = missing.features,
...)
names(x = updated.features) <- missing.features
for (miss in names(x = updated.features)) {
index <- which(x == miss)
x[index] <- updated.features[miss]
}
}, error = function(...) {
warning("Could not reach HGNC's gene names database",
call. = FALSE, immediate. = TRUE)
})
missing.features <- setdiff(x = x, y = rownames(x = object))
if (length(x = missing.features) > 0) {
warning("The following features are still not present in the object: ",
paste(missing.features, collapse = ", "),
call. = FALSE, immediate. = TRUE)
}
}
}
return(intersect(x = x, y = rownames(x = object)))
})
cluster.length <- length(x = features)
}
if (!all(LengthCheck(values = features))) {
warning(paste("Could not find enough features in the object from the following feature lists:",
paste(names(x = which(x = !LengthCheck(values = features)))),
"Attempting to match case..."))
features <- lapply(X = features.old, FUN = CaseMatch,
match = rownames(x = object))
}
if (!all(LengthCheck(values = features))) {
stop(paste("The following feature lists do not have enough features present in the object:",
paste(names(x = which(x = !LengthCheck(values = features)))),
"exiting..."))
}
pool <- pool %||% rownames(x = object)
data.avg <- Matrix::rowMeans(x = assay.data[pool, ])
data.avg <- data.avg[order(data.avg)]
data.cut <- cut_number(x = data.avg + rnorm(n = length(data.avg))/1e+30,
n = nbin, labels = FALSE, right = FALSE)
names(x = data.cut) <- names(x = data.avg)
ctrl.use <- vector(mode = "list", length = cluster.length)
for (i in 1:cluster.length) {
features.use <- features[[i]]
for (j in 1:length(x = features.use)) {
ctrl.use[[i]] <- c(ctrl.use[[i]], names(x = sample(x = data.cut[which(x = data.cut ==
data.cut[features.use[j]])], size = ctrl, replace = FALSE)))
}
}
ctrl.use <- lapply(X = ctrl.use, FUN = unique)
ctrl.scores <- matrix(data = numeric(length = 1L), nrow = length(x = ctrl.use),
ncol = ncol(x = object))
for (i in 1:length(ctrl.use)) {
features.use <- ctrl.use[[i]]
ctrl.scores[i, ] <- Matrix::colMeans(x = assay.data[features.use,
])
}
features.scores <- matrix(data = numeric(length = 1L), nrow = cluster.length,
ncol = ncol(x = object))
for (i in 1:cluster.length) {
features.use <- features[[i]]
data.use <- assay.data[features.use, , drop = FALSE]
features.scores[i, ] <- Matrix::colMeans(x = data.use)
}
features.scores.use <- features.scores - ctrl.scores
rownames(x = features.scores.use) <- paste0(name, 1:cluster.length)
features.scores.use <- as.data.frame(x = t(x = features.scores.use))
rownames(x = features.scores.use) <- colnames(x = object)
object[[colnames(x = features.scores.use)]] <- features.scores.use
CheckGC()
DefaultAssay(object = object) <- assay.old
return(object)
}
<bytecode: 0x7fc8df13f008>
<environment: namespace:Seurat>
5. 参考文献：
网友评论