Integration by Harmony_on_SCTransform_DATA

1. load libraries

2. Load Seurat Object

  All_samples_Merged

An object of class Seurat 
62625 features across 46976 samples within 6 assays 
Active assay: SCT (25901 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: pca, umap, integrated_dr, ref.umap

3. Data PREPERATION

ElbowPlot(All_samples_Merged)

# plot
before <- DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)



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

NA
NA
NA

4. Harmony-integration

# run Harmony -----------
All_samples_Merged_Harmony_Integrated <- All_samples_Merged %>%
  RunHarmony(group.by.vars = 'orig.ident', plot_convergence = FALSE)

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

All_samples_Merged_Harmony_Integrated@reductions

$pca
A dimensional reduction object with key PC_ 
 Number of dimensions: 50 
 Number of cells: 46976 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

$umap
A dimensional reduction object with key umap_ 
 Number of dimensions: 2 
 Number of cells: 46976 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

$integrated_dr
A dimensional reduction object with key integrateddr_ 
 Number of dimensions: 50 
 Number of cells: 46976 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

$ref.umap
A dimensional reduction object with key UMAP_ 
 Number of dimensions: 2 
 Number of cells: 46976 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$harmony
A dimensional reduction object with key harmony_ 
 Number of dimensions: 50 
 Number of cells: 46976 
 Projected dimensional reduction calculated:  TRUE 
 Jackstraw run: FALSE 
 Computed using assay: SCT 

All_samples_Merged_Harmony_Integrated.embed <- Embeddings(All_samples_Merged_Harmony_Integrated, "harmony")
All_samples_Merged_Harmony_Integrated.embed[1:10,1:10]

                       harmony_1  harmony_2 harmony_3 harmony_4  harmony_5     harmony_6  harmony_7 harmony_8  harmony_9
L1_AAACCTGAGGGCTTCC-1   3.930011  0.7522262 -2.110055  4.773869  2.1750463   0.817792877  1.5898612 -1.657756 -6.4216000
L1_AAACCTGGTGCAGGTA-1 -10.425288 -0.7127626 -2.571288  9.321215 -0.7627203   9.776243840 -0.1976164  1.364135  1.3972471
L1_AAACCTGGTTAAAGTG-1 -13.570800 -8.8580263  1.544664  6.615586  1.2933943   3.989043556  1.0894727 -2.996915  4.5499349
L1_AAACCTGTCAGGTAAA-1  -1.827420 -5.5824516  1.981894  2.077667  2.0384860   0.002298915 -2.8482661 -1.327224  4.5527668
L1_AAACCTGTCCCTGACT-1   5.005414  1.9898605 -6.825308  2.266979 -1.0807969   1.235699777  1.0048150 -1.139788 -3.4462636
L1_AAACCTGTCCTTCAAT-1 -12.327984 -5.0089178 -1.818867  6.351098 -1.1000114   5.851202967  3.6618221 -1.357876  0.8172072
L1_AAACCTGTCTTGCAAG-1 -13.943001 -3.1867878 -2.059097  6.016835 -2.2830143   7.949084787  4.0782511 -2.923581  4.0621949
L1_AAACGGGAGGCTAGAC-1  -3.347339 -5.2998819  1.413381  5.067610  0.8416488  -0.382886526 -1.6293074 -4.091839  3.4044487
L1_AAACGGGAGGGTATCG-1  -3.941779 -1.4028403  1.345789  5.063883 -1.1777797   0.477552100 -0.1966721 -3.171414  4.1788198
L1_AAACGGGAGGGTTCCC-1  12.852076  0.5849897  2.823838 -4.883418  0.9214000 -10.244227045 -2.6604167  1.136342 -2.5039421
                      harmony_10
L1_AAACCTGAGGGCTTCC-1  1.4913999
L1_AAACCTGGTGCAGGTA-1 -0.9712920
L1_AAACCTGGTTAAAGTG-1  3.1960910
L1_AAACCTGTCAGGTAAA-1  1.0732609
L1_AAACCTGTCCCTGACT-1  0.9597765
L1_AAACCTGTCCTTCAAT-1  2.4723096
L1_AAACCTGTCTTGCAAG-1  2.9507128
L1_AAACGGGAGGCTAGAC-1  2.5150597
L1_AAACGGGAGGGTATCG-1  3.7960851
L1_AAACGGGAGGGTTCCC-1 -0.0541667

# Do UMAP and clustering using ** Harmony embeddings instead of PCA **
All_samples_Merged_Harmony_Integrated <- All_samples_Merged_Harmony_Integrated %>%
   RunUMAP(reduction = 'harmony', dims = 1:13) 

15:13:57 UMAP embedding parameters a = 0.9922 b = 1.112
15:13:57 Read 46976 rows and found 13 numeric columns
15:13:57 Using Annoy for neighbor search, n_neighbors = 30
15:13:57 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:14:00 Writing NN index file to temp file /tmp/Rtmp2eQfbV/filef36db6c1e4956
15:14:00 Searching Annoy index using 1 thread, search_k = 3000
15:14:11 Annoy recall = 100%
15:14:12 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:14:14 Initializing from normalized Laplacian + noise (using RSpectra)
15:14:15 Commencing optimization for 200 epochs, with 1955416 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:14:30 Optimization finished

  

# visualize 


after <- DimPlot(All_samples_Merged_Harmony_Integrated, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)

before|after

DimPlot(All_samples_Merged_Harmony_Integrated,reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)

DimPlot(All_samples_Merged_Harmony_Integrated,reduction = "umap", group.by = "SCT_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)

DimPlot(All_samples_Merged_Harmony_Integrated, reduction = "umap", group.by = "predicted.celltype.l2", label = TRUE, label.box = TRUE, repel = TRUE)

FeaturePlot

myfeatures <- c("CD3E", "CD4", "CD8A", "CD8B", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FOXP3", "TIGIT", "KIR3DL2", "GZMA", "CCL17", "CCL5", "CD52", "CD7", "CD26")

FeaturePlot(All_samples_Merged_Harmony_Integrated, reduction = "umap", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()

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

5. Save the Seurat object as an Robj file

# save(All_samples_Merged_Harmony_Integrated, file = "All_samples_Merged_Harmony_Integrated.Robj")