---
title: "GSE158422: DNA methylation in Lung squamous cell carcinoma - gene methylation enrichment analysis"
author: "The GMEA team"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    fig_width: 7
    fig_height: 7
theme: cosmo
---

## Introduction

Here we are splitting the dataset and testing concordance between LA and AL approaches.

```{r,packages}

suppressPackageStartupMessages({
  library("limma")
  library("eulerr")
  library("IlluminaHumanMethylation450kmanifest")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  library("tictoc")
  library("mitch")
  library("kableExtra")
  library("beeswarm")
})

CORES=28

```

## Load data

* annotations

* probe sets

* gene sets

* design matrix

* mval matrix

```{r,annotation}

anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","Regulatory_Feature_Group","Islands_Name","Relation_to_Island")])

gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
sets <- split(rep(names(gp2), lengths(gp2)), unlist(gp2))

summary(unlist(lapply(sets,length)))

genesets <- gmt_import("https://ziemann-lab.net/public/gmea_prototype/ReactomePathways.gmt")

if (!file.exists("GSE158422_design.rds")) {
  download.file("https://ziemann-lab.net/public/gmea_prototype/GSE158422_design.rds", "GSE158422_design.rds")
}
design <- readRDS("GSE158422_design.rds")

if (!file.exists("GSE158422_design.rds")) {
 download.file("https://ziemann-lab.net/public/gmea_prototype/GSE158422_mx.rds","GSE158422_mx.rds")
}
mval <- readRDS("GSE158422_mx.rds")

boxplot(list("normal"=matrix(colMeans(mval),ncol=2)[,2],"tumor"=matrix(colMeans(mval),ncol=2)[,1]),
  main="mean probe methylation mval")

```

Reactome pathways were downloaded on the 14th Sept 2023 from MsigDB.

```{r,gsdl}

gs_entrez <- gmt_import("c2.cp.reactome.v2023.1.Hs.entrez.gmt")

gs_symbols <- gmt_import("c2.cp.reactome.v2023.1.Hs.symbols.gmt")

```

# GMEA FUNCTIONS

There are three approaches.
First, limma can be done first, followed by aggregation by gene and then enrichment
test (called LA).
Secondly, probe values can be aggregated by gene, followed by limma, followed by
enrichment test (AL).
Third, probe values can be aggregated to gene level, and then aggregated to gene set
level followed by limma test (AA).

## LA

This process runs limma first and then aggregates the results before doing an enrichment test.

```{r,lafunc}

# limma
runlimma <- function(mval,design,myann) {
  fit.reduced <- lmFit(mval,design)
  fit.reduced <- eBayes(fit.reduced)
  summary(decideTests(fit.reduced))
  dm <- topTable(fit.reduced,coef=4, number = Inf)
  dm <- merge(myann,dm,by=0)
  dm <- dm[order(dm$P.Value),]
  rownames(dm) <- dm$Row.names
  dm$Row.names=NULL
  return(dm)
}

# aggregate
agg <- function(dm,cores=1) {
  gn <- unique(unlist(strsplit( dm$UCSC_RefGene_Name ,";")))
  gnl <- strsplit( dm$UCSC_RefGene_Name ,";")
  gnl <- mclapply(gnl,unique,mc.cores=cores)
  dm$UCSC_RefGene_Name <- gnl
  l <- mclapply(1:nrow(dm), function(i) {
    a <- dm[i,]
    len <- length(a[[1]][[1]])
    tvals <- as.numeric(rep(a["t"],len))
    genes <- a[[1]][[1]]
    data.frame(genes,tvals)
  },mc.cores=cores)
  df <- do.call(rbind,l)
  keep <- names(which(table(df$genes)>1))
  df <- df[df$genes %in% keep,]
  gn <- unique(df$genes)
  gme_res <- lapply( 1:length(gn), function(i) {
    g <- gn[i]
    tstats <- df[which(df$genes==g),"tvals"]
    myn <- length(tstats)
    mymean <- mean(tstats)
    mymedian <- median(tstats)
    if ( length(tstats) > 2 ) {
      ttest <- t.test(tstats)
      pval <- ttest$p.value
    } else {
      pval = 1
    }
    res <- c("gene"=g,"nprobes"=myn,"mean"=mymean,
      "median"=mymedian, pval=pval)
  } )
  gme_res_df <- do.call(rbind, gme_res)
  rownames(gme_res_df) <- gme_res_df[,1]
  gme_res_df <- gme_res_df[,-1]
  tmp <- apply(gme_res_df,2,as.numeric)
  rownames(tmp) <- rownames(gme_res_df)
  gme_res_df <- as.data.frame(tmp)
  gme_res_df$sig <- -log10(gme_res_df[,4])
  gme_res_df <- gme_res_df[order(-gme_res_df$sig),]
  gme_res_df$fdr <- p.adjust(gme_res_df$pval)
  out <- list("df"=df,"gme_res_df"=gme_res_df)
  return(out)
}

# enrich
ttenrich <- function(m,genesets,cores=1) {
  res <- mclapply( 1:length(genesets), function(i) {
    gs <- genesets[i]
    name <- names(gs)
    n_members <- length(which(rownames(m) %in% gs[[1]]))
    if ( n_members > 4 ) {
      tstats <- m[which(rownames(m) %in% gs[[1]]),]
      myn <- length(tstats)
      mymean <- mean(tstats)
      mymedian <- median(tstats)
      wt <- t.test(tstats)
      res <- c(name,myn,mymean,mymedian,wt$p.value)
    }
  } , mc.cores = cores)
  res_df <- do.call(rbind, res)
  rownames(res_df) <- res_df[,1]
  res_df <- res_df[,-1]
  colnames(res_df) <- c("n_genes","t_mean","t_median","pval")
  tmp <- apply(res_df,2,as.numeric)
  rownames(tmp) <- rownames(res_df)
  res_df <- tmp
  res_df <- as.data.frame(res_df)
  res_df <- res_df[order(res_df$pval),]
  res_df$logp <- -log10(res_df$pval )
  res_df$fdr <- p.adjust(res_df$pval,method="fdr")
  res_df[order(abs(res_df$pval)),]
  return(res_df)
}

```

## AL: Aggregate limma functions

Functions for aggregate-limma-enrich approach.

```{r,alfunc}

# chromosome by chromosome will be much faster
magg <- function(mval,myann,cores=1){
  gn <- unique(unlist(strsplit( myann$UCSC_RefGene_Name ,";")))
  gnl <- strsplit( myann$UCSC_RefGene_Name ,";")
  gnl <- mclapply(gnl,unique,mc.cores=cores)
  myann$gnl <- gnl
  keep <- rownames(subset(myann,UCSC_RefGene_Name!=""))
  mx <- mval[rownames(mval) %in% keep,]
  mymed <- function(g) {
    probes <- rownames(myann[grep(g,myann$gnl),])
    rows <- which(rownames(mx) %in% probes)
    if ( length(rows) > 1 ) {
      b <- mx[rows,]
      med <- apply(b,2,median)
      med <- matrix(med,nrow=1)
      colnames(med) <- colnames(b)
      rownames(med) <- g
      return(med)
    }
  }
  med <- mclapply(gn,mymed,mc.cores=cores)
  med <- med[lapply(med,length)>0]
  medf <- do.call(rbind,med)
  return(medf)
}

chragg <- function(mval,myann,cores=1){
  annodf <- as.data.frame(anno)
  keep <- rownames(subset(myann,UCSC_RefGene_Name!=""))
  mx <- mval[rownames(mval) %in% keep,]
  chrs <- unique(anno$chr)
  myorder <- unlist(lapply(chrs,function(mychr) { nrow( annodf[annodf$chr==mychr,] ) } ))
  chrs <- chrs[order(-myorder)]
  leadercores <- floor(sqrt(cores))
  workercores <- ceiling(sqrt(cores))
  chrmedf <- mclapply(chrs,function(chr) {
    chrfrag <- annodf[annodf$chr==chr,]
    chrprobes <-rownames(chrfrag)
    chrmx <- mx[rownames(mx) %in% chrprobes,]
    chranno <- myann[rownames(myann) %in% chrprobes,]
    chrmedf <- magg(mval=chrmx,myann=chranno,cores=workercores)
    return(chrmedf)
  },mc.cores=leadercores)
  medf <- do.call(rbind, chrmedf)
  return(medf)
}

agglimma <- function(medf,design) {
  fit.reduced <- lmFit(medf,design)
  fit.reduced <- eBayes(fit.reduced)
  dmagg <- topTable(fit.reduced,coef=ncol(design), number = Inf)
  nondup <- !duplicated(dmagg$ID)
  dmagg <- dmagg[nondup,]
  rownames(dmagg) <- dmagg$ID
  dmagg$ID = NULL
  return(dmagg)
}

```

## AA Aggregate-aggregate-limma functions

```{r,aafunc}

gsagg <- function(x,genesets,cores=1) {
  meds <- mclapply(1:length(genesets), function(i) {
    gs = genesets[[i]]
    xx <- x[rownames(x) %in% gs,]
    med <- apply(xx,2,median)
  },mc.cores=cores)
  mymed <- do.call(rbind,meds)
  rownames(mymed) <- names(genesets)
  as.data.frame(mymed)
}

aalimma <- function(agag,design) {
  fit.reduced <- lmFit(agag,design)
  fit.reduced <- eBayes(fit.reduced)
  dmagg <- topTable(fit.reduced,coef=ncol(design), number = Inf)
  return(dmagg)
}

aal <- function(mval,myann,genesets,design,cores=1) {
  medf <- chragg(mval,myann,cores=CORES)
  agag <- gsagg(x=medf,genesets=gs_symbols,cores=CORES)
  aal <- aalimma(agag=agag,design=design)
  return(aal)
}

aalres <- aal(mval=mval,myann=myann,genesets=gs_symbols,design=design,cores=CORES)

```

## Prep

```{r,prep}

# prep
sex <- as.data.frame(design)$sex
tumor <- as.data.frame(design)$tumor
patient <- as.character(unlist(lapply(1:ncol(mval),function(i) {c(i,i)})))
patient <- head(patient,ncol(mval))
design <- model.matrix(~ patient + tumor )
rownames(design) <- colnames(mval)

```

## LA workflow split and analyse

split 74 samples 37 patients split into 18 + 18
Conduct LA workflow analysis.
Then make overlap.

```{r,laflow}

splitanalyze_la <- function(seed){
  set.seed(seed)
  a <- sample(1:37,18)
  a <- a[order(a)]
  nona <- setdiff(1:37,a)
  set.seed(seed)
  b <- sample(nona,18)
  b <- b[order(b)]
  a <- c( (a*2)-1 ,a*2 )
  a <- a[order(a)]
  b <- c( (b*2)-1 ,b*2 )
  b <- b[order(b)]

  design_a <- design[a,]
  design_a <- design_a[,colSums(design_a)>0]
  mval_a <- mval[,a]
  design_b <- design[b,]
  design_b <- design_b[,colSums(design_b)>0]
  mval_b <- mval[,b]

  # la1
  dm1 <- runlimma(mval_a,design_a,myann)
  dmagg1 <- agg(dm1,cores=CORES)
  m1 <- dmagg1$gme_res_df[,"mean",drop=FALSE]
  lares1 <- ttenrich(m=m1,genesets=gs_symbols,cores=CORES)

  # la2
  dm2 <- runlimma(mval_b,design_b,myann)
  dmagg2 <- agg(dm2,cores=CORES)
  m2 <- dmagg2$gme_res_df[,"mean",drop=FALSE]
  lares2 <- ttenrich(m=m2,genesets=gs_symbols,cores=CORES)

  list("lares1"=lares1,"lares2"=lares2)

}

lares <- lapply( seq(100,5000,100) ,splitanalyze_la)

```

## AL workflow split and analyse

```{r,alflow}

splitanalyze_al <- function(seed){
  set.seed(seed)
  a <- sample(1:37,18)
  a <- a[order(a)]
  nona <- setdiff(1:37,a)
  set.seed(seed)
  b <- sample(nona,18)
  b <- b[order(b)]
  a <- c( (a*2)-1 ,a*2 )
  a <- a[order(a)]
  b <- c( (b*2)-1 ,b*2 )
  b <- b[order(b)]

  design_a <- design[a,]
  design_a <- design_a[,colSums(design_a)>0]
  mval_a <- mval[,a]
  design_b <- design[b,]
  design_b <- design_b[,colSums(design_b)>0]
  mval_b <- mval[,b]

  #al
  medf1 <- chragg(mval_a,myann,cores=CORES)
  magg1 <- agglimma(medf1,design_a)
  m1 <- as.data.frame(magg1$t)
  rownames(m1) <- rownames(magg1)
  colnames(m1) <- "t"
  alres1 <- ttenrich(m=m1,genesets=gs_symbols,cores=CORES)

  #al
  medf2 <- chragg(mval_b,myann,cores=CORES)
  magg2 <- agglimma(medf2,design_b)
  m2 <- as.data.frame(magg2$t)
  rownames(m2) <- rownames(magg2)
  colnames(m2) <- "t"
  alres2 <- ttenrich(m=m2,genesets=gs_symbols,cores=CORES)

  list("alres1"=alres1,"alres2"=alres2)
}

alres <- lapply( seq(100,5000,100) ,splitanalyze_al)

```

## AA workflow split and analyse

```{r,aaflow}

splitanalyze_aa <- function(seed){
  set.seed(seed)
  a <- sample(1:37,18)
  a <- a[order(a)]
  nona <- setdiff(1:37,a)
  set.seed(seed)
  b <- sample(nona,18)
  b <- b[order(b)]
  a <- c( (a*2)-1 ,a*2 )
  a <- a[order(a)]
  b <- c( (b*2)-1 ,b*2 )
  b <- b[order(b)]

  design_a <- design[a,]
  design_a <- design_a[,colSums(design_a)>0]
  mval_a <- mval[,a]
  design_b <- design[b,]
  design_b <- design_b[,colSums(design_b)>0]
  mval_b <- mval[,b]

  #aa1
  aares1 <- aal(mval=mval_a,myann=myann,genesets=gs_symbols,
    design=design_a,cores=CORES)

  #aa2
  aares2 <- aal(mval=mval_b,myann=myann,genesets=gs_symbols,
    design=design_b,cores=CORES)

  list("aares1"=aares1,"aares2"=aares2)
}

aares <- lapply( seq(100,5000,100) ,splitanalyze_aa)

```

## Compare

```{r,summarise1}

lacompare <- function(lares){
  lares1 <- lares[[1]]
  lares2 <- lares[[2]]
  lasig1 <- subset(lares1,fdr<0.05)
  lasig2 <- subset(lares2,fdr<0.05)
  lasig1_up <- rownames(subset(lasig1,t_median>0))
  lasig1_dn <- rownames(subset(lasig1,t_median<0))
  lasig2_up <- rownames(subset(lasig2,t_median>0))
  lasig2_dn <- rownames(subset(lasig2,t_median<0))

  total <- length(unique(c(lasig1_up,lasig2_up,lasig1_dn,lasig2_dn)))
  common <- length(intersect(lasig1_up, lasig2_up)) + length(intersect(lasig1_dn, lasig2_dn))
  discordant <- length(intersect(lasig1_up, lasig2_dn)) + length(intersect(lasig1_dn, lasig2_up))
  uncommon = total - common - discordant

  laout <- data.frame(total,common,discordant,uncommon)
  laout$p_comm <-laout$common / laout$total
  laout$p_disc <-laout$discordant / laout$total
  laout$p_uncom <-laout$uncommon / laout$total
  laout
}

alcompare <-function(alres){
  alres1 <- alres[[1]]
  alres2 <- alres[[2]]
  alsig1 <- subset(alres1,fdr<0.05)
  alsig2 <- subset(alres2,fdr<0.05)

  alsig1_up <- rownames(subset(alsig1,t_median>0))
  alsig1_dn <- rownames(subset(alsig1,t_median<0))
  alsig2_up <- rownames(subset(alsig2,t_median>0))
  alsig2_dn <- rownames(subset(alsig2,t_median<0))

  total <- length(unique(c(alsig1_up,alsig2_up,alsig1_dn,alsig2_dn)))
  common <- length(intersect(alsig1_up, alsig2_up)) + length(intersect(alsig1_dn, alsig2_dn))
  discordant <- length(intersect(alsig1_up, alsig2_dn)) + length(intersect(alsig1_dn, alsig2_up))
  uncommon = total - common - discordant
  result <- c("total"=total,"common"=common,"discordant"=discordant,"uncommon"=uncommon)

  alout <- data.frame(total,common,discordant,uncommon)
  alout$p_comm <-alout$common / alout$total
  alout$p_disc <-alout$discordant / alout$total
  alout$p_uncom <-alout$uncommon / alout$total
  alout
}

aacompare <-function(aares){
  aares1 <- aares[[1]]
  aares2 <- aares[[2]]
  aasig1 <- subset(aares1,adj.P.Val<0.05)
  aasig2 <- subset(aares2,adj.P.Val<0.05)

  aasig1_up <- rownames(subset(aasig1,logFC>0))
  aasig1_dn <- rownames(subset(aasig1,logFC<0))
  aasig2_up <- rownames(subset(aasig2,logFC>0))
  aasig2_dn <- rownames(subset(aasig2,logFC<0))

  total <- length(unique(c(aasig1_up,aasig2_up,aasig1_dn,aasig2_dn)))
  common <- length(intersect(aasig1_up, aasig2_up)) + length(intersect(aasig1_dn, aasig2_dn))
  discordant <- length(intersect(aasig1_up, aasig2_dn)) + length(intersect(aasig1_dn, aasig2_up))
  uncommon = total - common - discordant
  result <- c("total"=total,"common"=common,"discordant"=discordant,"uncommon"=uncommon)

  aaout <- data.frame(total,common,discordant,uncommon)
  aaout$p_comm <- aaout$common / aaout$total
  aaout$p_disc <- aaout$discordant / aaout$total
  aaout$p_uncom <-aaout$uncommon / aaout$total
  aaout
}

```

```{r,compare1}

laout <- do.call(rbind,lapply(lares,lacompare))

alout <- do.call(rbind,lapply(alres,alcompare))

aaout <- do.call(rbind,lapply(aares,aacompare))

laout

alout

aaout

```

## Session information

```{r,save}

save.image("GSE158422_split.Rdata")

sessionInfo()

```
