Integration by CCA_by_K

1. load libraries

2. Load Seurat Object

# load("AllSample_corrected_AzimuthAnnotated_L1.Robj")

# merged_seurat_filtered <- All_samples_Merged

merged_seurat_filtered

An object of class Seurat 
36724 features across 49193 samples within 5 assays 
Active assay: RNA (36601 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 4 other assays present: ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 4 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap

3. Data PREPERATION

# perform standard workflow steps to figure out if we see any batch effects --------
# merged_seurat_filtered <- NormalizeData(object = merged_seurat_filtered, verbose = FALSE)
# merged_seurat_filtered <- FindVariableFeatures(object = merged_seurat_filtered , selection.method = "vst", nfeatures = 2000,verbose = FALSE)
#                                                
# merged_seurat_filtered <- ScaleData(object = merged_seurat_filtered)
# 
# 
# Variable_genes <- merged_seurat_filtered@assays$RNA@var.features
# 
# Variables_genes_after_exclusion <- Variable_genes[!grepl("^HLA-|^Xist", Variable_genes)]
#   
# merged_seurat_filtered <- RunPCA(merged_seurat_filtered,
#                          features = Variables_genes_after_exclusion,
#                          do.print = TRUE, 
#                          pcs.print = 1:5, 
#                          genes.print = 15)

ElbowPlot(merged_seurat_filtered)

# merged_seurat_filtered <- FindNeighbors(object = merged_seurat_filtered, dims = 1:20, verbose = FALSE)
# merged_seurat_filtered <- FindClusters(object = merged_seurat_filtered, resolution = 1.2)
# merged_seurat_filtered <- RunUMAP(object = merged_seurat_filtered, dims = 1:20)


# plot
p1 <- DimPlot(merged_seurat_filtered, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)
p2 <- DimPlot(merged_seurat_filtered, reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)

grid.arrange(p1, p2, ncol = 2, nrow = 2)

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

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

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

4. CCA-integration

# perform integration to correct for batch effects ------
# obj.list <- SplitObject(merged_seurat_filtered, split.by = "orig.ident")
# for(i in 1:length(obj.list)){
#   obj.list[[i]] <- NormalizeData(object = obj.list[[i]], verbose = FALSE)
#   obj.list[[i]] <- FindVariableFeatures(object = obj.list[[i]], selection.method = "vst", nfeatures = 2000,verbose = FALSE)
# }



# select integration features
# features <- SelectIntegrationFeatures(object.list = obj.list)



# find integration anchors (CCA)
# # anchors <- FindIntegrationAnchors(object.list = obj.list,
#                        anchor.features = features)

# integrate data
# seurat.integrated <- IntegrateData(anchorset = anchors)


# Scale data, run PCA and UMAP and visualize integrated data
# seurat.integrated <- ScaleData(object = seurat.integrated)
# 
# seurat.integrated <-  RunPCA(merged_seurat_filtered,
#                          features = Variables_genes_after_exclusion,
#                          do.print = TRUE, 
#                          pcs.print = 1:5, 
#                          genes.print = 15)
# 
 seurat.integrated <- RunUMAP(object = seurat.integrated, dims = 1:12)

16:44:49 UMAP embedding parameters a = 0.9922 b = 1.112
16:44:49 Read 49193 rows and found 12 numeric columns
16:44:49 Using Annoy for neighbor search, n_neighbors = 30
16:44:49 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:44:52 Writing NN index file to temp file /tmp/Rtmpo9cMvG/filec773d7ed006a3
16:44:52 Searching Annoy index using 1 thread, search_k = 3000
16:45:03 Annoy recall = 100%
16:45:04 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
16:45:06 Initializing from normalized Laplacian + noise (using RSpectra)
16:46:29 Commencing optimization for 200 epochs, with 2051184 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:46:45 Optimization finished

p3 <- DimPlot(seurat.integrated, reduction = "umap", group.by = "cell_line", label = TRUE, label.box = TRUE, repel = TRUE)
p4 <- DimPlot(seurat.integrated,reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)


grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)

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

DimPlot(seurat.integrated, reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE, repel = TRUE)

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

FeaturePlot

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")

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

6. Save the Seurat object as an Robj file

save(seurat.integrated, file = "seurat.integrated_K_code_CCA.Robj")