module scores (eg from genescorematrix, like seurat addModuleScore)

[juliabelk]

Would it be possible to add a function to combine multiple gene scores? Such as the addModuleScore function in Seurat?

https://rdrr.io/github/satijalab/seurat/man/AddModuleScore.html

Thanks!

[rcorces]

+1 I like this idea

[jgranja24]

Hi @juliabelk,

Sorry for the delays. I am now adding an addModuleScore to the new release. To use this you can use the following little tutorial below-

#Load Hidden Helper Functions
fn <- unclass(lsf.str(envir = asNamespace("ArchR"), all = TRUE))
for (i in seq_along(fn)) {
    tryCatch({
        eval(parse(text = paste0(fn[i], "<-ArchR:::", fn[i])))
    }, error = function(x) {
    })
}

#New Function
addModuleScore <- function(
	ArchRProj = NULL,
	useMatrix = NULL,
	features = NULL,
	groupBy = "Clusters",
	nBin = 25,
	nBgd = 100,
	name = "Module",
	seed = 1,
	threads = getArchRThreads(),
	logFile = createLogFile("addModuleScore")
	){

	if(!is.null(seed)) set.seed(seed)

	#Get Feature DF
	featureDF <- ArchR:::.getFeatureDF(head(getArrowFiles(ArchRProj),2), subGroup=useMatrix)
    rownames(featureDF) <- paste0(featureDF$seqnames, ":", featureDF$idx)
    featureDF$Match <- seq_len(nrow(featureDF))

    if(useMatrix %ni% getAvailableMatrices(ArchRProj)){
    	stop("useMatrix not in available matrices! See getAvailableMatrices!")
    }

	matrixClass <- h5read(getArrowFiles(ArchRProj)[1], paste0(useMatrix, "/Info/Class"))

	if(matrixClass == "Sparse.Assays.Matrix"){
		if(!all(unlist(lapply(unlist(features), function(x) grepl(":",x))))){
		  .logMessage("When accessing features from a matrix of class Sparse.Assays.Matrix it requires seqnames\n(denoted by seqnames:name) specifying to which assay to pull the feature from.\nIf confused, try getFeatures(ArchRProj, useMatrix) to list out available formats for input!", logFile = logFile)
		  stop("When accessing features from a matrix of class Sparse.Assays.Matrix it requires seqnames\n(denoted by seqnames:name) specifying to which assay to pull the feature from.\nIf confused, try getFeatures(ArchRProj, useMatrix) to list out available formats for input!")
		}
	}

	if(grepl(":",unlist(features)[1])){

		sname <- stringr::str_split(unlist(features),pattern=":",simplify=TRUE)[,1]
		name <- stringr::str_split(unlist(features),pattern=":",simplify=TRUE)[,2]

		idx <- lapply(seq_along(name), function(x){
		  ix <- intersect(which(tolower(name[x]) == tolower(featureDF$name)), BiocGenerics::which(tolower(sname[x]) == tolower(featureDF$seqnames)))
		  if(length(ix)==0){
		    .logStop(sprintf("FeatureName (%s) does not exist! See getFeatures", name[x]), logFile = logFile)
		  }
		  ix
		}) %>% unlist

	}else{

		idx <- lapply(seq_along(unlist(features)), function(x){
		  ix <- which(tolower(unlist(features)[x]) == tolower(featureDF$name))[1]
		  if(length(ix)==0){
		    .logStop(sprintf("FeatureName (%s) does not exist! See getFeatures", unlist(features)[x]), logFile = logFile)
		  }
		  ix
		}) %>% unlist

	}

	if(is.null(names(features))){
		names(features) <- paste0(name, seq_along(features))
	}else{
		names(features) <- paste0(name, ".", names(features))
	}

	featuresUse <- featureDF[idx,]
	featuresUse$Module <- Rle(stack(features)[,2])

	#Get Averages
	rS <- ArchR:::.getRowSums(ArrowFiles = getArrowFiles(ArchRProj), useMatrix = useMatrix)
	rS <- rS[order(rS[,3]), ]
	rS$Bins <- Rle(ggplot2::cut_number(x = rS[,3] + rnorm(length(rS[,3]))/1e30, n = nBin, labels = FALSE, right = FALSE))
	rS$Match <- match(paste0(rS$seqnames, ":", rS$idx), rownames(featureDF))
	
	if(nBgd > min(rS$Bins@lengths)){
		stop("nBgd must be lower than ", min(rS$Bins@lengths), "!")
	}

	idxMatch <- match(paste0(featuresUse$seqnames, ":", featuresUse$idx), paste0(rS$seqnames, ":", rS$idx))
	featuresUse$Bins <- as.vector(rS$Bins[idxMatch])
	
	#MakeLists
	featureList <- split(featuresUse$Match, featuresUse$Module)
	moduleList <- split(featuresUse$Bins, featuresUse$Module)
	binList <- split(rS$Match, rS$Bins)

	dfM <- lapply(seq_along(featureList), function(x){
		message("Computing Module ",x, " of ", length(featureList))
		binx <- binList[moduleList[[x]]]
		idxFgd <- featureList[[x]]
		idxBgd <- unlist(lapply(binx, function(x) sample(x, nBgd)), use.names=FALSE)
		m <- ArchR:::.getPartialMatrix(
			ArrowFiles = getArrowFiles(ArchRProj),
			featureDF = featureDF[c(idxFgd, idxBgd), ],
			useMatrix = useMatrix,
			cellNames = ArchRProj$cellNames,
			threads = threads,
			verbose = FALSE,
			doSampleCells = FALSE
		)
		Matrix::colMeans(m[seq_along(idxFgd), ]) - Matrix::colMeans(m[-seq_along(idxFgd), ])
	}) %>% Reduce("cbind", .)

	for(x in seq_len(ncol(dfM))){
		ArchRProj <- addCellColData(ArchRProj, data = dfM[,x], name=names(featureList)[x], cells=rownames(dfM), force = TRUE)
	}

	ArchRProj

}

#Test Function
projHeme3 <- addImputeWeights(projHeme3)
features <- list(
	BScore = c("MS4A1", "CD79A", "CD74"),
	TScore = c("CD3D", "CD8A", "GZMB", "CCR7", "LEF1")
)
projHeme3 <- addModuleScore(projHeme3, features = features, useMatrix = "GeneScoreMatrix")
p1 <- plotEmbedding(projHeme3, name="Module.TScore", imputeWeights = getImputeWeights(projHeme3))
p2 <- plotEmbedding(projHeme3, name="Module.BScore", imputeWeights = getImputeWeights(projHeme3))
plotPDF(ggAlignPlots(p1,p2,draw=F,type="h"))

Hope this helps

Jeff

[juliabelk]

Just a follow up to this - is it possible to do a ModuleScore of peaks? I have been trying but

getFeatures(proj,useMatrix = "PeakMatrix")

seems to return NULL

[rcorces]

getFeatures() returns NULL because the PeakMatrix features dont have a column called name. Not sure if this is by design.
@jgranja24 ?

[dagarfield]

Hey all -- I think I encountered an issue with this function in cases where you have only one set of genes. That is to say, this works as a feature set:
features <- list( earlyPMC = earlyPMC, germMarkers = germMarkers )

But this doesn't:
features <- list( earlyPMC = earlyPMC)

e.g

ArchRProj <- addModuleScore(ArchRProj, features = features, useMatrix = "GeneScoreMatrix")

Computing Module 1 of 1

Warning message in seq_len(ncol(dfM)):
“first element used of 'length.out' argument”
Error in seq_len(ncol(dfM)): argument must be coercible to non-negative integer
Traceback:

1. addModuleScore(projSpurp_doubletRemoved_filtered2, features = features,
   . useMatrix = "GeneScoreMatrix")

The issue seems to come from the very last part of the function
for (x in seq_len(ncol(dfM))) { ArchRProj <- addCellColData(ArchRProj, data = dfM[, x], name = names(featureList)[x], cells = rownames(dfM), force = TRUE) }

This is because you get back a vector for dfM rather than a matrix/dataframe if your input list for features contains only a single set of genes. It is easy to fix by also adding a dummy feature set, but is not intuitive to the user.

[rcorces]

Thanks @dagarfield - @jgranja24 will take a look.

Also, you can use permalinks to point us to the exact line of code you're referring to. Like this:

ArchR/R/ModuleScore.R

Line 97 in 2f022a4

dfM <- lapply(seq_along(featureList), function(x){

[NoemieL]

Hi. This is very nice and it works well with GeneScoreMatrix. I wanted to do the same with GeneIntegrationMatrix but it is not adapted and returns this error:

Error in addModuleScore(proj_all_3, features = c("TNFRSF17", "SDC1", "CD19", :
nBgd must be lower than 1!

Do you think it is possible to manage it?
Thanks

[rcorces]

Looking at the code, it looks like addModuleScore() may not have been coded to handle gene expression? But i agree that this is an important feature to have. @jgranja24 ?

[j-andrews7]

Note that the documentation for this function appears to be for the addImputeWeights() function currently:
https://www.archrproject.com/reference/addModuleScore.html

Would also really like to see this implemented for gene expression.

[rcorces]

@j-andrews7 - thanks for pointing that out. I fixed the documentation via 9366af8 in release_1.0.2 branch

Still looking to carve out time to implement this for gene expression.