-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path5-Eno_compute.R
230 lines (224 loc) · 7.35 KB
/
5-Eno_compute.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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
# metacell -----
setwd("/home/pingxr/Atherosis_0723/20230204_Atherosis/code/Endo/")
library(Seurat)
library(tidyverse)
source("/home/tutorial/20220802_scrna-m6A/code/computing/custom_function.R")
source("/home/tutorial/20220802_scrna-m6A/code/visualization/custom_plot_function.R")
seurat_obj_path <-
"/home/pingxr/Atherosis_0723/20230204_Atherosis/result/seurat/seurat_integration_anno.rds"
selected_cell_type <- "Endothelial"
k <- 20
outdir <- "/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/"
outdir <-
file.path(outdir, selected_cell_type, str_c("k_", k), "data")
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = T)
}
seurat_obj <- seurat_obj_path %>% readRDS()
sub_seurat_obj <- seurat_obj %>%
subset(cell_type == selected_cell_type)
sub_seurat_obj <- sub_seurat_obj %>%
cat_construct_metacells(k = k, name = selected_cell_type)
metacell_obj <- hdWGCNA::GetMetacellObject(sub_seurat_obj)
table(metacell_obj$donor)
sub_seurat_obj %>% saveRDS(file.path(outdir, "seurat_obj.rds"))
head(sub_seurat_obj)
# cor--------
setwd("/home/pingxr/Atherosis_0723/20230204_Atherosis/code/Endo/")
library(Seurat)
library(tidyverse)
source("/home/tutorial/20220802_scrna-m6A/code/computing/custom_function.R")
seurat_obj_path <-
"/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/Endothelial/k_20/data/seurat_obj.rds"
selected_genes <- c("METTL3","METTL14","WTAP","FTO","ALKBH5","IGF2BP2","IGF2BP3","YTHDF1","YTHDF2")
selected_cell_type <- "Endothelial"
k <- 20
outdir <- "/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/"
outdir <-
file.path(outdir, selected_cell_type, str_c("k_", k), "correlation")
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = T)
}
seurat_obj <- seurat_obj_path %>% readRDS()
metacell_obj <- seurat_obj %>% hdWGCNA::GetMetacellObject()
metacell_obj %>% dim()
cell_types <- names(table(metacell_obj$cell_type))
selected_genes <-
selected_genes[selected_genes %in% rownames(metacell_obj)]
all_genes <- rownames(metacell_obj)
genes_cor_res <- metacell_obj %>%
cat_correlation(feature_x = selected_genes,
feature_y = all_genes)
write_csv(genes_cor_res,
file = file.path(outdir,
"gene.cor.csv"))
category <- "H"
pathway_cor_res <- metacell_obj %>%
cat_correlation(feature_x = selected_genes,
feature_y = category,
outdir = outdir
)
write_csv(pathway_cor_res,
file = file.path(outdir,
str_c(category, ".cor.csv")))
category <- "GO:BP"
pathway_cor_res1 <- metacell_obj %>%
cat_correlation(feature_x = selected_genes,
feature_y = category,
outdir = outdir
)
write_csv(pathway_cor_res1,
file = file.path(outdir,
str_c(category, ".cor.csv")))
category <- "CP:KEGG"
pathway_cor_res2 <- metacell_obj %>%
cat_correlation(feature_x = selected_genes,
feature_y = category,
outdir = outdir
)
write_csv(pathway_cor_res2,
file = file.path(outdir,
str_c(category, ".cor.csv")))
# wgcna -------
setwd("/home/pingxr")
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
source("/home/tutorial/20220802_scrna-m6A/code/computing/custom_function.R")
selected_cell_type <- "Endothelial"
selected_genes <- c("METTL3","METTL14","WTAP","FTO","ALKBH5","IGF2BP2",
"IGF2BP3","YTHDF1","YTHDF2")
k <- 20
outdir <- "./Atherosis_0723/20230204_Atherosis/result/Endo" #注意设置路径问题
outdir <-
file.path(outdir, selected_cell_type, str_c("k_", k), "wgcna")
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = T)
}
seurat_obj <- readRDS("/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/Endothelial/k_20/data/seurat_obj.rds")
seurat_obj <- SetDatExpr(
seurat_obj,
group_name = selected_cell_type,
group.by = "cell_type",
assay = "RNA"
)
seurat_obj <- TestSoftPowers(seurat_obj,
setDatExpr = FALSE)
plot_list <-
PlotSoftPowers(seurat_obj)
p <- wrap_plots(plot_list, ncol = 2)
p
ggsave(
file.path(outdir, "test_soft_powers.pdf"),
plot = p,
height = 6,
width = 10
)
power_table <-
GetPowerTable(seurat_obj)
head(power_table)
seurat_obj <-
ConstructNetwork(seurat_obj,
soft_power = 8,
setDatExpr = FALSE,
tom_outdir = outdir,
nThreads = 10)
pdf(file.path(outdir, "plot_dendrogram.pdf"))
PlotDendrogram(seurat_obj, main = str_c(selected_cell_type, " hdWGCNA Dendrogram"))
dev.off()
TOM <- GetTOM(seurat_obj)
seurat_obj <-
ModuleEigengenes(seurat_obj,
group.by.vars = "donor")
hMEs <-
GetMEs(seurat_obj)
MEs <-
GetMEs(seurat_obj, harmonized = FALSE)
seurat_obj <-
ModuleConnectivity(seurat_obj,
group.by = "cell_type",
group_name = selected_cell_type)
p <-
PlotKMEs(seurat_obj, ncol = 5)
p
ggsave(
file.path(outdir, "plot_kmes.pdf"),
plot = p,
height = 6,
width = 10
)
modules <-
GetModules(seurat_obj)
modules %>%
write_csv(file = file.path(outdir, "modules.csv"))
head(modules[, 1:6])
names(table(modules$module))
modules %>%
dplyr::filter(gene_name %in% selected_genes) %>%
group_by(module) %>%
dplyr::count(gene_name) %>%
write_csv(file = file.path(outdir, "modules_genes.csv"))
library(enrichR)
library(GeneOverlap)
dbs <- c("GO_Biological_Process_2021",
"KEGG_2021_Human",
"MSigDB_Hallmark_2020")
seurat_obj <-
RunEnrichr(seurat_obj,
dbs = dbs,
max_genes = 100)
saveRDS(seurat_obj, file = file.path(outdir, "wgcna.rds"))
enrichrtable <- hdWGCNA::GetEnrichrTable(seurat_obj)
enrichrtable %>%
write_csv(file.path(outdir, "enrichr_table.csv"))
# deg -----
setwd("/home/pingxr/Atherosis_0723/20230204_Atherosis/code/Endo/")
library(Seurat)
library(tidyverse)
source("/home/tutorial/20220802_scrna-m6A/code/computing/custom_function.R")
seurat_obj_path <-
"/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/Endothelial/k_20/data/seurat_obj.rds"
selected_genes <- c("ALKBH5")
selected_cell_type <- "Endothelial"
k <- 20
outdir <- "/home/pingxr/Atherosis_0723/20230204_Atherosis/result/Endo/"
outdir <-
file.path(outdir, selected_cell_type, str_c("k_", k), "deg")
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = T)
}
seurat_obj <- seurat_obj_path %>% readRDS()
seurat_obj
DimPlot(seurat_obj,group.by = "cell_type")
metacell_obj <- hdWGCNA::GetMetacellObject(seurat_obj)
deg_metacell_obj <- metacell_obj |> cat_deg(
group_by = selected_genes,
cell_types = selected_cell_type,
min_pct = 0,
mode = 'median',
logfc_threshold = 0
)
table(metacell_obj@assays$RNA@counts['ALKBH5', ] > median(metacell_obj@assays$RNA@counts['ALKBH5', ]))
write_csv(deg_metacell_obj, file.path(outdir, "deg_metacell_obj.csv"))
enriched <- deg_metacell_obj |>
filter(group == "ALKBH5") |>
cat_enrich()
saveRDS(enriched, file.path(outdir, "ALKBH5_enriched.rds"))
gsea_res <- deg_metacell_obj |>
filter(group == "ALKBH5") |>
cat_gsea(category = "H"
)
saveRDS(gsea_res, file.path(outdir, "h_gsea.rds"))
gsea_res1 <- deg_metacell_obj |>
filter(group == "ALKBH5") |>
cat_gsea(category = "C2"
)
saveRDS(gsea_res1, file.path(outdir, "C2_gsea.rds"))
gsea_res2 <- deg_metacell_obj |>
filter(group == "ALKBH5") |>
cat_gsea(category = "C5"
)
saveRDS(gsea_res2, file.path(outdir, "C5_gsea.rds"))