satijalab / seurat

R toolkit for single cell genomics
http://www.satijalab.org/seurat
Other
2.29k stars 914 forks source link

Integrating multiple objects with IntegrateEmbeddings() erases all previous reductions #8917

Open hannanw opened 5 months ago

hannanw commented 5 months ago

I have four different 10x multiome seurat object (one for each experiment). I am able to integrate the RNA assay without much issue.

merge_obj <- merge(expt1, c(expt2,expt3,expt4))
DefaultAssay(merge_obj) <- "SCT"
obj.list <- SplitObject(merge_obj, split.by = "orig.ident")
obj.features <- SelectIntegrationFeatures(obj.list, nfeatures = 3000)
VariableFeatures(merge_obj[["SCT"]]) <- obj.features
merge_obj <- RunPCA(merge_obj, dims=1:50, verbose = FALSE)
merge_obj_integration <- IntegrateLayers(object = merge_obj, method = CCAIntegration, normalization.method = "SCT", verbose = F)
merge_obj_integration <- RunPCA(merge_obj_integration, dims=1:50, verbose = FALSE)
merge_obj_integration <- FindNeighbors(merge_obj_integration, reduction = "integrated.dr", dims = 1:16)
merge_obj_integration <- FindClusters(merge_obj_integration, res=0.2,method='igraph',algorithm=4)

>merge_obj_integration
An object of class Seurat
243934 features across 23666 samples within 3 assays
Active assay: SCT (23167 features, 2960 variable features)
 3 layers present: counts, data, scale.data
 2 other assays present: RNA, ATAC
 2 dimensional reductions calculated: pca, integrated.dr

However, when I try to do the same integration for the ATAC assay IntegrateEmbeddings() removes all prior reductions. I followed the step from the Signac scATAC-seq data integration here.

DefaultAssay(merge_obj_integration) <- "ATAC"
merge_obj_integration <- FindTopFeatures(merge_obj_integration, min.cutoff = 10)
merge_obj_integration <- RunTFIDF(merge_obj_integration)
merge_obj_integration <- RunSVD(merge_obj_integration)
merge_obj_integration <- RunUMAP(merge_obj_integration, reduction = "lsi", dims = 2:50)

> merge_obj_integration
An object of class Seurat
243934 features across 23666 samples within 3 assays
Active assay: ATAC (184166 features, 183851 variable features)
 2 layers present: counts, data
 2 other assays present: RNA, SCT
 5 dimensional reductions calculated: pca, integrated.dr, integrated_RNA_umap, lsi, umap

atac.obj.list <- SplitObject(merge_obj_integration, split.by = "orig.ident")

integration.anchors <- FindIntegrationAnchors(
  object.list = atac.obj.list,
  anchor.features = rownames(merge_obj_integration),
  reduction = "rlsi",
  dims = 2:50
)

merge_obj_integration_ATAC <- IntegrateEmbeddings(
  anchorset = integration.anchors,
  reductions = merge_obj_integration[["lsi"]],
  k.weight = 50,
  new.reduction.name = "integrated_lsi",
  dims.to.integrate = 2:50
)

> merge_obj_integration_ATAC
An object of class Seurat
243934 features across 23666 samples within 3 assays
Active assay: ATAC (184166 features, 0 variable features)
 2 layers present: counts, data
 2 other assays present: RNA, SCT
 1 dimensional reduction calculated: integrated_lsi

Am I missing or doing something wrong? Any help/advice would be greatly appreciated!

> sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 18.04.6 LTS

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] ROCR_1.0-11                       KernSmooth_2.23-22
 [3] sctransform_0.4.1                 fields_15.2
 [5] viridisLite_0.4.2                 spam_2.10-0
 [7] leiden_0.4.3.1                    reticulate_1.35.0
 [9] scDblFinder_1.16.0                SingleCellExperiment_1.24.0
[11] SummarizedExperiment_1.32.0       MatrixGenerics_1.14.0
[13] matrixStats_1.2.0                 DoubletFinder_2.0.4
[15] glmGamPoi_1.14.3                  hdf5r_1.3.9
[17] BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.70.2
[19] rtracklayer_1.62.0                BiocIO_1.12.0
[21] Biostrings_2.70.2                 XVector_0.42.0
[23] EnsDb.Hsapiens.v86_2.99.0         ensembldb_2.26.0
[25] AnnotationFilter_1.26.0           GenomicFeatures_1.54.3
[27] AnnotationDbi_1.64.1              Biobase_2.62.0
[29] GenomicRanges_1.54.1              GenomeInfoDb_1.38.6
[31] IRanges_2.36.0                    S4Vectors_0.40.2
[33] BiocGenerics_0.48.1               ggforce_0.4.2
[35] qlcMatrix_0.9.7                   sparsesvd_0.2-2
[37] slam_0.1-50                       Matrix_1.6-5
[39] Signac_1.12.0                     ggplot2_3.4.4
[41] dplyr_1.1.4                       SeuratWrappers_0.3.4
[43] Seurat_5.0.1                      SeuratObject_5.0.1
[45] sp_2.1-3
IzzyCoales commented 2 months ago

Hi! Has this been resolved? I am currently coming across the same issue.