Integration by Harmony_by_K

1. load libraries

2. Load Seurat Object

 load("AllSample_corrected_AzimuthAnnotated_L1.Robj")

 merged_seurat_filtered <- All_samples_Merged

merged_seurat_filtered

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

# # perform standard workflow steps to figure out if we see any batch effects --------
# merged_seurat_filtered <- NormalizeData(object = merged_seurat_filtered, verbose = FALSE)
# merged_seurat_filtered <- FindVariableFeatures(object = merged_seurat_filtered , selection.method = "vst", nfeatures = 2000,verbose = FALSE)
# 
# merged_seurat_filtered <- ScaleData(object = merged_seurat_filtered)
# 
# 
# Variable_genes <- merged_seurat_filtered@assays$RNA@var.features
# 
# Variables_genes_after_exclusion <- Variable_genes[!grepl("^HLA-|^Xist", Variable_genes)]
# 
# merged_seurat_filtered <- RunPCA(merged_seurat_filtered,
#                          features = Variables_genes_after_exclusion,
#                          do.print = TRUE,
#                          pcs.print = 1:5,
#                          genes.print = 15)

ElbowPlot(merged_seurat_filtered)

# merged_seurat_filtered <- RunUMAP(object = merged_seurat_filtered, dims = 1:20)


# plot
before <- DimPlot(merged_seurat_filtered, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)



DimPlot(merged_seurat_filtered, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)

NA
NA
NA

4. Harmony-integration

# run Harmony -----------
# merged_seurat_filtered.harmony <- merged_seurat_filtered %>%
#   RunHarmony(group.by.vars = 'orig.ident', plot_convergence = FALSE)
# 
# merged_seurat_filtered.harmony@reductions
# 
# merged_seurat_filtered.harmony.embed <- Embeddings(merged_seurat_filtered.harmony, "harmony")
merged_seurat_filtered.harmony.embed[1:10,1:10]

                      harmony_1 harmony_2  harmony_3  harmony_4  harmony_5 harmony_6  harmony_7  harmony_8
L1_AAACCTGAGGGCTTCC-1 -6.427262  8.254969  0.7715261  1.3011146 -1.9436480 1.0627125  0.6935481 -0.3321792
L1_AAACCTGGTGCAGGTA-1 -6.873649  9.065062  1.2588632 -0.6417722 -3.6111313 0.1986480 -1.7125518 -1.3283037
L1_AAACCTGGTTAAAGTG-1 -5.077434  7.937975  1.9739861 -3.4979441  0.4035189 0.1413428  0.6675758  0.3393447
L1_AAACCTGTCAGGTAAA-1 -5.625381  8.585351 -0.2443243 -1.1042725  0.5691545 0.1484833  1.4564738  0.9102239
L1_AAACCTGTCCCTGACT-1 -5.430924  6.315944 -0.1937682  2.4432789  0.0980270 1.7624673  2.5231999  0.0290877
L1_AAACCTGTCCTTCAAT-1 -5.183905  7.760936  0.8473519 -1.8797225 -0.3723902 0.7151263  0.3032509  0.3750337
L1_AAACCTGTCTTGCAAG-1 -6.166992  8.118478  1.7496105 -1.8267948 -1.2020396 0.4086726 -0.1776375 -0.3375075
L1_AAACGGGAGGCTAGAC-1 -5.334131  7.449693  0.4832964 -0.2057412  0.6959762 0.9464938  2.0398931 -0.1466619
L1_AAACGGGAGGGTATCG-1 -3.965787  8.008437  1.7013528 -3.2969653 -0.1033832 0.5942966  0.1964398  0.1416826
L1_AAACGGGAGGGTTCCC-1 -3.599167  4.826824 -1.3444358  3.8005011  1.3107254 0.1699336  2.5707134  0.7357237
                        harmony_9   harmony_10
L1_AAACCTGAGGGCTTCC-1  0.30567983 -1.459491708
L1_AAACCTGGTGCAGGTA-1  0.16344334 -3.122505580
L1_AAACCTGGTTAAAGTG-1  0.62933951  1.222807346
L1_AAACCTGTCAGGTAAA-1  1.21593701 -0.005167528
L1_AAACCTGTCCCTGACT-1  0.01445315 -0.797898497
L1_AAACCTGTCCTTCAAT-1  0.46542327 -0.611959137
L1_AAACCTGTCTTGCAAG-1 -0.78707360  0.894576523
L1_AAACGGGAGGCTAGAC-1  0.41925378  0.026021925
L1_AAACGGGAGGGTATCG-1 -0.42878922  5.479253362
L1_AAACGGGAGGGTTCCC-1  0.99484956 -1.548533529

# Do UMAP and clustering using ** Harmony embeddings instead of PCA **
# merged_seurat_filtered.harmony <- merged_seurat_filtered.harmony %>%
#   RunUMAP(reduction = 'harmony', dims = 1:20) %>%
#   FindNeighbors(reduction = "harmony", dims = 1:20) %>%
#   FindClusters(resolution = 0.5)

# visualize 


after <- DimPlot(merged_seurat_filtered.harmony, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)

before|after

DimPlot(merged_seurat_filtered.harmony,reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)

DimPlot(merged_seurat_filtered.harmony, reduction = "umap", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)

FeaturePlot

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")

FeaturePlot(merged_seurat_filtered.harmony, reduction = "umap", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()

6. Save the Seurat object as an Robj file

# save(seurat.integrated, file = "seurat.integrated_K_code_Harmony.Robj")