-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_dge.R
132 lines (100 loc) · 2.61 KB
/
plot_dge.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
library(tidyverse)
library(janitor)
library(pheatmap)
library(limma)
library(fs)
dir_create("figures")
# load data ---------------------------------------------------------------
data_matrix <-
read_delim(
"data/proteinGroups_log2_substractMedian_Batchcorrected.csv",
delim = ";"
) %>%
column_to_rownames(var = "Accession") %>%
as.matrix()
colnames(data_matrix) <- str_replace_all(
colnames(data_matrix),
c(Hpyl = "Hp", Alwof = "Aci", uninduced = "Uninduced")
)
df_data_unind_vs_treated <-
read_csv("analysis/results_limma.csv") %>%
clean_names()
# data exploration PCA & correlations -------------------------------------
plot(density(data_matrix))
boxplot(data_matrix,
las = 2,
ylab = "log2 intensity")
pheatmap(cor(data_matrix, method = "spearman"))
#plot PCA
mds <- plotMDS(
data_matrix,
gene.selection = "common",
var.explained = TRUE
)
tibble(
pc = 1:7,
variance = mds$var.explained[1:7] * 100
) %>%
ggplot(aes(pc, variance)) +
geom_col() +
geom_text(aes(label = round(variance, digits = 2)), vjust = -0.2) +
xlab("principal component number") +
ylab("variance (%)")
# visualisation of DGE results --------------------------------------------
## heatmaps ---------------------------------------------------------------
plot_comparison <- function(logfc_col, p_col, file) {
df_filtered <-
df_data_unind_vs_treated %>%
filter({{p_col}} < 0.05)
ann_row <-
df_filtered %>%
column_to_rownames("accession") %>%
mutate(
log2FC = if_else({{logfc_col}} > 0, "positive", "negative"),
.keep = "none"
)
png(
file,
height = 100,
width = 100,
unit = "mm",
res = 600
)
pheatmap(
data_matrix[df_filtered$accession,],
cluster_cols = TRUE,
show_rownames = FALSE,
annotation_colors = list(log2FC = c(negative = "blue", positive = "red")),
annotation_row = ann_row,
scale = "row"
)
dev.off()
}
dev.off()
# aci vs control
plot_comparison(
coef_t3alwofvs_control,
p_value_adj_t3alwofvs_control,
"figures/heatmap_aci_vs_cntrl_600dpi.png"
)
# hp vs control
plot_comparison(
coef_t2hpylvs_control,
p_value_adj_t2hpylvs_control,
"figures/heatmap_hp_vs_cntrl_600dpi.png"
)
# lps vs control
plot_comparison(
coef_t1lp_svs_control,
p_value_adj_t1lp_svs_control,
"figures/heatmap_lps_vs_cntrl_600dpi.png"
)
## Figure 2B volcano plot -------------------------------------------------
ggplot(
df_data_unind_vs_treated,
aes(coef_t2hpylvs_control, -log10(p_value_adj_t2hpylvs_control))
) +
geom_point() +
xlim(-3, 4) +
ylab("-log10(padj)") +
xlab("log2 fold change")