## Ref. /ddn/gs1/home/jinj2/LondonData/PACE/asthma/analysis/reesese/forPaper/forest_plots_update.r
library(forestplot)
library(compare)
library(metafor)

## Read in the meta data and the CpG lists of interest
path <- c("/ddn/gs1/home/jinj2/LondonData/START/AdultSmoking/MetaAnalysis/MetaAnalysis_ALHS_Results_Adj_For_Neu_CORRECT/Datasets_For_MetaAnalysis/")
adultSmkStart <- read.csv(paste0(path,"START_AdultSmk.csv"))
adultSmkAlhs <- read.csv(paste0(path,"ALHS_AdultSmk.csv"))
adultSmkShs <- read.csv(paste0(path,"SHS_AdultSmk.csv"))
inUteroSmkExpStart <- read.csv(paste0(path,"START_InUteroSmkExposure.csv"))
inUteroSmkExpAlhs <- read.csv(paste0(path,"ALHS_InUteroSmkExposure.csv"))
etsStart <- read.csv(paste0(path,"START_ETS.csv"))
etsAlhs <- read.csv(paste0(path,"ALHS_ETS.csv"))
adultSmkSexIntStart <- read.csv(paste0(path,"START_AdultSmk_SexInteraction.csv"))
adultSmkSexIntAlhs <- read.csv(paste0(path,"ALHS_AdultSmk_SexInteraction.csv"))
lists <- read.csv("/ddn/gs1/group/london/START/AdultSmoking/Forest_Plots/FDR_Sig_CpGs_Across_4_MetaAnalyses_For_ForestPlots.csv")

## Take the meta data for a given CpG
adultSmkStart1 <- adultSmkStart[,c("Probe","Beta","SE","P")]
adultSmkAlhs1 <- adultSmkAlhs[,c("Probe","Beta","SE","P")]
adultSmkShs1 <- adultSmkShs[,c("Probe","Beta","SE","P")]
cg05575921_adultSmkStart <- data.frame(rbind(adultSmkStart1[adultSmkStart1$Probe=="cg05575921",],
    adultSmkAlhs1[adultSmkAlhs1$Probe=="cg05575921",],adultSmkShs1[adultSmkShs1$Probe=="cg05575921",]))
rownames(cg05575921_adultSmkStart) <- c("adultSmkStart","adultSmkALHS","adultSmkSHS")

## Meta analysis using RMA (robust multiple average)
meta <- data.frame(cg05575921_adultSmkStart$Beta, cg05575921_adultSmkStart$SE)
rownames(meta) <- rownames(cg05575921_adultSmkStart)
colnames(meta) <- c("coef","se")
res <- rma(yi=meta$coef, sei=meta$se,weighted=TRUE)
Meta <- data.frame(cg05575921_adultSmkStart$Probe[1], res$beta, res$se, res$pval)
colnames(Meta) <- c("Probe","Beta","SE","P")
Mat <- rbind(cg05575921_adultSmkStart,Meta)
rownames(Mat)[4] <- "Meta_Analysis"

## Calculate the mean, ci of odds_ratio 1 per
source("/ddn/gs1/group/london/norway/methylation/reesese/functions/f.conf_int_methy.r")	#f.format.perc, f.conf_int_methy.r
f.odds_ratio_1per<-function(x){exp(x/100)}
f.conf_int_methy <- function(results, level = 0.95, per1perOR=TRUE, ...) {
        cf <- results$Beta 
        se <- results$SE
        parm <- results$Probe
        a <- (1 - level)/2
        a <- c(a, 1 - a)
        pct <- f.format.perc(a, 3)
        fac <- qnorm(a)
        ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,pct))
        ci[] <- cf + se %o% fac
        if (per1perOR){
                exp(ci/100)
        }else{
                ci
        }
}
Mat$mean <- f.odds_ratio_1per(Mat$Beta)
Mat$CIlower <- data.frame(f.conf_int_methy(Mat))[,1]
Mat$CIupper <- data.frame(f.conf_int_methy(Mat))[,2]

## ploting out forest plot
## Generate table text
##Data
meta.dat <- data.frame(
mean  = c(NA, Mat$Beta),
lower = c(NA, Mat$CIlower),
upper = c(NA, Mat$CIupper))
meta.dat

##For making table text
meta <- data.frame(
mean  = Mat$Beta,
lower = Mat$CIlower,
upper = Mat$CIupper)
meta

cohorts <- rownames(Mat)
#tabletext<-cbind(c("Study", paste0(cohorts,"(N=",ntotal,")")),
tabletext<-cbind(c("Study", cohorts), c("Mean",round(meta$mean,6)),  
                 c("CI2.5", round(meta$lower,6)),c("CI97.5", round(meta$upper,6)))
pdf("cg05575921_adultSmkStart.pdf")
#forestplot(tabletext, meta.dat, new_page = TRUE,
forestplot(tabletext, meta.dat, new_page = TRUE,
           title = paste0("Start",".","ForestPlot of cg055759212"),
           ##boxsize=0.25, ref.vline.at = numeric(0),   #for fixed size of box
           #txt_gp = fpTxtGp(ticks = gpar(cex=0.7), xlab  = gpar(cex = 0.9)), xlab ="Regression Coefficient",
           txt_gp = fpTxtGp(ticks = gpar(cex=0.7), xlab  = gpar(cex = 0.9)), xlab ="OR",
           #fn.ci_norm=c("fpDrawDiamondCI",rep("fpDrawNormalCI",times=11),"fpDrawDiamondCI"),
           #is.summary=c(TRUE,rep(FALSE,16),TRUE),
           is.summary=c(TRUE,rep(FALSE,4),TRUE),
           col=fpColors(box="royalblue",line="darkblue", summary="royalblue")) #
dev.off()