SCTransform layer by layer - unexpected results

[RobinPranter]

Thank you for an excellent package!

I have a question about SCTransform that has been brewing for a while but I have not been able to find an answer.

I typically investigate if any normalization makes sense by plotting the normalized vs un-normalized and expecting something like this:

This particular plot was generated using this code:

P01 <- cbind(Matrix::rowSums(green.SO[["RNA"]]$counts)[
                                        rownames(green.SO[["RNA"]]$counts) %in% 
                                          rownames(green.SO[["SCT"]]$data)],
             Matrix::rowSums(green.SO[["SCT"]]$data)) %>%
  as.data.frame() %>%
  dplyr::rename(raw = V1, norm = V2) %>%
  ggplot(aes(raw, norm)) +
  geom_point() +
  labs(title = "Normalized data vs raw counts (each dot is a gene)") +
  theme(text = element_text(size = 6))
P01

Where green.SO is a seurat object that was normalized using SCTransform.

However, I have come to understand that if I want to do integration (eg harmony) across my different samples it is preferable to do the normalization separately sample by sample rather than to normalize across the entire dataset. For example like this:

# Normalize each sample separately
full.SO[["RNA"]] <- split(full.SO[["RNA"]], f = full.SO$orig.ident)
DefaultAssay(full.SO) <- "RNA"
full.SO <- SCTransform(full.SO,
                       assay = "RNA",
                       new.assay.name = "SCT_ByPlate",
                       vars.to.regress =c("percent.mt", "G2M.Score", "S.Score"),
                       vst.flavor = "v2",
                       variable.features.n = 3000)
# Calculate PCA
full.SO <- RunPCA(full.SO, assay = "SCT_ByPlate", 
                  reduction.name = "pca_SCT_ByPlate")

# Integrate data sets
full.SO <- IntegrateLayers(object = full.SO, 
                           method = HarmonyIntegration,
                           orig.reduction = "pca_SCT_ByPlate", 
                           new.reduction = 'harmony',
                           assay = "SCT_ByPlate", 
                           verbose = FALSE)

# Calc umap
full.SO <- RunUMAP(full.SO, dims = 1:30, reduction = "harmony", 
                   reduction.name = "umap_harmony")

When I plot the same kind of plot it look nothing like what I would expect:

layer.list <- full.SO@meta.data$orig.ident %>% unique()

layer <- subset(full.SO, subset = orig.ident == layer.list[3])

cbind(Matrix::rowSums(layer[["RNA"]]$counts)[
                                    rownames(layer[["RNA"]]$counts) %in% 
                                      rownames(layer[["SCT_ByPlate"]]$data)],
             Matrix::rowSums(layer[["SCT_ByPlate"]]$data)) %>%
  as.data.frame() %>%
  dplyr::rename(raw = V1, norm = V2) %>%
  ggplot(aes(raw, norm)) +
  geom_point() +
  labs(title = "Normalized data vs raw counts (each dot is a gene)") +
  theme(text = element_text(size = 6))

How can this be? Has any one else encountered this? Am I perhaps implementing the calls incorrectly, just failing in my attempt to visualize the outcome, or is this potentially a problem? Perhaps I am missing the point somehow?

I would be very grateful for any input from the community or the developers!

Best,
Robin