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)
library(xtable)
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.min.variance2 <- subset(merged.mouse.eQTL.min.variance2, neg_log_lung_pvalue <= 30 & liver_residual_variance <=2.5)
library(ggplot2)
# ggplot(merged.mouse.eQTL.min.variance2, aes(x=abs_lung.beta, y=neg_log_lung_pvalue)) +geom_point()+geom_smooth(method=lm)
lm(abs_lung.beta ~ neg_log_lung_pvalue, data=merged.mouse.eQTL.min.variance2)
##
## Call:
## lm(formula = abs_lung.beta ~ neg_log_lung_pvalue, data = merged.mouse.eQTL.min.variance2)
##
## Coefficients:
## (Intercept) neg_log_lung_pvalue
## 0.004804 0.017509
cor.test(merged.mouse.eQTL.min.variance2$abs_lung.beta, merged.mouse.eQTL.min.variance2$neg_log_lung_pvalue, method = "spearman")
## Warning in
## cor.test.default(merged.mouse.eQTL.min.variance2$abs_lung.beta, : Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: merged.mouse.eQTL.min.variance2$abs_lung.beta and merged.mouse.eQTL.min.variance2$neg_log_lung_pvalue
## S = 2.8698e+10, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8707747
cor(merged.mouse.eQTL.min.variance2$abs_lung.beta, merged.mouse.eQTL.min.variance2$neg_log_lung_pvalue, method = "spearman")
## [1] 0.8707747
# fit1<-summary(lm(abs_liver.beta ~ abs_lung.beta, data=merged.mouse.eQTL.min.variance2))
# fit1
# tau<-fit1$sigma**2
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)
# liver.betas=Z*gama+T^2
lmsummary<-summary(lm(abs_liver.beta~-1+Z, data=merged.mouse.eQTL))
lmsummary
##
## Call:
## lm(formula = abs_liver.beta ~ -1 + Z, data = merged.mouse.eQTL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59950 -0.03173 -0.01859 0.00630 2.04058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ZX1 0.026447 0.001116 23.70 <2e-16 ***
## ZZ 0.490621 0.010731 45.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09817 on 11002 degrees of freedom
## Multiple R-squared: 0.3313, Adjusted R-squared: 0.3312
## F-statistic: 2726 on 2 and 11002 DF, p-value: < 2.2e-16
tau<-lmsummary$sigma**2
tau
## [1] 0.009637149
outputVector<-c(lmsummary$coefficients[,1],lmsummary$coefficients[,2])
write.table(matrix(outputVector,length(Z[1,])),file="2016-04-01_hm_tau_coefficients with intercept.txt",col.names=FALSE,row.names=FALSE,quote=FALSE)
hm_tau_coefficients<-read.table(file="2016-04-01_hm_tau_coefficients with intercept.txt", header=T)
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 )
summary(omega.diag)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09185 0.36550 0.52670 0.54510 0.70940 0.99540
# 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
summary(regbeta)
## V1
## Min. :0.02645
## 1st Qu.:0.03486
## Median :0.04241
## Mean :0.05425
## 3rd Qu.:0.05718
## Max. :0.93746
markers1<-as.character(markers)
outputVector<-c(markers1,betas.hat,betas.tieda)
write.table(matrix(outputVector,rowLength),file="2016-04-01_hm_tau_hmresults.txt",col.names=FALSE,row.names=FALSE,quote=FALSE)
liver.mouse.eQTL.bayesian<-read.table(file="2016-04-01_hm_tau_hmresults.txt")
colnames(liver.mouse.eQTL.bayesian)<-c( "ensembl_id", "betas.hat","betas.tieda")
liver.mouse.eQTL.bayesian.all<- merge(liver.mouse.eQTL.bayesian, merged.mouse.eQTL.min.variance2, by = "ensembl_id")
write.table(liver.mouse.eQTL.bayesian.all,file="2016-04-01_liver.mouse.eQTL.bayesian.all.txt")
liver.mouse.eQTL.bayesian<-read.table(file="2016-04-01_liver.mouse.eQTL.bayesian.all.txt")
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat betas.tieda liver.beta liver.beta_se
## 1 ENSMUSG00000000001 0.015986111 0.03353994 -0.015986111 0.025624670
## 2 ENSMUSG00000000003 0.016880682 0.02363075 0.016880682 0.012653279
## 3 ENSMUSG00000000037 0.003865079 0.01407009 -0.003865079 0.013629289
## 4 ENSMUSG00000000049 0.042233333 0.04158436 0.042233333 0.007808461
## 5 ENSMUSG00000000056 0.063604072 0.03999481 0.063604072 0.051481198
## 6 ENSMUSG00000000058 0.025611111 0.03063312 -0.025611111 0.012982317
## liver_residual_variance lung_pvalue lung.beta lung.beta_se
## 1 0.018910763 0.12293526 -0.03269000 0.02077547
## 2 0.003757142 0.06295461 0.02954960 0.01547936
## 3 0.004681090 0.30681836 -0.01759477 0.01701272
## 4 0.001829162 0.79280526 -0.02388393 0.09036397
## 5 0.078095912 0.21831618 0.02167500 0.01734983
## 6 0.004853968 0.29250998 0.02885470 0.02707600
## liver_pvalue abs_liver.beta abs_lung.beta neg_log_lung_pvalue
## 1 5.377718e-01 0.015986111 0.03269000 2.0960974
## 2 1.929225e-01 0.016880682 0.02954960 2.7653413
## 3 7.788138e-01 0.003865079 0.01759477 1.1814994
## 4 9.084753e-06 0.042233333 0.02388393 0.2321777
## 5 2.269175e-01 0.063604072 0.02167500 1.5218109
## 6 5.846922e-02 0.025611111 0.02885470 1.2292565
liver.mouse.eQTL.bayesian<-subset(liver.mouse.eQTL.bayesian, select = c("ensembl_id", "betas.hat",
"liver.beta_se", "betas.tieda",
"liver_residual_variance", "liver_pvalue", "abs_lung.beta",
"abs_lung.beta", "neg_log_lung_pvalue"))
Tau_invert<-diag.inverse(Tau)
V_invert<-diag.inverse(V)
PS_invert<-Tau_invert+V_invert%*% Z %*% Z_transpose
PS<-diag.inverse(PS_invert)
ps<-diag(PS)
range(ps)
## [1] 0.0007348675 0.0095926506
library(reshape)
ps.long <- melt(ps)
head(ps.long)
## value
## 1 0.006381557
## 2 0.002701554
## 3 0.003150036
## 4 0.001536628
## 5 0.008578102
## 6 0.003226289
ps.long$betas.tieda.se<-(ps.long$value)^0.5
liver.mouse.eQTL.bayesian<-cbind(liver.mouse.eQTL.bayesian,ps.long$betas.tieda.se)
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat liver.beta_se betas.tieda
## 1 ENSMUSG00000000001 0.015986111 0.025624670 0.03353994
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02363075
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01407009
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.04158436
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.03999481
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.03063312
## liver_residual_variance liver_pvalue abs_lung.beta abs_lung.beta.1
## 1 0.018910763 5.377718e-01 0.03269000 0.03269000
## 2 0.003757142 1.929225e-01 0.02954960 0.02954960
## 3 0.004681090 7.788138e-01 0.01759477 0.01759477
## 4 0.001829162 9.084753e-06 0.02388393 0.02388393
## 5 0.078095912 2.269175e-01 0.02167500 0.02167500
## 6 0.004853968 5.846922e-02 0.02885470 0.02885470
## neg_log_lung_pvalue ps.long$betas.tieda.se
## 1 2.0960974 0.07988465
## 2 2.7653413 0.05197648
## 3 1.1814994 0.05612518
## 4 0.2321777 0.03919985
## 5 1.5218109 0.09261804
## 6 1.2292565 0.05680043
liver.mouse.eQTL.bayesian<-rename(liver.mouse.eQTL.bayesian, c("ps.long$betas.tieda.se"="betas.tieda.se", "liver.beta_se"="betas.hat.se"))
liver.mouse.eQTL.bayesian<-subset(liver.mouse.eQTL.bayesian, select = c("ensembl_id", "betas.hat", "betas.hat.se",
"betas.tieda", "betas.tieda.se", "liver_residual_variance",
"liver_pvalue", "abs_lung.beta", "neg_log_lung_pvalue"))
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat betas.hat.se betas.tieda betas.tieda.se
## 1 ENSMUSG00000000001 0.015986111 0.025624670 0.03353994 0.07988465
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02363075 0.05197648
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01407009 0.05612518
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.04158436 0.03919985
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.03999481 0.09261804
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.03063312 0.05680043
## 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
# library(tigerstats)
# pnormGC(0, region="below", mean=0.002352829,sd=0.09972950)
liver.mouse.eQTL.bayesian$p.below.0<-pnorm(0,liver.mouse.eQTL.bayesian$betas.tieda, liver.mouse.eQTL.bayesian$betas.tieda.se)
range(liver.mouse.eQTL.bayesian$p.below.0)
## [1] 2.170068e-77 4.534196e-01
write.table(liver.mouse.eQTL.bayesian,file="2016-04-01_liver.mouse.eQTL.bayesian with beta.txt")
liver.mouse.eQTL.bayesian <- read.table(file="2016-04-01_liver.mouse.eQTL.bayesian with beta.txt")
liver.mouse.eQTL.bayesian.4tau <- liver.mouse.eQTL.bayesian
rowLength<-nrow(liver.mouse.eQTL.bayesian.4tau)
min_lv <- log(min(liver.mouse.eQTL.bayesian.4tau$liver_residual_variance)/tau)
max_lv<-log(max(liver.mouse.eQTL.bayesian.4tau$liver_residual_variance)/tau)
min_nlogp <- min(liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue)
max_nlogp<-max(liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue)
liver.mouse.eQTL.bayesian.4tau$Tmm <- exp((max_lv-min_lv)/(min_nlogp - max_nlogp)*(liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue-min_nlogp)+max_lv)
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.03353994 0.07988465
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02363075 0.05197648
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01407009 0.05612518
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.04158436 0.03919985
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.03999481 0.09261804
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.03063312 0.05680043
## 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 Tmm
## 1 0.3372958 125.9912
## 2 0.3246830 106.0681
## 3 0.4010264 159.4049
## 4 0.1443837 203.4895
## 5 0.3329342 146.0456
## 6 0.2948360 157.4589
qplot(neg_log_lung_pvalue,Tmm, data = liver.mouse.eQTL.bayesian.4tau)
liver.mouse.eQTL.bayesian.4tau$nTau <- liver.mouse.eQTL.bayesian.4tau$Tmm*tau
summary(liver.mouse.eQTL.bayesian.4tau$nTau)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0009747 0.8371000 1.3470000 1.2330000 1.6920000 2.0820000
qplot(neg_log_lung_pvalue, nTau, data = liver.mouse.eQTL.bayesian.4tau)
nTau<- diag(liver.mouse.eQTL.bayesian.4tau$nTau, rowLength, rowLength)
ns<-V + nTau
nS <- diag.inverse(ns)
nomega<-diag.multi(nS, V)
liver.mouse.eQTL.bayesian.4tau$n.omega <- diag(nomega )
liver.mouse.eQTL.bayesian.4tau$n.beta_tieda <-nomega %*% Z %*% gamma + (identity-nomega) %*% betas.hat
nTau_invert<-diag.inverse(nTau)
V_invert<-diag.inverse(V)
nPS_invert<-nTau_invert + (V_invert%*% Z%*%Z_transpose)
nPS<-diag.inverse(nPS_invert)
nps<-diag(nPS)
nps.long <- melt(nps)
liver.mouse.eQTL.bayesian.4tau$n.betas.tieda.se <- (nps.long$value)^0.5
liver.mouse.eQTL.bayesian.4tau$n.p.below.0 <- pnorm(0, liver.mouse.eQTL.bayesian.4tau$n.beta_tieda, liver.mouse.eQTL.bayesian.4tau$n.betas.tieda.se)
summary(liver.mouse.eQTL.bayesian.4tau$liver_residual_variance)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0009747 0.0055510 0.0107300 0.0283400 0.0235300 2.0820000
summary(liver.mouse.eQTL.bayesian.4tau$nTau)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0009747 0.8371000 1.3470000 1.2330000 1.6920000 2.0820000
nrow(liver.mouse.eQTL.bayesian.4tau)
## [1] 11004
sorted.liver.mouse.eQTL.bayesian.4tau <- liver.mouse.eQTL.bayesian.4tau[order(liver.mouse.eQTL.bayesian.4tau$betas.hat),]
#head(sorted.liver.mouse.eQTL.bayesian.4tau)
#tail(sorted.liver.mouse.eQTL.bayesian.4tau)
summary(sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.01715 0.03254 0.05668 0.06264 1.85700
sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan[sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta < 0.01715] <- "a"
sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan[sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta >= 0.01715 & sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta < 0.03257 ] <- "b"
sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan[sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta >= 0.03257 & sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta < 0.06265 ] <- "c"
sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan[sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta >= 0.06265 ] <- "d"
sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan <- as.factor(sorted.liver.mouse.eQTL.bayesian.4tau$abs_lung.beta_quan)
# ggplot(sorted.liver.mouse.eQTL.bayesian.4tau, aes(x=betas.hat, y=betas.tieda, color=abs_lung.beta_quan)) +geom_point()+geom_smooth(method=lm)
options(stringsAsFactors=FALSE)
library(EBcoexpress)
## Loading required package: EBarrays
## Loading required package: Biobase
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, as.vector, cbind,
## colnames, do.call, duplicated, eval, evalq, Filter, Find, get,
## grep, grepl, intersect, is.unsorted, lapply, lengths, Map,
## mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
## setdiff, sort, table, tapply, union, unique, unlist, unsplit
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: lattice
## Loading required package: mclust
## Package 'mclust' version 5.1
## Type 'citation("mclust")' for citing this R package in publications.
## Loading required package: minqa
##
## Attaching package: 'EBcoexpress'
## The following object is masked from 'package:EBarrays':
##
## crit.fun
liver.mouse.eQTL.bayesian.4tau$posteriorprobability <- liver.mouse.eQTL.bayesian.4tau$p.below.0
threshold1 <- crit.fun(liver.mouse.eQTL.bayesian.4tau$posteriorprobability, 0.05)
threshold1
## [1] 0.8883265
newdata1 <- subset(liver.mouse.eQTL.bayesian.4tau, 1-posteriorprobability >= threshold1)
dim(newdata1)
## [1] 795 17
threshold2 <- crit.fun(liver.mouse.eQTL.bayesian.4tau$posteriorprobability, 0.1)
threshold2
## [1] 0.8118636
newdata2 <- subset(liver.mouse.eQTL.bayesian.4tau, 1-posteriorprobability >= threshold2)
dim(newdata2)
## [1] 1544 17
threshold3 <- crit.fun(liver.mouse.eQTL.bayesian.4tau$posteriorprobability, 0.2)
threshold3
## [1] 0.7062924
newdata3 <- subset(liver.mouse.eQTL.bayesian.4tau, 1-posteriorprobability >= threshold3)
dim(newdata3)
## [1] 4522 17
liver.ASE <- read.csv(file= "ASE.genetics.113.153882-6.csv")
liver.ASE.symbol <- unique(liver.ASE$geneID)
length(liver.ASE.symbol)
## [1] 440
head(liver.ASE.symbol)
## [1] "Rrs1" "6330578E17Rik" "Aox3" "Cps1"
## [5] "Ugt1a10" "Heatr7b1"
#liver.ASE <-liver.ASE[liver.ASE$geneID %in% names(table(liver.ASE$geneID))[table(liver.ASE$geneID) >=2], ]
sub.liver.ASE <-liver.ASE[which(liver.ASE$pvalBH.DxB7 < 0.05 & liver.ASE$pvalBH.BxD7 < 0.05), ]
dim(sub.liver.ASE)
## [1] 1044 19
liver.ASE.symbol <- unique(liver.ASE$geneID)
liver.ASE.symbol <- unique(sub.liver.ASE$geneID)
length(liver.ASE.symbol)
## [1] 367
head(liver.ASE.symbol)
## [1] "Rrs1" "6330578E17Rik" "Aox3" "Cps1"
## [5] "Ugt1a10" "Heatr7b1"
library(biomaRt)
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
liver.ASE.ensembl <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol" , values =liver.ASE.symbol ,mart=mouse)
dim(liver.ASE.ensembl)
## [1] 327 2
length(unique(liver.ASE.ensembl$ensembl_gene_id))
## [1] 327
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.03353994 0.07988465
## 2 ENSMUSG00000000003 0.016880682 0.012653279 0.02363075 0.05197648
## 3 ENSMUSG00000000037 0.003865079 0.013629289 0.01407009 0.05612518
## 4 ENSMUSG00000000049 0.042233333 0.007808461 0.04158436 0.03919985
## 5 ENSMUSG00000000056 0.063604072 0.051481198 0.03999481 0.09261804
## 6 ENSMUSG00000000058 0.025611111 0.012982317 0.03063312 0.05680043
## 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 Tmm nTau n.omega n.beta_tieda n.betas.tieda.se
## 1 0.3372958 125.9912 1.214196 0.0153358734 0.016392504 0.13638613
## 2 0.3246830 106.0681 1.022194 0.0036621059 0.016968807 0.06115659
## 3 0.4010264 159.4049 1.536209 0.0030379125 0.003959906 0.06830395
## 4 0.1443837 203.4895 1.961059 0.0009318726 0.042229542 0.04273659
## 5 0.3329342 146.0456 1.407463 0.0525700369 0.062209774 0.27195126
## 6 0.2948360 157.4589 1.517455 0.0031885566 0.025658916 0.06953042
## n.p.below.0 posteriorprobability
## 1 0.4521656 0.3372958
## 2 0.3907116 0.3246830
## 3 0.4768844 0.4010264
## 4 0.1615432 0.1443837
## 5 0.4095303 0.3329342
## 6 0.3560521 0.2948360
dim(liver.mouse.eQTL.bayesian.4tau)
## [1] 11004 17
liver.mouse.eQTL.bayesian.4tau$rank.p.below.0 <- rank(liver.mouse.eQTL.bayesian.4tau$p.below.0)
liver.mouse.eQTL.bayesian.4tau$rank.n.p.below.0 <- rank(liver.mouse.eQTL.bayesian.4tau$n.p.below.0)
liver.mouse.eQTL.bayesian.4tau$rank.liver.pvalue <- rank(liver.mouse.eQTL.bayesian.4tau$liver_pvalue)
sorted.liver.mouse.eQTL.bayesian.4tau <- liver.mouse.eQTL.bayesian.4tau[order(liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue), ]
tail(sorted.liver.mouse.eQTL.bayesian.4tau)
## ensembl_id betas.hat betas.hat.se betas.tieda
## 10430 ENSMUSG00000070583 0.10587054 0.01454742 0.1526573
## 815 ENSMUSG00000007656 0.19846429 0.05320415 0.1016293
## 2602 ENSMUSG00000022620 0.12533333 0.05152972 0.1114037
## 4865 ENSMUSG00000028656 1.06462698 0.02702678 0.7892947
## 8585 ENSMUSG00000044827 0.04306109 0.01035170 0.1117917
## 8654 ENSMUSG00000045594 0.27253968 0.03469755 0.1811110
## betas.tieda.se liver_residual_variance liver_pvalue abs_lung.beta
## 10430 0.05896365 0.006320602 6.359505e-08 0.4026468
## 815 0.09292958 0.084543033 8.618941e-04 0.1307400
## 2602 0.09257664 0.079659376 2.165185e-02 0.1697263
## 4865 0.06395509 0.018407266 4.552986e-26 1.2610486
## 8585 0.04323172 0.003157578 2.733728e-04 0.6015134
## 8654 0.08485335 0.030338790 1.481601e-08 0.2560458
## neg_log_lung_pvalue p.below.0 Tmm nTau n.omega
## 10430 27.69862 4.812648e-03 0.1739836 0.0016767064 0.7903412
## 815 27.84248 1.370617e-01 0.1676635 0.0016157983 0.9812463
## 2602 28.89764 1.144172e-01 0.1278132 0.0012317551 0.9847727
## 4865 29.46768 2.710558e-35 0.1103830 0.0010637771 0.9453662
## 8585 29.63010 4.856650e-03 0.1058666 0.0010202523 0.7557938
## 8654 29.80771 1.640501e-02 0.1011394 0.0009746952 0.9688730
## n.beta_tieda n.betas.tieda.se n.p.below.0 posteriorprobability
## 10430 0.1992285 0.03579954 1.309915e-08 4.812648e-03
## 815 0.0926140 0.03981191 1.000157e-02 1.370617e-01
## 2602 0.1099562 0.03482050 7.948168e-04 1.144172e-01
## 4865 0.6680621 0.03041826 3.286415e-107 2.710558e-35
## 8585 0.2535505 0.02661764 8.199868e-22 4.856650e-03
## 8654 0.1558185 0.03069907 1.930601e-07 1.640501e-02
## rank.p.below.0 rank.n.p.below.0 rank.liver.pvalue
## 10430 130 20 390
## 815 1005 153 1238
## 2602 817 79 2336
## 4865 3 1 5
## 8585 131 6 1060
## 8654 219 25 323
result <- matrix(, nrow(liver.mouse.eQTL.bayesian.4tau), 8)
colnames(result)<-c("bayrank","bay_PPV","bay_TPR","bay_FDR", "orirank","ori_PPV","ori_TPR","ori_FPR" )
for (i in 1:nrow(liver.mouse.eQTL.bayesian.4tau))
{
newdata1 <- subset(liver.mouse.eQTL.bayesian.4tau, rank.n.p.below.0 <= i)
overlap.newdata1 <- newdata1[newdata1$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
result[i, 1] <- i
result[i, 2] <- nrow(overlap.newdata1)/nrow(newdata1)
result[i, 3] <- nrow(overlap.newdata1)/nrow(liver.ASE.ensembl)
newdata2 <- subset(liver.mouse.eQTL.bayesian.4tau, rank.liver.pvalue <= i)
overlap.newdata2 <- newdata2[newdata2$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
result[i, 5] <- i
result[i, 6] <- nrow(overlap.newdata2)/nrow(newdata2)
result[i, 7] <- nrow(overlap.newdata2)/nrow(liver.ASE.ensembl)
}
head(result)
## bayrank bay_PPV bay_TPR bay_FDR orirank ori_PPV ori_TPR
## [1,] 1 0.0000000 0.000000000 NA 1 0.0000000 0.000000000
## [2,] 2 0.0000000 0.000000000 NA 2 0.5000000 0.003058104
## [3,] 3 0.3333333 0.003058104 NA 3 0.6666667 0.006116208
## [4,] 4 0.2500000 0.003058104 NA 4 0.5000000 0.006116208
## [5,] 5 0.2000000 0.003058104 NA 5 0.4000000 0.006116208
## [6,] 6 0.1666667 0.003058104 NA 6 0.5000000 0.009174312
## ori_FPR
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NA
tail(result)
## bayrank bay_PPV bay_TPR bay_FDR orirank ori_PPV ori_TPR
## [10999,] 10999 0.02409310 0.8103976 NA 10999 0.02409310 0.8103976
## [11000,] 11000 0.02409091 0.8103976 NA 11000 0.02409091 0.8103976
## [11001,] 11001 0.02408872 0.8103976 NA 11001 0.02408872 0.8103976
## [11002,] 11002 0.02408653 0.8103976 NA 11002 0.02408653 0.8103976
## [11003,] 11003 0.02408434 0.8103976 NA 11003 0.02408434 0.8103976
## [11004,] 11004 0.02408215 0.8103976 NA 11004 0.02408215 0.8103976
## ori_FPR
## [10999,] NA
## [11000,] NA
## [11001,] NA
## [11002,] NA
## [11003,] NA
## [11004,] NA
head(result)
## bayrank bay_PPV bay_TPR bay_FDR orirank ori_PPV ori_TPR
## [1,] 1 0.0000000 0.000000000 NA 1 0.0000000 0.000000000
## [2,] 2 0.0000000 0.000000000 NA 2 0.5000000 0.003058104
## [3,] 3 0.3333333 0.003058104 NA 3 0.6666667 0.006116208
## [4,] 4 0.2500000 0.003058104 NA 4 0.5000000 0.006116208
## [5,] 5 0.2000000 0.003058104 NA 5 0.4000000 0.006116208
## [6,] 6 0.1666667 0.003058104 NA 6 0.5000000 0.009174312
## ori_FPR
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NA
plot(result[, 1], result[, 3], type="l", col="red", xlab="Ranking", ylab="TPR", ylim=c(0, 1) )
par(new=TRUE)
plot( result[, 1], result[, 7], type="l", col="green", xlab="Ranking", ylab="TPR", ylim=c(0, 1) )
legend("topright", inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
plot(result[, 1], result[, 2], type="l", col="red", xlab="Ranking", ylab="PPV", ylim=c(0, 1))
par(new=TRUE)
plot(result[, 5], result[, 6], type="l", col="green", xlab="Ranking", ylab="PPV", ylim=c(0, 1))
legend("topright", inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
plot(result[, 1], result[, 2], type="l", col="red", xlab="Ranking", ylab="PPV", ylim=c(0, 1), xlim=c(0, 300))
par(new=TRUE)
plot(result[, 5], result[, 6], type="l", col="green", xlab="Ranking", ylab="PPV",
ylim=c(0, 1), xlim=c(0, 300))
legend("topright", inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
