harmony integration on SCT

1. load libraries

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
SS_All_samples_Merged <- load("SS_All_Sample_Merged_Azimuth_ProjectTils_singleR_ANNOTATION_on_My_UMAP0.7_HPC.Robj")

All_samples_Merged

An object of class Seurat 
63193 features across 49193 samples within 6 assays 
Active assay: SCT (26469 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. QC

Idents(All_samples_Merged) <- "cell_line"
VlnPlot(All_samples_Merged, features = c("nFeature_RNA", 
                                    "nCount_RNA", 
                                    "percent.mito"), 
                                      ncol = 3)

FeatureScatter(All_samples_Merged, 
               feature1 = "nCount_RNA", 
               feature2 = "nFeature_RNA") +
  geom_smooth(method = 'lm')

##FeatureScatter is typically used to visualize feature-feature relationships ##for anything calculated by the object, ##i.e. columns in object metadata, PC scores etc.

FeatureScatter(All_samples_Merged, 
               feature1 = "nCount_RNA", 
               feature2 = "percent.mito")+
  geom_smooth(method = 'lm')

FeatureScatter(All_samples_Merged, 
               feature1 = "nCount_RNA", 
               feature2 = "nFeature_RNA")+
  geom_smooth(method = 'lm')

4. Perform PCA

ElbowPlot(All_samples_Merged)

NA
NA

5. Harmony TEST

P1 <- DimPlot(object = All_samples_Merged, group.by = "cell_line") + 
  labs(title = 'Colored by cellline')
P1

P2 <- DimPlot(object = All_samples_Merged, group.by = "predicted.celltype.l2") + 
  labs(title = 'Colored by celltype')
P2

P3 <- DimPlot(object = All_samples_Merged, group.by = "cell_line", reduction = "pca") + 
  labs(title = 'Colored by cellline')
P3

P4 <- DimPlot(object = All_samples_Merged, group.by = "predicted.celltype.l2", reduction = "pca") + 
  labs(title = 'Colored by celltype')
P4

cowplot::plot_grid(P1, P2, P3, P4, nrow = 2)

6. Apply Harmony

All_samples_Merged_integrated <- RunUMAP(All_samples_Merged_integrated, reduction = "harmony", dims = 1:8)

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

DimPlot(object = All_samples_Merged_integrated, group.by = "predicted.celltype.l2", label = T, label.box = T, repel = T, reduction = "umap")

DimPlot(object = All_samples_Merged_integrated, group.by = "predicted.celltype.l2", label = T, label.box = T, repel = T, reduction = "harmony")

DimPlot(object = All_samples_Merged_integrated, group.by = "predicted.celltype.l2", label = T, label.box = T, repel = T, reduction = "integrated_dr")

DimPlot(object = All_samples_Merged_integrated, group.by = "predicted.celltype.l2", label = T, label.box = T, repel = T, reduction = "pca")

DimPlot(object = All_samples_Merged_integrated, group.by = "predicted.celltype.l2", label = T, label.box = T, repel = T, reduction = "ref.umap")

DimPlot(object = All_samples_Merged_integrated, group.by = "cell_line", label = T, label.box = T, repel = T, reduction = "umap")

DimPlot(object = All_samples_Merged_integrated, group.by = "seurat_clusters", label = T, label.box = T, repel = T, reduction = "umap")

7. Cell Distribution

cell_distribution_table <- table(All_samples_Merged_integrated$cell_line, All_samples_Merged_integrated$seurat_clusters)

cell_distribution_df <- as.data.frame.matrix(cell_distribution_table)


print(cell_distribution_df)

write.csv(cell_distribution_df, file = "15-2-integration_table_HARMONY-TEST2.csv", row.names = TRUE)

Save the Seurat object as an Robj file———————————————————————-

# save
#save(All_samples_Merged_integrated,file = "15-2_All_samples_Merged_integrated.Robj")

```