library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
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
rho <-c(0.1, 1:10)*tau
rho
## [1] 0.0009699749 0.0096997493 0.0193994987 0.0290992480 0.0387989974
## [6] 0.0484987467 0.0581984961 0.0678982454 0.0775979948 0.0872977441
## [11] 0.0969974934
liver.mouse.eQTL.bayesian.4tau <- read.table(file="2016-03-10_liver.mouse.eQTL.bayesian.4tau.txt", header=T)
head(liver.mouse.eQTL.bayesian.4tau)
## ensembl_id betas.hat betas.hat.se betas.tieda betas.tieda.se
## 1 ENSMUSG00000000001 0.015986111 0.025624670 0.03356494 0.09848523
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02362737 0.09847664
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01406380 0.09847794
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.04160507 0.09846453
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.04011077 0.09848676
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.03064532 0.09847901
## liver_residual_variance liver_pvalue abs_lung.beta neg_log_lung_pvalue
## 1 0.018910763 5.377718e-01 0.03269000 2.0960974
## 2 0.003757142 1.929225e-01 0.02954960 2.7653413
## 3 0.004681090 7.788138e-01 0.01759477 1.1814994
## 4 0.001829162 9.084753e-06 0.02388393 0.2321777
## 5 0.078095912 2.269175e-01 0.02167500 1.5218109
## 6 0.004853968 5.846922e-02 0.02885470 1.2292565
## p.below.0 betas.tieda2m betas.tieda3rd betas.tieda4max regbeta
## 1 0.3666226 0.03283457 0.027701239 0.016111099 0.04258151
## 2 0.4051928 0.02305804 0.020149177 0.016903330 0.04104519
## 3 0.4432194 0.01325407 0.008975207 0.003901657 0.03519672
## 4 0.3363160 0.04166609 0.041953121 0.042231526 0.03827346
## 5 0.3419046 0.04043802 0.043405355 0.063098914 0.03719283
## 6 0.3778295 0.03025000 0.028148566 0.025629383 0.04070523
## posteriorprobability regbeta.1 regbeta.2 omega.diag rank.p.below.0
## 1 0.3666226 0.04258151 0.04258151 0.6609725 5928
## 2 0.4051928 0.04104519 0.04104519 0.2791983 8575
## 3 0.4432194 0.03519672 0.03519672 0.3255088 10600
## 4 0.3363160 0.03827346 0.03827346 0.1586587 4135
## 5 0.3419046 0.03719283 0.03719283 0.8895190 4439
## 6 0.3778295 0.04070523 0.04070523 0.3335208 6650
## rank.liver.pvalue fzm ratio_fzm tmm10 nTau
## 1 8485 -1.0896765 -0.07132899 0.8485374 0.008230600
## 2 5210 -0.7990276 -0.05230344 0.8865364 0.008599181
## 3 10100 -1.4868813 -0.09732957 0.7992275 0.007752307
## 4 697 -1.8991665 -0.12431730 0.7510740 0.007285229
## 5 5615 -1.3390859 -0.08765505 0.8172312 0.007926938
## 6 3182 -1.4661407 -0.09597192 0.8017299 0.007776579
## n.betas.tieda nps n.betas.tieda.se n.p.below.0
## 1 0.03451648 0.008230260 0.09072078 0.3517985
## 2 0.02422829 0.008597280 0.09272152 0.3969299
## 3 0.01566123 0.007750977 0.08803964 0.4294053
## 4 0.04143863 0.007282337 0.08533661 0.3136288
## 5 0.03962660 0.007926857 0.08903290 0.3281316
## 6 0.03141184 0.007775371 0.08817807 0.3608334
liver.mouse.eQTL.bayesian.4tau$fzm <- -log10(exp(-liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue))-2
liver.mouse.eQTL.bayesian.4tau$ratio_fzm <- liver.mouse.eQTL.bayesian.4tau$fzm/max(liver.mouse.eQTL.bayesian.4tau$fzm)
range(liver.mouse.eQTL.bayesian.4tau$fzm)
## [1] -2.00000 15.27677
range(liver.mouse.eQTL.bayesian.4tau$ratio_fzm)
## [1] -0.1309177 1.0000000
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" )
rho <- c(0.1, 1:10)*tau
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 <-qr.solve (ns)
nomega<-nS %*% V
rho.optimization[ ,4] <- diag(nomega )
rho.optimization[ ,5] <-nomega %*% Z %*% gamma + (identity-nomega) %*% betas.hat
nTau_invert<-qr.solve(nTau)
V_invert<-qr.solve(V)
nPS_invert<-nTau_invert+V_invert%*% Z %*% Z_transpose
nPS<-qr.solve(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-14_liver.mouse.eQTL.bayesian.result.txt",col.names=TRUE,row.names=FALSE,quote=FALSE)
liver.mouse.eQTL.bayesian.result <- read.table(file="2016-03-14_liver.mouse.eQTL.bayesian.result.txt", header=T)
head(liver.mouse.eQTL.bayesian.result)
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 1 0.0009699749 1.178498 0.01143114 0.6232557 0.03256185 0.1069142
## 2 0.0009699749 1.127985 0.01094117 0.2473711 0.02285828 0.1069049
## 3 0.0009699749 1.251208 0.01213641 0.2905302 0.01296787 0.1069063
## 4 0.0009699749 1.331427 0.01291451 0.1379427 0.04168710 0.1068918
## 5 0.0009699749 1.223644 0.01186904 0.8723164 0.04056511 0.1069159
## 6 0.0009699749 1.247303 0.01209853 0.2980618 0.03011009 0.1069075
## p.below.0
## 1 0.3803505
## 2 0.4153442
## 3 0.4517262
## 4 0.3482707
## 5 0.3521912
## 6 0.3891074
tail(liver.mouse.eQTL.bayesian.result)
## rho tmm tau omega beta_tieda
## 121094 0.09699749 0.8254897 0.008007043 0.7371938 0.04571824
## 121095 0.09699749 0.9558005 0.009271026 0.4634656 0.04020653
## 121096 0.09699749 0.7865097 0.007628946 0.4225144 0.02233009
## 121097 0.09699749 0.8067374 0.007825151 0.6686728 0.04100612
## 121098 0.09699749 0.8830012 0.008564890 0.5855247 0.03630458
## 121099 0.09699749 0.8786000 0.008522200 0.7491483 0.03951004
## n.betas.tieda.se p.below.0
## 121094 0.08772680 0.3011334
## 121095 0.08772353 0.3233568
## 121096 0.08772205 0.3995334
## 121097 0.08772620 0.3200951
## 121098 0.08772544 0.3394945
## 121099 0.08772686 0.3262198
result.df <-liver.mouse.eQTL.bayesian.result
result.df$rho.class <- factor(result.df$rho/tau)
new.result.df <- subset(result.df, rho.class == 0.1 | rho.class == 1 | rho.class == 2 |rho.class == 4 |rho.class == 6 |rho.class == 10)
new.result.df <- droplevels(new.result.df)
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, new.result.df)
head(liver.mouse.eQTL.bayesian.4tau[ , c(1:3, 8, 18, 4, 5, 10)])
## ensembl_id betas.hat betas.hat.se abs_lung.beta omega.diag
## 1 ENSMUSG00000000001 0.015986111 0.025624670 0.03269000 0.6609725
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02954960 0.2791983
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01759477 0.3255088
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.02388393 0.1586587
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.02167500 0.8895190
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.02885470 0.3335208
## betas.tieda betas.tieda.se p.below.0
## 1 0.03356494 0.09848523 0.3666226
## 2 0.02362737 0.09847664 0.4051928
## 3 0.01406380 0.09847794 0.4432194
## 4 0.04160507 0.09846453 0.3363160
## 5 0.04011077 0.09848676 0.3419046
## 6 0.03064532 0.09847901 0.3778295
head(new.result.df)
## ensembl_id betas.hat rho tmm tau
## 1 ENSMUSG00000000001 0.015986111 0.0009699749 1.178498 0.01143114
## 2 ENSMUSG00000000003 0.016880682 0.0009699749 1.127985 0.01094117
## 3 ENSMUSG00000000037 0.003865079 0.0009699749 1.251208 0.01213641
## 4 ENSMUSG00000000049 0.042233333 0.0009699749 1.331427 0.01291451
## 5 ENSMUSG00000000056 0.063604072 0.0009699749 1.223644 0.01186904
## 6 ENSMUSG00000000058 0.025611111 0.0009699749 1.247303 0.01209853
## omega beta_tieda n.betas.tieda.se p.below.0 rho.class
## 1 0.6232557 0.03256185 0.1069142 0.3803505 0.1
## 2 0.2473711 0.02285828 0.1069049 0.4153442 0.1
## 3 0.2905302 0.01296787 0.1069063 0.4517262 0.1
## 4 0.1379427 0.04168710 0.1068918 0.3482707 0.1
## 5 0.8723164 0.04056511 0.1069159 0.3521912 0.1
## 6 0.2980618 0.03011009 0.1069075 0.3891074 0.1
by(new.result.df[ , c(4, 6, 7, 9)], new.result.df[, "rho.class"], head)
## new.result.df[, "rho.class"]: 0.1
## tmm omega beta_tieda p.below.0
## 1 1.178498 0.6232557 0.03256185 0.3803505
## 2 1.127985 0.2473711 0.02285828 0.4153442
## 3 1.251208 0.2905302 0.01296787 0.4517262
## 4 1.331427 0.1379427 0.04168710 0.3482707
## 5 1.223644 0.8723164 0.04056511 0.3521912
## 6 1.247303 0.2980618 0.03011009 0.3891074
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 1
## tmm omega beta_tieda p.below.0
## 11010 1 0.6609725 0.03356494 0.3666226
## 11011 1 0.2791983 0.02362737 0.4051928
## 11012 1 0.3255088 0.01406380 0.4432194
## 11013 1 0.1586587 0.04160507 0.3363160
## 11014 1 0.8895190 0.04011077 0.3419046
## 11015 1 0.3335208 0.03064532 0.3778295
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 2
## tmm omega beta_tieda p.below.0
## 22019 0.9517608 0.6759223 0.03396254 0.3606657
## 22020 0.9643953 0.2929752 0.02396028 0.4006533
## 22021 0.9347616 0.3404913 0.01453323 0.4393405
## 22022 0.9174381 0.1678728 0.04156858 0.3311821
## 22023 0.9410511 0.8959769 0.03994021 0.3374430
## 22024 0.9356417 0.3486815 0.03087415 0.3728685
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 4
## tmm omega beta_tieda p.below.0
## 33028 0.9058487 0.6876481 0.03427439 0.3557599
## 33029 0.9300583 0.3042955 0.02423383 0.3968526
## 33030 0.8737793 0.3527316 0.01491674 0.4360619
## 33031 0.8416926 0.1755596 0.04153814 0.3269738
## 33032 0.8855771 0.9009079 0.03980997 0.3337662
## 33033 0.8754254 0.3610551 0.03106092 0.3687507
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 6
## tmm omega beta_tieda p.below.0
## 44037 0.8800255 0.7049307 0.03473403 0.3481167
## 44038 0.9105421 0.3218717 0.02465855 0.3908187
## 44039 0.8399682 0.3716111 0.01550827 0.4308048
## 44040 0.8003175 0.1877062 0.04149004 0.3204483
## 44041 0.8546555 0.9079699 0.03962346 0.3280297
## 44042 0.8420140 0.3801184 0.03134867 0.3622763
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 10
## tmm omega beta_tieda p.below.0
## 55046 0.8485374 0.7107426 0.03488860 0.3454269
## 55047 0.8865364 0.3280365 0.02480752 0.3886625
## 55048 0.7992275 0.3781973 0.01571462 0.4289119
## 55049 0.7510740 0.1920290 0.04147292 0.3181597
## 55050 0.8172312 0.9102915 0.03956214 0.3260079
## 55051 0.8017299 0.3867626 0.03144895 0.3599806
by(new.result.df[ , c(4, 6, 7, 9)], new.result.df[, "rho.class"], summary)
## new.result.df[, "rho.class"]: 0.1
## tmm omega beta_tieda p.below.0
## Min. :0.100 Min. :0.07857 Min. :0.004036 Min. :0.0000
## 1st Qu.:1.072 1st Qu.:0.32689 1st Qu.:0.023895 1st Qu.:0.3158
## Median :1.210 Median :0.48420 Median :0.034619 Median :0.3730
## Mean :1.139 Mean :0.51019 Mean :0.045569 Mean :0.3487
## 3rd Qu.:1.282 3rd Qu.:0.67320 3rd Qu.:0.051254 3rd Qu.:0.4116
## Max. :1.352 Max. :0.99715 Max. :1.265591 Max. :0.4849
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 1
## tmm omega beta_tieda p.below.0
## Min. :1 Min. :0.09131 Min. :0.004582 Min. :0.0000
## 1st Qu.:1 1st Qu.:0.36400 1st Qu.:0.024829 1st Qu.:0.3011
## Median :1 Median :0.52523 Median :0.035162 Median :0.3605
## Mean :1 Mean :0.54385 Mean :0.045850 Mean :0.3361
## 3rd Qu.:1 3rd Qu.:0.70825 3rd Qu.:0.051346 3rd Qu.:0.4005
## Max. :1 Max. :0.99758 Max. :1.223416 Max. :0.4814
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 2
## tmm omega beta_tieda p.below.0
## Min. :0.9133 Min. :0.09707 Min. :0.004825 Min. :0.0000
## 1st Qu.:0.9279 1st Qu.:0.37976 1st Qu.:0.025222 1st Qu.:0.2947
## Median :0.9442 Median :0.54202 Median :0.035407 Median :0.3550
## Mean :0.9710 Mean :0.55764 Mean :0.045983 Mean :0.3307
## 3rd Qu.:0.9793 3rd Qu.:0.72200 3rd Qu.:0.051398 3rd Qu.:0.3955
## Max. :2.0000 Max. :0.99774 Max. :1.207403 Max. :0.4798
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 4
## tmm omega beta_tieda p.below.0
## Min. :0.8340 Min. :0.1019 Min. :0.005028 Min. :0.0000
## 1st Qu.:0.8610 1st Qu.:0.3926 1st Qu.:0.025507 1st Qu.:0.2893
## Median :0.8916 Median :0.5554 Median :0.035647 Median :0.3503
## Mean :0.9489 Mean :0.5687 Mean :0.046095 Mean :0.3262
## 3rd Qu.:0.9590 3rd Qu.:0.7327 3rd Qu.:0.051461 3rd Qu.:0.3916
## Max. :4.0000 Max. :0.9979 Max. :1.195109 Max. :0.4784
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 6
## tmm omega beta_tieda p.below.0
## Min. :0.7909 Min. :0.1096 Min. :0.005348 Min. :0.0000
## 1st Qu.:0.8241 1st Qu.:0.4122 1st Qu.:0.026017 1st Qu.:0.2810
## Median :0.8622 Median :0.5755 Median :0.035895 Median :0.3433
## Mean :0.9396 Mean :0.5853 Mean :0.046275 Mean :0.3192
## 3rd Qu.:0.9473 3rd Qu.:0.7484 3rd Qu.:0.051595 3rd Qu.:0.3850
## Max. :6.0000 Max. :0.9980 Max. :1.177399 Max. :0.4760
## --------------------------------------------------------
## new.result.df[, "rho.class"]: 10
## tmm omega beta_tieda p.below.0
## Min. : 0.7397 Min. :0.1124 Min. :0.005462 Min. :0.0000
## 1st Qu.: 0.7798 1st Qu.:0.4191 1st Qu.:0.026165 1st Qu.:0.2781
## Median : 0.8265 Median :0.5823 Median :0.036004 Median :0.3407
## Mean : 0.9321 Mean :0.5909 Mean :0.046340 Mean :0.3168
## 3rd Qu.: 0.9328 3rd Qu.:0.7537 3rd Qu.:0.051635 3rd Qu.:0.3828
## Max. :10.0000 Max. :0.9981 Max. :1.171550 Max. :0.4752
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)
