Integration by CCA-Harmony

1. load libraries

2. Load Seurat Object

load("AllSample_corrected_AzimuthAnnotated_L1.Robj")

All_samples_Merged

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 = "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.0.5")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_CCA", group.by = "CCA_snn_res.0.5")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.0.5")+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.0.5")+NoAxes()+ggtitle("PCA integrated")

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

    DimPlot(alldata.int, reduction = "umap_CCA", group.by = "CCA_snn_res.0.5")+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 = "tsne_CCA", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)+NoAxes()+ggtitle("tSNE 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")

5. Harmony-integration

alldata.int@active.assay = "RNA"

VariableFeatures(alldata.int) = hvgs_all_final
alldata.int = ScaleData(alldata.int, verbose = TRUE)

Centering and scaling data matrix

  |                                                                                                            
  |                                                                                                      |   0%
  |                                                                                                            
  |===================================================                                                   |  50%
  |                                                                                                            
  |======================================================================================================| 100%

alldata.int = RunPCA(alldata.int, reduction.name = "pca_harmony")

PC_ 1 
Positive:  PNRC1, CD52, MALAT1, GIMAP5, TCF7, YPEL3, BTG1, BTG2, FCMR, LEF1 
       ZFP36, TXNIP, SELL, SH3BGRL3, ZFP36L2, FYB1, IL7R, CD44, PBXIP1, IKZF3 
       KLF2, LBH, GZMM, LAPTM5, LTB, SYNE2, CD27, NCF1, GAS5, CD7 
Negative:  TUBA1B, FABP5, H2AFZ, JPT1, MTHFD2, UBE2S, LGALS1, HMGB2, HSPD1, TXN 
       TYMS, PRDX1, PRELID1, PTTG1, GAPDH, IL2RA, CD70, NPM1, TUBB, CCND2 
       PKM, ENO1, PCLAF, HNRNPAB, NME1, CLIC1, ATP5MC3, BATF3, MTDH, RRM2 
PC_ 2 
Positive:  CST7, XCL1, MATK, ZBTB16, GZMM, XCL2, C1QBP, KIR2DL4, ESYT2, KRT86 
       TRBV23-1, KRT81, CXCR3, SRRT, TOX, NKG7, EIF4A1, PFN1, GYPC, TRBV20-1 
       NME2, IL32, LTB, CD7, CORO1B, IKZF2, MYO1E, LSR, MRPL12, IL12RB2 
Negative:  SQSTM1, ANXA1, NFKBIA, CD74, FAM107B, FTH1, KYNU, SAT1, LGALS3, ZEB2 
       EMP3, CFLAR, TYMP, SLC7A11, CXCL8, MARCKS, DOCK4, LMNA, TYROBP, C15orf48 
       FCER1G, CTSZ, NCF2, TIMP1, BCL3, BIRC3, CCR7, RBM47, CD2, AHNAK 
PC_ 3 
Positive:  RBPMS, TRAV17, TENM3, PPBP, PPP2R2B, LMNA, RPL22L1, PLD1, IQCG, NCR3 
       VAMP5, SRGN, SPOCK1, CTAG2, SLC7A11-AS1, STMN1, STAT1, ANXA2, KIF2A, FAM241A 
       TNFSF10, MACROD2, AC010967.1, MT2A, TK1, PLAAT3, TUBA4A, RYR2, CREB3L3, IGFBP3 
Negative:  TNFRSF4, SYT4, BACE2, CYBA, PHLDA2, CA2, CCL17, PTP4A3, ADGRB3, GRIA4 
       LY6E, ONECUT2, EEF1A2, IGHE, PON2, PRG4, MCTP2, KRT7, AC112770.1, THY1 
       MSC, ACTN1, SMIM3, FGGY, IGKC, PLPP1, CFI, MAP4K4, RHOC, NME4 
PC_ 4 
Positive:  TYROBP, DOCK4, FCER1G, C15orf48, THBS1, KYNU, CXCL8, MMP9, EPB41L3, MS4A7 
       SLC43A2, LYZ, TNFAIP2, CXCL5, SLC7A11, GLIS3, GAS7, CCR1, RBM47, LILRB4 
       RNF130, TIMP2, NCF2, LUCAT1, ZEB2, ZMIZ1, VCAN, PLEK, EMP1, MT1G 
Negative:  RPL7, CCR7, NPM1, CD74, EEF2, MIR155HG, GIMAP5, SELL, ALOX5AP, ARID5A 
       LTA, YBX3, SNHG8, CCL17, S1PR1, SLC35F3, TCF7, HSP90AB1, FCER2, TSHZ2 
       RPL22L1, RXFP1, BATF3, SOCS1, GAS5, DANCR, CCR4, AL590550.1, CA2, CD2 
PC_ 5 
Positive:  S100A4, DUSP4, ENTPD1, F2R, TP73, AL136456.1, RAB37, LINC02694, CD2, EEF1A2 
       S100A6, GPAT3, CRIP1, WNT5B, RPS6KA5, TMEM163, ITM2A, QPRT, TNFRSF18, TMSB4X 
       LIME1, HIST1H1C, PTCHD1-AS, LINC01934, PALLD, GAS2L1, IRX3, AHNAK, MYO1G, LAPTM5 
Negative:  FCER2, LTA, MIR155HG, DNAJC12, RPL23, C1QBP, PPBP, PPID, RXFP1, KYNU 
       SLC7A11-AS1, COL19A1, NFIB, MIR3945HG, TXNDC17, ZNF804A, CSMD1, AC068672.2, SLC35F3, CCL17 
       RBM47, ODC1, MPP1, UNC13C, EIF5A, DOCK4, ACSL1, ARL4A, IFNGR2, FCER1G 
Warning: Key 'PC_' taken, using 'pcaharmony_' instead

alldata.int <- RunHarmony(
  alldata.int,
  group.by.vars = "orig.ident",
  reduction.use = "pca_harmony",
  dims.use = 1:15,
  assay.use = "RNA")

Transposing data matrix
Initializing state using k-means centroids initialization
Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations

alldata.int <- RunUMAP(alldata.int, dims = 1:15, reduction = "harmony", reduction.name = "umap_harmony")

11:59:08 UMAP embedding parameters a = 0.9922 b = 1.112
11:59:08 Read 49193 rows and found 15 numeric columns
11:59:08 Using Annoy for neighbor search, n_neighbors = 30
11:59:08 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:59:11 Writing NN index file to temp file /tmp/RtmpDstfLq/fileac4b85dd5c26a
11:59:11 Searching Annoy index using 1 thread, search_k = 3000
11:59:23 Annoy recall = 100%
11:59:23 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
11:59:25 Initializing from normalized Laplacian + noise (using RSpectra)
11:59:27 Commencing optimization for 200 epochs, with 2060712 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:59:43 Optimization finished

alldata.int <- RunTSNE(alldata.int, dims = 1:15, reduction = "harmony", reduction.name = "tsne_harmony")

Warning: Key 'tSNE_' taken, using 'tsneharmony_' instead

alldata.int <- FindNeighbors(alldata.int, reduction = "pca_harmony", dims = 1:15, verbose = FALSE)
alldata.int <- FindClusters(alldata.int, resolution = 0.5, verbose = FALSE)


 
 
 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.0.5")+NoAxes()+ggtitle("UMAP integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.0.5", 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 = "RNA_snn_res.0.5")+NoAxes()+ggtitle("PCA integrated"),
    DimPlot(alldata.int, reduction = "tsne_harmony", group.by = "RNA_snn_res.0.5")+NoAxes()+ggtitle("tSNE integrated"),
    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.0.5")+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.0.5")+NoAxes()+ggtitle("UMAP integrated")

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.0.5", 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 = "RNA_snn_res.0.5")+NoAxes()+ggtitle("PCA integrated")

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

    DimPlot(alldata.int, reduction = "umap_harmony", group.by = "RNA_snn_res.0.5")+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")

NA
NA

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