Multiple Harmony integrations of PBMC10x

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
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 PREPARATION and Harmony Integration-1

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "cell_line",  # Replace with the metadata column specifying batch or cell line
  dims.use = 1:22,  # Use the same dimensions as PCA
  assay.use = "SCT",
plot_convergence = TRUE
  )

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

# Run UMAP on the new Harmony reduction
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22, 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 session20:52:39 UMAP embedding parameters a = 0.9922 b = 1.112
20:52:39 Read 59355 rows and found 22 numeric columns
20:52:39 Using Annoy for neighbor search, n_neighbors = 30
20:52:39 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:52:47 Writing NN index file to temp file /tmp/RtmpXZmaNf/file3adae7a1148d
20:52:47 Searching Annoy index using 1 thread, search_k = 3000
20:53:16 Annoy recall = 100%
20:53:17 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:53:21 Initializing from normalized Laplacian + noise (using RSpectra)
20:53:24 Commencing optimization for 200 epochs, with 2569038 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:54:02 Optimization finished

# Find neighbors and clusters using the Harmony reduction
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)

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

Number of nodes: 59355
Number of edges: 1769146

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8990
Number of communities: 18
Elapsed time: 20 seconds

# Visualize results
p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line") + 
  ggtitle("Harmony Integration - By Cell Line")
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters") + 
  ggtitle("Harmony Integration - By Clusters")

# Compare with original UMAP
p3 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line") + 
  ggtitle("Original Integration - By Cell Line")
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters") + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p1 + p2)

print(p3 + p4)

NA
NA

Harmony Visualization-1

# Visualize results
p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")

# Compare with original UMAP
p3 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p1 + p2)

print(p3 + p4)

4. Data PREPARATION and Harmony Integration-2

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "cell_line",  # Replace with the metadata column specifying batch or cell line
  dims.use = 1:22,  # Use the same dimensions as PCA
  assay.use = "SCT",
  Theta = 0.5,
  lambda =0.5,
  plot_convergence = TRUE
  )

Transposing data matrix
Initializing state using k-means centroids initialization
Avis : Les étapes de transfer (quick-TRANSfer stage) ont dépassé le maximum (= 2967750)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

# Run UMAP on the new Harmony reduction
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22, reduction.name = "umap.harmony")

21:23:30 UMAP embedding parameters a = 0.9922 b = 1.112
21:23:30 Read 59355 rows and found 22 numeric columns
21:23:30 Using Annoy for neighbor search, n_neighbors = 30
21:23:30 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
21:23:37 Writing NN index file to temp file /tmp/RtmpXZmaNf/file3adae591aca7d
21:23:37 Searching Annoy index using 1 thread, search_k = 3000
21:24:06 Annoy recall = 100%
21:24:07 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
21:24:11 Initializing from normalized Laplacian + noise (using RSpectra)
21:24:14 Commencing optimization for 200 epochs, with 2575888 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
21:24:52 Optimization finished

# Find neighbors and clusters using the Harmony reduction
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)

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

Number of nodes: 59355
Number of edges: 1756087

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8918
Number of communities: 17
Elapsed time: 26 seconds

Harmony Visualization-2

# Visualize results
p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")

# Compare with original UMAP
p3 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p1 + p2)

print(p3 + p4)

5. Data PREPARATION and Harmony Integration-3

# Create a new metadata column for grouping
All_samples_Merged$sample_group <- case_when(
  grepl("^L", All_samples_Merged$cell_line) ~ "Cell_Line",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"
)

# Create the cell line grouping correctly
All_samples_Merged$cell_line_group <- case_when(
  grepl("^L[1-2]", All_samples_Merged$cell_line) ~ "P1",
  grepl("^L[3-4]", All_samples_Merged$cell_line) ~ "P2",
  grepl("^L[5-7]", All_samples_Merged$cell_line) ~ "P3",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"  # This catches any unexpected values
)

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "sample_group",
  dims.use = 1:22,  # Increased to capture more variation
  theta = c(0.5),  # Adjust these values as needed
  plot_convergence = TRUE
)

Transposing data matrix
Initializing state using k-means centroids initialization
Avis : Les étapes de transfer (quick-TRANSfer stage) ont dépassé le maximum (= 2967750)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 converged after 2 iterations

# Run UMAP on the new Harmony reduction
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22, reduction.name = "umap.harmony")

22:10:52 UMAP embedding parameters a = 0.9922 b = 1.112
22:10:52 Read 59355 rows and found 22 numeric columns
22:10:52 Using Annoy for neighbor search, n_neighbors = 30
22:10:52 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:11:00 Writing NN index file to temp file /tmp/RtmpXZmaNf/file3adae5b8271ac
22:11:22 Searching Annoy index using 1 thread, search_k = 3000
22:11:47 Annoy recall = 100%
22:11:49 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
22:12:36 Initializing from normalized Laplacian + noise (using RSpectra)
22:12:46 Commencing optimization for 200 epochs, with 2535858 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:13:22 Optimization finished

# Find neighbors and clusters using the Harmony reduction
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)

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

Number of nodes: 59355
Number of edges: 1956184

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9485
Number of communities: 21
Elapsed time: 46 seconds

Harmony Visualization-3

p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group", label = T, label.box = T)
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)
p3 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

p1 + p2 + p3

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

# Compare with original UMAP
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")
p5 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p4 + p5)

DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")

DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Visualize results
p6 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")
p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")
# Print the plots
print(p6 + p7)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")

p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")

6. Data PREPARATION and Harmony Integration-4

# Create a new metadata column for grouping
All_samples_Merged$sample_group <- case_when(
  grepl("^L", All_samples_Merged$cell_line) ~ "Cell_Line",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"
)

# Create the cell line grouping correctly
All_samples_Merged$cell_line_group <- case_when(
  grepl("^L[1-2]", All_samples_Merged$cell_line) ~ "P1",
  grepl("^L[3-4]", All_samples_Merged$cell_line) ~ "P2",
  grepl("^L[5-7]", All_samples_Merged$cell_line) ~ "P3",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"  # This catches any unexpected values
)

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "cell_line_group",
  dims.use = 1:22,  # Increased to capture more variation
  theta = c(0.5),  # Adjust these values as needed
  plot_convergence = TRUE
)

Transposing data matrix
Initializing state using k-means centroids initialization
Avis : Les étapes de transfer (quick-TRANSfer stage) ont dépassé le maximum (= 2967750)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 converged after 2 iterations

# Run UMAP on the new Harmony reduction
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22, reduction.name = "umap.harmony")

22:21:18 UMAP embedding parameters a = 0.9922 b = 1.112
22:21:18 Read 59355 rows and found 22 numeric columns
22:21:18 Using Annoy for neighbor search, n_neighbors = 30
22:21:18 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:21:26 Writing NN index file to temp file /tmp/RtmpXZmaNf/file3adae14396554
22:21:26 Searching Annoy index using 1 thread, search_k = 3000
22:21:51 Annoy recall = 100%
22:21:52 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
22:21:57 Initializing from normalized Laplacian + noise (using RSpectra)
22:22:03 Commencing optimization for 200 epochs, with 2541136 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:22:41 Optimization finished

# Find neighbors and clusters using the Harmony reduction
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)

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

Number of nodes: 59355
Number of edges: 1951462

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9474
Number of communities: 21
Elapsed time: 24 seconds

Harmony Visualization-4

p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group", label = T, label.box = T)
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)
p3 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

p1 + p2 + p3

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

# Compare with original UMAP
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")
p5 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p4 + p5)

DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")

DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Visualize results
p6 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")
p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")
# Print the plots
print(p6 + p7)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")

p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")

7. Data PREPARATION and Harmony Integration-5

# Create a new metadata column for grouping
All_samples_Merged$sample_group <- case_when(
  grepl("^L", All_samples_Merged$cell_line) ~ "Cell_Line",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"
)

# Create the cell line grouping correctly
All_samples_Merged$cell_line_group <- case_when(
  grepl("^L[1-2]", All_samples_Merged$cell_line) ~ "P1",
  grepl("^L[3-4]", All_samples_Merged$cell_line) ~ "P2",
  grepl("^L[5-7]", All_samples_Merged$cell_line) ~ "P3",
  All_samples_Merged$cell_line == "PBMC" ~ "PBMC",
  All_samples_Merged$cell_line == "PBMC_10x" ~ "PBMC_10x",
  TRUE ~ "Other"  # This catches any unexpected values
)

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "cell_line",
  dims.use = 1:22,  # Increased to capture more variation
  theta = c(0.5),  # Adjust these values as needed
  plot_convergence = TRUE
)

Transposing data matrix
Initializing state using k-means centroids initialization
Avis : Les étapes de transfer (quick-TRANSfer stage) ont dépassé le maximum (= 2967750)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 converged after 2 iterations

# Run UMAP on the new Harmony reduction
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22, reduction.name = "umap.harmony")

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

# Find neighbors and clusters using the Harmony reduction
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)

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

Number of nodes: 59355
Number of edges: 1882456

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9226
Number of communities: 18
Elapsed time: 24 seconds

Harmony Visualization-5

p1 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group", label = T, label.box = T)
p2 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)
p3 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

p1 + p2 + p3

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "sample_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line_group",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

# Compare with original UMAP
p4 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")
p5 <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Print the plots
print(p4 + p5)

DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line",label = T, label.box = T) + 
  ggtitle("Original Integration - By Cell Line")

DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Original Integration - By Clusters")

# Visualize results
p6 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")
p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")
# Print the plots
print(p6 + p7)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line", label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Cell Line")

p7 <- DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "seurat_clusters",label = T, label.box = T) + 
  ggtitle("Harmony Integration - By Clusters")

Marker Gene Visualization

# Set marker genes specific to requested immune cell types
myfeatures <- c("CD19", "CD79A", "MS4A1", # B cells
                "CD14", "LYZ", "FCGR3A", # Monocytes
                "CSF1R", "CD68", # Macrophages
                "NKG7", "GNLY", "KIR3DL1", # NK cells
                "MKI67", # Proliferating NK cells
                "CD34", "KIT", # HSPCs
                "CD3E", "CCR7", # T cells
                "SELL", "CD45RO", # Tnaive, Tcm
                "CD44", "CD45RA") # Tem, Temra
# Visualize marker genes for Harmony
FeaturePlot(All_samples_Merged, features = myfeatures, reduction = "umap.harmony", ncol = 4) + 
  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

#save(All_samples_Merged, file = "integrated_All_samples_Merged_with_PBMC10x.Robj")