Integration by Different Methods-rpca

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
load("../5-SS_ScRNA_Data_Analysis/4-ScSS_MyAnalysis_on_SS/2-Marie-Objects-Merged-into-Seurat-Object/1-Merged_using_marie_obj/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% ~09s          
  |++++++++++++++++++++++++++++++++                  | 62% ~07s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~04s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=18s  

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 23s      

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 54s      

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

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

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

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

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

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

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

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

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

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

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

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

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

  |++++++++++++++++++++++++++++++++++++              | 71% ~59s          

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

  |++++++++++++++++++++++++++++++++++++++            | 75% ~51s          

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

  |++++++++++++++++++++++++++++++++++++++++          | 79% ~44s          

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

  |++++++++++++++++++++++++++++++++++++++++++        | 82% ~37s          

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% ~22s          

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 33s

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   13482045   720.1   19769137  1055.8   19769137  1055.8
Vcells 1958609699 14943.1 5019866774 38298.6 5006253492 38194.7

FeaturePlot

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

Avis : Found the following features in more than one assay, excluding the default. We will not include these in the final data frame: CD4Avis : Could not find CD3E in the default search locations, found in 'RNA' assay insteadAvis : Could not find CD14 in the default search locations, found in 'RNA' assay insteadAvis : Could not find FCGR3A in the default search locations, found in 'RNA' assay insteadAvis : Could not find FCER1A in the default search locations, found in 'RNA' assay insteadAvis : The following requested variables were not found: CD4

6. Save the Seurat object as an Robj file