-
Notifications
You must be signed in to change notification settings - Fork 0
/
funcomics_dea.R
132 lines (79 loc) · 2.84 KB
/
funcomics_dea.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
# Copyright (c) [2020] [Ricardo O. Ramirez Flores]
#'Makes dotplots of differential expression analysis
#'
library(tidyverse)
library(Seurat)
library(clustree)
library(progeny)
library(dorothea)
library(cowplot)
source("visiumtools/reading2clustering.R")
source("visiumtools/funcomics.R")
source("visiumtools/differential_test.R")
source("visiumtools/misty_pipelines.R")
source("visiumtools/misty_utils.R")
sample_ids = set_names(c("V19S23-097_A1","V19S23-097_B1"))
# Fix dea results
for(slide in sample_ids){
print(slide)
dea_file = sprintf("results/single_slide/%s/%s_diff_features_ldvgann.rds",
slide,slide)
dea_res = readRDS(dea_file)
dea_res_fixed = lapply(dea_res,function(x){
new_res = x %>% dplyr::mutate(cluster = as.numeric(as.character(cluster))) %>%
arrange(cluster,p_val_adj)
return(new_res)
})
saveRDS(dea_res_fixed, file = dea_file)
}
# Plot general features
for(slide in sample_ids){
print(slide)
slide_file = sprintf("./results/single_slide/%s/%s.rds",
slide,slide)
visium_slide = readRDS(slide_file)
plot_file = sprintf("./results/single_slide/%s/%s_funcomics_feat.pdf",
slide,slide)
func_assays = c("progeny","ctscores")
pdf(height = 12, width = 15,
file = plot_file)
for(a in Assays(visium_slide)[Assays(visium_slide) %in% func_assays]){
DefaultAssay(visium_slide) = a
plot(SpatialPlot(visium_slide,image.alpha = 0,
features = rownames(visium_slide),
pt.size.factor = 6))
}
dev.off()
}
# First generate excel file with summaries:
for(slide in sample_ids){
print(slide)
dea_file = sprintf("results/single_slide/%s/%s_diff_features_ldvgann.rds",
slide,slide)
dea_excel = sprintf("results/single_slide/%s/%s_diff_features_ldvgann.xlsx",
slide,slide)
dea_res = readRDS(dea_file)
diff_expr_xlsx(differential_features = dea_res,
only_pos = T,
excel_out = dea_excel)
}
# Generate dotplots
for(slide in sample_ids){
print(slide)
slide_file = sprintf("./results/single_slide/%s/%s.rds",
slide,slide)
visium_slide = readRDS(slide_file)
dea_file = sprintf("results/single_slide/%s/%s_diff_features_ldvgann.rds",
slide,slide)
dea_res = readRDS(dea_file)
dea_plots = sprintf("results/single_slide/%s/%s_diff_features_dplots_ldvgann.pdf",
slide,slide)
dot_plots = get_dot_list(visium_slide = visium_slide,
dea_res = dea_res,
top_genes = 5,
p_val_thrsh = 0.05)
pdf(file = dea_plots, height = 15, width = 10)
lapply(dot_plots,plot)
dev.off()
}