Integration by Different Methods-Harmony

1. load libraries

2. Load Seurat Object

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


#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
load("../5-SS_ScRNA_Data_Analysis/4-ScSS_MyAnalysis_on_SS/5-HPC-Run_on_Robj/Integration-by-different-Methods-2.Robj")

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

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

alldata.int <- RunHarmony(
  alldata.int,
  group.by.vars = "orig.ident",
  reduction.use = "pca_harmony",
  dims.use = 1:12,
  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 converged after 4 iterations

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

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 session17:14:56 UMAP embedding parameters a = 0.9922 b = 1.112
17:14:56 Read 49193 rows and found 12 numeric columns
17:14:56 Using Annoy for neighbor search, n_neighbors = 30
17:14:56 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:15:02 Writing NN index file to temp file /tmp/RtmpmTbV8R/file1494bf319569
17:15:02 Searching Annoy index using 1 thread, search_k = 3000
17:15:26 Annoy recall = 100%
17:15:27 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
17:15:31 Initializing from normalized Laplacian + noise (using RSpectra)
17:15:34 Commencing optimization for 200 epochs, with 2017658 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:16:05 Optimization finished

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

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

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


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

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

5. Save the Seurat object as an Robj file

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