Integration by Different Methods-rpca-CCA-Harmony

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
load("SS_merged_marie_obj-2-Variation.Robj")

All_samples_Merged

An object of class Seurat 
36629 features across 49193 samples within 2 assays 
Active assay: RNA (36601 features, 0 variable features)
 2 layers present: counts, data
 1 other assay present: ADT

3. Data PREPARATION

alldata <- All_samples_Merged

alldata.list <- SplitObject(alldata, split.by = "orig.ident")

for (i in 1:length(alldata.list)) {
    alldata.list[[i]] <- NormalizeData(alldata.list[[i]], verbose = FALSE)
    alldata.list[[i]] <- FindVariableFeatures(alldata.list[[i]], selection.method = "vst", nfeatures = 2000,verbose = FALSE)
}

# get the variable genes from all the datasets.
hvgs_per_dataset <- lapply(alldata.list, function(x) { x@assays$RNA@var.features })

# also add in the variable genes that was selected on the whole dataset
hvgs_per_dataset$all = VariableFeatures(alldata)

temp <- unique(unlist(hvgs_per_dataset))
overlap <- sapply( hvgs_per_dataset , function(x) { temp %in% x } )
pheatmap::pheatmap(t(overlap*1),cluster_rows = F ,
                   color = c("grey90","grey20"))

hvgs_all = SelectIntegrationFeatures(alldata.list)
hvgs_per_dataset$all_ranks = hvgs_all

temp <- unique(unlist(hvgs_per_dataset))
overlap <- sapply( hvgs_per_dataset , function(x) { temp %in% x } )
pheatmap::pheatmap(t(overlap*1),cluster_rows = F ,
                   color = c("grey90","grey20"))

alldata.list <- lapply(X = alldata.list, FUN = function(x) {
    x <- ScaleData(x, features = hvgs_all, verbose = FALSE)
    x <- RunPCA(x, features = hvgs_all, verbose = FALSE)
})

4. rpca-integration

alldata.anchors <- FindIntegrationAnchors(object.list = alldata.list, dims = 1:12, reduction = "rpca", anchor.features = hvgs_all)

Scaling features for provided objects

  |                                                  | 0 % ~calculating  
  |+++++++                                           | 12% ~05s          
  |+++++++++++++                                     | 25% ~04s          
  |+++++++++++++++++++                               | 38% ~04s          
  |+++++++++++++++++++++++++                         | 50% ~03s          
  |++++++++++++++++++++++++++++++++                  | 62% ~02s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06s  

Computing within dataset neighborhoods

  |                                                  | 0 % ~calculating  
  |+++++++                                           | 12% ~17s          
  |+++++++++++++                                     | 25% ~13s          
  |+++++++++++++++++++                               | 38% ~11s          
  |+++++++++++++++++++++++++                         | 50% ~08s          
  |++++++++++++++++++++++++++++++++                  | 62% ~06s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~04s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=17s  

Finding all pairwise anchors

  |                                                  | 0 % ~calculating  

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 798 anchors

  |++                                                | 4 % ~03m 29s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1024 anchors

  |++++                                              | 7 % ~03m 43s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1251 anchors

  |++++++                                            | 11% ~03m 26s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1784 anchors

  |++++++++                                          | 14% ~03m 17s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 628 anchors

  |+++++++++                                         | 18% ~03m 05s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 878 anchors

  |+++++++++++                                       | 21% ~02m 55s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1614 anchors

  |+++++++++++++                                     | 25% ~02m 46s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 554 anchors

  |+++++++++++++++                                   | 29% ~02m 37s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 920 anchors

  |+++++++++++++++++                                 | 32% ~02m 29s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 2612 anchors

  |++++++++++++++++++                                | 36% ~02m 21s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 2016 anchors

  |++++++++++++++++++++                              | 39% ~02m 12s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 757 anchors

  |++++++++++++++++++++++                            | 43% ~02m 03s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 924 anchors

  |++++++++++++++++++++++++                          | 46% ~01m 55s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1104 anchors

  |+++++++++++++++++++++++++                         | 50% ~01m 48s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1352 anchors

  |+++++++++++++++++++++++++++                       | 54% ~01m 40s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1053 anchors

  |+++++++++++++++++++++++++++++                     | 57% ~01m 32s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 644 anchors

  |+++++++++++++++++++++++++++++++                   | 61% ~01m 24s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 889 anchors

  |+++++++++++++++++++++++++++++++++                 | 64% ~01m 16s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 2061 anchors

  |++++++++++++++++++++++++++++++++++                | 68% ~01m 08s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 2628 anchors

  |++++++++++++++++++++++++++++++++++++              | 71% ~01m 00s      

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 1536 anchors

  |++++++++++++++++++++++++++++++++++++++            | 75% ~52s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 283 anchors

  |++++++++++++++++++++++++++++++++++++++++          | 79% ~45s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 294 anchors

  |++++++++++++++++++++++++++++++++++++++++++        | 82% ~38s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 281 anchors

  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~30s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 262 anchors

  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~23s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 255 anchors

  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~15s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 378 anchors

  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~08s          

Projecting new data onto SVD
Projecting new data onto SVD
Finding neighborhoods
Finding anchors
    Found 269 anchors

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 35s

alldata.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:12, new.assay.name = "rpca")

Merging dataset 7 into 5
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
Merging dataset 6 into 1
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
Merging dataset 4 into 5 7
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
Merging dataset 2 into 3
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
Merging dataset 1 6 into 5 7 4
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
Merging dataset 3 2 into 5 7 4 1 6
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
Merging dataset 8 into 5 7 4 1 6 3 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

names(alldata.int@assays)

[1] "RNA"  "ADT"  "rpca"

alldata.int@active.assay

[1] "rpca"

DefaultAssay(alldata.int) <- "rpca"

#Run Dimensionality reduction on integrated space
alldata.int <- ScaleData(alldata.int, verbose = FALSE)
alldata.int <- RunPCA(alldata.int, features = hvgs_all, reduction.name = "pca_rpca", do.print = TRUE, pcs.print = 1:5, genes.print = 15, verbose = FALSE)
alldata.int <- RunUMAP(alldata.int, reduction = "pca_rpca", reduction.name = "umap_rpca", dims = 1:12, verbose = FALSE)

Avis : The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

alldata.int <- RunTSNE(alldata.int, reduction = "pca_rpca", reduction.name = "tsne_rpca", dims = 1:12, verbose = FALSE)
alldata.int <- FindNeighbors(alldata.int, reduction = "pca_rpca", dims = 1:12, verbose = FALSE)
alldata.int <- FindClusters(alldata.int, resolution = 1.2, verbose = FALSE)



wrap_plots(

    DimPlot(alldata.int, reduction = "pca_rpca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_rpca", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_rpca", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),

    DimPlot(alldata.int, reduction = "pca_rpca", group.by = "rpca_snn_res.1.2")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_rpca", group.by = "rpca_snn_res.1.2")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_rpca", group.by = "rpca_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated"),
    ncol = 3) + plot_layout(guides = "collect")

NA
NA
NA

clean memory

# remove all objects that will not be used
rm(alldata, alldata.anchors)

gc()

             used    (Mb) gc trigger    (Mb)   max used    (Mb)
Ncells   13372812   714.2   22952785  1225.9   22952785  1225.9
Vcells 1958246802 14940.3 5018916606 38291.3 5018904699 38291.3

5. CCA-integration

clean memory

# remove all objects that will not be used
rm(alldata.anchors)

gc()

FeaturePlot

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")
FeaturePlot(alldata.int, reduction = "umap_CCA", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()

NA
NA

6. Harmony-integration

FeaturePlot

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")
FeaturePlot(alldata.int, reduction = "umap_harmony", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()

NA
NA

7. Save the Seurat object as an Robj file

#save(alldata.int, file = "Integration-by-different-Methods-1.Robj")