/
pleiotropy_rmvmr.R
179 lines (138 loc) · 5.22 KB
/
pleiotropy_rmvmr.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
#' pleiotropy_rmvmr
#'
#' Generates Q-statistics quantifying the degree of heterogeneity in univariate Radial MR analyses applying a correction using the
#' output from [`ivw_rmvmr`]. The function returns two data frames. The first data frame includes the global Q-statistic for each exposure after applying
#' a correction, as well as a corresponding p-value. The second data frame contains the individual Q-statistic for each SNP in the corrected univariate
#' analyses, relative to the exposure given in column \code{exposure}.
#'
#' @param r_input A formatted data frame using the [`format_rmvmr`] function or an object of class `MRMVInput` from [`MendelianRandomization::mr_mvinput`]
#' @param rmvmr An object containing the output from the [`ivw_rmvmr`] function of class \code{IVW_RMVMR}.
#'
#' @return An object of class \code{"RMVMR_Q"} containing the following components:\describe{
#' \item{\code{gq}}{A data frame containing the global Q-statistic and p-value after applying a correction for each exposure}
#' \item{\code{qdat}}{A data frame containing the individual Q-statistic and p-value for each SNP after applying a correction for each exposure}
#' }
#'
#' @author Wes Spiller; Eleanor Sanderson; Jack Bowden.
#' @references Spiller, W., et al., Estimating and visualising multivariable Mendelian randomization analyses within a radial framework. Forthcoming.
#' @export
#' @examples
#' f.data <- format_rmvmr(
#' BXGs = rawdat_rmvmr[,c("ldl_beta","hdl_beta","tg_beta")],
#' BYG = rawdat_rmvmr$sbp_beta,
#' seBXGs = rawdat_rmvmr[,c("ldl_se","hdl_se","tg_se")],
#' seBYG = rawdat_rmvmr$sbp_se,
#' RSID = rawdat_rmvmr$snp)
#' rmvmr_output <- ivw_rmvmr(f.data, FALSE)
#' q_object <- pleiotropy_rmvmr(f.data, rmvmr_output)
#' q_object$gq
#' head(q_object$qdat)
pleiotropy_rmvmr <- function(r_input, rmvmr){
# convert MRMVInput object to mvmr_format
if ("MRMVInput" %in% class(r_input)) {
r_input <- mrmvinput_to_rmvmr_format(r_input)
}
# Perform check that r_input has been formatted using format_rmvmr function
if(!("rmvmr_format" %in%
class(r_input))) {
stop('The class of the data object must be "rmvmr_format", please resave the object with the output of format_rmvmr().')
}
# Extract MVMR estimates
rmvmr<-rmvmr$coef
#Define number of exposures included in MVMR model
exp.number<-length(names(r_input)[-c(1,2,3)])/2
#Define matrix for identifying IV1 satisfying variants using F>10.
f.vec<-matrix(0L, nrow = length(r_input[,1]), ncol = exp.number)
for(i in 1:exp.number){
f.vec[,i]<- r_input[,3+i]^2/r_input[,3 + exp.number + i]^2
for(j in 1:length(r_input[,1])){
if(f.vec[j,i] < 10){
f.vec[j,i]<-0
}else{
f.vec[j,i]<-1
}
}
}
#Define null variable for univariate MR data
Xlist<-NULL
#Obtain univariate MR data for each exposure
for(i in 1:exp.number){
Xsub<-r_input[f.vec[,i] == 1,]
Xrad.dat <- RadialMR::format_radial(Xsub[,3+i],Xsub[,2],Xsub[,3 + exp.number + i],Xsub[,3],Xsub[,1])
X.res <- RadialMR::ivw_radial(Xrad.dat,0.05/nrow(Xrad.dat),1,0.0001,F)
if(is.null(Xlist)){
Xlist<-X.res
}else{
Xlist<-append(Xlist,X.res)
}
}
#Create combined data frame of univariate values
p.dat<-NULL
for(i in 1:exp.number){
Xdat<-data.frame(Xlist[5 + ((i-1)* 13)])
Xdat$Group<-rep(i,nrow(Xdat))
names(Xdat)<-c("SNP","Wj","BetaWj","Qj","Qj_Chi","Outliers","Group")
if(is.null(p.dat)){
p.dat<-Xdat
}else{
p.dat<-rbind(p.dat,Xdat)
}
}
p.dat[,7]<-as.factor(p.dat[,7])
for(i in 1:exp.number){
levels(p.dat[,7])[i] <- paste0("Exposure_",i,collapse="")
}
##AL
Ratios<-NULL
for(i in 1:exp.number){
tdat<-r_input[r_input$SNP %in% p.dat[p.dat$Group==levels(p.dat$Group)[i],]$SNP,]
Ratio_temp<-tdat[,2] / tdat[,(3+i)]
for(j in 1:exp.number){
if(i == j){
Ratio_temp<-Ratio_temp
}else{
Ratio_temp<- Ratio_temp - ((tdat[,(3+j)]*rmvmr[j,1]) / tdat[,(3+i)])
}
}
if(is.null(Ratios)){
Ratios<-Ratio_temp
}else{
Ratios<-c(Ratios,Ratio_temp)
}
}
#QJ calculations
p.dat$coratios<-Ratios
Qjvec<-NULL
for(i in 1:exp.number){
tdat<-p.dat[p.dat$Group==levels(p.dat$Group)[i],]
Qj<-tdat$Wj * (tdat$coratios - rmvmr[i,1])^2
if(is.null(Qjvec)){
Qjvec<-Qj
}else{
Qjvec<-c(Qjvec,Qj)
}
}
p.dat$Qjcor<-Qjvec
#Define matrix for recording total Q statistics.
Qj_out<-matrix(0L, nrow = exp.number, ncol = 2)
indqj<-rep(0,length(p.dat[,1]))
for(i in 1:exp.number){
Qj_out[i,1]<-sum(p.dat[p.dat$Group==levels(p.dat$Group)[i],]$Qjcor)
Qj_out[i,2]<-stats::pchisq(Qj_out[i,1],nrow(p.dat[p.dat$Group==levels(p.dat$Group)[i],])-exp.number,lower.tail = FALSE)
}
TotalQs<-data.frame(Qj_out)
names(TotalQs)<- c("q_statistic","p_value")
row.names(TotalQs)<-levels(p.dat$Group)
for(i in 1:length(p.dat[,1])){
indqj[i]<- stats::pchisq(p.dat$Qjcor[i],1,lower.tail = FALSE)
}
p.dat$corQjchi<-indqj
out_data<-p.dat[,c(1,2,8,9,10,7)]
names(out_data)<-c("snp","wj","corrected_beta","qj","qj_p","ref_exposure")
multi_return <- function() {
Out_list <- list("gq" = TotalQs, "qdat" = out_data)
class(Out_list)<-"RMVMR_Q"
return(Out_list)
}
OUT<-multi_return()
}