Harmony integrations of PBMC10x by cell_line-theta-0.5 both

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control after filtration
load("../22-Seurat_Integrate/0-R_Objects/CD4Tcells_SCTnormalized_done_on_HPC_inluding_Patient_origin.robj")



# Visualize before Harmony integration
DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "Patient_origin",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Cell Line")

before <- 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 = "cell_line",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By Cell Line")

DimPlot(All_samples_Merged, 
              reduction = "umap", 
              group.by = "SCT_snn_res.0.5",
              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")

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.5)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11   12   13
  B intermediate       0    3    0    0    0    0    2    0    0    0    0    0    2    0
  B memory             8    6    1    0   85    0   30    2    0  115    4    0    1    0
  CD14 Mono            0    1    0    0    0    0    4    0    0    7    0    0    0    0
  CD4 CTL              0    0    0    0    0   12    0    0    0    0    0    0    0    1
  CD4 Naive            0    8    0    0    0  517    0    0 1479    0    0   37    0    1
  CD4 Proliferating 5448 2474 5388 2852 3954    0 3256 2863    6 1270 1407    0   93    0
  CD4 TCM            871 3414  522  269  536 4214  106   29 1838  457   46  425   49   54
  CD4 TEM              0    1    0    0    0   61    0    0   21    0    0    1    0    0
  CD8 Proliferating    0    0    0    0    1    0    0    0    0    1    0    0    0    0
  CD8 TCM              0    1    0   16    0    0    0    0    0    0    0    0    0    0
  CD8 TEM              0    1    0    8    3    0    2    0    0    1    0    0    0    0
  cDC1                 0    0    0    0    5    0    2    0    0    0    0    0    1    0
  cDC2                 0    1    2    0    3    0   10    0    0   36    0    0    0    1
  dnT                  0    3    1    1    1    0    2    0    0    3    0    1    3    0
  HSPC                57   10    1    0  211    0  678  483    0    5  358    0    2    0
  NK Proliferating     4   40   23 2785  237    0   10   12    0   22    1    0   27    0
  Treg                15   14    1    0    1    0    0    0    0    0    0    1   13    0

3. Perform Harmony Integration

# Perform Harmony integration
All_samples_Merged <- RunHarmony(All_samples_Merged, 
                                 group.by.vars = c("cell_line"), 
                                 reduction.use = "pca", 
                                 dim.use = 1:15,
                                 theta = c(0.5),
                                 assay.use = "SCT")

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

# Check if Harmony integration ran successfully
print(names(All_samples_Merged@reductions))  # Should include "harmony"

[1] "integrated_dr" "ref.umap"      "pca"           "umap"          "harmony"      

# Find neighbors using the Harmony reduction and explicitly name the graph
All_samples_Merged <- FindNeighbors(All_samples_Merged,
                                    reduction = "harmony",   # Harmony reduction used
                                    dims = 1:15,             # Use first 15 dimensions of the Harmony reduction
                                    graph.name = "harmony_snn")  # Explicitly name the graph

Computing nearest neighbor graph
Computing SNN
Only one graph name supplied, storing nearest-neighbor graph only

# Check if the "harmony_snn" graph is present
print(names(All_samples_Merged@graphs))  # Should now include "harmony_snn"

[1] "SCT_nn"      "SCT_snn"     "harmony_snn"

# Find clusters for each resolution and store them
resolutions <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2)
for (res in resolutions) {
  cluster_name <- paste0("harmony_res_", res)  # Dynamic cluster name
  All_samples_Merged <- FindClusters(
    object = All_samples_Merged,
    graph.name = "harmony_snn",               # Graph created in FindNeighbors
    resolution = res,                         # Resolution for clustering
    verbose = FALSE
  )
  # Add cluster identities to metadata
  All_samples_Merged[[cluster_name]] <- Idents(All_samples_Merged)
}

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

Warning: 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 session18:06:20 UMAP embedding parameters a = 0.9922 b = 1.112
18:06:20 Read 49372 rows and found 15 numeric columns
18:06:20 Using Annoy for neighbor search, n_neighbors = 30
18:06:20 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
18:06:23 Writing NN index file to temp file /tmp/RtmpW4aDAH/file1480677f3ac2f
18:06:23 Searching Annoy index using 1 thread, search_k = 3000
18:06:35 Annoy recall = 100%
18:06:36 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
18:06:37 Initializing from normalized Laplacian + noise (using RSpectra)
18:06:39 Commencing optimization for 200 epochs, with 2025940 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
18:06:55 Optimization finished

4. Visualize Harmony Integrated Data

# Visualization after Harmony

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

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

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

# Print comparison plots
p3 + p4

print(p5)

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

before|after

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

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

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

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$cell_line)

                   
                      L1   L2   L3   L4   L5   L6   L7 PBMC PBMC_10x
  B intermediate       0    0    2    1    2    2    0    0        0
  B memory             0    0   11    1   38   82  120    0        0
  CD14 Mono            0    0    1    0    5    0    6    0        0
  CD4 CTL              0    0    0    0    0    0    0   12        1
  CD4 Naive            0    0    0    7    0    0    0  523     1512
  CD4 Proliferating 2461 2852 5452 5391 4732 4002 4115    0        6
  CD4 TCM           3320  270  887  562  178  557  517 4576     1963
  CD4 TEM              1    0    0    0    0    0    0   60       23
  CD8 Proliferating    0    0    0    0    0    1    1    0        0
  CD8 TCM              1   16    0    0    0    0    0    0        0
  CD8 TEM              1    8    0    0    2    3    1    0        0
  cDC1                 0    0    0    0    2    6    0    0        0
  cDC2                 0    0    0    4   11    3   35    0        0
  dnT                  2    3    0    1    2    5    2    0        0
  HSPC                 0    0   60    7 1035  213  490    0        0
  NK Proliferating    38 2785    6   24   11  259   38    0        0
  Treg                 1    1    9    9    4   15    6    0        0

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$harmony_res_0.7)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11   12   13
  B intermediate       0    0    0    0    0    0    0    2    1    1    0    3    0    0
  B memory             3    0    6    0    0    6    0   19    8  151    0   59    0    0
  CD14 Mono            0    0    0    0    0    0    0    0    0   11    0    0    0    1
  CD4 CTL              0    8    0    2    0    0    0    1    1    0    0    0    0    1
  CD4 Naive            0 1317    0    2    0    0    0    8   32    0  660    0   22    1
  CD4 Proliferating 5149   22 5156 2527 5181 2848 3341  880 1875 1020    0  886  126    0
  CD4 TCM            134 4870  944  109   63   40   17 2455 1384 1458  917  207  178   54
  CD4 TEM              0   79    0    1    0    0    0    0    4    0    0    0    0    0
  CD8 Proliferating    0    0    0    0    0    0    0    0    0    1    0    1    0    0
  CD8 TCM              0    0    0    4    0    0    0   13    0    0    0    0    0    0
  CD8 TEM              0    0    0    1    0    0    0    4    1    6    0    0    3    0
  cDC1                 0    0    0    0    0    0    1    0    0    2    0    4    1    0
  cDC2                 2    0    0    0    0    1    0    0    0   44    0    4    0    2
  dnT                  0    0    0    0    0    0    0    0    2    7    0    0    6    0
  HSPC               986    0   32    0   24  660   10   17    7   28    0   38    3    0
  NK Proliferating   115    4    0 2748   25    7  144    6   66    4    0   11   31    0
  Treg                 0    0    0    0    1    0    0    2   17    1    0    0   24    0

Visualize Harmony Integrated Data distribution

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$cell_line)

                   
                      L1   L2   L3   L4   L5   L6   L7 PBMC PBMC_10x
  B intermediate       0    0    2    1    2    2    0    0        0
  B memory             0    0   11    1   38   82  120    0        0
  CD14 Mono            0    0    1    0    5    0    6    0        0
  CD4 CTL              0    0    0    0    0    0    0   12        1
  CD4 Naive            0    0    0    7    0    0    0  523     1512
  CD4 Proliferating 2461 2852 5452 5391 4732 4002 4115    0        6
  CD4 TCM           3320  270  887  562  178  557  517 4576     1963
  CD4 TEM              1    0    0    0    0    0    0   60       23
  CD8 Proliferating    0    0    0    0    0    1    1    0        0
  CD8 TCM              1   16    0    0    0    0    0    0        0
  CD8 TEM              1    8    0    0    2    3    1    0        0
  cDC1                 0    0    0    0    2    6    0    0        0
  cDC2                 0    0    0    4   11    3   35    0        0
  dnT                  2    3    0    1    2    5    2    0        0
  HSPC                 0    0   60    7 1035  213  490    0        0
  NK Proliferating    38 2785    6   24   11  259   38    0        0
  Treg                 1    1    9    9    4   15    6    0        0

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$harmony_res_0.7)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11   12   13
  B intermediate       0    0    0    0    0    0    0    2    1    1    0    3    0    0
  B memory             3    0    6    0    0    6    0   19    8  151    0   59    0    0
  CD14 Mono            0    0    0    0    0    0    0    0    0   11    0    0    0    1
  CD4 CTL              0    8    0    2    0    0    0    1    1    0    0    0    0    1
  CD4 Naive            0 1317    0    2    0    0    0    8   32    0  660    0   22    1
  CD4 Proliferating 5149   22 5156 2527 5181 2848 3341  880 1875 1020    0  886  126    0
  CD4 TCM            134 4870  944  109   63   40   17 2455 1384 1458  917  207  178   54
  CD4 TEM              0   79    0    1    0    0    0    0    4    0    0    0    0    0
  CD8 Proliferating    0    0    0    0    0    0    0    0    0    1    0    1    0    0
  CD8 TCM              0    0    0    4    0    0    0   13    0    0    0    0    0    0
  CD8 TEM              0    0    0    1    0    0    0    4    1    6    0    0    3    0
  cDC1                 0    0    0    0    0    0    1    0    0    2    0    4    1    0
  cDC2                 2    0    0    0    0    1    0    0    0   44    0    4    0    2
  dnT                  0    0    0    0    0    0    0    0    2    7    0    0    6    0
  HSPC               986    0   32    0   24  660   10   17    7   28    0   38    3    0
  NK Proliferating   115    4    0 2748   25    7  144    6   66    4    0   11   31    0
  Treg                 0    0    0    0    1    0    0    2   17    1    0    0   24    0

table(All_samples_Merged$cell_line, All_samples_Merged$harmony_res_0.7)

          
              0    1    2    3    4    5    6    7    8    9   10   11   12   13
  L1          0   14   84   12   29    1  342 2333 2614  372    0    2   22    0
  L2          0    5    0 5362    0    0   25  476   55    2    0    0   10    0
  L3          0    0 3695    0 2379    0    5  181   20   15    0   77   55    1
  L4          0    0 2348    0 2885    0   12  100   12  385    0  251   12    2
  L5       1758    0    7    0    0 3224  340  104   35  240    0  304   10    0
  L6       1981    1    1    0    0  109 1626   73   61  714    0  459  123    0
  L7       2649    0    3    0    1  228 1163  100   70  992    0  105   20    0
  PBMC        0 3928    0   20    0    0    0   34  458   12  619   15   51   34
  PBMC_10x    1 2352    0    0    0    0    0    6   73    2  958    0   91   22

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

cd4_feature_plot1 <- FeaturePlot(
  All_samples_Merged, 
  features = myfeatures1, 
  reduction = "umap", 
  ncol = 4
) + 
  ggtitle("CD4 T Cell Marker Expression - Harmony Integration") +
  NoLegend()

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

# Display the plot
print(cd4_feature_plot1)

# 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_plot2 <- FeaturePlot(
  All_samples_Merged, 
  features = cd4_markers, 
  reduction = "umap", 
  ncol = 4
) + 
  ggtitle("CD4 T Cell Marker Expression - Harmony Integration") +
  NoLegend()

# Display the plot
print(cd4_feature_plot2)

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", 
            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:20], # First 8 markers
        stack = TRUE,
        flip = TRUE) +
        ggtitle("CD4 T Cell Marker Distribution Across Clusters")

NA
NA

6. Save the Seurat object as an Robj file

#save(All_samples_Merged, file = "../0-R_Objects/CD4Tcells_harmony_integrated_0.5_theta_patientorigin_cell_line.Robj")