Harmony Integration of PBMC10x-part5

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
SS_All_samples_Merged <- load("../../../0-IMP-OBJECTS/All_Samples_Merged_with_10x_Azitmuth_Annotated_SCT_HPC_without_harmony_integration.robj")

All_samples_Merged

An object of class Seurat 
64169 features across 59355 samples within 6 assays 
Active assay: SCT (27417 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 4 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap

3. Data PREPERATION

ElbowPlot(All_samples_Merged)

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



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

NA
NA
NA

4. Harmony-integration

# run Harmony -----------
All_samples_Merged_Harmony_Integrated <- All_samples_Merged %>%
  RunHarmony(group.by.vars = 'orig.ident', plot_convergence = FALSE)

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 converged after 3 iterations

All_samples_Merged_Harmony_Integrated@reductions

$integrated_dr
A dimensional reduction object with key integrateddr_ 
 Number of dimensions: 50 
 Number of cells: 59355 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$ref.umap
A dimensional reduction object with key UMAP_ 
 Number of dimensions: 2 
 Number of cells: 59355 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$pca
A dimensional reduction object with key PC_ 
 Number of dimensions: 50 
 Number of cells: 59355 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

$umap
A dimensional reduction object with key umap_ 
 Number of dimensions: 2 
 Number of cells: 59355 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

$harmony
A dimensional reduction object with key harmony_ 
 Number of dimensions: 50 
 Number of cells: 59355 
 Projected dimensional reduction calculated:  TRUE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

All_samples_Merged_Harmony_Integrated.embed <- Embeddings(All_samples_Merged_Harmony_Integrated, "harmony")
All_samples_Merged_Harmony_Integrated.embed[1:10,1:10]

                        harmony_1 harmony_2 harmony_3 harmony_4 harmony_5   harmony_6 harmony_7  harmony_8  harmony_9 harmony_10
L1_AAACCTGAGGGCTTCC-1 -24.9425010  9.729485 -2.603145 -4.801961  9.290611  0.70200038  3.527128 -2.0336592 -2.9965266 -0.6076247
L1_AAACCTGGTGCAGGTA-1 -13.1935667  1.809197 -7.899042 -2.298168  5.933104  0.73286641  5.266604 -3.7656839 -1.5253091  3.3371281
L1_AAACCTGGTTAAAGTG-1   3.5891781 -2.936323 -7.792930  2.040950 -2.604144  0.94453792  4.399951  0.4034785 -1.8202440  0.9328230
L1_AAACCTGTCAGGTAAA-1   7.0089393  4.336626 -3.333358 -1.719958  1.846025  0.08728845  1.298744  1.5439122  0.1570706  1.4365544
L1_AAACCTGTCCCTGACT-1 -21.9519404  5.354331 -1.923206 -4.699053  9.830525  1.97438886 -1.175053 -2.0316367 -1.4417260 -1.0526408
L1_AAACCTGTCCTTCAAT-1  -1.6364858 -3.559325 -6.498843  1.742341 -0.893657  1.01139857  3.132217 -3.0447033 -3.1399042  1.1422886
L1_AAACCTGTCTTGCAAG-1  -0.8480336 -4.533099 -6.374012  1.146192 -1.427099  0.81612901  1.022530 -2.3285230 -2.2636863  0.3815743
L1_AAACGGGAGGCTAGAC-1   1.5774096  2.718044 -3.403277 -1.919959  3.054539  1.45739699  2.152590  1.9496453  0.7156881 -0.3798055
L1_AAACGGGAGGGTATCG-1   4.4011655 -2.366682 -2.105957 -1.753435 -1.048252  2.21873581  3.606342 -0.1197602  0.4755433 -1.2568843
L1_AAACGGGAGGGTTCCC-1  -6.2110019 10.223965  3.783698 -6.717456  6.398826 -0.62214007 -1.842301  3.3429751  1.5598661 -5.0025031

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

15:56:49 UMAP embedding parameters a = 0.9922 b = 1.112
15:56:49 Read 59355 rows and found 22 numeric columns
15:56:49 Using Annoy for neighbor search, n_neighbors = 30
15:56:49 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:56:58 Writing NN index file to temp file /tmp/Rtmp1rPu2p/file37d5229bcc0a3
15:56:58 Searching Annoy index using 1 thread, search_k = 3000
15:57:26 Annoy recall = 100%
15:57:28 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:57:32 Initializing from normalized Laplacian + noise (using RSpectra)
15:57:36 Commencing optimization for 200 epochs, with 2570420 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:58:15 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: 59355
Number of edges: 1775376

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9022
Number of communities: 16
Elapsed time: 29 seconds

# visualize 


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

before|after

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

DimPlot(All_samples_Merged_Harmony_Integrated,reduction = "umap", group.by = "SCT_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)

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

5. Marker Gene Visualization

# Visualize marker genes for Harmony
FeaturePlot(All_samples_Merged_Harmony_Integrated, features = myfeatures, reduction = "umap", ncol = 5) + 
  ggtitle("Marker Gene Expression - Harmony Integration") +
  NoLegend()

Avis : Could not find CD45RO in the default search locations, found in 'ADT' assay insteadAvis : Could not find CD45RA in the default search locations, found in 'ADT' assay instead

6. Save the Seurat object as an Robj file

```