-
Notifications
You must be signed in to change notification settings - Fork 0
/
03_reprocess_resolutions.R
61 lines (52 loc) · 1.75 KB
/
03_reprocess_resolutions.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
library(Seurat)
library(dplyr)
library(ggplot2)
source("process_seu.R")
source("explore_seu.R")
source("enrichment.R")
datadir <- "data/"
outdir <- "outs/"
seu <- readRDS(paste0(datadir, "seu_filtered.RDS")) %>%
process_seu(
feature_type = "var",
nfeatures = 3000,
ndims = 30,
n_neighbours = 30,
resolution = 1
)
head(seu)
seu <- FindClusters(seu, resolution = 0.6)
colnames([email protected])[match("seurat_clusters", colnames([email protected]))] <- "clusters_lowres"
seu <- FindClusters(seu, resolution = 1.4)
colnames([email protected])[match("seurat_clusters", colnames([email protected]))] <- "clusters_highres"
rename(clusters_medres = clusters_umap) %>%
select(1:7, clusters_lowres, clusters_medres, clusters_highres)
saveRDS(seu, paste0(datadir, "seu_processed.RDS"))
explore_seu(seu = seu, shiny.dir = paste0(outdir, "Shota_Jul2022_resolutions"))
markers_enrichment <- function(seu, idents, id, dir_markers, dir_enrichment, organism){
Idents(seu) <- idents
message("Finding markers for ", idents)
markers <- FindAllMarkers(seu)
write.csv(markers, paste0(dir_markers, "/markers_", id, ".csv"))
for(cl in levels(Idents(seu))){
message("Running enrichment for cluster ", cl)
res <- enrichment(
genes = markers$gene[markers$cluster == cl],
background = rownames(seu)[rowSums(seu@assays$RNA@counts[, Idents(seu) == cl]) > 0],
file = paste0(dir_enrichment, "/enrichment_", id, "_cl_", cl, ".csv"),
organism = organism
)
}
invisible(seu)
}
for(id in c("lowres", "medres", "highres")){
markers_enrichment(
seu = seu,
idents = paste0("clusters_", id),
id = id,
dir_markers = outdir,
dir_enrichment = paste0(outdir, "/enrichment/", id),
organism = "hsapiens"
)
}