forest.R

setwd("/Users/zd821/Documents/project2/TCGA")

library(forestplot)
library(dplyr)
# Cochrane data from the 'rmeta'-package
cochrane_from_rmeta <- structure(list(mean  = c(1.354, 0.358, 1.357, 0.722, 1.536, 1.293, 0.439, 0.768, 1.315, 1.247, 1.377, 1.631, 1.228, 1.293,  2.792, 0.142, 4.162, 0.207, 31.640, 4.143, 0.205, 0.325, 3.559, 10.730, 1.802, 15.150, 1.941, 7.081), 
                                      lower = c(1.145, 0.203, 1.138, 0.597, 1.187, 1.106, 0.262, 0.646, 1.089, 1.070, 1.089, 1.126, 1.026, 1.005,  1.783, 0.046, 1.666, 0.085, 1.450, 2.048, 0.085, 0.201, 1.713, 3.625, 1.220, 2.990, 1.326, 1.839),
                                      upper = c(1.602, 0.634, 1.619, 0.873, 1.989, 1.511, 0.737, 0.914, 1.587, 1.454, 1.741, 2.363, 1.470, 1.663,  4.371, 0.431, 10.396, 0.507, 690.379, 8.379, 0.491, 0.525, 7.392, 31.755, 2.662, 76.759, 2.844, 27.271),
                                      group = c(rep("TCGA-PRAD",14,),rep("MSKCC",14))),
                                 .Names = c("mean", "lower", "upper","group"), 
                                 row.names = c(1:28), 
                                 class = "data.frame")

tabletext <- cbind(Gene = c("COL1A1","FREM2","APLP1","DPT","MATN3","THBS2","ANGPT1","FBLN1","MDK","COL10A1","LRRC15","PDGFB","VCAN","LAMC3"), 
                  "HR(MSKCC)" = c("2.792", "0.142", "4.162", "0.207", "31.640", "4.143", "0.205", "0.325", "3.559", "10.730", "1.802", "15.150", "1.941", "7.081"), 
                  "p-value" = c("< 0.001", "< 0.001", "0.002", "< 0.001", "0.028", "< 0.001", "< 0.001", "< 0.001", "< 0.001", "< 0.001", "0.003", "0.001", "< 0.001", "0.004"), 
                  "HR(TCGA-PRAD)" = c("1.354", "0.358", "1.357", "0.722", "1.536", "1.293", "0.439", "0.768", "1.315", "1.247", "1.377", "1.631", "1.228", "1.293"), 
                  "p-value" = c("< 0.001", "< 0.001", "< 0.001", "< 0.001", "0.001", "0.001", "0.002", "0.003", "0.004", "0.005", "0.008", "0.010", "0.025", "0.046"))


pdf("forest_plot.pdf")

cochrane_from_rmeta %>% group_by(group) %>%
  forestplot(legend = c("MSKCC", "TCGA-PRAD"), fn.ci_norm = fpDrawDiamondCI, boxsize = .25,ci.vertices = TRUE,
             ci.vertices.height = 0.1,line.margin = .2, labeltext = tabletext, clip=c(0.05, 50), xticks = c(0.05,0.1,0.2,0.5,1,2,5,10,20,50),xlog = TRUE, col = fpColors(box = c("#5070D4", "#D76045")))

dev.off()