Hi, I've been using
seurat-wrappers (https://github.com/satijalab/seurat-wrappers) to align my datasets using MNN but I have some issues with downstream analysis. I'm having problem with projecting my dimension reduction onto the full dataset. I realized that the output of
RunFastMNN is a merged
Seurat object with the MNN dimensional reduction information but without the feature-level metadata such as variable feature information. (see issue here: https://github.com/satijalab/seurat-wrappers/issues/15)
However, it is nice to examine the projection of the data with this reduction and I would like to retain the information in my MNN-aligned
Seurat object. I worked around the issue by first performing MNN correction with
batchelor then converting it into a
Seurat object, then save the highly variable genes list into the variable feature slot in the
I wonder if what I did is correct & logical. Did anyone experience the same issue?
My code is as below. Do my steps seem logical?
Thank you very much.
## Load Seurat object so <- readRDS(file = paste0(output, "/PBMC/SO_merge.Rds")) ## Create SingleCellExperiment object sce <- as.SingleCellExperiment(so) rowData(sce) <- NULL reducedDim(sce) <- NULL reducedDim(sce, type = "UMAP") <- NULL ## Correct by sample ID s11 <- sce[ , grepl("S11", sce$orig.ident)] s12 <- sce[ , grepl("S12", sce$orig.ident)] s13 <- sce[ , grepl("S13", sce$orig.ident)] s14 <- sce[ , grepl("S14", sce$orig.ident)] s15 <- sce[ , grepl("S15", sce$orig.ident)] s16 <- sce[ , grepl("S16", sce$orig.ident)] all.sce <- list(S11 = s11, S12 = s12, S13 = s13, S14 = s14, S15 = s15, S16 = s16) ## Subset all batches to common universe of genes universe <- Reduce(intersect, lapply(all.sce, rownames)) all.sce <- lapply(all.sce, "[", i = universe,) ## Adjust scaling to equalize sequencing coverage normed.sce <- do.call(multiBatchNorm, all.sce) ## Find highly variable genes all.var <- lapply(all.sce, modelGeneVar) combined.var <- do.call(combineVar, all.var) hvg.list <- rownames(combined.var)[combined.var$bio > 0] ## Correct batch effect set.seed(920101) mnn.sce <- do.call(fastMNN, c(normed.sce, list(subset.row = hvg.list))) ## Save computed MNN into SCE object, then convert to Seurat object reducedDim(sce, "MNN") <- reducedDim(mnn.sce, "corrected") so.fastmnn <- as.Seurat(sce) ## Keep highly variable genes list into Seurat object so.fastmnn@assays$RNA@var.features <- hvg.list ## Scale data & project loadings so.fastmnn <- ScaleData(so.fastmnn) ProjectDim(so.fastmnn, reduction = "mnn", dims.print = 1:2, nfeatures.print = 5)