Understanding scores and significance

[alexwskh]

Hello and thank you for your hard work developing a great tool.

I see there are a few other posts examining this question but I still haven't found my answer, though apologies if I have missed it. I am trying to understand the context of these scores. I understand that you can look at the relative enrichment of scores between clusters in your dataset, so you can easily make comparative observations (the score for x signature is statistically significantly higher in cluster A then cluster B).

I am trying to understand the absolute measure of the score though... In standard GSEA, the metric I normally look at is the FDR.. That is to say, can I be confident that this gene set is actually enriched in the top set of genes as opposed to random chance. I'm struggling to understand the equivalent here. I can say that cluster A has an average of 0.10 versus cluster B having an average of 0.05 for UCell score. Does an average UCell score of 0.10 represent a significant (defined as majorly not due to chance) enrichment?

Thanks for any assistance here.

[mass-a]

Hi @alexwskh!
thanks for the question.

As you said, the most common use of these signature scores is for comparison between groups of cells or clusters. But I agree that in some cases an absolute interpretation of UCell scores "a-la-GSEA" could be useful.

I would address this question in a numerical way: for a given signature S of length N, generate 1,000 random signatures of the same length N. Then draw the distribution of average scores over the 1,000 background sets, and ask where the observed UCell score sits in this distribution. The fraction of random scores that are higher than the observed score can be thought of as a p-value for the observation not being part of the random distribution.

Note that this would be a sample-wise (or cluster-wise) test for a signature's UCell score, rather than a cell-wise test. Does it make sense? otherwise I can try to make some sample code.

Best
-m

[alexwskh]

That does make sense, thanks for taking the time to answer. I would happily take some sample code if you were interested, since you're certainly much more proficient then I am! Just to parrot and make sure I understand...

You have your cluster of interest. You have your mean UCell score for that cluster. Since UCell scoring is independent-per-cell, if you wanted to test a specific cluster, you could just subset those cells and run UCell for your 1000 random permutations of length N genesets. If I find that 95% of randomly generated gene sets at that size are have a mean UCell score lower than my test set, I can equate that to a p value of 0.05 (roughly).

Is UCell particularly sensitive to gene set size? I know this is a metric covered in some comparative reviews of bioinformatic tools, but I'm not really illuminated in my understanding there yet. I have fairly heterogenous integrated scRNA set from 7 different conditions, and I'm looking at UCell over something like ssGSEA/GSVA since from my reading it'll handle dropout events better...

[mass-a]

Sure, here's some example code that may help.

library(UCell)
library(Seurat)

#get data
options(timeout = 3000)
url <- "https://figshare.com/ndownloader/files/54054143"
download.file(url, destfile = "data/pbmc_multimodal.Tcells_ds9k.rds")
pbmc.Tcell <- readRDS("data/pbmc_multimodal.Tcells_ds9k.rds")

#define a signature to test
sig.obs <- list(CD4_Tcell = c("CD4","CD40LG"))

#generate 1000 random signatures of sample size
nrep <- 1000
lgt <- length(sig.obs$CD4_Tcell)

genes_all <- rownames(pbmc.Tcell)
rand <- lapply(seq(1,nrep), function(x) {
  sample(genes_all, lgt, replace = FALSE)
})
names(rand) <- paste0("r", seq(1, length(rand)))

Calculate UCell scores for target and background

sig.all <- c(sig.obs, rand)
pbmc.Tcell  <- AddModuleScore_UCell(pbmc.Tcell, features = sig.all, ncores = 4, name = NULL)

Is the target significantly higher than BG?

means <- apply(pbmc.Tcell@meta.data[,names(rand)], 2, mean)
obs <- mean(pbmc.Tcell@meta.data[,names(sig.obs)])
hist(means, breaks=100)
abline(v=obs, col="red")

#estimate p-value
p <- 1 - sum(obs > means) / length(means)
p

[1] 0.11

It looks like this signature for CD4 T cells is not significantly higher than other random 2-gene signatures, at least on the whole dataset.
What if we do the test in a subset of CD4 T cells only?

#same test, but for a particular cluster
cd4t <- subset(pbmc.Tcell, subset=celltype.l2=="CD4 TEM")
means <- apply(cd4t@meta.data[,names(rand)], 2, mean)
obs <- mean(cd4t@meta.data[,names(sig.obs)])

hist(means, breaks=100)
abline(v=obs, col="red")

p <- 1 - sum(obs > means) / length(means)
p

[1] 0.052

UCell scores shouldn't be too sensitive to gene set size. However, small gene sets are bound to be noisier, and that could especially affect the background distributions (in the example above, it's more likely to strike on a random 2-gene combinations where both genes are highly expressed, compared to e.g. a 10-gene signature).

I hope this helps!
-massimo

[alexwskh]

Helps a lot! Thanks Massimo, that's really kind of you to take time out of your day and write this up for me step by step. I look forward to trying this out on my own data.

I have another question for you, should you feel inclined to answer. Given an integrated dataset where the feature space for each contributing member can vary quite heavily, should I be considering an alternative value for rmax other than 1500? Below I have the detected genes for each of my sets, and as you can see some are very sparse, while some are quite deep.

[mass-a]

No problem! this question came up before, so perhaps the example can be useful to others.

Yes you could consider lower values for the maxRank parameter, seeing that some of your datasets seem to have low sequencing depth. One rule of thumb could be the median number of detected genes of the shallowest dataset?
That said, it would also be important to have a good QC on the data and set some minimal criteria for including cells (min number of genes, counts per cells, mitochondrial content, etc.).

[alexwskh]

Thanks you for your suggestions, I'll play around with them. Yeah, always nice to stumble across the forum post addressing your current problem.

If you don't mind me asking, what kind of QC metrics do you usually go with? I've kind of started from nothing and have been learning as I go along, but it seems like QC parameters vary pretty wildly between labs. My current threshold is quite generous I suppose... I keep cells with >1000 counts <50000 counts, and < 10% mitochondrial content. I didn't have a hard gene threshold, but I only have a single cell in this subset that has <200 genes. I am also running ambient RNA removal with decontX and then doublet finder. I do all of this on each individual dataset before SCTransforming the data and integrating with Seurat.

[GWW]

This approach is really interesting. Would a similar cell-wise approach be viable? One issue I have when quantifying signatures across samples is deciding if the signature is expressed in that sample or what would be a reasonable cutoff for real expression versus random noise.

I was thinking along the lines of combining what Seurat does by binning genes by their expression and then for each signature:

Pick a random set of genes from the same bins, calculate the ucell score for each cell, and then repeat this a bunch of times and then calculate a p-value similar to what you suggested above.