-
Notifications
You must be signed in to change notification settings - Fork 0
/
2_Cite-seq_Analysis_R.R
76 lines (58 loc) · 3.36 KB
/
2_Cite-seq_Analysis_R.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
# Load Seurat file
load("Swarbrick_metastatic_lymphnode.seurat.Rdata")
# Preprocessing: Remove the error cell (last one)
citeCells <- subset(citeCells, cells =colnames(citeCells)[1:3411] )
# Load Autoencoder file
ae_data = read.csv("autoencoder_50dim_test.csv")
row.names(ae_data) = row.names([email protected])
ae_data = as.matrix(ae_data)
# Import AE to Seurat object
citeCells[["ae"]] <- CreateDimReducObject(embeddings = ae_data, key = "AE_", assay = DefaultAssay(citeCells))
# Clustering
citeCells <- FindNeighbors(citeCells,reduction="ae",dims = 1:50)
citeCells <- FindClusters(citeCells, resolution = 0.2)
# Run Umap and Visualize
citeCells = RunUMAP(object = citeCells,reduction.key = "aeumap_",
reduction = "ae",dims = 1:50,n.components = 2)
DimPlot(citeCells, reduction = "umap",label = T, pt.size = 1) + NoLegend()
################################ DE Analysis #########################################
library(tidyverse)
# Find top 10 and visualize RNA markers
citeCells.rna.markers <- FindAllMarkers(citeCells, max.cells.per.ident = 100, only.pos = TRUE)
top10 <- citeCells.rna.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(citeCells,features = top10$gene)
# Find top 10 and visualize ADT markers
DefaultAssay(citeCells) <- "ADT"
citeCells.adt.markers <- FindAllMarkers(citeCells, max.cells.per.ident = 100, only.pos = TRUE)
top10_adt <- citeCells.adt.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(citeCells,features = top10_adt$gene)
################################ Annotation #########################################
library(SingleR)
counts <- GetAssayData(citeCells, slot = "data")
edit_genename <- function(x){
x <- sub("GRCh38-","",x)
x <- sub("mm10---","",x)
return(x)
}
row.names(counts) = edit_genename(row.names(counts))
singler = CreateSinglerObject(counts, annot = NULL, "OZ", min.genes = 0,
technology = "10X", species = "Human", citation = "",
ref.list = list(), normalize.gene.length = F, variable.genes = "de",
fine.tune = F, do.signatures = F, clusters = [email protected], do.main.types = T,
reduce.file.size = T, numCores = SingleR.numCores)
human <- singler$singler[1][[1]]$SingleR.clusters.main$labels
singler = CreateSinglerObject(counts, annot = NULL, "OZ", min.genes = 0,
technology = "10X", species = "Mouse", citation = "",
ref.list = list(), normalize.gene.length = F, variable.genes = "de",
fine.tune = F, do.signatures = F, clusters = [email protected], do.main.types = T,
reduce.file.size = T, numCores = SingleR.numCores)
mouse <- singler$singler[1][[1]]$SingleR.clusters.main$labels
# Look at the heatmap, I manually annotate based on the RNA markers start with mm_ is
# Mouse cell, GRCh38_ is Human cell
# Then we got these labels:
new.cluster.ids <- c("HG_T_cells_1","MM_Macrophages","HG_B_cell","HG_T_cells_2","MM_Monocytes","HG_Endothelial_cells",
"HG_Monocyte", "MM_Epithelial cells", "HG_NK_cell","HG_Pre-B_cell_CD34-")
names(new.cluster.ids) <- levels(citeCells)
citeCells <- RenameIdents(citeCells, new.cluster.ids)
# Visualize it again
DimPlot(citeCells, reduction = "umap",label = T, pt.size = 1) + NoLegend()