-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathScript_04_Sensitivity_analysis_limma.Rmd
156 lines (127 loc) · 6.15 KB
/
Script_04_Sensitivity_analysis_limma.Rmd
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
---
title: "Script_04_Sensitivity_analysis_limma"
author: "Yunfeng Liu"
date: "1/16/2023"
output: html_document
---
# Script information
To see if dglm results are robust and convincing, we used another R package called "limma" to fit the linear model and compared the effect size of aDMCs from limma and dglm.
```{r}
#set library path.
.libPaths("~/researchdrive/yliu/Rlibs")
```
```{r}
#Load data.
load(file ="./Output/Output_02/Output_02_sex-separate_analysis.RData")
load(file ="./Output/Output_02/Output_02_aDMCs_effectsize_BIOS.RData")
```
```{r}
#Load all necessary libraries.
library(limma)
library(bacon)
library(tidyverse)
```
```{r}
#### Run actual models with limma
## Xfemales
fit_Xfemales <- lmFit(RIN.mvalues_Xfemales, design_sv_Xfemales)
se_Xfemales<-fit_Xfemales$stdev.unscaled*fit_Xfemales$sigma
tstat_Xfemales <- fit_Xfemales$coef/se_Xfemales
pval_Xfemales<-2*pt(-abs(tstat_Xfemales),fit_Xfemales$df.residual)
table(pval_Xfemales[,2]<0.05) # 4423
padj_Xfemales <- p.adjust(pval_Xfemales[,2], method="bonf")
table(padj_Xfemales<0.05) # 1109
# Extract significant statistics and save them as a dataframe
dat_Xfemales<- as.data.frame(cbind(fit_Xfemales$coefficients[,2],se_Xfemales[,2],tstat_Xfemales[,2],pval_Xfemales[,2],padj_Xfemales))
colnames(dat_Xfemales)<-c("effectsize","standard.error","t-satistics","p.value","p.adj")
## Xmales
fit_Xmales <- lmFit(RIN.mvalues_Xmales, design_sv_Xmales)
se_Xmales<-fit_Xmales$stdev.unscaled*fit_Xmales$sigma
tstat_Xmales <- fit_Xmales$coef/se_Xmales
pval_Xmales<-2*pt(-abs(tstat_Xmales),fit_Xmales$df.residual)
table(pval_Xmales[,2]<0.05) # 4083
padj_Xmales <- p.adjust(pval_Xmales[,2], method="bonf")
table(padj_Xmales<0.05) # 2245
# Extract significant statistics and save them as a dataframe
dat_Xmales<- as.data.frame(cbind(fit_Xmales$coefficients[,2],se_Xmales[,2],tstat_Xmales[,2],pval_Xmales[,2],padj_Xmales))
colnames(dat_Xmales)<-c("effectsize","standard.error","t-satistics","p.value","p.adj")
```
```{r}
#### Run bacon
## Correct bias and inflation for Xfemales
set.seed(1)
bc_Xfemales<- bacon(teststatistics=dat_Xfemales[,3])
# Calculate bias and inflation
bias(bc_Xfemales) # 0.08
inflation(bc_Xfemales) # 2.3
#Extract the BACON-adjusted t-statistics and p-values
tstats_Xfemales<-tstat(bc_Xfemales)
pvals_Xfemales <-pval(bc_Xfemales)
table(pvals_Xfemales<0.05) # 1104
padjs_Xfemales <- as.matrix(p.adjust(pvals_Xfemales, method="bonf"))
table(padjs_Xfemales<0.05) # 78
# Extract BACON-adjusted effectsize.
set.seed(1)
bc_Xfemales_es<-bacon(NULL,dat_Xfemales[,1],dat_Xfemales[,2])
es_Xfemales<-es(bc_Xfemales_es)
# Summary Bacon statistics
dat_Xfemales_bacon<-as.data.frame(cbind(es_Xfemales,tstats_Xfemales,pvals_Xfemales,padjs_Xfemales))
colnames(dat_Xfemales_bacon)<-c("bacon_effectsize","bacon_tstatistics","bacon_pval","bacon_p.adj")
rownames(dat_Xfemales_bacon)<-make.names(rownames(dat_Xfemales),unique = TRUE)
## Correct bias and inflation for Xmales
set.seed(1)
bc_Xmales<- bacon(teststatistics =dat_Xmales[,3])
# Calculate bias and inflation
bias(bc_Xmales) # -0.06
inflation(bc_Xmales) # 1.2
#Extract the BACON-adjusted t-statistics and p-values.
tstats_Xmales<-tstat(bc_Xmales)
pvals_Xmales <- pval(bc_Xmales)
table(pvals_Xmales<0.05) # 3569
padjs_Xmales <- as.matrix(p.adjust(pvals_Xmales, method="bonf"))
table(padjs_Xmales<0.05) # 1834
# Extract BACON-adjusted effectsize
set.seed(1)
bc_Xmales_es<-bacon(NULL,dat_Xmales[,1],dat_Xmales[,2])
es_Xmales<-es(bc_Xmales_es)
# Summary Bacon statistics
dat_Xmales_bacon<-as.data.frame(cbind(es_Xmales,tstats_Xmales,pvals_Xmales,padjs_Xmales))
colnames(dat_Xmales_bacon)<-c("bacon_effectsize","bacon_tstatistics","bacon_pval","bacon_p.adj")
rownames(dat_Xmales_bacon)<-make.names(rownames(dat_Xmales),unique = TRUE)
save(dat_Xfemales_bacon,dat_Xmales_bacon,file = "Output_04_aDMCs_limma_results.RData")
```
```{r}
## Compare aDMCs effectsize between dglm and limma in females
limma_aDMCs_females<-dat_Xfemales_bacon[rownames(df_aDMCs_effectsize_BIOS),]
df_effectsize_aDMCs_females<-data.frame(dglm=df_aDMCs_effectsize_BIOS[,1],limma=limma_aDMCs_females[,1])
rownames(df_effectsize_aDMCs_females)<-make.names(rownames(df_aDMCs_effectsize_BIOS),unique = TRUE)
df_effectsize_aDMCs_females$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
## Compare aDMCs effectsize between dglm and limma in males
limma_aDMCs_males<-dat_Xmales_bacon[rownames(df_aDMCs_effectsize_BIOS),]
df_effectsize_aDMCs_males<-data.frame(dglm=df_aDMCs_effectsize_BIOS[,2],limma=limma_aDMCs_males[,1])
rownames(df_effectsize_aDMCs_males)<-make.names(rownames(df_aDMCs_effectsize_BIOS),unique = TRUE)
df_effectsize_aDMCs_males$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
save(df_effectsize_aDMCs_females,df_effectsize_aDMCs_males,file="Output_04_Sensitivity analysis.RData")
# Scatter plot of aDMCs effect size between DGLM and limma
df_effectsize_aDMCs_females$sex<-"Females"
df_effectsize_aDMCs_males$sex<-"Males"
df_effectsize_aDMCs<-rbind(df_effectsize_aDMCs_females,df_effectsize_aDMCs_males)
df_effectsize_aDMCs$aDMCs<-factor(df_effectsize_aDMCs$aDMCs,levels = unique(df_effectsize_aDMCs$aDMCs))
df_effectsize_aDMCs$sex<-factor(df_effectsize_aDMCs$sex,levels = c("Males","Females"))
tiff("Figure_S2_Quality check of aDMCs.tiff", units="in", width=14, height=10, res=600,compression = 'lzw')
ggplot(df_effectsize_aDMCs) +
aes(x = dglm,y = limma,color=aDMCs) +
theme_grey(base_size = 18) +
facet_grid(~sex,switch ="y",scales = "fixed") +
geom_point(size = 3, alpha = 0.5) +
geom_vline(aes(xintercept=0)) +
geom_hline(aes(yintercept=0)) +
geom_abline(slope=1, intercept = 0)+
scale_x_continuous(limits=c(-0.04, 0.04))+
scale_y_continuous(limits=c(-0.04, 0.04))+
scale_color_manual(values = c("#F08080", "#56B4E9", "#32CD32"))+
theme_bw()+
theme(panel.spacing.x = unit(.5, "cm"),legend.position = "bottom",legend.text =element_text(size = 30),legend.title = element_text(size = 30),axis.title.x = element_text(size = 30),axis.title.y = element_text(size = 30),axis.text.x = element_text(size = 30,angle = 45,hjust = 1),axis.text.y = element_text(size = 30),strip.text.x = element_text(size = 23))+
labs(x="Effect size in dglm",y="Effect size in limma")
dev.off()
```