Harmony integrations of PBMC10x-part4

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

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

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 session09:31:02 UMAP embedding parameters a = 0.9922 b = 1.112
09:31:02 Read 59355 rows and found 22 numeric columns
09:31:02 Using Annoy for neighbor search, n_neighbors = 30
09:31:02 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:31:10 Writing NN index file to temp file /tmp/Rtmpcrhtj7/file3c27e46c0ca39
09:31:10 Searching Annoy index using 1 thread, search_k = 3000
09:31:36 Annoy recall = 100%
09:31:37 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
09:31:41 Initializing from normalized Laplacian + noise (using RSpectra)
09:31:50 Commencing optimization for 200 epochs, with 2508680 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:32:27 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: 1961526

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9505
Number of communities: 22
Elapsed time: 20 seconds

Harmony Visualization-1

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

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

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "predicted.celltype.l2",label = T, label.box = T) + 
  ggtitle("Harmony Integration - Annotations")

4. Data PREPARATION and Harmony Integration-2

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

09:43:42 UMAP embedding parameters a = 0.9922 b = 1.112
09:43:42 Read 59355 rows and found 22 numeric columns
09:43:42 Using Annoy for neighbor search, n_neighbors = 30
09:43:42 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:43:50 Writing NN index file to temp file /tmp/Rtmpcrhtj7/file3c27e72d964be
09:43:50 Searching Annoy index using 1 thread, search_k = 3000
09:44:15 Annoy recall = 100%
09:44:16 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
09:44:20 Initializing from normalized Laplacian + noise (using RSpectra)
09:44:27 Commencing optimization for 200 epochs, with 2541136 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:45:05 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: 25 seconds

Harmony Visualization-2

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

Registered S3 method overwritten by 'rmarkdown':
  method         from
  print.paged_df     

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

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

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "predicted.celltype.l2",label = T, label.box = T) + 
  ggtitle("Harmony Integration - Annotations")

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

09:53:56 UMAP embedding parameters a = 0.9922 b = 1.112
09:53:56 Read 59355 rows and found 22 numeric columns
09:53:56 Using Annoy for neighbor search, n_neighbors = 30
09:53:56 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:54:03 Writing NN index file to temp file /tmp/Rtmpcrhtj7/file3c27e24cfb76b
09:54:03 Searching Annoy index using 1 thread, search_k = 3000
09:54:29 Annoy recall = 100%
09:54:31 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
09:54:35 Initializing from normalized Laplacian + noise (using RSpectra)
09:54:40 Commencing optimization for 200 epochs, with 2543152 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:55:19 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: 25 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")

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

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "predicted.celltype.l2",label = T, label.box = T) + 
  ggtitle("Harmony Integration - Annotations")

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 = "../../../0-IMP-OBJECTS/Harmony_integrated_All_samples_Merged_with_PBMC10x.Robj")