library(GWASbrewer)
library(tidyr)
library(dplyr)
library(ggplot2)
library(mrclust)
Background
data("ld_mat_list")
data("AF")
<- function(n, ld_mat_list) {
make_ld_mat
<- tibble(ind = 1:nrow(ld_mat_list[[1]]), lind = 1:nrow(ld_mat_list[[1]]), i = 1)
ind for(i in 2:length(ld_mat_list)) {
<- bind_rows(ind, tibble(ind = max(ind$ind) + 1:nrow(ld_mat_list[[i]]), lind = 1:nrow(ld_mat_list[[i]]), i = i))
ind
}<- ind[ind$ind <= n,]
ind <- matrix(0, n, n)
mat for(i1 in 1:length(unique(ind$i))) {
<- subset(ind, i == i1)
d $ind[1]:max(d$ind), d$ind[1]:max(d$ind)] <- as.matrix(ld_mat_list[[i1]][d$lind[1]:max(d$lind), d$lind[1]:max(d$lind)])
mat[d
}return(mat)
}<- make_ld_mat(1000, ld_mat_list)
ldmat dim(ldmat)
set.seed(1)
<- matrix(c(0, sqrt(0.25), 0, sqrt(0.15),
G 0, 0, 0, sqrt(0.1),
sqrt(0.2), 0, 0, -sqrt(0.3),
0, 0, 0, 0), nrow = 4, byrow = TRUE)
colnames(G) <- row.names(G) <- c("X", "Y", "Z", "W")
G
<- sim_mv(G = G,
sim_dat1_LD J = 1000,
N = 50000,
h2 = c(0.3, 0.3, 0.5, 0.4),
pi = 1/1000,
R_LD = list(ldmat=Matrix(ldmat)),
af = AF)
<- inner_join(
b $beta_hat %>% as_tibble() %>% rename(X=V1, Y=V2, Z=V3, W=V4) %>% mutate(pos=1:n()) %>%
sim_dat1_LDpivot_longer(c("X", "Y", "Z", "W")) %>% rename(bhat=value),
$se_beta_hat %>% as_tibble() %>% rename(X=V1, Y=V2, Z=V3, W=V4) %>% mutate(pos=1:n()) %>%
sim_dat1_LDpivot_longer(c("X", "Y", "Z", "W")) %>% rename(se=value)
)
b
<- inner_join(
bm $beta_joint %>% as_tibble() %>% mutate(pos=1:n()) %>%
sim_dat1_LDpivot_longer(c("X", "Y", "Z", "W")) %>% rename(bhat=value),
$se_beta_hat %>% as_tibble() %>% rename(X=V1, Y=V2, Z=V3, W=V4) %>% mutate(pos=1:n()) %>%
sim_dat1_LDpivot_longer(c("X", "Y", "Z", "W")) %>% rename(se=value)
)
bm
%>% ggplot(., aes(x=pos, y=bhat)) +
b geom_point() +
facet_grid(name ~ .)
%>% ggplot(., aes(x=pos, y=bhat)) +
bm geom_point() +
facet_grid(name ~ .)
<- b %>% select(pos, name, bhat) %>% pivot_wider(names_from=c(name), values_from=bhat)
bw
<- b %>% select(pos, name, se) %>% pivot_wider(names_from=c(name), values_from=se)
sew
pairs(bw[,-1])
<- function(BetaIV.in, seBetaIV.in, phi)
beta
{#Bandwidth rule - modified Silverman's rule proposed by Bickel (2002)
<- 0.9*(min(stats::sd(BetaIV.in), stats::mad(BetaIV.in)))/length(BetaIV.in)^(1/5)
s
#Standardised weights
<- seBetaIV.in^-2/sum(seBetaIV.in^-2)
weights
<- NULL
beta
for(cur_phi in phi)
{#Define the actual bandwidth
<- max(0.00000001, s*cur_phi)
h #Compute the smoothed empirical density function
<- stats::density(BetaIV.in, weights=weights, bw=h)
densityIV #Extract the point with the highest density as the point estimate
length(beta)+1] <- densityIV$x[densityIV$y==max(densityIV$y)]
beta[
}return(beta)
}
<- bw$Z / bw$W
biv abs(bw$Z) < 0.01 | abs(bw$Y) < 0.01] <- NA
biv[plot(biv)
<- bw$W / bw$Z
biv abs(bw$Z) < 0.01 | abs(bw$Y) < 0.01] <- NA
biv[plot(biv)
<- bw$W / bw$X
biv abs(bw$Z) < 0.01 | abs(bw$Y) < 0.01] <- NA
biv[plot(biv)
<- abs(bw$W) > 0.05 | abs(bw$X) > 0.05
ind pairs(bw[ind,-1])
<- bw[ind,]
bw <- sew[ind,]
sew = mr_clust_em(theta = bw$W/bw$X, theta_se = sew$W/abs(bw$X), bx = bw$X, by = bw$W, bxse = sew$X, byse = sew$Y, obs_names = bw$pos) res_em
head(res_em$results$best)
= res_em$plots$two_stage + ggplot2::xlim(0, max(abs(bx)+2*bxse)) + ggplot2::xlab("Genetic association with SBP") + ggplot2::ylab("Genetic association with CAD") + ggplot2::ggtitle("")
plot.sbp.best
plot.sbp.best
build graph
steiger dir x -> y
- get the causal variants across all traits
- get conditional on all variants
get_conditional
<- t(ldmat) %*% ldmat
xpx 1:10,1:10]
xpx[
<- solve(xpx)
xpxi 1:10,1:10]
xpxi[
<- princomp(xpx)
pc
names(pc)
$loadings[1:10,1:10]
pc
<- eigen(ldmat)
xe
class(xe)
2]][1:10,1:10]
xe[[class(xe)
<- tcrossprod(xe$vectors, tcrossprod(xe$vectors, diag(xe$values)))
x
1:10,1:10]
x[1:10,1:10]
ldmat[
= t(v)l
rho
<- xe$vectors %*% t(xe$vectors %*% diag(xe$values))
rho dim(rho)
1:10,1:10]
rho[
<- solve(xe$vectors) vi
sessionInfo()