R version 3.6.1 (2019-07-05) -- "Action of the Toes"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #################################################################################################################################################################
> # Purpose: Make forest plot using the EWAS results from multiple cohorts and their meta result (by rma). 
> #          Note: CIs are calculated from the se values of the EWAS results and numbers of samples are used just for plot text
> # Request: Stephanie London, Mi Kyeong, 11/18/2019, Forest plots using R. 
> # Input: /ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/ARIC|GTP|Inchianti.....
> #        /ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/off-methylation-smoke-all-models123.txt 
> # Output: PACE-Charge_forestplot.pdf 
> # Ref: SampleCode_ForestPlot_forJP_MKL_11182019.R; modified_to_Sinjini_rma.r 
> # WDir: /ddn/gs1/home/jinj2/LondonData/PACE/PACE-CHARGE_smoking/Check/Jianping/Charge_cohort_forestPlot/ 
> # Script: forestPlot_cg05575921.r
> # Dec. 27, 2019, Jianping Jin
> #################################################################################################################################################################
> 
> library(data.table)
> library(metafor)
> library(forestplot)
> 
> ## Read in data (p-value, coef, se) from cohorts
> fhs <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/off-methylation-smoke-all-models123.txt"))
> fhs_red <- fhs[,which(colnames(fhs)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SD_SmokerVsNeverSmokerModel2Fave"))]
> colnames(fhs_red) = c(".id","pvalue","coef","se")
> 
> aric <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/ARIC_smoking_EWAS_model2-combined.csv"))
> aric_red <- aric[,which(colnames(aric)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(aric_red) <- c(".id","coef","se","pvalue")
> 
> gtp <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/GTP_smoking_model2-combined.csv"))
> gtp_red <- gtp[,which(colnames(gtp)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(gtp_red) <- c(".id","pvalue","coef","se")
> 
> inchianti <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/Inchianti_smoking_EWAS_model2-combined.csv"))
> inchianti_red <- inchianti[,which(colnames(inchianti)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(inchianti_red) <- c(".id","pvalue","coef","se")
> 
> kora <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/KORA_smoking_EWAS_model2-combined.csv"))
> kora_red <- kora[,which(colnames(kora)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(kora_red) <- c(".id","pvalue","coef","se")
> 
> lbc1921 <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/LBC1921_smoking_EWAS_model2-combined.csv"))
> lbc1921_red <- lbc1921[,which(colnames(lbc1921)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(lbc1921_red) <- c(".id","pvalue","coef","se")
> 
> lbc1936 <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/LBC1936_smoking_EWAS_model2-combined.csv"))
> lbc1936_red <- lbc1936[,which(colnames(lbc1936)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(lbc1936_red) <- c(".id","pvalue","coef","se")
> 
> nas <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/NAS_smoking_EWAS_model2-combined.csv"))
> nas_red <- nas[,which(colnames(nas)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(nas_red) <- c(".id","pvalue","coef","se")
> 
> rs <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/RS_smoking_EWAS_model2-combined.csv"))
> rs_red <- rs[,which(colnames(rs)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(rs_red) <- c(".id","pvalue","coef","se")
> 
> goldn <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/GOLDN_smoking_EWAS_model2-combined.csv"))
> goldn_red <- goldn[,which(colnames(goldn)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(goldn_red) <- c(".id","coef","se","pvalue")
> 
> 
> mesa <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/MESA_smoking_model2-combined.csv"))
> mesa_red <- mesa[,which(colnames(mesa)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(mesa_red) <- c(".id","coef","se","pvalue")
> 
> epic <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/EPIC_smoking_model2-combined.csv"))
> epic_red <- epic[,which(colnames(epic)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(epic_red) <- c(".id","pvalue","coef","se")
> 
> norfolk <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/EPICNorfolk_smoking_model2-combined.csv"))
> norfolk_red <- norfolk[,which(colnames(norfolk)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(norfolk_red) <- c(".id","pvalue","coef","se")
> 
> 
> chsea <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/CHS-EA_smoking_model2-combined.csv"))
> chsea_red <- chsea[,which(colnames(chsea)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(chsea_red) <- c(".id","pvalue","coef","se")
> 
> chsaa <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/CHS-AA_smoking_model2-combined.csv"))
> chsaa_red <- chsaa[,which(colnames(chsaa)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(chsaa_red) <- c(".id","pvalue","coef","se")
> 
> genoa <- data.frame(fread("/ddn/gs1/group/london/CHARGE_Methylation/Smokers-methylation/other-cohorts/GENOA_smoking_model2-combined.csv"))
> genoa_red <- genoa[,which(colnames(genoa)%in%c("V1","Fx_SmokerVsNeverSmokerModel2Fave","P_SmokerVsNeverSmokerModel2Fave","SE_SmokerVsNeverSmokerModel2Fave"))]
> colnames(genoa_red) <- c(".id","pvalue","coef","se")
> 
> ########################################################################################################################
> ## generate all cpgs
> options(digits=4)
> 
> all_cpgs <- Reduce(union, list(fhs_red$.id, aric_red$.id, gtp_red$.id, inchianti_red$.id, kora_red$.id, lbc1921_red$.id,
+ lbc1936_red$.id, nas_red$.id, rs_red$.id, goldn_red$.id, mesa_red$.id, epic_red$.id, norfolk_red$.id, chsea_red$.id,
+ chsaa_red$.id, genoa_red$.id));
> length(all_cpgs)
[1] 485577
> 
> ## Take target CpG(s)
> i <- which(all_cpgs=="cg05575921")
> cpg <- all_cpgs[i]    
> coef_cpg <- c(fhs_red[which(fhs_red$.id==cpg),"coef"],
+ aric_red[which(aric_red$.id==cpg),"coef"],
+ gtp_red[which(gtp_red$.id==cpg),"coef"],
+ inchianti_red[which(inchianti_red$.id==cpg),"coef"],
+ kora_red[which(kora_red$.id==cpg),"coef"],
+ lbc1921_red[which(lbc1921_red$.id==cpg),"coef"],
+ lbc1936_red[which(lbc1936_red$.id==cpg),"coef"],
+ nas_red[which(nas_red$.id==cpg),"coef"],
+ rs_red[which(rs_red$.id==cpg),"coef"],
+ goldn_red[which(goldn_red$.id==cpg),"coef"],
+ mesa_red[which(mesa_red$.id==cpg),"coef"],
+ epic_red[which(epic_red$.id==cpg),"coef"],
+ norfolk_red[which(norfolk_red$.id==cpg),"coef"],
+ chsea_red[which(chsea_red$.id==cpg),"coef"],
+ chsaa_red[which(chsaa_red$.id==cpg),"coef"],
+ genoa_red[which(genoa_red$.id==cpg),"coef"])
> 
> se_cpg <- c(fhs_red[which(fhs_red$.id==cpg),"se"],
+ aric_red[which(aric_red$.id==cpg),"se"],
+ gtp_red[which(gtp_red$.id==cpg),"se"],
+ inchianti_red[which(inchianti_red$.id==cpg),"se"],
+ kora_red[which(kora_red$.id==cpg),"se"],
+ lbc1921_red[which(lbc1921_red$.id==cpg),"se"],
+ lbc1936_red[which(lbc1936_red$.id==cpg),"se"],
+ nas_red[which(nas_red$.id==cpg),"se"],
+ rs_red[which(rs_red$.id==cpg),"se"],
+ goldn_red[which(goldn_red$.id==cpg),"se"],
+ mesa_red[which(mesa_red$.id==cpg),"se"],
+ epic_red[which(epic_red$.id==cpg),"se"],
+ norfolk_red[which(norfolk_red$.id==cpg),"se"],
+ chsea_red[which(chsea_red$.id==cpg),"se"],
+ chsaa_red[which(chsaa_red$.id==cpg),"se"],
+ genoa_red[which(genoa_red$.id==cpg),"se"])
> 
> coef_se <- data.frame(coef_cpg, se_cpg)
> rownames(coef_se) <- c("fhs","aric","gtp","inchianti","kora","lbc1921","lbc1936","nas","rs","goldn","mesa","epic","norfolk","chsea","chsaa","genoa")
> 
> ## Metal analysis using RMA (robust multiple average) 
> res_fe <- rma(yi=coef_se$coef_cpg, sei=coef_se$se_cpg, method="FE",weighted=TRUE)
> res_dl <- rma(yi=coef_se$coef_cpg, sei=coef_se$se_cpg, method="DL",weighted=TRUE)
> #options(digits=4)
> 
> Mat <- rbind(coef_se, c(res_fe$beta, res_fe$se), c(res_dl$beta, res_dl$se))
> rownames(Mat)[17,18] <- c("Meta_FE", "Meta_DL")