library(Matrix)
library(qpgraph)

load("param.RData")
load("eqtlnet.q0.fdr.nrr.sel.RData")

## remove cutoffs of the final selected eQTL network
eqtlnet.q0.nrr <- resetCutoffs(eqtlnet.q0.fdr.nrr.sel)
## store the p-values of the marginal LRT for every gene-marker pair
pvalmat <- eqtlnet.q0.nrr@pvaluesG0[geneNames(param), markerNames(param)]

## select a qp-graph with a conservative NRR cutoff epsilon < 0.1
eqtlnet.q0.nrr.eps01 <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.nrr, epsilon=0.1)

## get all cis eQTLs at cis radius <= 500bp and 10kb
ciseqtl.q0.nrr.eps01.cisr500 <- ciseQTL(eqtlnet.q0.nrr.eps01, cisr=500)
ciseqtl.q0.nrr.eps01.cisr10000 <- ciseQTL(eqtlnet.q0.nrr.eps01, cisr=10000)

## figure out what the raw p-value cutoff is for the number of selected eQTLs at epsilon < 0.1
neQTLeps01 <- nrow(alleQTL(eqtlnet.q0.nrr.eps01))
pvalcutoffeps01 <- sort(as.vector(pvalmat))[neQTLeps01]
eqtlnet.q0.nrr.pvaleps01 <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.nrr,
                                                p.value=pvalcutoffeps01, method="none")

## get all cis eQTLs at cis radius <= 500bp and 10kb
ciseqtl.q0.nrr.pvaleps01.cisr500 <- ciseQTL(eqtlnet.q0.nrr.pvaleps01, cisr=500)
ciseqtl.q0.nrr.pvaleps01.cisr10000 <- ciseQTL(eqtlnet.q0.nrr.pvaleps01, cisr=10000)

## calculate cis enrichment with a conservative NRR cutoff
cisenrichmentConservativeTable <- array(c(nrow(ciseqtl.q0.nrr.eps01.cisr500),
                                          nrow(ciseqtl.q0.nrr.pvaleps01.cisr500),
                                          nrow(ciseqtl.q0.nrr.eps01.cisr500)/nrow(ciseqtl.q0.nrr.pvaleps01.cisr500),
                                          nrow(ciseqtl.q0.nrr.eps01.cisr10000),
                                          nrow(ciseqtl.q0.nrr.pvaleps01.cisr10000),
                                          nrow(ciseqtl.q0.nrr.eps01.cisr10000)/nrow(ciseqtl.q0.nrr.pvaleps01.cisr10000)),
                                        dim=c(3, 2), dimnames=list(Method=c("qpgraph", "marginal", "enrichment"),
                                                                   cisRadius=c("500bp", "10kb")))
cisenrichmentConservativeTable
##             cisRadius
## Method             500bp         10kb
##   qpgraph    104.0000000 1469.0000000
##   marginal    56.0000000  784.0000000
##   enrichment   1.8571429    1.8737245

## select a qp-graph with a liberal NRR cutoff epsilon < 0.5
eqtlnet.q0.nrr.eps05 <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.nrr, epsilon=0.5)

## get all cis eQTLs at cis radius <= 500bp and 10kb
ciseqtl.q0.nrr.eps05.cisr500 <- ciseQTL(eqtlnet.q0.nrr.eps05, cisr=500)
ciseqtl.q0.nrr.eps05.cisr10000 <- ciseQTL(eqtlnet.q0.nrr.eps05, cisr=10000)

## figure out what the raw p-value cutoff is for the number of selected eQTLs at epsilon < 0.5
neQTLeps05 <- nrow(alleQTL(eqtlnet.q0.nrr.eps05))
pvalcutoffeps05 <- sort(as.vector(pvalmat))[neQTLeps05]
eqtlnet.q0.nrr.pvaleps05 <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.nrr,
                                                p.value=pvalcutoffeps05, method="none")

## get all cis eQTLs at cis radius <= 500bp and 10kb
ciseqtl.q0.nrr.pvaleps05.cisr500 <- ciseQTL(eqtlnet.q0.nrr.pvaleps05, cisr=500)
ciseqtl.q0.nrr.pvaleps05.cisr10000 <- ciseQTL(eqtlnet.q0.nrr.pvaleps05, cisr=10000)

## calculate cis enrichment with a liberal NRR cutoff
cisenrichmentLiberalTable <- array(c(nrow(ciseqtl.q0.nrr.eps05.cisr500),
                                     nrow(ciseqtl.q0.nrr.pvaleps05.cisr500),
                                     nrow(ciseqtl.q0.nrr.eps05.cisr500)/nrow(ciseqtl.q0.nrr.pvaleps05.cisr500),
                                     nrow(ciseqtl.q0.nrr.eps05.cisr10000),
                                     nrow(ciseqtl.q0.nrr.pvaleps05.cisr10000),
                                     nrow(ciseqtl.q0.nrr.eps05.cisr10000)/nrow(ciseqtl.q0.nrr.pvaleps05.cisr10000)),
                                   dim=c(3, 2), dimnames=list(Method=c("qpgraph", "marginal", "enrichment"),
                                                                   cisRadius=c("500bp", "10kb")))
cisenrichmentLiberalTable
##             cisRadius
## Method        500bp        10kb
##   qpgraph    369.00 5878.000000
##   marginal   225.00 3390.000000
##   enrichment   1.64    1.733923


## build the dot plot showing the cis enrichment
labelPanel <- function(lab, font=2, cex=2, offsetx=0.05, offsety=0.05) {
  par(xpd=TRUE)
  w <- par("usr")[2] - par("usr")[1]
  h <- par("usr")[4] - par("usr")[3]
  text(par("usr")[1]-w*offsetx, par("usr")[4]+h*offsety, lab, font=font, cex=cex)
  par(xpd=FALSE)
}

png(file="Figure3.png", height=9, width=9, units="in", res=300)
par(mfrow=c(2, 2), mar=c(3, 3, 1, 2), mgp=c(1, 1, 0))
plot(eqtlnet.q0.nrr.eps05, xlab="Ordered markers", ylab="Ordered genes", axes=FALSE, cex.lab=2)
labelPanel(LETTERS[1], offsetx=0.05, offsety=0.01)
plot(eqtlnet.q0.nrr.pvaleps05, xlab="Ordered markers", ylab="Ordered genes", axes=FALSE, cex.lab=2)
labelPanel(LETTERS[2], offsetx=0.05, offsety=0.01)
plot(eqtlnet.q0.nrr.eps01, xlab="Ordered markers", ylab="Ordered genes", axes=FALSE, cex.lab=2)
labelPanel(LETTERS[3], offsetx=0.05, offsety=0.01)
plot(eqtlnet.q0.nrr.pvaleps01, xlab="Ordered markers", ylab="Ordered genes", axes=FALSE, cex.lab=2)
labelPanel(LETTERS[4], offsetx=0.05, offsety=0.01)
dev.off()