---
title: "Applying GMEA to meta-regression results"
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

Import meta-regression results to run mitch.

## Requirements

**IMPORTANT** Remove the mitch package and reinstall with devtools directly from GitHub.

```
devtools::install_github("markziemann/mitch")
```

Load packages.

```{r,packages}

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

CORES=8

```

## Load data

Reactome pathways were downloaded on the 11th January.

```{r,load1}

if(!file.exists("ReactomePathways.gmt")) {
  download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
    destfile="ReactomePathways.gmt.zip")
  unzip("ReactomePathways.gmt.zip")
}

gs <- gmt_import("ReactomePathways.gmt")
head(gs)

```

## Curate the annotation

Use all probes on the chip.

Updating gene names using the HGNC package because a lot of them have been updated over time.

```{r,anno1}

tic()
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt <- stack(gp2)
colnames(gt) <- c("gene","probe")
gt$probe <- as.character(gt$probe)
dim(gt)
str(gt)
toc() #9.0s

tic()
#new.hgnc.table <- getCurrentHumanMap() Dec 2023
if(!file.exists("new.hgnc.table.rds")) {
  download.file("https://github.com/markziemann/gmea/raw/main/figs/new.hgnc.table.rds",
    destfile="new.hgnc.table.rds")
}
new.hgnc.table <- readRDS("new.hgnc.table.rds")
fix <- checkGeneSymbols(gt$gene,map=new.hgnc.table)
fix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))
gt$gene <- fix$Suggested.Symbol
toc()

```

## Load meta-regression data

```{r,loaddata}

x <- readRDS("models_sex_all_blood_need_cleaning>=10_5000.rds")
df <- as.data.frame(as.matrix(unlist(lapply(x,function(y) { y$model$zval[2] }))))
colnames(df) <- "zval"
hist(df$zval)
which(is.na(df$zval))
rm(x)

```

## Mitch import

Ensure that mitch version has aggregate by mean, not sum.
Delete the bioconductor version and install with devtools::install_github

```{r,mg1}

# mitch aggregate function needs to be changes to `mean`
tic()
m <- mitch_import(x=df,DEtype="prescored",geneTable=gt)
toc() # 7.25 sec

hist(m$x)
dim(m)

head(m) %>%
  kbl(caption = "Example of differential gene methylation scores used for pathway analysis") %>%
  kable_paper("hover", full_width = F)

```

## Run mitch

```{r,mitch1}

genenames <- unique(unname(unlist(gs)))
length(which(rownames(m) %in% genenames))
length(which(!rownames(m) %in% genenames))

tic()
mres <- mitch_calc(x=m,genesets=gs,minsetsize=1,cores=16, priority="significance")
toc() #82.639 sec

mtable <- mres$enrichment_result
up <- subset(mtable,s.dist>0 & p.adjustANOVA<0.05)
dn <- subset(mtable,s.dist<0 & p.adjustANOVA<0.05)
nrow(up)
nrow(dn)

head(up,10) %>%
  kbl(caption = "Top significant pathways with higher methylation") %>%
  kable_paper("hover", full_width = F)

head(dn,10) %>%
  kbl(caption = "Top significant pathways with lower methylation") %>%
  kable_paper("hover", full_width = F)

top <- rbind(up[1:10,],dn[1:10,])
top <- top[order(top$s.dist),]

barnames <- gsub("_"," ",gsub("REACTOME_","",top$set))
cols <- as.character(((sign(top$s.dist)+1)/2)+1)
cols <- gsub("1","blue",gsub("2","red",cols))

par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(abs(top$s.dist),horiz=TRUE,las=1,names.arg=barnames, cex.names=0.8,
  cex.axis=0.8,col=cols, xlab="Enrichment score",main="Reactomes")
mtext("whole gene")
grid()

par( mar = c(5.1, 4.1, 4.1, 2.1) )

```

Make a html report and some charts.

```{r,report1}

mitch_report(res=mres,outfile="mandy_mitchreport.html",overwrite=TRUE)
mitch_plots(res=mres,outfile="mandy_mitchcharts.pdf")

```

## Session information

```{r,save}

save.image("mandy_gmea.Rdata")

sessionInfo()

```
