working with Seurat v5 layers?

[dpearton]

Hello,

I am working with Seurat v5 with my samples in different layers to test the different integration methods. Unfortunately the version of STACAS I am working with does not work with V5 layers. I get the following error when I try run FindAnchors.STACAS

stacas_anchors <- FindAnchors.STACAS(obj.list, 
+                                      anchor.features = nfeatures,
+                                      dims = 1:ndim)

Computing 2000 integration features
Error in `GetAssayData()`:
! GetAssayData doesn't work for multiple layers in v5 assay.
Run `rlang::last_trace()` to see where the error occurred.

`

rlang::last_trace()`

<error/ You can run 'object <- JoinLayers(object = object, layers = layer)'.>
Error in `GetAssayData()`:
! GetAssayData doesn't work for multiple layers in v5 assay.
---
Backtrace:
     ▆
  1. ├─... %>% RunUMAP(dims = 1:ndim)
  2. ├─Seurat::RunUMAP(., dims = 1:ndim)
  3. └─STACAS::Run.STACAS(., dims = 1:ndim, anchor.features = nfeatures)
  4.   └─STACAS::FindAnchors.STACAS(...)
  5.     └─base::lapply(...)
  6.       └─STACAS (local) FUN(X[[i]], ...)
  7.         └─STACAS::FindVariableFeatures.STACAS(x, nfeat = n.this, genesBlockList = genesBlockList)
  8.           ├─base::apply(...)
  9.           ├─SeuratObject::GetAssayData(obj, assay = assay, slot = "data")
 10.           └─SeuratObject:::GetAssayData.Seurat(obj, assay = assay, slot = "data")
 11.             ├─SeuratObject::GetAssayData(object = object[[assay]], layer = layer)
 12.             └─SeuratObject:::GetAssayData.StdAssay(object = object[[assay]], layer = layer)
Run rlang::last_trace(drop = FALSE) to see 1 hidden frame.

Is there any chance of adding support for layers in STACAS?

Thank you

[mass-a]

Hello!
STACAS is made to work on a list of objects, where each object contains one sample. It should be sufficient to split out your object by sample (maybe after merging the layers, if your samples are organized in layers). For example:

pbmcsca.integrated <- NormalizeData(pbmcsca) |>
    SplitObject(split.by = "Method")|>
    Run.STACAS()

I understand that Seurat have changed their integration methods to use single objects in a layer structure, but you will need a list of objects in STACAS - at least for now.

[dpearton]

Hi, Thank you for the very prompt reply. The issue appears to be that I was trying to run STACAS based on combined samples (i.e. trying to integrate two different sequencing runs each comprised of multiple samples) rather than individual samples. If I split by sample (orig.ident) it works. I am getting an error when I FindAnchors.STACAS (for each sample). Warning: Different features in new layer data than already exists for scale.data Will this be an issue? - the integration does proceed after these warnings. Thanks again for your time. Yours, David

…

________________________________ From: Massimo Andreatta ***@***.***> Sent: 08 December 2023 10:33 To: carmonalab/STACAS ***@***.***> Cc: David Pearton ***@***.***>; Author ***@***.***> Subject: Re: [carmonalab/STACAS] working with Seurat v5 layers? (Issue #30) You don't often get email from ***@***.*** Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> Hello! STACAS is made to work on a list of objects, where each object contains one sample. It should be sufficient to split out your object by sample (maybe after merging the layers, if your samples are organized in layers). For example: pbmcsca.integrated <- NormalizeData(pbmcsca) |> SplitObject(split.by = "Method")|> Run.STACAS() I understand that Seurat have changed their integration methods to use single objects in a layer structure, but you will need a list of objects in STACAS - at least for now. — Reply to this email directly, view it on GitHub<#30 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BCYA63N2JBMKOHYSZVCK5EDYILUI7AVCNFSM6AAAAABAMNGQBWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNBWHEZTENRSG4>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[mass-a]

Good that you got it to work.
As for the warning you posted, it's probably safe to ignore it. My guess is that there is already a scale.data layer in the object, and it's being overwritten by using a different set of features (i.e. the highly variable genes calculated by STACAS).

[domi84]

Hello,
thanks for this great tool!
Just to update that now Run.STACAS is getting a warning with Seurat 5.3, that probably will be transformed in an error in the future:

WST.semisup <- list.WST %>% Run.STACAS(dims = 1:ndim, anchor.features = 1500, cell.labels = "clusters")

Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
ℹ Please use the `layer` argument instead.
ℹ The deprecated feature was likely used in the Seurat package.
  Please report the issue at <https://github.com/satijalab/seurat/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

[leonardHerault]

Hi,

I also ran into this issue recently and I wrote this function that can be used with IntegrateLayers seurat v5 funtion.
It's based on the code of the RPCAIntegration function in the seurat v5 package:

STACASIntegration <- function(object = NULL, assay = NULL, layers = NULL, orig = NULL, 
                               new.reduction = "integrated.dr", reference = NULL, features = NULL, 
                               normalization.method = c("LogNormalize"), dims = 1:30, scale.data = FALSE,
                               k.filter = NA, scale.layer = "scale.data", dims.to.integrate = NULL, 
                               k.weight = 100, weight.reduction = NULL, sd.weight = 1, sample.tree = NULL, 
                               preserve.order = FALSE, verbose = TRUE,labels.df = NULL,anchor.filtering.mtx=NULL, ...) 
{
  
  op <- options(Seurat.object.assay.version = "v3", Seurat.object.assay.calcn = FALSE)
  on.exit(expr = options(op), add = TRUE)
  normalization.method <- match.arg(arg = normalization.method)
  features <- features %||% SelectIntegrationFeatures5(object = object)
  assay <- assay %||% "RNA"
  layers <- layers %||% Layers(object = object, search = "data")
  ncells <- sapply(X = layers, FUN = function(x) {
    ncell <- dim(object[x])[2]
    return(ncell)
  })
  if (min(ncells) < max(dims)) {
    abort(message = "At least one layer has fewer cells than dimensions specified, please lower 'dims' accordingly.")
  }
  
  object.list <- list()
  for (i in seq_along(along.with = layers)) {
    object.list[[i]] <- suppressMessages(suppressWarnings(CreateSeuratObject(counts = NULL, 
                                                                             data = object[layers[i]][features, ])))
    VariableFeatures(object = object.list[[i]]) <- features
    object.list[[i]] <- AddMetaData(object.list[[i]],metadata = labels.df[colnames(object.list[[i]]),],col.name = colnames(labels.df))
  }
  
  
  message("Finding anchors using STACAS...")
  
  stacas.anchors <- FindAnchors.STACAS(object.list = object.list, assay = rep(assay,length(object.list)), 
                                       anchor.features = features, reference = reference, 
                                       min.sample.size = 0, max.seed.objects = length(object.list), 
                                       genesBlockList = NULL, dims = dims, 
                                       cell.labels = colnames(labels.df)[1], 
                                       scale.data = scale.data, 
                                       verbose = verbose,...)
  
  if (verbose) {
    print("initial anchors:")
    print((stacas.anchors))
  }
  
  slot(object = stacas.anchors, name = "object.list") <- lapply(X = slot(object = stacas.anchors, 
                                                                         name = "object.list"), FUN = function(x) {
                                                                           suppressWarnings(expr = x <- DietSeurat(x, features = features[1:2]))
                                                                           return(x)
                                                                         })
  
  unsup.sample.tree <- SampleTree.STACAS(
    anchorset = stacas.anchors,
    obj.names = names(object.list),
    semisupervised = F
  )
  
  if (verbose) {
    print("Sample tree:")
    print((unsup.sample.tree))
  }
  
  object_merged <- IntegrateEmbeddings(anchorset = stacas.anchors,
                                       reductions = orig, new.reduction.name = paste0("unsup_",new.reduction),
                                       dims.to.integrate = dims.to.integrate, k.weight = k.weight,
                                       weight.reduction = weight.reduction, sd.weight = sd.weight,
                                       sample.tree = unsup.sample.tree, preserve.order = preserve.order,
                                       verbose = verbose)
  
  output.list <- list(object_merged[[paste0("unsup_",new.reduction)]])
  names(output.list)<- paste0("unsup_",new.reduction) 
  
  remove(object_merged)
  gc(verbose = F)
  
  message("unsupervised integration done")
  
  ## semi sup integration with anchor filtering
  
  if (is.null(anchor.filtering.mtx)) {
    anchor.filtering.mtx <- matrix(F,ncol = length(unique(labels.df[,1])),
                                   nrow =  length(unique(labels.df[,1])))
    colnames(anchor.filtering.mtx) <- unique(labels.df[,1])
    rownames(anchor.filtering.mtx) <- unique(labels.df[,1])
    
    diag(anchor.filtering.mtx) <- T
    anchor.filtering.mtx["unknown",] <- T
    anchor.filtering.mtx[,"unknown"] <- T
  }
  
  if (verbose) {
    message("anchor filtering matrix:")
    print(anchor.filtering.mtx)
  }
  
  filtered.stacas.anchors <- subset(stacas.anchors, group.by = colnames(labels.df)[1], ident.matrix = anchor.filtering.mtx)
  
  if (verbose) {
    print("filtered anchors:")
    print((filtered.stacas.anchors))
  }
  
  remove(stacas.anchors)
  gc(verbose = F)
  
  
  
  slot(object = filtered.stacas.anchors, name = "object.list") <- lapply(X = slot(object = filtered.stacas.anchors,
                                                                                  name = "object.list"), FUN = function(x) {
                                                                                    suppressWarnings(expr = x <- DietSeurat(x, features = features[1:2]))
                                                                                    return(x)
                                                                                  })
  
  semisup.sample.tree <- SampleTree.STACAS(
    anchorset = filtered.stacas.anchors,
    obj.names = names(object.list),
    semisupervised = T
  )
  
  object_merged <- IntegrateEmbeddings(anchorset = filtered.stacas.anchors,
                                       reductions = orig, new.reduction.name = paste0("semisup_",new.reduction),
                                       dims.to.integrate = dims.to.integrate, k.weight = k.weight,
                                       weight.reduction = weight.reduction, sd.weight = sd.weight,
                                       sample.tree = semisup.sample.tree, preserve.order = preserve.order,
                                       verbose = verbose)
  
  output.list[[paste0("semisup_",new.reduction)]] <- object_merged[[paste0("semisup_",new.reduction)]]
  
  
  return(output.list)
}

Then sobj being a splitted seurat v5 object, it can be integrated with STACAS as follow:

sobj[["RNA"]] <- split(sobj[["RNA"]], f = sobj[["your_sample_batch_col"][,1]) # Be carefull it doesn't preserve sample "appearing" order in the metadata compared to SplitObject

sobj <- NormalizeData(sobj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()

sobj <- IntegrateLayers(
  object = sobj, 
  method = STACASIntegration,
  orig.reduction = "pca", 
  new.reduction = "STACAS", 
  dims = 1:50,
  verbose = T,
  labels.df = sobj[["your_guiding_label_metadata_column"]]     #  "unknown" for unlabelled cells free to pair with every cells types
)

Note that with this function anchor filtering is performed using the seurat subset method for class 'AnchorSet' using the keyword "unknown" for unlabelled cells.

Another important difference with standard STACAS workflow is the use of seurat IntegrateEmbeddings function correcting directly an unintegrated dimension reduction instead of IntegrateData function correcting the scale count matrix.
Coupled with seurat v5 layer management it allows a significant reduction in cpu time and memory usage.

It''s performing unsupervised and (semi)supervised integration in a row to save the time of anchor findings but an option to perform only one or the other could easily be added.

If you want to use reference samples be careful that Seurat split function doesn't not conserve the order compared to the previous SplitObject function.

Thus the correct indices to use for the reference argument in IntegrateLayers can be retrieved like this:

layer.data.names <- stringr::str_split_fixed(SeuratObject::Layers(object = sobj, search = "data"),pattern = "\\.",Inf)[,2]
ref.index.layer <- match(reference_sample_name_vector, layer.data.names)

sobj <- IntegrateLayers(
  object = sobj, 
  method = STACASIntegration,
  orig.reduction = "pca", 
  new.reduction = "STACAS", 
  dims = 1:50,
  verbose = T,
  reference = ref.index.layer,      
  labels.df = sobj[["your_guiding_label_metadata_column"]]     #  "unknown" for unlabelled cells free to pair with every cells
)

@mass-a, @sjcarmona let me know what you think about it and if this could be more formally added to the STACAS package!

Best,

Léonard

[mass-a]

Hey @leonardHerault
thanks a lot!

I will look into how we can integrate these ideas into STACAS. Of course, if you are motivated, feel free to do a pull request with your suggested changes.

-m