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]



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

17:41:37 UMAP embedding parameters a = 0.9922 b = 1.112
17:41:37 Read 49193 rows and found 30 numeric columns
17:41:37 Using Annoy for neighbor search, n_neighbors = 30
17:41:37 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:41:41 Writing NN index file to temp file /tmp/Rtmph6UMoo/filec8ace5144214a
17:41:41 Searching Annoy index using 1 thread, search_k = 3000
17:41:52 Annoy recall = 100%
17:41:53 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
17:41:55 Initializing from normalized Laplacian + noise (using RSpectra)
17:41:57 Commencing optimization for 200 epochs, with 2110034 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:42:13 Optimization finished
Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49193
Number of edges: 1704400

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9163
Number of communities: 18
Elapsed time: 9 seconds

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

5. Save the Seurat object as an Robj file

 save(merged_seurat_filtered.harmony, file = "seurat.integrated_K_code_Harmony.Robj")