Difference between vst_out$y and get_residuals()

[frederikziebell]

Hi Christoph,

what is the difference between what is stored in the y slot of the return from vst() and what get_residuals() returns? From the documentation, vst_out$y contains pearson residuals but get_residuals() internally builds an additional regression model?

I noted that vst() and Seurat::SCTransform() give different results, potentially because SCTransform() internally calls get_residuals(). Here's an example:

pbmc_transf <- CreateSeuratObject(pbmc)
pbmc_transf <- SCTransform(pbmc_transf,do.center=F)

as.matrix(pbmc_transf@assays$SCT@data) %>% rowMeans() %>% qplot(bins=100) + xlim(-1,5) + ggtitle("Seurat data")
as.matrix(pbmc_transf@assays$SCT@scale.data) %>% rowMeans() %>% qplot(bins=100) + xlim(-1,5) + ggtitle("Seurat scale.data")
as.matrix(sctransform::vst(pbmc)$y) %>% rowMeans() %>% qplot(bins=100) + xlim(-1,5) + ggtitle("vst_out$y")

Could you please clarify (1) why get_residuals() is not simply a call to the y slot of the vst() output? (2) if get_residuals() can be used to reproduce the scale.data slot in Seurat (3) and maybe why the data slot in Seurat has only positive results? I actually thought that this would already contain pearson residuals.

Thank you!

[ChristophH]

There is no difference between the y slot of the return list of vst() and the output of get_residuals()

> vst_out <- vst(pbmc, return_cell_attr = TRUE)
> vst_residuals <- get_residuals(vst_out)
> identical(vst$y, vst_residuals)
[1] TRUE

Seurat::SCTransform applies a different clipping to the residuals and also centers the genes. Check the documentation of Seurat::SCTransform for details.

For more information what the slots in a Seurat Assay are, have a look at the Assay documentation (?Seurat::Assay)

[frederikziebell]

Thank you for the clarification, I haven't spotted the different clipping.