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_embed[1:10,1:10]

                         harmony_1 harmony_2  harmony_3 harmony_4  harmony_5  harmony_6   harmony_7 harmony_8 harmony_9 harmony_10
L3_B_AAACCTGAGGGATACC-1 -3.9095998  13.78526  -4.940704 3.4989230  2.9649714 -13.216809   3.3915887 -2.297467  8.096964 0.34628277
L3_B_AAACCTGCACATAACC-1 -1.3611464  21.34146 -17.044432 5.9969237 -2.9998383   3.949790  -5.6249972 -1.898653  3.561800 6.84990233
L3_B_AAACCTGCATCATCCC-1 -3.5550039  23.78250 -15.272387 4.8177437 -1.0650071   4.664416 -10.8897892 -4.944097  3.308021 7.16126840
L3_B_AAACCTGCATCCAACA-1  1.3002780  24.19582 -15.231390 2.8545064 -1.1484979   7.368370  -5.9550965 -2.142644  4.005343 1.78313977
L3_B_AAACCTGCATTCACTT-1 -0.6022814  21.95116 -12.323902 4.0463080  0.5432117   4.035378  -9.1942187 -3.512917  2.885230 3.07000318
L3_B_AAACCTGGTTCCATGA-1  1.1153356  22.02548 -15.093286 2.6848739 -1.6559538   4.798130  -4.1466364 -1.199545  3.149105 2.03709964
L3_B_AAACCTGTCAATACCG-1 -2.7448097  22.04932 -14.532976 7.4362953  0.5027631   2.688179  -5.3077889 -1.784279  8.101161 7.56392449
L3_B_AAACCTGTCACGACTA-1  1.6666778  28.62823 -20.204502 5.2013266  1.9786183   2.723314 -11.6712529 -3.229143  3.335289 0.07957117
L3_B_AAACCTGTCGTATCAG-1  1.4759493  26.66196 -16.969509 3.3140030 -0.8816951   4.965293  -4.7040881 -1.791957  1.530372 2.00108530
L3_B_AAACCTGTCGTGACAT-1  1.0025095  15.91482  -5.394434 0.1306539 -1.4207220  -0.490923   0.5903668 -1.469661  1.656719 1.74736733

#Visualize
Before1 <- DimPlot(object = All_samples_Merged, group.by = "cell_line", label = T, label.box = T, repel = T, reduction = "umap")
After1 <-  DimPlot(object = All_samples_Merged_integrated, group.by = "cell_line", label = T, label.box = T, repel = T, reduction = "umap")

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

Before1 | After1

Before2 | After2


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

NA
NA

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

```