# Meta-analysis for LL vs RR

# Datasets: GSE125943 and GSE74481

library("MetaDE")
library("readxl")
library("org.Hs.eg.db")
library("AnnotationDbi")
library("WriteXLS")

# Reading Expression Matrix from GSE125943

GSE125943 <- read.csv(".../GSE125943/GSE125943_LL_vs_RR.csv", header = TRUE)
head(GSE125943, 2)

# Formatting

table(duplicated(GSE125943$ENTREZID))
row.names(GSE125943) <- as.character(GSE125943$ENTREZID)
GSE125943 <- GSE125943[,c(13:24)]
head(GSE125943, 2)
label_GSE125943 <- ifelse(grepl("LL", colnames(GSE125943)) == TRUE, '1', '0')
label_GSE125943

# Reading Expression Matrix from GSE74481

GSE74481 <- read.csv(".../GSE74481/GSE74481_LL_vs_RR.csv", header = TRUE)
head(GSE74481, 2)

# Formatting 

table(duplicated(GSE74481$ENTREZID))
row.names(GSE74481) <- as.character(GSE74481$ENTREZID)
GSE74481 <- GSE74481[, c(2:18)]
head(GSE74481, 2)
label_GSE74481 <- ifelse(grepl("LL", colnames(GSE74481)) == TRUE, '1', '0') 
label_GSE74481

#Constructing list

i <- Reduce(intersect, x = list(row.names(GSE125943), rownames(GSE74481)))

# exclude genes from ribosomal proteins
ribosome <- read.table("../ribosomal_proteins.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
head(ribosome, 3)
i <- setdiff(i, ribosome$gene_id)

LL_vs_RR <- list(data.matrix(GSE125943[i, ]), data.matrix(GSE74481[i, ]))
rm(i)

names(LL_vs_RR) <- c("GSE125943", "GSE74481")
K <- length(LL_vs_RR)
label <- list(GSE125943 = label_GSE125943, GSE74481 = label_GSE74481)

clin.data <- lapply(label, function(x) {data.frame(x)} )
for (k in 1:length(clin.data)){
        colnames(clin.data[[k]]) <- "label"
}

# meta analysis using REM
meta.res <- MetaDE(data = LL_vs_RR, clin.data = clin.data, data.type = "continuous", resp.type = "twoclass", response = 'label', ind.method = rep('limma', 4), meta.method = "REM", select.group = c('0', '1'), ref.level = c('0'), paired = rep(FALSE, length(LL_vs_RR)), REM.type = "HO", tail = 'abs')

# saving main holder object 

saveRDS(meta.res, file = "obj/LL_vs_RR_meta.rds", compress = "bzip2")

# Constructing a matrix to store results

LL_vs_RR.rem2 <- matrix(data = c(meta.res$meta.analysis$mu.hat, meta.res$meta.analysis$mu.var, meta.res$meta.analysis$tau2, meta.res$meta.analysis$FDR), nrow = length(meta.res$meta.analysis$mu.hat), ncol = 4, dimnames = list(row.names(meta.res$meta.analysis$FDR), c("muhat", "muvar", "tau2", "FDR")), byrow = FALSE)
LL_vs_RR.rem2 <- cbind(meta.res$ind.ES[row.names(LL_vs_RR.rem2), ], LL_vs_RR.rem2)
LL_vs_RR.rem2 <- as.data.frame(LL_vs_RR.rem2)
LL_vs_RR.rem2$Symbol <- select(org.Hs.eg.db, keys = row.names(LL_vs_RR.rem2),columns = c("SYMBOL"), keytype = "ENTREZID")$SYMBOL 
colnames(LL_vs_RR.rem2)[1:4] <- names(LL_vs_RR)
LL_vs_RR.rem2 <- merge(LL_vs_RR.rem2, meta.res$ind.Var, by = "row.names")
colnames(LL_vs_RR.rem2)[11:14] <- paste("var", names(LL_vs_RR), sep = "_")

WriteXLS(LL_vs_RR.rem2, ExcelFileName = "results/LL_vs_RR_ES.xls")

rm(list = setdiff(ls(), c("meta.res", 'draw.DEnumber_all')))

# Generating P-values using limma

p_GSE125943 <- read.csv("GSE125943_LL_vs_RR.csv", header = TRUE)
p_GSE125943 <- p_GSE125943[,c("ENTREZID", "P_value", "adjusted_P_value")]
p_GSE125943 <- p_GSE125943[!is.na(p_GSE125943$ENTREZID),]
row.names(p_GSE125943) <- p_GSE125943$ENTREZID

p_GSE74481 <- as.data.frame(read_xls("../../GSE74481/results/results.LL_vs_RR.xls"))
p_GSE74481 <- p_GSE74481[,c("ENTREZID", "P_value", "adjusted_P_value")]
p_GSE74481 <- p_GSE74481[!is.na(p_GSE74481$ENTREZID),]
row.names(p_GSE74481) <- p_GSE74481$ENTREZID

## Genes common to both the studies
i <- Reduce(intersect, list(row.names(p_GSE125943), row.names(p_GSE74481)))

## exclude genes from ribosomal proteins
ribosome <- read.table("../ribosomal_proteins.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
i <- setdiff(i, ribosome$gene_id)
tmp <- cbind(p_GSE125943[i,]$P_value, p_GSE74481[i,]$P_value)
row.names(tmp) <- i
colnames(tmp) <- c("GSE125943", "GSE74481")
rm(i)

# Storing P_values from both studies
LL_vs_RR_pvalues <- list(p = tmp)
rm(tmp)

# Meta-analysis using other methods

LL_vs_RR<- MetaDE.pvalue(LL_vs_RR_pvalues, meta.method = c("fisher","maxP", "SR"), rth = 3, parametric = FALSE)

#Merging results

res <- cbind(LL_vs_RR$ind.p, LL_vs_RR$meta.analysis$pval)
tmp <- cbind(meta.res$meta.analysis$pval, meta.res$meta.analysis$FDR)
res <- cbind(res, tmp[row.names(res),][,1])
colnames(res)[8] <- "REM" 
rm(tmp)

i <- Reduce(intersect, list(row.names(p_GSE125943), row.names(p_GSE74481)))
tmp <- cbind(p_GSE125943[i,]$adjusted_P_val, p_GSE74481[i,]$adjusted_P_val)
row.names(tmp) <- i
colnames(tmp) <- c("GSE125943","GSE74481")
res <- cbind(tmp[i, ], LL_vs_RR$meta.analysis$FDR[i, ])
res <- cbind(res[i, ], meta.res$meta.analysis$FDR[i,])
colnames(res)[8] <- "REM" 
rm(tmp)
res <- as.data.frame(res)
res$Symbol <- select(org.Hs.eg.db, keys = row.names(res), columns = "SYMBOL", "ENTREZID")$SYMBOL 

# Saving results from other methods
WriteXLS(res, ExcelFileName = ".../LL_vs_RR_allMethods.xls", row.names = TRUE)

# Intersection between all methods 
intersection <- data.frame(ENTREZID = Reduce(intersect, list(sdef, fisher, maxp, sr)))
intersection$SYMBOL <- select(org.Hs.eg.db, keys = as.character(intersection$ENTREZID), columns = "SYMBOL", "ENTREZID")$SYMBOL
WriteXLS(intersection, ExcelFileName = ".../Intersection.xls")