代码改的比较粗糙,自用。
AddModuleScore_bulk <- function (object, features, nbin = 24, ctrl = 100,
k = FALSE, assay = NULL, name = "Cluster", seed = 1)
{
set.seed(seed = 1)
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 = FALSE, yes = ", attempting to find updated synonyms",
no = ", not searching for symbol synonyms"),
call. = FALSE, immediate. = TRUE)
}
return(intersect(x = x, y = rownames(x = object)))
})
cluster.length <- length(x = features)
pool <- rownames(x = object)
data.avg <- apply(object,1, mean)
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(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, ] <- apply(object[features.use,], 2, mean)
}
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 <- object[features.use, , drop = FALSE]
features.scores[i, ] <- apply(data.use, 2, mean)
}
features.scores.use <- features.scores - ctrl.scores
colnames(features.scores.use) <- colnames(object)
rownames(features.scores.use) <- name
return(features.scores.use)
}
输入的是一个N×M的矩阵,features是用于打分的基因list。一般情况下输入一个list,代码中的cluster.length=1。
score <- AddModuleScore_bulk(expr, features, name='score')