Harmony integrations of PBMC10x by cell_line-theta-0

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
)

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

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

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

NA
NA
NA

Harmony Visualization-1

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = c("cell_line"),
  theta = c(0),
  dims.use = 1:22,  # Increased to capture more variation
  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 session11:32:46 UMAP embedding parameters a = 0.9922 b = 1.112
11:32:46 Read 59355 rows and found 22 numeric columns
11:32:46 Using Annoy for neighbor search, n_neighbors = 30
11:32:46 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:32:54 Writing NN index file to temp file /tmp/RtmpyqnRv8/file46c3e342a84ee
11:32:54 Searching Annoy index using 1 thread, search_k = 3000
11:33:19 Annoy recall = 100%
11:33:21 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
11:33:25 Initializing from normalized Laplacian + noise (using RSpectra)
11:33:30 Commencing optimization for 200 epochs, with 2528458 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:34:08 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: 1902411

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9303
Number of communities: 19
Elapsed time: 18 seconds

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

4. Save the Seurat object as an Robj file

#save(All_samples_Merged, file = "../../../0-IMP-OBJECTS/Harmony_integrated_All_samples_Merged_with_PBMC10x.Robj")