Merging Poglio and Gaydosik Data and harmony integration

1. load libraries

1. Load your objects

# Initialize new column with NA
combined_merged$PatientOrigin <- NA

# Assign groups for your in-house data ("ourData")
combined_merged$PatientOrigin[combined_merged$batch %in% c("L1", "L2")] <- "P1_ourData"
combined_merged$PatientOrigin[combined_merged$batch %in% c("L3", "L4")] <- "P2_ourData"
combined_merged$PatientOrigin[combined_merged$batch %in% c("L5", "L6", "L7")] <- "P3_ourData"
combined_merged$PatientOrigin[combined_merged$batch == "PBMC"] <- "PBMC_ourData"

# For GaydosikData cell lines
combined_merged$PatientOrigin[combined_merged$batch == "SS1"] <- "SS1_GaydosikData"
combined_merged$PatientOrigin[combined_merged$batch == "SS2"] <- "SS2_GaydosikData"
combined_merged$PatientOrigin[combined_merged$batch == "SS3"] <- "SS3_GaydosikData"


# For GaydosikData healthy controls
combined_merged$PatientOrigin[combined_merged$batch == "HC1_2"] <- "HC1_2_GaydosikData"
combined_merged$PatientOrigin[combined_merged$batch == "HC3_4"] <- "HC3_4_GaydosikData"


# View table to check annotations
table(combined_merged$PatientOrigin)

HC1_2_GaydosikData HC3_4_GaydosikData         P1_ourData         P2_ourData         P3_ourData 
              2035               3161              11543              12389              16300 
      PBMC_ourData   SS1_GaydosikData   SS2_GaydosikData   SS3_GaydosikData 
              8349               3106              12626              19751 

Azimuth Visualization

# Visualize before Harmony integration


before <- DimPlot(combined_merged, 
              reduction = "umap", 
              group.by = "batch",
              label = TRUE, 
              label.box = TRUE) + 
      ggtitle("Before Harmony - By batch")


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

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

NA
NA
NA

2. Harmony Integration

# Load required libraries
library(Seurat)
library(harmony)
library(ggplot2)

combined_merged <- RunHarmony(
  combined_merged,
  group.by.vars = "batch",
  theta = 0.5  # adjust theta to control the strength of batch correction
)

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 results in harmony embeddings
harmony_embeddings <- Embeddings(combined_merged, reduction = "harmony")
head(harmony_embeddings)

                       harmony_1  harmony_2  harmony_3 harmony_4   harmony_5  harmony_6   harmony_7
L1_AAACCTGAGGGCTTCC-1 -8.4750456  0.8854753  0.1401754 2.1990896  1.99165045 -1.1816483  0.44677463
L1_AAACCTGGTGCAGGTA-1 -8.3379548 -2.6029428 -2.6011599 1.7754279  0.13045781  0.4137099  1.27850107
L1_AAACCTGGTTAAAGTG-1 -0.5410056 -3.5228316 -2.3804017 2.8374817 -0.47797340  4.3211876 -1.18451733
L1_AAACCTGTCAGGTAAA-1 -2.7911759 -0.1853321 -1.1712509 2.8107926  4.41329695  1.6508389 -1.56309915
L1_AAACCTGTCCCTGACT-1 -6.3302279  2.0987431  0.2207280 0.5620687 -0.05523098  0.4901045  0.03762563
L1_AAACCTGTCCTTCAAT-1 -5.6239711 -4.3537904 -3.1391900 1.3893146 -1.44556143  1.8877147  0.21300686
                        harmony_8 harmony_9 harmony_10 harmony_11 harmony_12 harmony_13 harmony_14
L1_AAACCTGAGGGCTTCC-1  2.36514011  3.127279  0.4602793 -1.6114981  0.5061754  2.1628814 -0.6128945
L1_AAACCTGGTGCAGGTA-1  4.27880033  5.095448  0.2312765 -4.4421756  1.7672432  1.2133324  0.8234847
L1_AAACCTGGTTAAAGTG-1  2.16975771  1.108583  0.5675762 -1.6467962 -0.3732625  1.8597098 -0.1031082
L1_AAACCTGTCAGGTAAA-1 -0.03076737  1.109587 -0.1892057 -0.1130248 -0.2365182  2.5942269  0.8414373
L1_AAACCTGTCCCTGACT-1  2.43285974  3.402622  0.1324123  0.3903623 -0.1792789  0.4252663  3.1965223
L1_AAACCTGTCCTTCAAT-1  2.01587453  2.970901  0.3737679 -1.1943186  0.8593797  1.0170968  0.9335179
                      harmony_15 harmony_16 harmony_17 harmony_18 harmony_19 harmony_20 harmony_21
L1_AAACCTGAGGGCTTCC-1  1.0398135 -0.4683957 -0.2801826 -0.8264834 -1.0171784  0.8460426  1.9405317
L1_AAACCTGGTGCAGGTA-1 -0.9251038  0.7258993 -1.4106062  0.6334174  2.3871142  1.4532488 -1.9922143
L1_AAACCTGGTTAAAGTG-1 -0.4050778 -0.2678150 -0.2848796 -1.3021338  2.0668257 -1.9703928 -2.1866791
L1_AAACCTGTCAGGTAAA-1  0.2462399 -0.6749231  0.7968253 -1.0619897 -0.3084888 -1.7052556 -3.8229159
L1_AAACCTGTCCCTGACT-1  0.6785505 -0.9361458  0.3016969 -1.6185393 -1.1868071  1.4461868  0.8042642
L1_AAACCTGTCCTTCAAT-1 -0.2026303 -0.4765092 -1.0070694 -0.5452248  0.9656883  0.2480860 -0.7882930
                      harmony_22 harmony_23 harmony_24 harmony_25 harmony_26 harmony_27  harmony_28
L1_AAACCTGAGGGCTTCC-1   4.028553  -1.033822  1.4017062  0.1117544 -1.7026521  1.1877495 -1.91034950
L1_AAACCTGGTGCAGGTA-1   1.703688  -1.541049  1.6010351 -1.2287615 -1.9013886 -0.2203379  0.65638076
L1_AAACCTGGTTAAAGTG-1   2.916394  -2.852085 -1.8670028 -3.2926887  1.7438544  3.5040060 -1.17858868
L1_AAACCTGTCAGGTAAA-1   1.662924  -1.756008 -0.1004504 -1.3040248 -0.6504728  1.5533556  0.06531409
L1_AAACCTGTCCCTGACT-1   3.572124  -0.719536  1.5383190 -0.5764302 -1.2224538  2.3430671 -0.99892629
L1_AAACCTGTCCTTCAAT-1   2.646898  -0.102884  0.8827323 -2.8045128  0.5609568  1.4242606 -0.33447445
                      harmony_29    harmony_30 harmony_31  harmony_32  harmony_33    harmony_34
L1_AAACCTGAGGGCTTCC-1 -1.0240970  2.1396213449 -1.8407529  0.35841191 -0.03769201  0.5194146142
L1_AAACCTGGTGCAGGTA-1 -1.9240000 -0.1481913470  0.3251502  1.41704379 -0.28269844  0.3608422995
L1_AAACCTGGTTAAAGTG-1  0.9065098  0.0009278782  1.7338590 -0.66634902 -0.24217725 -0.0001117634
L1_AAACCTGTCAGGTAAA-1  0.2089973  0.4462016717  1.0603453  1.72509074  0.65397057 -1.9838686622
L1_AAACCTGTCCCTGACT-1 -1.4908581  3.2544934226 -1.2669499 -0.07941294  0.58782739 -2.1388657895
L1_AAACCTGTCCTTCAAT-1 -1.7956210  0.5080979042 -1.3404014  1.05968868  2.31584181 -1.0675585158
                        harmony_35 harmony_36 harmony_37  harmony_38   harmony_39 harmony_40  harmony_41
L1_AAACCTGAGGGCTTCC-1  0.006824832  1.3664203 -0.9878709 -0.81792774  0.265949650 -0.7566405 -0.23343638
L1_AAACCTGGTGCAGGTA-1  0.148262761  0.7322868 -2.4405793 -0.08209858  0.001702334  1.1200969 -1.11369348
L1_AAACCTGGTTAAAGTG-1 -2.195343222 -0.1930475 -1.0234834  0.38330762  0.331622009  1.8097380 -0.05100157
L1_AAACCTGTCAGGTAAA-1  0.781373047  0.4907636 -0.9565876  0.73475198  0.242065893  2.9998273 -1.02473623
L1_AAACCTGTCCCTGACT-1  0.494478600  2.5959648 -0.3639646  0.62847950 -0.290949226 -0.9177849  2.60081226
L1_AAACCTGTCCTTCAAT-1  1.300231174  0.1685702 -1.3123533  0.10410045  0.657947797  2.5046247 -0.36929743
                      harmony_42  harmony_43 harmony_44 harmony_45 harmony_46 harmony_47 harmony_48
L1_AAACCTGAGGGCTTCC-1  1.0381313  1.12147875 -0.6561348  -2.245302 -0.6367506  0.8362259  0.8817999
L1_AAACCTGGTGCAGGTA-1  2.1672494  0.87381871  2.0184628  -2.299327 -0.8510006  1.3979526 -0.6793890
L1_AAACCTGGTTAAAGTG-1  0.9418567  1.38147237 -2.1894980  -2.464035 -0.6785400  1.8934146 -0.3517092
L1_AAACCTGTCAGGTAAA-1 -0.9384953 -0.09418887 -1.7604112  -3.323799 -1.0587938  2.8751137 -0.8379009
L1_AAACCTGTCCCTGACT-1  0.5233654  0.42657846 -3.1843582  -2.796636 -0.2257227  0.7560893  2.5573819
L1_AAACCTGTCCTTCAAT-1 -0.7350375  0.29228211 -0.6675299  -2.539335 -1.8033000  2.2041533 -2.2827512
                      harmony_49 harmony_50
L1_AAACCTGAGGGCTTCC-1 -1.7889981  3.7244008
L1_AAACCTGGTGCAGGTA-1 -2.0276452  1.0149594
L1_AAACCTGGTTAAAGTG-1 -3.5868594 -0.0403791
L1_AAACCTGTCAGGTAAA-1  2.6972577 -0.6643848
L1_AAACCTGTCCCTGACT-1  0.5913205  2.5870745
L1_AAACCTGTCCTTCAAT-1 -1.2440532  2.8360701

# Set the seed for clustering steps
set.seed(123)

# Run UMAP on Harmony embeddings
combined_merged <- RunUMAP(combined_merged, reduction = "harmony", dims = 1:20)

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

# Set the seed for clustering steps
set.seed(123)

# Optionally, find neighbors and clusters (if you plan to do clustering analysis)
combined_merged <- FindNeighbors(combined_merged, reduction = "harmony", dims = 1:20)

Computing nearest neighbor graph
Computing SNN

combined_merged <- FindClusters(combined_merged, reduction = "harmony", resolution = 0.5)  # Adjust resolution as needed

Warning: The following arguments are not used: reductionWarning: The following arguments are not used: reduction

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

Number of nodes: 89260
Number of edges: 2802622

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9300
Number of communities: 25
Elapsed time: 27 seconds

# Visualize UMAP
DimPlot(combined_merged, reduction = "umap", group.by = "batch", label = TRUE, label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP of Harmony-Integrated Data")

# Visualize UMAP with batch/cell line information
DimPlot(combined_merged, reduction = "umap", group.by = "batch", label = TRUE,label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP - Colored by Cell Line (After Harmony Integration)")

# Visualize UMAP with clusters
DimPlot(combined_merged, reduction = "umap", group.by = "seurat_clusters", label = TRUE,label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP - Clustered Data (After Harmony Integration)")

# # Visualize specific cell types or other metadata
# DimPlot(combined_merged, reduction = "umap", group.by = "predicted.celltype.l2", label = TRUE, pt.size = 0.5) +
#     ggtitle("UMAP - Cell Types After Harmony Integration")

3. Harmony Integration

# Load required libraries
library(Seurat)
library(harmony)
library(ggplot2)

combined_merged <- RunHarmony(
  combined_merged,
  group.by.vars = "batch",
  theta = 2  # adjust theta to control the strength of batch correction
)

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

# Check results in harmony embeddings
harmony_embeddings <- Embeddings(combined_merged, reduction = "harmony")
head(harmony_embeddings)

                       harmony_1  harmony_2  harmony_3   harmony_4  harmony_5  harmony_6   harmony_7
L1_AAACCTGAGGGCTTCC-1 -13.130652  4.1857397  2.3096055  1.07713536  1.3392787 -1.4616402 -0.74013843
L1_AAACCTGGTGCAGGTA-1  -6.331099 -4.2364777 -1.6578999  0.48072923 -0.7166837 -0.3836477  0.38212220
L1_AAACCTGGTTAAAGTG-1   4.109785 -2.8770329 -1.1426544  0.36082214 -1.9318698  4.4097104 -0.62552946
L1_AAACCTGTCAGGTAAA-1   2.513392 -0.4499775 -0.1713337  1.45979095  4.0064590  0.1245493  0.05403058
L1_AAACCTGTCCCTGACT-1 -12.293540  6.5390882  2.2526562 -0.02760549 -0.3065715  0.4991373 -1.34916146
L1_AAACCTGTCCTTCAAT-1  -1.808765 -6.5980648 -2.3571369 -0.35828323 -2.6418057  0.7235441  0.11692685
                      harmony_8  harmony_9    harmony_10 harmony_11 harmony_12 harmony_13 harmony_14
L1_AAACCTGAGGGCTTCC-1 1.9472457 -0.4227151  0.2128791057 -2.7259327  1.9450612  0.5933810 -0.1443880
L1_AAACCTGGTGCAGGTA-1 3.6129529  4.9911543 -0.5897096557 -3.6067919  0.9350275  0.7344714  0.7331883
L1_AAACCTGGTTAAAGTG-1 2.1171233 -0.6665914  0.4487200811 -2.2874269 -0.8514458  0.9915625 -0.6806862
L1_AAACCTGTCAGGTAAA-1 0.3264429  1.2188703  0.0003380007 -1.4258733 -0.3262691  1.6033884  0.1027395
L1_AAACCTGTCCCTGACT-1 2.2573817 -0.3906132 -0.2333447377 -1.1084644  1.9335322 -0.9326940  3.6146485
L1_AAACCTGTCCTTCAAT-1 1.3719350  2.1931744 -0.1742342957 -0.5659776  0.3963054  0.6789831  0.5644654
                       harmony_15  harmony_16 harmony_17  harmony_18 harmony_19 harmony_20 harmony_21
L1_AAACCTGAGGGCTTCC-1  1.73293897  0.06350125 -1.1063943 -0.42192822 -0.1404072  0.4444815  2.3206268
L1_AAACCTGGTGCAGGTA-1  0.03942636  0.27946263 -1.1957612  0.22830289  2.0420645  0.8597665 -0.2960559
L1_AAACCTGGTTAAAGTG-1 -0.10804865  0.01562314 -0.5758438 -0.57211062  2.1264347 -2.5512170 -0.5799774
L1_AAACCTGTCAGGTAAA-1  0.25607102  0.12354290  0.5567859  0.09994409 -0.1620333 -1.6793873 -1.8209371
L1_AAACCTGTCCCTGACT-1  1.54193203 -0.41550355 -0.5218507 -1.11487643 -0.3631852  1.2429335  1.3787488
L1_AAACCTGTCCTTCAAT-1  0.52099469 -0.80969897 -0.8025478 -1.14820517  0.4580647  0.2365674  1.2494401
                      harmony_22 harmony_23 harmony_24  harmony_25 harmony_26  harmony_27  harmony_28
L1_AAACCTGAGGGCTTCC-1  2.0337336  1.2278846 -0.6271286 -0.08668679 -1.2298754  0.82804153 -2.12150210
L1_AAACCTGGTGCAGGTA-1 -0.7529617 -0.9837455  0.4042462 -1.07887551 -1.4955840 -0.95293706 -0.45924109
L1_AAACCTGGTTAAAGTG-1  0.1384171 -1.7055452 -2.8815506 -2.55894150  2.3717785  2.68043547 -1.20529560
L1_AAACCTGTCAGGTAAA-1 -0.8192792 -0.3301827 -0.7682123 -0.41712867 -0.1551393  0.06414529 -0.02009605
L1_AAACCTGTCCCTGACT-1  1.3557148  2.0830590 -0.7741903 -0.81586603 -1.0934556  2.09961228 -1.10067749
L1_AAACCTGTCCTTCAAT-1  0.2603825  0.5911880 -0.0290027 -2.01524234  1.3209301  0.36924078 -1.38275285
                      harmony_29 harmony_30 harmony_31  harmony_32  harmony_33 harmony_34 harmony_35
L1_AAACCTGAGGGCTTCC-1 -0.5862336  1.2366802 -1.1071854 -0.05219346  0.13570691  1.2411198 -0.1835257
L1_AAACCTGGTGCAGGTA-1 -2.0355548  0.6448050 -0.4438547  1.71043175 -0.02151297  0.8432863 -0.6603375
L1_AAACCTGGTTAAAGTG-1  0.9791214  0.1086857  1.1965701 -0.85148608  0.05849890  0.7831148 -1.7198419
L1_AAACCTGTCAGGTAAA-1  0.2992763  0.3611554  0.1857759  1.07625925  0.49353014 -1.2399473  0.9538114
L1_AAACCTGTCCCTGACT-1 -1.0307396  2.2204170 -0.6199576 -0.44925799  0.78477521 -1.2603782  0.4220003
L1_AAACCTGTCCTTCAAT-1 -1.8477668  1.3231075 -2.2390734  0.77450113  2.33294445 -0.5910081  0.4325602
                       harmony_36 harmony_37 harmony_38 harmony_39 harmony_40  harmony_41 harmony_42
L1_AAACCTGAGGGCTTCC-1  0.02024424  0.2371840 -0.4296925  0.9888164 -0.8026483 -0.26420295  0.2784906
L1_AAACCTGGTGCAGGTA-1 -0.01272612 -0.7254534 -0.4219414  0.5047916  0.1832884 -1.26597816  1.7957999
L1_AAACCTGGTTAAAGTG-1 -1.49796704  0.5820428  0.0863170  0.7951917  0.4429106  0.06694545  0.7394944
L1_AAACCTGTCAGGTAAA-1 -0.49562046  0.3483260  0.4442329  0.1339161  0.9872057 -0.70573467 -0.9363063
L1_AAACCTGTCCCTGACT-1  1.26443861  0.9613657  0.9763342  0.3935418 -1.2540818  2.64288636 -0.2501256
L1_AAACCTGTCCTTCAAT-1 -0.65436973  0.1795634 -0.2694841  0.9906574  1.1198927 -0.29581276 -0.8238996
                       harmony_43 harmony_44  harmony_45  harmony_46  harmony_47 harmony_48 harmony_49
L1_AAACCTGAGGGCTTCC-1 -0.64679826  0.3931976 -1.06107350 -0.24485586 -0.03933205  0.4558822 -1.5062035
L1_AAACCTGGTGCAGGTA-1 -0.03624431  2.9812669 -0.25175302 -0.31867680 -0.59210736 -0.9055993 -1.4921237
L1_AAACCTGGTTAAAGTG-1  0.40631615 -0.8727721 -0.09756242 -0.40529807 -0.63758395 -0.7903837 -2.7724919
L1_AAACCTGTCAGGTAAA-1 -0.77411932 -0.3857263 -0.63462059 -0.01057857  0.87039065 -0.7590134  2.9882936
L1_AAACCTGTCCCTGACT-1 -1.40098651 -2.0307238 -1.37990812  0.25510703 -0.36354803  2.1960797  0.9042908
L1_AAACCTGTCCTTCAAT-1 -0.44909258  0.5524946 -0.16608751 -1.16350928 -0.16016030 -2.3758927 -0.7817324
                      harmony_50
L1_AAACCTGAGGGCTTCC-1  4.0122642
L1_AAACCTGGTGCAGGTA-1  1.4399945
L1_AAACCTGGTTAAAGTG-1  0.5588393
L1_AAACCTGTCAGGTAAA-1 -0.2926596
L1_AAACCTGTCCCTGACT-1  2.8513963
L1_AAACCTGTCCTTCAAT-1  3.2179203

# Set the seed for clustering steps
set.seed(123)

# Run UMAP on Harmony embeddings
combined_merged <- RunUMAP(combined_merged, reduction = "harmony", dims = 1:20)

14:49:35 UMAP embedding parameters a = 0.9922 b = 1.112
14:49:35 Read 89260 rows and found 20 numeric columns
14:49:35 Using Annoy for neighbor search, n_neighbors = 30
14:49:35 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:49:42 Writing NN index file to temp file /tmp/RtmpF9MIW5/file851d12301c54
14:49:42 Searching Annoy index using 1 thread, search_k = 3000
14:50:08 Annoy recall = 100%
14:50:09 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
14:50:11 Initializing from normalized Laplacian + noise (using RSpectra)
14:50:14 Commencing optimization for 200 epochs, with 3927432 positive edges
14:50:14 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:50:44 Optimization finished

# Set the seed for clustering steps
set.seed(123)

# Optionally, find neighbors and clusters (if you plan to do clustering analysis)
combined_merged <- FindNeighbors(combined_merged, reduction = "harmony", dims = 1:20)

Computing nearest neighbor graph
Computing SNN

combined_merged <- FindClusters(combined_merged, reduction = "harmony", resolution = 0.5)  # Adjust resolution as needed

Warning: The following arguments are not used: reductionWarning: The following arguments are not used: reduction

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

Number of nodes: 89260
Number of edges: 2671945

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9084
Number of communities: 24
Elapsed time: 29 seconds

# Visualize UMAP
DimPlot(combined_merged, reduction = "umap", group.by = "batch", label = TRUE, label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP of Harmony-Integrated Data")

# Visualize UMAP with batch/cell line information
DimPlot(combined_merged, reduction = "umap", group.by = "batch", label = TRUE,label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP - Colored by Cell Line (After Harmony Integration)")

# Visualize UMAP with clusters
DimPlot(combined_merged, reduction = "umap", group.by = "seurat_clusters", label = TRUE,label.box = T, repel = T, pt.size = 0.5) +
    ggtitle("UMAP - Clustered Data (After Harmony Integration)")

# # Visualize specific cell types or other metadata
# DimPlot(combined_merged, reduction = "umap", group.by = "predicted.celltype.l2", label = TRUE, pt.size = 0.5) +
#     ggtitle("UMAP - Cell Types After Harmony Integration")

Save the Seurat object as an Robj file

saveRDS(combined_merged, file = "Poglio_Gaydosik_combined_merged_Harmonized.rds")

3. Azimuth Annotation

library(Seurat)
library(SeuratData)
library(Azimuth)

# The RunAzimuth function can take a Seurat object as input
combined_merged <- RunAzimuth(combined_merged, reference = "pbmcref")

Warning: Overwriting miscellanous data for modelWarning: Adding a dimensional reduction (refUMAP) without the associated assay being presentWarning: Adding a dimensional reduction (refUMAP) without the associated assay being presentdetected inputs from HUMAN with id type Gene.name
reference rownames detected HUMAN with id type Gene.name
Normalizing query using reference SCT model
Warning: 113 features of the features specified were not present in both the reference query assays. 
Continuing with remaining 4887 features.Projecting cell embeddings
Error: Cannot add new cells with [[<-