Addmodulescore is a built-in functions in Seurat.

How it works?

The first step is to cut all genes in ==nbin== groups.

1
2
3
4
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))/1e30, n = nbin, labels = FALSE, right = FALSE)

The second step is to build background genes.

Firstly, to subset the ctrl number genes from same features group;

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
        ))
      )
    }
}

Secondly, to calculate each control case's gene average expression.

1
2
3
4
5
6
7
8
9
10
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, ])
	}

The third step is to calculate each geneset's features mean expression scores. The final geneset's score is from features.scores minus ctrl.scores.

1
2
3
4
5
6
7
8
9
10
11
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

Reference