IntegrateLayers

[Gimenagomez2230]

Hi, I have 2 cells Lines (run in different sequencing runs), exposed to 8 different compounds each, so 16 samples.
I am interested in the effect of the compounds (exposure).

My steps:

1. filter each of the 16 samples by nfeatures, mt and doublets, as independent Seurat objects.
2. merge() in one object, added one column for Sample, one for Line and one for Exposure, and then
   merged_seurat <- NormalizeData(object = merged_seurat)
   merged_seurat <- FindVariableFeatures(object =merged_seurat)
   merged_seurat <- ScaleData(merged_seurat)
   merged_seurat <- RunPCA(merged_seurat)
   merged_seurat <- merged_seurat %>%
   RunUMAP(dims = 1:20) %>%
   FindNeighbors( dims = 1:20) %>%
   FindClusters(resolution = 0.1, cluster.name = "unintegrated_clusters_0.1")

when I plot DimPlot(merged_seurat, reduction = "umap", split.by = "Line") I see that Line1 and Line 2 do not overlap much, while if I split.by = exposure they do.

Then I think I need to Integrate by Line (to remove the effect of Line/batch).
But I don't think I should integrate by Exposure (as they still share clusters, and I am afraid of overcorrecting the effect of the exposure which is what I am trying to catch). * Is my thinking correct??????

3. ways to integrate:

A - IntegrateLayers():

 I am supposed to join-split but if I split by Sample I get again the 16 layers of counts + 16 data + scale.data. ** is this redundant?
                obj <- merged_seurat
                obj[["RNA"]] <- JoinLayers(obj[["RNA"]])
                obj[["RNA"]] <- split(obj[["RNA"]], f = obj$Sample)
                obj <- NormalizeData(obj)
                obj <- FindVariableFeatures(obj)
                obj <- ScaleData(obj)
                obj <- RunPCA(obj)
                obj <- IntegrateLayers(object = obj, method = CCAIntegration,
                                       orig.reduction = "pca", new.reduction = "integrated.cca",
                                       verbose = FALSE)
  
 Then when I IntegrateLayers, as I am considering the 16 layers (samples) I think I am integrating all the samples (not just by Line but also exposure.).. which I thought I didnt want it if I just wanted to remove or integrate Line/batch effect.

Another option would be to split by Line and get 2 layers of counts + 2 data + scaledata

                obj <- merged_seurat
                obj[["RNA"]] <- JoinLayers(obj[["RNA"]])
                obj[["RNA"]] <- split(obj[["RNA"]], f = obj$Line)
                obj <- NormalizeData(obj)
                obj <- FindVariableFeatures(obj)
                obj <- ScaleData(obj)
                obj <- RunPCA(obj)
                obj <- IntegrateLayers(object = obj, method = CCAIntegration,
                   orig.reduction = "pca", new.reduction = "integrated.cca",
                   verbose = FALSE)

Then I would be integrating by Line only. which I think is more appropriate for my design but the normalization and find variables is not being done by sample but by Line. *** Is this correct?

I can do this with Harmony as well... but in any case the selection of layers to integrate depends on how I split.

B- RunHarmony():

               merged_seurat.harmony <- merged_seurat %>%
               RunHarmony(group.by.vars = 'Line', plot_convergence = FALSE)

Here I am doing the normalization and find variables in my merged_seurat which is the result of step 2 (and therefore the normalization and findvariables etc is being done by sample as they are independent layers in the merged_seurat object, as if I has joined-split by sample), but I am running Harmony group.by.vars = 'LINE' (without having split by Line the layers). So is not the same as running Harmony in with integrateLayers() after splitting by line as in approach A.

**** which Harmony approach would be correct?

C- FindIntegrationAnchors() CCA

            obj.list <- SplitObject(merged_seurat, split.by = 'Line')
            for(i in 1:length(obj.list)){
              obj.list[[i]] <- NormalizeData(object = obj.list[[i]])
              obj.list[[i]] <- FindVariableFeatures(object = obj.list[[i]])
            }
                            
            features <- SelectIntegrationFeatures(object.list = obj.list)
            anchors <- FindIntegrationAnchors(object.list = obj.list, anchor.features = features) 
            seurat.integrated <- IntegrateData(anchorset = anchors)
                            
            seurat.integrated <- ScaleData(object = seurat.integrated)
            seurat.integrated <- RunPCA(object = seurat.integrated)
            seurat.integrated <- RunUMAP(object = seurat.integrated, dims = 1:20)

Here, it seems I am splitting by Line for the integration, but the normalizeData and FindVariablesFeatures is being done for each sample.

***** would this be correct?
******this is failing in the step
FindIntegrationAnchors(object.list = obj.list, anchor.features = features)

            Error in getGlobalsAndPackages(expr, envir = envir, globals = globals) :
              The total size of the 7 globals exported for future expression (‘FUN()’) is 639.44 MiB.. This exceeds the maximum allowed size of 500.00 MiB (option 'future.globals.maxSize'). The three largest globals are ‘FUN’ (433.91 MiB of class ‘function’), ‘query’ (195.55 MiB of class ‘numeric’) and ‘index’ (9.98 MiB of class ‘S4’)

..........................................

In summary, I am very confused in how I should split the layers, if the normalization and find variables should be done by sample or by "Line-group of samples of the same batch", and how to tell each function the specific variable I want to consider for integration (e.g. Line, Exposure, all samples).

Final question for Harmony, either with RunHarmony or IntegrateLayers(HarmonyIntegration) I am getting several warnings:
9: Quick-TRANSfer stage steps exceeded maximum (= 12420850)
and with IntegrateLayers(HarmonyIntegration) also many deprecated as for example
1: HarmonyMatrix is deprecated and will be removed in the future from the API in the future
2: Warning: The parameters do_pca and npcs are deprecated. ...
... 9: Quick-TRANSfer stage steps exceeded maximum (= 12420850)

Is it ok to ignore all these warnings?

I am sorry it is a very long question... but found many forums with same concerns and no straight answer so I wanted to be clear.
Thank you so much!

Gimena

[longmanz]

Hi,
A. If you have 2 different cell lines, then it is better Not to integrate based on "cell lines" since they are biologically different from each other. But if they are the same cell lines but just with different sequencing runs, then it is fine to merge them based on "cell lines".
B / C. You do not need to run normalizeData for each sample individually. It will be the same when running it for all samples simultaneously. As for FindVariablesFeatures, you might not need to run it when considering all samples together.

As for the warning messages, I do not think they are of particular concern.

Since this is more of a analysis question, I will move it to the Discussion channel to get more people from the community to engage.