Metafeatures of a seurat object: cannot recompute the exact values of vst.variance.expected and vst.variance.standardized

[gael-millot]

Hello.

Using the pbmc seurat object, available here :
https://dl.pasteur.fr/fop/bRtuYjMS/pbmc.rds

I have no problem recomputing vst.mean and vst.variance, but, concerning vst.variance.expected and vst.variance.standardized, this is the closest but not identical I could find:

seu_obj <- readRDS("C:\\Users\\gmillot\\Documents\\Git_projects\\seurat_b_cells\\dataset\\pbmc\\pbmc.rds")
if(base::all(base::class(seu_obj@assays$RNA) == "Assay")) {
  metafeat_df <- as.data.frame(seu_obj@[email protected])
  counts_mat <- as.matrix(seu_obj@assays$RNA@counts) # [rownames(seurat_scale), colnames(seurat_scale)]
}else if (base::all(base::class(seu_obj@assays$RNA) == "Assay5")) {
  metafeat_df <- as.data.frame(seu_obj@[email protected])
  counts_mat <- as.matrix(seu_obj@assays$RNA@layers$counts) # [rownames(seurat_scale), colnames(seurat_scale)]
}

vst.mean <- rowMeans(counts_mat)
vst.variance <- apply(counts_mat, 1, var)

valid <- vst.mean > 0
fit <- loess(log10(vst.variance[valid]) ~ log10(vst.mean[valid])
)

# Predict expected variance for all genes (including those with mean == 0, will be NA)
predicted_log10 <- predict(fit, log10(vst.mean))
vst.variance.expected <- 10^predicted_log10

# Standardize (as ratio, not z-score)
vst.variance.standardized <- vst.variance / vst.variance.expected

# Assemble
metafeat_df_recomputed <- data.frame(
  vst.mean = vst.mean,
  vst.variance = vst.variance,
  vst.variance.expected = vst.variance.expected,
  vst.variance.standardized = vst.variance.standardized
)
head(metafeat_df_recomputed)
head(metafeat_df)

I really need to understand how the values used in differential analysis are computed. The documentation I found and ChatGPT are not efficient.
Any idea ?

Thanks !

Gael.