Integration by CCA-Harmony

1. load libraries

2. Load Seurat Object

load("AllSample_corrected_AzimuthAnnotated_L1.Robj")

All_samples_Merged

An object of class Seurat 
36724 features across 49193 samples within 5 assays 
Active assay: RNA (36601 features, 0 variable features)
 2 layers present: counts, data
 4 other assays present: ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 2 dimensional reductions calculated: integrated_dr, ref.umap

3. Data PREPERATION

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. CCA-integration

wrap_plots(

    
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "orig.ident", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.1.2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    
    ncol = 2) + plot_layout(guides = "collect")


wrap_plots(

    DimPlot(alldata.int, reduction = "pca_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),
  
    DimPlot(alldata.int, reduction = "pca_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "pca_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated"),
    
    ncol = 3) + plot_layout(guides = "collect")


DimPlot(alldata.int, reduction = "pca_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")

  
    DimPlot(alldata.int, reduction = "pca_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated")

    
    DimPlot(alldata.int, reduction = "pca_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated")

5. Harmony-integration

wrap_plots(

    
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.1.2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated"),
     DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    ncol = 2) + plot_layout(guides = "collect")

wrap_plots(

    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),
  
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated"),
    
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated"),
    ncol = 3) + plot_layout(guides = "collect")


DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated")

    
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.1.2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated")

    
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated")

     DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated")

   

DimPlot(alldata.int, reduction = "pca_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")

  
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "CCA_snn_res.1.2")+NoAxes()+ggtitle("UMAP integrated")

    
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2")+NoAxes()+ggtitle("UMAP integrated")

    
    DimPlot(alldata.int, reduction = "pca_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("PCA integrated")

    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("tSNE integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("UMAP integrated")

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()

FeaturePlot(alldata.int, reduction = "umap_CCA", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()

6. Save the Seurat object as an Robj file

save(alldata.int, file = "Integrated_by_CCA_Harmony.Robj")