Problem with running IntegrateData -- failing to select a single reference model

[tcikezu]

I ran the integration pipeline, clustered, and then found a subset that I wanted to analyze further. I re-ran the integration pipeline with this subset (beginning from the RNA assay, recomputed SCT assay, ...), and found this error while calling IntegrateData:

Error in model.list[[which(reference.model)]] :
  attempt to select less than one element in get1index
Calls: IntegrateData
Execution halted

The error is produced if we set reference=1 while calling IntegrateData on a previously integrated dataset. Relevant lines are from integration.R (lines 1094 - 1119)

  if (normalization.method == "SCT") {
    if (is.null(x = Tool(object = reference.integrated, slot = "Integration"))) {
      reference.sample <- slot(object = anchorset, name = "reference.objects")
    } else {
      reference.sample <- SampleIntegrationOrder(
        tree = slot(
          object = reference.integrated,
          name = "tools"
        )$Integration@sample.tree
      )[1]
    }
    reference.cells <- Cells(x = object.list[[reference.sample]])
    reference.model <- NULL
    if (length(x = model.list) > 0) {
      reference.model <- sapply(X = model.list, FUN = function(model) {
        reference.check <- FALSE
        model.cells <- Cells(x = model)
        if (length(x = model.cells) > 0 &
            length(x = setdiff(x = model.cells, y = reference.cells)) == 0) {
          reference.check <- TRUE
        }
        return(reference.check)
        }
      )
      reference.model <- model.list[[which(reference.model)]]
    }
  }

The issue is Tool(object = reference.integrated, slot = "Integration")) does not return NULL, even if reference.integrated is the result of running PairwiseIntegrateReference on an anchorset withreference.object slot of length 1. (The Integration Tool slot probably remains from the previous integration.) Is there an easy way to remove this slot from just the reference model in case it exists? If that can be done in PairwiseIntegrateReference then I think this issue can be resolved.

[yuhanH]

HI, @tcikezu
Would you mind to show the scripts you used for the integration, subset and re-integration?
Thanks.

[tcikezu]

Hi @yuhanH, here is a script to replicate the error:

library(Seurat)
library(sctransform)
library(future.apply)
library(SeuratData)

plan(multicore, workers=4)

data(pbmc3k)
pbmc3k <- pbmc3k[sample(1:nrow(pbmc3k), 6000), sample(1:ncol(pbmc3k), 1000)]
groups <- sample(c("g1", "g2", "g3", "g4"), size = ncol(pbmc3k), replace = TRUE)
names(groups) <- colnames(pbmc3k)
pbmc3k <- AddMetaData(object = pbmc3k, metadata = groups, col.name = "group")
pbmc3k.list <- SplitObject(pbmc3k, split.by = "group")

pbmc3k.sc.list <- future_lapply(pbmc3k.list, FUN=function(x) SCTransform(x))
pbmc3k.features <- SelectIntegrationFeatures(pbmc3k.sc.list, nfeatures=2000)
pbmc3k.prepared.list <- PrepSCTIntegration(pbmc3k.sc.list, anchor.features = pbmc3k.features)
pbmc3k.anchors <- FindIntegrationAnchors(pbmc3k.prepared.list,
                                         anchor.features = pbmc3k.features,
                                         normalization.method='SCT',
                                         dims=1:30,
                                         reference=1)

pbmc3k.int <- IntegrateData(anchorset = pbmc3k.anchors, normalization.method="SCT", verbose=T)
pbmc3k.int.subset <- subset(pbmc3k.int, subset = group %in% c('g1', 'g2', 'g3'))
pbmc3k.int.subset.list <- SplitObject(pbmc3k.int.subset, split.by = 'group')
pbmc3k.int.subset.sc.list <- future_lapply(pbmc3k.int.subset.list, FUN = function(x) SCTransform(x))
pbmc3k.int.subset.features <- SelectIntegrationFeatures(pbmc3k.int.subset.sc.list, nfeatures=2000)
pbmc3k.int.subset.prepared.list <- PrepSCTIntegration(pbmc3k.int.subset.sc.list, anchor.features = pbmc3k.int.subset.features)
pbmc3k.int.subset.anchors <- FindIntegrationAnchors(pbmc3k.int.subset.prepared.list, anchor.features = pbmc3k.int.subset.features, normalization.method = 'SCT', dims=1:30, reference = 1)
pbmc3k.int.subset.int <- IntegrateData(anchorset=pbmc3k.int.subset.anchors, normalization.method="SCT", verbose=T)

Many thanks for looking into this.

[yuhanH]

Hi,
What is the default assay for pbmc3k.int.subset? If it is not RNA, it may cause some error in the downstream.

[tcikezu]

Hi @yuhanH -- default assay used was 'integrated'. I re-tried again, this time adding this line:

DefaultAssay(pbmc3k.int.subset) <- 'RNA'

before running SplitObject. Same error still occurs.

[yuhanH]

A safest strategy is that you create an object only containing RNA assay of your subset cells.
But this is a bug. Thanks for reporting this. We will fix soon.

[tcikezu]

Okay, thank you!

[yuhanH]

@tcikezu It is fixed in the develop branch. Thanks for reporting this.

[BenjaminDEMAILLE]

Hi ! it's not fixed. I had the same error with Seurat 4.2

features <- SelectIntegrationFeatures(object.list = Plasticity_spatial_immune_split, 
                                      nfeatures = 3000)

Plasticity_spatial_immune_split <- PrepSCTIntegration(object.list = Plasticity_spatial_immune_split, 
                                                      anchor.features = features)

Plasticity_spatial_immune_split <- lapply(X = Plasticity_spatial_immune_split, 
                                          FUN = RunPCA, 
                                          features = features)

anchors <- FindIntegrationAnchors(object.list = Plasticity_spatial_immune_split, 
                                  normalization.method = "SCT",
                                  anchor.features = features, 
                                  dims = 1:30, reduction = "rpca", 
                                  k.anchor = 20)

Plasticity_spatial.combined.sct<- IntegrateData(anchorset = anchors, 
                                                normalization.method = "SCT", 
                                                dims = 1:30)

Merging dataset 1 into 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data
Error in model.list[[which(reference.model)]] :
recursive indexing failed at level 2