Figure 2. Weight effect of lung P value
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
library(reshape)
library(knitr)
merged.mouse.eQTL.min.variance2<-read.table(file="2016-02-22-merged.mouse.eQTL.min.variance.txt", header=T)
merged.mouse.eQTL.min.variance2$abs_liver.beta<-abs(merged.mouse.eQTL.min.variance2$liver.beta)
merged.mouse.eQTL.min.variance2$abs_lung.beta<-abs(merged.mouse.eQTL.min.variance2$lung.beta)
merged.mouse.eQTL.min.variance2$neg_log_lung_pvalue<--log(merged.mouse.eQTL.min.variance2$lung_pvalue)
merged.mouse.eQTL<-merged.mouse.eQTL.min.variance2
markers<-merged.mouse.eQTL[, 1]
# Yg=Ag + Bg*Xsnp+V
betas.hat<-abs(merged.mouse.eQTL[,2])
se<-(merged.mouse.eQTL[,3])
liver_residual_variance<-merged.mouse.eQTL[,4]
Z<-as.matrix(merged.mouse.eQTL[,10])
Z<-replace(Z,is.na(Z),0)
Z<-data.frame(1,Z)
Z<-as.matrix(Z)
rowLength<-length(markers)
lmsummary<-summary(lm(abs_liver.beta~-1+Z, data=merged.mouse.eQTL))
tau<-lmsummary$sigma**2
tau
## [1] 0.009699749
gamma<-as.matrix(lmsummary$coefficients[,1])
Z_transpose<-t(Z)
identity<-diag(nrow=rowLength)
betas.hat<-as.matrix(betas.hat)
V <- matrix(0, rowLength, rowLength)
diag(V) <-liver_residual_variance
Tau<- diag(tau, rowLength, rowLength)
s<-V + Tau
diag.inverse <- function(x){diag(1/diag(x), nrow(x), ncol(x))}
diag.multi <- function(x,y){diag(diag(x)*diag(y), nrow(x), ncol(x))}
S <-diag.inverse(s)
omega<-diag.multi(S, V)
omega.diag<-diag(omega )
# betas.thea<- S %*% Z %*% gamma + (identity-S) %*% betas.hat
betas.tieda<- omega %*% Z %*% gamma + (identity-omega) %*% betas.hat
# crbetas.tieda<- cromega %*% Z %*% gamma + (identity-cromega) %*% betas.hat
regbeta <-Z %*% gamma
liver.mouse.eQTL.bayesian.4tau <- read.table(file="2016-03-15_liver.mouse.eQTL.bayesian.4tau.txt", head=TRUE)
liver.mouse.eQTL.bayesian.4tau$fzm <- -log10(exp(-liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue))
range(liver.mouse.eQTL.bayesian.4tau$fzm)
## [1] 8.678947e-16 1.727677e+01
liver.mouse.eQTL.bayesian.4tau$ratio_fzm <- max(liver.mouse.eQTL.bayesian.4tau$fzm)/liver.mouse.eQTL.bayesian.4tau$fzm
range(liver.mouse.eQTL.bayesian.4tau$ratio_fzm)
## [1] 1.000000e+00 1.990653e+16
rho.optimization <- matrix(0, nrow=nrow(liver.mouse.eQTL.bayesian.4tau), ncol=7)
colnames(rho.optimization)<-c("rho","tmm","tau", "omega","beta_tieda", "n.betas.tieda.se","p.below.0" )
nomega.diag<-diag(nomega )
rho <- seq(1,51, by=10)*tau
rho
result <- NULL
for (i in 1:length(rho)) {
rho.optimization[ ,1] <- rho[i]
rho.optimization[ ,2] <- (rho[i]/tau)^liver.mouse.eQTL.bayesian.4tau$ratio_fzm
rho.optimization[ ,3] <-tau*((rho[i]/tau)^liver.mouse.eQTL.bayesian.4tau$ratio_fzm)
nTau<- diag(rho.optimization[i,3], rowLength, rowLength)
ns<-V + nTau
nS <- diag.inverse(ns)
nomega<-diag.multi(nS, V)
rho.optimization[ ,4] <- diag(nomega )
rho.optimization[ ,5] <-nomega %*% Z %*% gamma + (identity-nomega) %*% betas.hat
nTau_invert<-diag.inverse(nTau)
V_invert<-diag.inverse(V)
nPS_invert<-nTau_invert+ diag.multi(diag.multi(V_invert, Z), Z_transpose)
nPS<-diag.inverse(nPS_invert)
nps<-diag(nPS)
nps.long <- melt(nps)
rho.optimization[ ,6] <-(nps.long$value)^0.5
rho.optimization[ ,7] <- pnorm(0, rho.optimization[ ,5], rho.optimization[ ,6])
result <- rbind(result,rho.optimization)
}
result$p.below.0 <- pnorm(0, rho.optimization[ ,5], rho.optimization[ ,6])
write.table(result, file="2016-03-22_liver.mouse.eQTL.bayesian.result.txt",col.names=TRUE,row.names=FALSE,quote=FALSE)
liver.mouse.eQTL.bayesian.result <- read.table(file="2016-03-22_liver.mouse.eQTL.bayesian.result.txt", header=T)
knitr::kable(head(liver.mouse.eQTL.bayesian.result))
| 0.0096997 |
1 |
0.0096997 |
0.6609725 |
0.0335649 |
0.0800704 |
0.3375376 |
| 0.0096997 |
1 |
0.0096997 |
0.2791983 |
0.0236274 |
0.0983765 |
0.4050981 |
| 0.0096997 |
1 |
0.0096997 |
0.3255088 |
0.0140638 |
0.0561903 |
0.4011821 |
| 0.0096997 |
1 |
0.0096997 |
0.1586587 |
0.0416051 |
0.0982601 |
0.3359952 |
| 0.0096997 |
1 |
0.0096997 |
0.8895190 |
0.0401108 |
0.0928876 |
0.3329360 |
| 0.0096997 |
1 |
0.0096997 |
0.3335208 |
0.0306453 |
0.0984015 |
0.3777363 |
knitr::kable(tail(liver.mouse.eQTL.bayesian.result))
| 66049 |
0.4946872 |
1.180091 |
0.0114466 |
0.6634209 |
0.0473521 |
0.1046325 |
0.3254342 |
| 66050 |
0.4946872 |
1.472411 |
0.0142820 |
0.3777155 |
0.0383391 |
0.0643198 |
0.2755646 |
| 66051 |
0.4946872 |
1.096974 |
0.0106404 |
0.3395461 |
0.0197970 |
0.1045669 |
0.4249194 |
| 66052 |
0.4946872 |
1.139850 |
0.0110563 |
0.5864554 |
0.0408051 |
0.0801457 |
0.3053284 |
| 66053 |
0.4946872 |
1.306410 |
0.0126718 |
0.4981594 |
0.0320585 |
0.1046096 |
0.3796280 |
| 66054 |
0.4946872 |
1.296591 |
0.0125766 |
0.6772621 |
0.0399182 |
0.0861273 |
0.3215104 |
result.df <-liver.mouse.eQTL.bayesian.result
result.df$rho.class <- factor(result.df$rho/tau)
a <-liver.mouse.eQTL.bayesian.4tau[, 1:2]
a <-rbind(a, a, a, a, a,a)
dim(a)
## [1] 66054 2
new.result.df<-cbind(a, result.df)
knitr::kable(head(liver.mouse.eQTL.bayesian.4tau)) #[ , c(1:3, 8, 18, 4, 5, 10)]
| ENSMUSG00000000001 |
0.0159861 |
0.0256247 |
0.0335649 |
0.0984852 |
0.0189108 |
0.5377718 |
0.0326900 |
2.0960974 |
0.3666226 |
0.0328346 |
0.0277012 |
0.0161111 |
0.0425815 |
0.3666226 |
0.0425815 |
0.0425815 |
0.6609725 |
0.9103235 |
18.97871 |
| ENSMUSG00000000003 |
0.0168807 |
0.0126533 |
0.0236274 |
0.0984766 |
0.0037571 |
0.1929225 |
0.0295496 |
2.7653413 |
0.4051928 |
0.0230580 |
0.0201492 |
0.0169033 |
0.0410452 |
0.4051928 |
0.0410452 |
0.0410452 |
0.2791983 |
1.2009724 |
14.38565 |
| ENSMUSG00000000037 |
0.0038651 |
0.0136293 |
0.0140638 |
0.0984779 |
0.0046811 |
0.7788138 |
0.0175948 |
1.1814994 |
0.4432194 |
0.0132541 |
0.0089752 |
0.0039017 |
0.0351967 |
0.4432194 |
0.0351967 |
0.0351967 |
0.3255088 |
0.5131187 |
33.67012 |
| ENSMUSG00000000049 |
0.0422333 |
0.0078085 |
0.0416051 |
0.0984645 |
0.0018292 |
0.0000091 |
0.0238839 |
0.2321777 |
0.3363160 |
0.0416661 |
0.0419531 |
0.0422315 |
0.0382735 |
0.3363160 |
0.0382735 |
0.0382735 |
0.1586587 |
0.1008335 |
171.33960 |
| ENSMUSG00000000056 |
0.0636041 |
0.0514812 |
0.0401108 |
0.0984868 |
0.0780959 |
0.2269175 |
0.0216750 |
1.5218109 |
0.3419046 |
0.0404380 |
0.0434054 |
0.0630989 |
0.0371928 |
0.3419046 |
0.0371928 |
0.0371928 |
0.8895190 |
0.6609141 |
26.14072 |
| ENSMUSG00000000058 |
0.0256111 |
0.0129823 |
0.0306453 |
0.0984790 |
0.0048540 |
0.0584692 |
0.0288547 |
1.2292565 |
0.3778295 |
0.0302500 |
0.0281486 |
0.0256294 |
0.0407052 |
0.3778295 |
0.0407052 |
0.0407052 |
0.3335208 |
0.5338593 |
32.36202 |
new.result.df2 <- new.result.df[, c(1,2,7,5,6,10)]
knitr::kable(head(new.result.df2))
| ENSMUSG00000000001 |
0.0159861 |
0.0335649 |
0.0096997 |
0.6609725 |
1 |
| ENSMUSG00000000003 |
0.0168807 |
0.0236274 |
0.0096997 |
0.2791983 |
1 |
| ENSMUSG00000000037 |
0.0038651 |
0.0140638 |
0.0096997 |
0.3255088 |
1 |
| ENSMUSG00000000049 |
0.0422333 |
0.0416051 |
0.0096997 |
0.1586587 |
1 |
| ENSMUSG00000000056 |
0.0636041 |
0.0401108 |
0.0096997 |
0.8895190 |
1 |
| ENSMUSG00000000058 |
0.0256111 |
0.0306453 |
0.0096997 |
0.3335208 |
1 |
by(new.result.df2[, c(1:5)], new.result.df2[, "rho.class"], head)
## new.result.df2[, "rho.class"]: 1
## ensembl_id betas.hat beta_tieda tau omega
## 1 ENSMUSG00000000001 0.015986111 0.03356494 0.009699749 0.6609725
## 2 ENSMUSG00000000003 0.016880682 0.02362737 0.009699749 0.2791983
## 3 ENSMUSG00000000037 0.003865079 0.01406380 0.009699749 0.3255088
## 4 ENSMUSG00000000049 0.042233333 0.04160507 0.009699749 0.1586587
## 5 ENSMUSG00000000056 0.063604072 0.04011077 0.009699749 0.8895190
## 6 ENSMUSG00000000058 0.025611111 0.03064532 0.009699749 0.3335208
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 11
## ensembl_id betas.hat beta_tieda tau omega
## 11010 ENSMUSG00000000001 0.015986111 0.03254657 0.011006067 0.6226813
## 11011 ENSMUSG00000000003 0.016880682 0.02284729 0.011459129 0.2469161
## 11012 ENSMUSG00000000037 0.003865079 0.01295208 0.010415732 0.2900264
## 11013 ENSMUSG00000000049 0.042233333 0.04168825 0.009836452 0.1376522
## 11014 ENSMUSG00000000056 0.063604072 0.04057231 0.010631596 0.8720437
## 11015 ENSMUSG00000000058 0.025611111 0.03010237 0.010445758 0.2975504
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 21
## ensembl_id betas.hat beta_tieda tau omega
## 22019 ENSMUSG00000000001 0.015986111 0.03301847 0.011387518 0.6404250
## 22020 ENSMUSG00000000003 0.016880682 0.02319654 0.011985962 0.2613693
## 22021 ENSMUSG00000000037 0.003865079 0.01345187 0.010617697 0.3059778
## 22022 ENSMUSG00000000049 0.042233333 0.04165140 0.009873644 0.1469577
## 22023 ENSMUSG00000000056 0.063604072 0.04035386 0.010897863 0.8803149
## 22024 ENSMUSG00000000058 0.025611111 0.03034663 0.010656575 0.3137328
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 31
## ensembl_id betas.hat beta_tieda tau omega
## 33028 ENSMUSG00000000001 0.015986111 0.03344512 0.011623617 0.6564670
## 33029 ENSMUSG00000000003 0.016880682 0.02353034 0.012314892 0.2751829
## 33030 ENSMUSG00000000037 0.003865079 0.01392642 0.010741226 0.3211240
## 33031 ENSMUSG00000000049 0.042233333 0.04161559 0.009896113 0.1560016
## 33032 ENSMUSG00000000056 0.063604072 0.04016319 0.011061443 0.8875340
## 33033 ENSMUSG00000000058 0.025611111 0.03057830 0.010785597 0.3290808
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 41
## ensembl_id betas.hat beta_tieda tau omega
## 44037 ENSMUSG00000000001 0.015986111 0.03269998 0.011796117 0.6284494
## 44038 ENSMUSG00000000003 0.016880682 0.02295863 0.012556573 0.2515237
## 44039 ENSMUSG00000000037 0.003865079 0.01311178 0.010830788 0.2951232
## 44040 ENSMUSG00000000049 0.042233333 0.04167657 0.009912274 0.1406015
## 44041 ENSMUSG00000000056 0.063604072 0.04050040 0.011180385 0.8747665
## 44042 ENSMUSG00000000058 0.025611111 0.03018045 0.010879181 0.3027231
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 51
## ensembl_id betas.hat beta_tieda tau omega
## 55046 ENSMUSG00000000001 0.015986111 0.03282735 0.011932555 0.6332387
## 55047 ENSMUSG00000000003 0.016880682 0.02305266 0.012748529 0.2554151
## 55048 ENSMUSG00000000037 0.003865079 0.01324638 0.010901223 0.2994194
## 55049 ENSMUSG00000000049 0.042233333 0.04166666 0.009924908 0.1431049
## 55050 ENSMUSG00000000056 0.063604072 0.04044135 0.011274123 0.8770022
## 55051 ENSMUSG00000000058 0.025611111 0.03024624 0.010952800 0.3070816
by(new.result.df2[, c(2:5)], new.result.df2[, "rho.class"], summary)
## new.result.df2[, "rho.class"]: 1
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004582 Min. :0.0097 Min. :0.09131
## 1st Qu.:0.01186 1st Qu.:0.024829 1st Qu.:0.0097 1st Qu.:0.36400
## Median :0.02430 Median :0.035162 Median :0.0097 Median :0.52523
## Mean :0.05436 Mean :0.045850 Mean :0.0097 Mean :0.54385
## 3rd Qu.:0.05360 3rd Qu.:0.051346 3rd Qu.:0.0097 3rd Qu.:0.70825
## Max. :2.76340 Max. :1.223416 Max. :0.0097 Max. :0.99758
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 11
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004028 Min. :0.00970 Min. :0.07839
## 1st Qu.:0.01186 1st Qu.:0.023889 1st Qu.:0.01018 1st Qu.:0.32635
## Median :0.02430 Median :0.034613 Median :0.01074 Median :0.48359
## Mean :0.05436 Mean :0.045565 Mean :0.01196 Mean :0.50969
## 3rd Qu.:0.05360 3rd Qu.:0.051252 3rd Qu.:0.01201 3rd Qu.:0.67266
## Max. :2.76340 Max. :1.266254 Max. :0.10670 Max. :0.99715
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 21
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004274 Min. :0.00970 Min. :0.08408
## 1st Qu.:0.01186 1st Qu.:0.024329 1st Qu.:0.01032 1st Qu.:0.34334
## Median :0.02430 Median :0.034861 Median :0.01104 Median :0.50265
## Mean :0.05436 Mean :0.045688 Mean :0.01283 Mean :0.52532
## 3rd Qu.:0.05360 3rd Qu.:0.051270 3rd Qu.:0.01272 3rd Qu.:0.68922
## Max. :2.76340 Max. :1.246068 Max. :0.20369 Max. :0.99736
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 31
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004512 Min. :0.00970 Min. :0.08966
## 1st Qu.:0.01186 1st Qu.:0.024718 1st Qu.:0.01040 1st Qu.:0.35938
## Median :0.02430 Median :0.035112 Median :0.01123 Median :0.52023
## Mean :0.05436 Mean :0.045813 Mean :0.01344 Mean :0.53974
## 3rd Qu.:0.05360 3rd Qu.:0.051296 3rd Qu.:0.01317 3rd Qu.:0.70410
## Max. :2.76340 Max. :1.228319 Max. :0.30069 Max. :0.99753
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 41
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004106 Min. :0.00970 Min. :0.08019
## 1st Qu.:0.01186 1st Qu.:0.024027 1st Qu.:0.01046 1st Qu.:0.33179
## Median :0.02430 Median :0.034684 Median :0.01136 Median :0.48974
## Mean :0.05436 Mean :0.045603 Mean :0.01393 Mean :0.51474
## 3rd Qu.:0.05360 3rd Qu.:0.051232 3rd Qu.:0.01350 3rd Qu.:0.67806
## Max. :2.76340 Max. :1.259626 Max. :0.39769 Max. :0.99722
## --------------------------------------------------------
## new.result.df2[, "rho.class"]: 51
## betas.hat beta_tieda tau omega
## Min. :0.00000 Min. :0.004172 Min. :0.00970 Min. :0.08172
## 1st Qu.:0.01186 1st Qu.:0.024136 1st Qu.:0.01050 1st Qu.:0.33637
## Median :0.02430 Median :0.034764 Median :0.01147 Median :0.49488
## Mean :0.05436 Mean :0.045636 Mean :0.01433 Mean :0.51895
## 3rd Qu.:0.05360 3rd Qu.:0.051227 3rd Qu.:0.01377 3rd Qu.:0.68253
## Max. :2.76340 Max. :1.254172 Max. :0.49469 Max. :0.99727
qplot(omega, 1-p.below.0, group=rho.class, colour = rho.class, shape=rho.class, data = new.result.df)
qplot(omega, beta_tieda, group=rho.class, colour = rho.class, shape=rho.class, data = new.result.df)
qplot(betas.hat, beta_tieda, group=rho.class, colour = rho.class, shape=rho.class, data = new.result.df)
