Harmony integrations of PBMC10x by cell_line-theta-0.5_removing non CD4Tcells and B cells from L4

1. load libraries

2. Load Seurat Object

# Check original UMAP before integration
p1 <- DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "cell_line",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Cell Line")

p2 <- DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Clusters")

# Print original plots
p1 + p2

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


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

DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "predicted.celltype.l1",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Annotation.l1")

DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "predicted.celltype.l2",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Annotation.l2")

DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "predicted.celltype.l3",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Annotation.l3")

Harmony Visualization-1

library(harmony)

All_samples_Merged <- RunHarmony(
  object = All_samples_Merged,
  group.by.vars = "cell_line",
  dims.use = 1:16,
  theta = 0.5,  # Lower theta to maintain biological differences
  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:16, reduction.name = "umap.harmony")

20:49:52 UMAP embedding parameters a = 0.9922 b = 1.112
20:49:52 Read 49388 rows and found 16 numeric columns
20:49:52 Using Annoy for neighbor search, n_neighbors = 30
20:49:52 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:49:58 Writing NN index file to temp file /tmp/Rtmpp3xsxj/file161ae333ffbabc
20:49:58 Searching Annoy index using 1 thread, search_k = 3000
20:50:19 Annoy recall = 100%
20:50:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:50:24 Initializing from normalized Laplacian + noise (using RSpectra)
20:50:27 Commencing optimization for 200 epochs, with 2089626 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:50:56 Optimization finished

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

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: 49388
Number of edges: 1563357

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9144
Number of communities: 15
Elapsed time: 15 seconds

# 7. Visualization after Harmony
# By cell line
p3 <- DimPlot(All_samples_Merged, 
              reduction = "umap.harmony", 
              group.by = "cell_line",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("After Harmony - By Cell Line")

# By clusters
p4 <- DimPlot(All_samples_Merged, 
              reduction = "umap.harmony", 
              group.by = "seurat_clusters",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("After Harmony - By Clusters")

# By cell type annotations
p5 <- DimPlot(All_samples_Merged, 
              reduction = "umap.harmony", 
              group.by = "predicted.celltype.l2",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("After Harmony - Cell Type Annotations")

# Print comparison plots
p3 + p4

print(p5)

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



# Define markers specific to CD4 T cells and their subsets
cd4_markers <- c(
  "CD4",          # General CD4 T cells
  "IL7R",         # Naive T cells
  "CCR7",         # T central memory (Tcm) cells
  "SELL",         # T naive cells
  "FOXP3",        # Regulatory T cells (Tregs)
  "IL2RA",        # Activated T cells
  "PDCD1",        # Exhausted T cells
  "LAG3",         # Exhausted T cells
  "TIGIT",        # Exhausted T cells
  "GATA3",        # Th2 cells
  "TBX21",        # Th1 cells
  "RORC",         # Th17 cells
  "BCL6"          # T follicular helper (Tfh) cells
)

# Visualize marker genes for CD4 T cells
cd4_feature_plot <- FeaturePlot(
  All_samples_Merged, 
  features = cd4_markers, 
  reduction = "umap.harmony", 
  ncol = 4
) + 
  ggtitle("CD4 T Cell Marker Expression - Harmony Integration") +
  NoLegend()

# Display the plot
print(cd4_feature_plot)

CD4 T Cell Marker Visualization

# Set marker genes specific to CD4 T cell biology and states
cd4_markers <- c(
    # Core T cell markers
    "CD3E",     # T cell marker
    "CD4",      # CD4 T cell marker
    
    # Naive/Memory markers
    "CCR7",     # Naive/Central memory
    "SELL",     # L-selectin, naive marker
    "CD27",     # Memory marker
    "IL7R",     # Naive/Memory marker
    
    # Activation/State markers
    "IL2RA",    # CD25, activation marker
    "CD69",     # Early activation
    "HLA-DRA",  # Activation marker
    
    # Exhaustion markers
    "PDCD1",    # PD-1
    "LAG3",     # Exhaustion marker
    "TIGIT",    # Exhaustion marker
    
    # Regulatory T cell markers
    "FOXP3",    # Treg marker
    "IL2RA",    # CD25, Treg marker
    "CTLA4",    # Treg/exhaustion marker
    
    # Effector/Function markers
    "IL2",      # T cell function
    "IFNG",     # Th1
    "IL4",      # Th2
    "IL13",     # Th2
    "IL17A"     # Th17
)

# Create feature plots with better visualization
FeaturePlot(All_samples_Merged, 
            features = cd4_markers, 
            reduction = "umap.harmony", 
            ncol = 4,
            pt.size = 0.1,           # Smaller point size for better resolution
            min.cutoff = "q1",       # Remove bottom 1% of expression
            max.cutoff = "q99",      # Remove top 1% of expression
            order = TRUE) +          # Plot highest expressing cells on top
    ggtitle("CD4 T Cell Marker Expression - Harmony Integration") +
    theme(plot.title = element_text(size = 16, face = "bold")) +
    NoLegend()

# Optional: Add violin plots to see expression distribution across clusters
VlnPlot(All_samples_Merged, 
        features = cd4_markers[1:8], # First 8 markers
        stack = TRUE,
        flip = TRUE) +
        ggtitle("CD4 T Cell Marker Distribution Across Clusters")

NA
NA

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