-
Notifications
You must be signed in to change notification settings - Fork 1
/
plot_simu.R
412 lines (336 loc) · 19 KB
/
plot_simu.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
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
setwd("/Users/cmx/Desktop/iGREX/simulation")
library(reshape2)
library(ggplot2)
library(gridExtra)
################################################ 012 mat: SNRy ~ n1 ###############################################
out1 <- read.table("geno_5prop23_SNRy0.3_n800_4000.txt",header = T)
out2 <- read.table("geno_5prop23_SNRy0.3_n1000_4000.txt",header = T)
out3 <- read.table("geno_5prop23_SNRy0.3_n2000_4000.txt",header = T)
out4 <- read.table("geno_5prop23_SNRy0.2_n800_4000.txt",header = T)
out5 <- read.table("geno_5prop23_SNRy0.2_n1000_4000.txt",header = T)
out6 <- read.table("geno_5prop23_SNRy0.2_n2000_4000.txt",header = T)
out7 <- read.table("geno_5prop23_SNRy0.1_n800_4000.txt",header = T)
out8 <- read.table("geno_5prop23_SNRy0.1_n1000_4000.txt",header = T)
out9 <- read.table("geno_5prop23_SNRy0.1_n2000_4000.txt",header = T)
n1 <- c(replicate(3,c(800,1000,2000)))
SNRy <- c(sapply(c(0.3,0.2,0.1),rep,3))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:9){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:10],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",REML0="REML0",MoM="MoM",MoM0="MoM0",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","GREX","Alternative"),n1=n1[i],SNRy=SNRy[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVEt","method","component","n1","PVEy")
out_all$method <- factor(out_all$method,levels = c("REML0","MoM0","REML","MoM","IGREX-s"))
Pa <- ggplot(data=out_all,aes(x=component,y=PVEt,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
facet_grid(PVEy~n1, labeller = label_bquote(cols="n"[r]*"="*.(n1),rows="PVE"[y]*"="*.(PVEy))) +
scale_fill_brewer(palette="RdBu",labels=c(REML="REML",REML0=bquote(REML[0]),Mom="MoM",MoM0=bquote(MoM[0]))) +
labs(title="a",y=expression(PVE))+
# scale_fill_discrete(labels=c(REML="REML",REML0=expression(REML[0]),Mom="MoM",MoM0=expression(MoM[0]),SummaryStat="SummaryStat"))+
theme(text = element_text(size=20))
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1_geno_PVEy_n1.pdf",width=13,height=10)
Pa
dev.off()
################################################ piB ###############################################
out1 <- read.table("sparsityB0.2_sparsityG0.2_SNRy0.3_n1000_4000.txt",header = T)
out2 <- read.table("sparsityB0.5_sparsityG0.2_SNRy0.3_n1000_4000.txt",header = T)
out3 <- read.table("sparsityB0.8_sparsityG0.2_SNRy0.3_n1000_4000.txt",header = T)
sparsityB <- c(0.2,0.5,0.8)
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:3){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:6],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",MoM="MoM",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","GREX","Alternative"),sparsityB=sparsityB[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVEt","method","component","sparsityB")
out_all$method <- factor(out_all$method,levels = c("REML","MoM","IGREX-s"))
Pb <- ggplot(data=out_all,aes(x=component,y=PVEt,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
labs(title="b",y=expression(PVE))+
facet_grid(.~sparsityB, labeller = label_bquote(cols=pi[beta]*"="*.(sparsityB))) +
scale_fill_brewer(palette="RdBu") +
theme(text = element_text(size=20))
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1_piB.pdf",width=13,height=5)
Pb
dev.off()
################################################ piG ###############################################
out1 <- read.table("sparsityB0.2_sparsityG0.2_SNRy0.3_n1000_4000.txt",header = T)
out2 <- read.table("sparsityB0.2_sparsityG0.5_SNRy0.3_n1000_4000.txt",header = T)
out3 <- read.table("sparsityB0.2_sparsityG0.8_SNRy0.3_n1000_4000.txt",header = T)
sparsityG <- c(0.2,0.5,0.8)
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:3){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:6],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",MoM="MoM",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","GREX","Alternative"),sparsityG=sparsityG[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVEt","method","component","sparsityG")
out_all$method <- factor(out_all$method,levels = c("REML","MoM","IGREX-s"))
Pc <- ggplot(data=out_all,aes(x=component,y=PVEt,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
labs(title="c",y=expression(PVE))+
facet_grid(.~sparsityG, labeller = label_bquote(cols=pi[alpha]*"="*.(sparsityG))) +
scale_fill_brewer(palette="RdBu") +
theme(text = element_text(size=20))
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1_piG.pdf",width=13,height=5)
Pc
dev.off()
############################################### 012 mat: auto-correlation ###############################################
out1 <- read.table("geno_rho0.1_5prop23_SNRy0.3_n1000_4000.txt",header = T)
out2 <- read.table("geno_rho0.3_5prop23_SNRy0.3_n1000_4000.txt",header = T)
out3 <- read.table("geno_rho0.5_5prop23_SNRy0.3_n1000_4000.txt",header = T)
out4 <- read.table("geno_rho0.8_5prop23_SNRy0.3_n1000_4000.txt",header = T)
# rho <- c(replicate(3,c(0.1,0.3,0.5,0.8)))
rho <- c(0.1,0.3,0.5,0.8)
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:4){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:10],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",REML0="REML0",MoM="MoM",MoM0="MoM0",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","GREX","Alternative"),rho=rho[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVEt","method","component","rho")
out_all$method <- factor(out_all$method,levels = c("REML0","MoM0","REML","MoM","IGREX-s"))
Pd <- ggplot(data=out_all,aes(x=component,y=PVEt,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
facet_grid(.~rho, labeller = label_bquote(cols=rho*"="*.(rho))) +
scale_fill_brewer(palette="RdBu",labels=c(REML="REML",REML0=bquote(REML[0]),Mom="MoM",MoM0=bquote(MoM[0]),SummaryStat="IGREX-s")) +
coord_cartesian(ylim = c(0.1,0.45)) +
labs(title="d",y=expression(PVE))+
theme(text = element_text(size=20))
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1_geno_autocorr.pdf",width=26,height=10)
Pd
dev.off()
################################################ 2 VC compare with RHOGE: n1 ###############################################
library(ggplot2)
library(reshape2)
out1 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.1_n800_4000.txt",header = T)
out2 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.3_n800_4000.txt",header = T)
out3 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.5_n800_4000.txt",header = T)
out4 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.7_n800_4000.txt",header = T)
out5 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.9_n800_4000.txt",header = T)
out6 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.1_n1000_4000.txt",header = T)
out7 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.3_n1000_4000.txt",header = T)
out8 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.5_n1000_4000.txt",header = T)
out9 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.7_n1000_4000.txt",header = T)
out10 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.9_n1000_4000.txt",header = T)
out11 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.1_n2000_4000.txt",header = T)
out12 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.3_n2000_4000.txt",header = T)
out13 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.5_n2000_4000.txt",header = T)
out14 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.7_n2000_4000.txt",header = T)
out15 <- read.table("iGREXvsRHOGE_SNRz0.2_SNRy0.9_n2000_4000.txt",header = T)
SNRy <- c(replicate(6,c(0.1,0.3,0.5,0.7,0.9)))
n1 <- c(replicate(2,c(sapply(c(800,1000,2000),rep,5))))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:15){
out <- get(paste("out",i,sep=""))
out <- melt(out[,c(1,2,3,5)],id.vars =NULL)
out <- cbind(out,SNRy=SNRy[i],n1=n1[i])
out_all <- rbind(out_all,out)
}
names(out_all) <- c("method","PVEy","PVE","n1")
out_all$PVE <- as.factor(out_all$PVE)
levels(out_all$method) <- c("REML","MoM","IGREX-s","RhoGE")
Pe <- ggplot(data=out_all,aes(x=PVE,y=PVEy,fill=method)) +
geom_boxplot() +
geom_hline(aes(yintercept = 0.2),linetype="dashed",color="blue") +
facet_grid(.~n1,scales="free", labeller = label_bquote(cols="n"[r]*"="*.(n1))) +
scale_fill_brewer(palette="RdBu") +
labs(title="e",x=expression(PVE[y]), y=expression(PVE[GREX]))+
theme(text = element_text(size=20))
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1_RHOGE_n1_PVEt.pdf",width=26,height=10)
Pe
dev.off()
gA <- ggplotGrob(Pa)
gA$layout$l[gA$layout$name == "title"] <- 1
gB <- ggplotGrob(Pb)
gB$layout$l[gB$layout$name == "title"] <- 1
gC <- ggplotGrob(Pc)
gC$layout$l[gC$layout$name == "title"] <- 1
gD <- ggplotGrob(Pd)
gD$layout$l[gD$layout$name == "title"] <- 1
gE <- ggplotGrob(Pe)
gE$layout$l[gE$layout$name == "title"] <- 1
Fig1 <- grid.arrange(gA,gB,gC,gD,gE,layout_matrix=cbind(c(1,1,4,5),c(2,3,4,5)))
ggsave("/Users/cmx/Desktop/Predixcan/iGREX_paper/Fig1.pdf",Fig1,width=20,height=20)
#############################################################################################################
################################################ Supplementary ##############################################
#############################################################################################################
################################################ PVEg ~ n1 ###############################################
out1 <- read.table("5prop14_SNRy0.3_n800_4000.txt",header = T)
out2 <- read.table("5prop14_SNRy0.3_n1000_4000.txt",header = T)
out3 <- read.table("5prop14_SNRy0.3_n2000_4000.txt",header = T)
out4 <- read.table("5prop23_SNRy0.3_n800_4000.txt",header = T)
out5 <- read.table("5prop23_SNRy0.3_n1000_4000.txt",header = T)
out6 <- read.table("5prop23_SNRy0.3_n2000_4000.txt",header = T)
out7 <- read.table("5prop32_SNRy0.3_n800_4000.txt",header = T)
out8 <- read.table("5prop32_SNRy0.3_n1000_4000.txt",header = T)
out9 <- read.table("5prop32_SNRy0.3_n2000_4000.txt",header = T)
out10 <- read.table("5prop41_SNRy0.3_n800_4000.txt",header = T)
out11 <- read.table("5prop41_SNRy0.3_n1000_4000.txt",header = T)
out12 <- read.table("5prop41_SNRy0.3_n2000_4000.txt",header = T)
n1 <- c(replicate(4,c(800,1000,2000)))
PVEg <- c(sapply(c(0.1,0.2,0.3,0.4),rep,3))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:12){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:10],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",REML0="REML0",MoM="MoM",MoM0="MoM0",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","gene","SNP"),n1=n1[i],PVEg=PVEg[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVE","method","component","n1","PVEg")
out_all$method <- factor(out_all$method,levels = c("REML0","MoM0","REML","MoM","IGREX-s"))
levels(out_all$component) <- c("GREX","Alternative")
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Supp_PVEg_n1.pdf",width=13,height=10)
P <- ggplot(data=out_all,aes(x=component,y=PVE,fill=method)) +
geom_boxplot() +
geom_hline(aes(yintercept = PVEg),linetype="dashed",color="blue") +
geom_hline(aes(yintercept = 0.5-PVEg),linetype="dashed",color="red") +
facet_grid(PVEg~n1, labeller = label_bquote(cols="n"[r]*"="*.(n1),rows="PVE"[GREX]*"="*.(PVEg))) +
scale_fill_brewer(palette="RdBu",labels=c(REML="REML",REML0=bquote(REML[0]),Mom="MoM",MoM0=bquote(MoM[0]))) +
labs(y=expression(PVE))+
theme(text = element_text(size=20),legend.position = "bottom")
P
dev.off()
############################################### piG ~ n1 ###############################################
out1 <- read.table("sparsityB0.2_sparsityG0.2_SNRy0.3_n800_4000.txt",header = T)
out2 <- read.table("sparsityB0.2_sparsityG0.5_SNRy0.3_n800_4000.txt",header = T)
out3 <- read.table("sparsityB0.2_sparsityG0.8_SNRy0.3_n800_4000.txt",header = T)
out4 <- read.table("sparsityB0.2_sparsityG0.2_SNRy0.3_n1000_4000.txt",header = T)
out5 <- read.table("sparsityB0.2_sparsityG0.5_SNRy0.3_n1000_4000.txt",header = T)
out6 <- read.table("sparsityB0.2_sparsityG0.8_SNRy0.3_n1000_4000.txt",header = T)
out7 <- read.table("sparsityB0.2_sparsityG0.2_SNRy0.3_n2000_4000.txt",header = T)
out8 <- read.table("sparsityB0.2_sparsityG0.5_SNRy0.3_n2000_4000.txt",header = T)
out9 <- read.table("sparsityB0.2_sparsityG0.8_SNRy0.3_n2000_4000.txt",header = T)
sparsityG <- c(replicate(4,c(0.2,0.5,0.8)))
n1 <- c(sapply(c(800,1000,2000),rep,3))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:9){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:6],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",MoM="MoM",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","gene","SNP"),sparsityG=sparsityG[i],n1=n1[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVE","method","component","sparsity","n1")
out_all$method <- factor(out_all$method,levels = c("REML","MoM","IGREX-s"))
levels(out_all$component) <- c("GREX","Alternative")
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Supp_piG_n1.pdf",width=13,height=10)
P <- ggplot(data=out_all,aes(x=component,y=PVE,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
facet_grid(sparsity~n1, labeller = label_bquote(cols="n"[r]*"="*.(n1),rows=pi[alpha]*"="*.(sparsity))) +
scale_fill_brewer(palette="RdBu",labels=c(REML="REML",REML0=bquote(REML[0]),Mom="MoM",MoM0=bquote(MoM[0]))) +
labs(y=expression(PVE))+
theme(text = element_text(size=20),legend.position = "bottom")
P
dev.off()
################################################ piB ~ SNRy ###############################################
out1 <- read.table("sparsityB0.2_SNRy0.1_n800_4000.txt",header = T)
out2 <- read.table("sparsityB0.5_SNRy0.1_n800_4000.txt",header = T)
out3 <- read.table("sparsityB0.8_SNRy0.1_n800_4000.txt",header = T)
out4 <- read.table("sparsityB0.2_SNRy0.2_n800_4000.txt",header = T)
out5 <- read.table("sparsityB0.5_SNRy0.2_n800_4000.txt",header = T)
out6 <- read.table("sparsityB0.8_SNRy0.2_n800_4000.txt",header = T)
out7 <- read.table("sparsityB0.2_SNRy0.3_n800_4000.txt",header = T)
out8 <- read.table("sparsityB0.5_SNRy0.3_n800_4000.txt",header = T)
out9 <- read.table("sparsityB0.8_SNRy0.3_n800_4000.txt",header = T)
out10 <- read.table("sparsityB0.2_SNRy0.5_n800_4000.txt",header = T)
out11 <- read.table("sparsityB0.5_SNRy0.5_n800_4000.txt",header = T)
out12 <- read.table("sparsityB0.8_SNRy0.5_n800_4000.txt",header = T)
sparsity <- c(replicate(4,c(0.2,0.5,0.8)))
SNRy <- c(sapply(c(0.1,0.2,0.3,0.5),rep,3))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:12){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:6],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",MoM="MoM",ss="IGREX-s"),
component=ifelse(substr(out[,1],4,4)=="g","gene","SNP"),sparsity=sparsity[i],SNRy=SNRy[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVE","method","component","sparsity","SNRy")
out_all$method <- factor(out_all$method,levels = c("REML","MoM","IGREX-s"))
levels(out_all$component) <- c("GREX","Alternative")
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Supp_piB_PVEy.pdf",width=13,height=10)
P <- ggplot(data=out_all,aes(x=component,y=PVE,fill=method)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
facet_grid(sparsity~SNRy, labeller = label_bquote(cols="PVE"[y]*"="*.(SNRy),rows=pi[beta ]*"="*.(sparsity))) +
scale_fill_brewer(palette="RdBu",labels=c(REML="REML",REML0=bquote(REML[0]),Mom="MoM",MoM0=bquote(MoM[0]))) +
labs(y=expression(PVE))+
theme(text = element_text(size=20),legend.position = "bottom")
P
dev.off()
################################################ 012 mat: m ~ n1 ###############################################
out1 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample100_prop23_SNRy0.3_n800_4000.txt",header = T)
out2 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample300_prop23_SNRy0.3_n800_4000.txt",header = T)
out3 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample500_prop23_SNRy0.3_n800_4000.txt",header = T)
out4 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample1000_prop23_SNRy0.3_n800_4000.txt",header = T)
out5 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample100_prop23_SNRy0.3_n1000_4000.txt",header = T)
out6 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample300_prop23_SNRy0.3_n1000_4000.txt",header = T)
out7 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample500_prop23_SNRy0.3_n1000_4000.txt",header = T)
out8 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample1000_prop23_SNRy0.3_n1000_4000.txt",header = T)
out9 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample100_prop23_SNRy0.3_n2000_4000.txt",header = T)
out10 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample300_prop23_SNRy0.3_n2000_4000.txt",header = T)
out11 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample500_prop23_SNRy0.3_n2000_4000.txt",header = T)
out12 <- read.table("/Users/cmx/Desktop/Predixcan/mHeritability/simulation_results/MoM_subsample1000_prop23_SNRy0.3_n2000_4000.txt",header = T)
subsample <- c(replicate(3,c(100,300,500,1000)))
n1 <- c(sapply(c(800,1000,2000),rep,4))
n_rep <- nrow(out1)
out_all <- data.frame()
for(i in 1:12){
out <- get(paste("out",i,sep=""))
out <- melt(out[,1:2],id.vars =NULL)
out <- cbind(out,method=sapply(substring(out[,1],6),switch,REML="REML",REML0="REML0",MoM="MoM",MoM0="MoM0",ss="SummaryStat"),
component=ifelse(substr(out[,1],4,4)=="g","gene","SNP"),subsample=subsample[i],n1=n1[i])
out <- out[,-1]
out_all <- rbind(out_all,out)
}
names(out_all) <- c("PVE","method","component","subsample","n1")
out_all$method <- factor(out_all$method,levels = c("REML0","MoM0","REML","MoM","SummaryStat"))
levels(out_all$component) <- c("GREX","Alternative")
pdf("/Users/cmx/Desktop/Predixcan/iGREX_paper/Supp_geno_m_n1.pdf",width=13,height=10)
P <- ggplot(data=out_all,aes(x=component,y=PVE)) +
geom_boxplot() +
geom_hline(yintercept = 0.2,linetype="dashed",color="blue") +
geom_hline(yintercept = 0.3,linetype="dashed",color="red") +
facet_grid(n1~subsample, labeller = label_bquote(cols="m"*"="*.(subsample),rows="n"[r]*"="*.(n1))) +
scale_fill_brewer(palette="RdBu") +
labs(y=expression(PVE))+
theme(text = element_text(size=20),legend.position = "bottom")
P
dev.off()