-
Notifications
You must be signed in to change notification settings - Fork 0
/
magma_gene_sets_plot.R
executable file
·146 lines (127 loc) · 5.58 KB
/
magma_gene_sets_plot.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
#!/usr/bin/env Rscript
library(dplyr)
library(ggplot2)
library(scales)
args = commandArgs(trailingOnly = TRUE)
prefix = "iPSYCH2015_EUR_"
suffix = "_C5_GO.gsa.out"
xdx = read.table(paste0(prefix, "xDx", suffix), header = T)
xdx = xdx %>%
filter(NGENES < 600) %>%
arrange(P) %>%
head(20) %>%
select(VARIABLE, FULL_NAME, BETA, BETA_STD, P) %>%
rename(GO = VARIABLE) %>%
mutate(TRAIT = "xDx") %>%
mutate(GWAS = "Case vs Cohort") %>%
mutate(BETA_HIGH = BETA + 1.96 * BETA_STD) %>%
mutate(BETA_LOW = BETA - 1.96 * BETA_STD)
xdx$GO = gsub("_.*", "", xdx$GO)
xdx$GO = gsub("^GO", "", xdx$GO)
xdx$FULL_NAME = substr(xdx$FULL_NAME, 6, nchar(xdx$FULL_NAME))
xdx$FULL_NAME = gsub("_", " ", xdx$FULL_NAME)
case_case = xdx[FALSE,]
pairwise = xdx[FALSE,]
for (trait in c("ADHD", "ANO", "AUT", "BIP", "MDD", "SCZ")) {
file = read.table(paste0(prefix, trait, suffix), header = T)
file_cc = read.table(paste0(prefix, trait, "_CC", suffix), header = T)
top_5 = file %>% arrange(P) %>% head(5) %>% select(FULL_NAME)
top_5_cc = file_cc %>% arrange(P) %>% head(5) %>% select(FULL_NAME)
to_subset = rbind(top_5, top_5_cc) %>% unique()
file = inner_join(to_subset, file, by = c("FULL_NAME")) %>%
filter(NGENES < 600) %>%
select(VARIABLE, FULL_NAME, BETA, BETA_STD, P) %>%
rename(GO = VARIABLE) %>%
mutate(GWAS = "Case vs Cohort") %>%
mutate(BETA_HIGH = BETA + 1.96 * BETA_STD) %>%
mutate(BETA_LOW = BETA - 1.96 * BETA_STD)
file_cc = inner_join(to_subset, file_cc, by = c("FULL_NAME")) %>%
filter(NGENES < 600) %>%
select(VARIABLE, FULL_NAME, BETA, BETA_STD, P) %>%
rename(GO = VARIABLE) %>%
mutate(GWAS = "Case vs Other Cases") %>%
mutate(BETA_HIGH = BETA + 1.96 * BETA_STD) %>%
mutate(BETA_LOW = BETA - 1.96 * BETA_STD)
merged = rbind(file, file_cc) %>% mutate(TRAIT = trait)
case_case = rbind(case_case, merged)
}
case_case$GO = gsub("_.*", "", case_case$GO)
case_case$GO = gsub("^GO", "", case_case$GO)
case_case$FULL_NAME = substr(case_case$FULL_NAME, 6, nchar(case_case$FULL_NAME))
case_case$FULL_NAME = gsub("_", " ", case_case$FULL_NAME)
for (trait in c("ADHD_AUT", "ADHD_ANO", "ADHD_BIP", "ADHD_MDD", "ADHD_SCZ",
"ANO_AUT", "ANO_BIP", "ANO_MDD", "ANO_SCZ",
"AUT_BIP", "AUT_MDD", "AUT_SCZ",
"BIP_MDD", "BIP_SCZ",
"MDD_SCZ")) {
file = read.table(paste0(prefix, trait, "_CC", suffix), header = T)
file = file %>%
filter(NGENES < 600) %>%
arrange(P) %>%
head(5) %>%
select(VARIABLE, FULL_NAME, BETA, BETA_STD, P) %>%
rename(GO = VARIABLE) %>%
mutate(TRAIT = trait) %>%
mutate(GWAS = "Case vs Case Pairwise") %>%
mutate(BETA_HIGH = BETA + 1.96 * BETA_STD) %>%
mutate(BETA_LOW = BETA - 1.96 * BETA_STD)
pairwise = rbind(pairwise, file)
}
pairwise$GO = gsub("_.*", "", pairwise$GO)
pairwise$GO = gsub("^GO", "", pairwise$GO)
pairwise$FULL_NAME = substr(pairwise$FULL_NAME, 6, nchar(pairwise$FULL_NAME))
pairwise$FULL_NAME = gsub("_", " ", pairwise$FULL_NAME)
png(paste0(args[1], "_xDx.png"), res = 300, width = 8, height = 8, units = "in")
ggplot(xdx, aes(y = FULL_NAME, x = BETA, shape = GO, fill = -log10(P))) +
geom_point(size = 3) +
geom_errorbarh(aes(xmin = BETA - 1.96 * BETA_STD,
xmax = BETA + 1.96 * BETA_STD),
height = 0.01, linewidth = 0.5) +
theme_classic() +
theme(axis.text.x = element_text(face = "bold"),
axis.text.y = element_text(face = "bold")) +
ylab("") +
scale_shape_manual(values = c(21, 22, 23)) +
scale_fill_gradient(low = "blue", high = "red")
dev.off()
png(paste0(args[1], "_Case_Case.png"), res = 300, width = 15, height = 15, units = "in")
ggplot(case_case, aes(y = FULL_NAME,
x = BETA,
shape = GO,
fill = -log10(P))) +
geom_point(size = 3, position = position_dodge(width = 0.9)) +
geom_errorbarh(aes(xmin = BETA - 1.96 * BETA_STD,
xmax = BETA + 1.96 * BETA_STD),
height = 0.01,
position = position_dodge(width = 0.9)) +
theme_bw() +
theme(axis.text.x = element_text(face = "bold"),
axis.text.y = element_text(size = 8,
hjust = 1,
face = "bold"),
legend.position = "bottom") +
scale_y_discrete(labels = label_wrap(40), guide = guide_axis(n.dodge = 2)) +
facet_grid(TRAIT ~ GWAS, scales = "free", space = "free") +
scale_fill_gradient(high = "red", low = "blue") +
ylab("") +
scale_shape_manual(values = c(21, 22, 23)) +
geom_vline(xintercept = 0, lty = 2)
dev.off()
png(paste0(args[1], "_Pairwise.png"), res = 300, width = 15, height = 15, units = "in")
ggplot(pairwise, aes(y = FULL_NAME, x = BETA, shape = GO, fill = -log10(P))) +
geom_point(size = 3) +
geom_errorbarh(aes(xmin = BETA - 1.96 * BETA_STD,
xmax = BETA + 1.96 * BETA_STD),
height = 0.01) +
theme_classic() +
theme(axis.text.x = element_text(face = "bold"),
axis.text.y = element_text(size = 8,
hjust = 1,
face = "bold"),
legend.position = "bottom") +
scale_y_discrete(labels = label_wrap(30)) +
ylab("") +
facet_wrap(TRAIT ~ ., scales = "free") +
scale_shape_manual(values = c(21, 22, 23)) +
scale_fill_gradient(low = "blue", high = "red")
dev.off()