Choufani2019.Rmd

---
title: "Choufani Reanalysis"
author: "The methylation in ART meta-analysis group"
date: "6/10/2020"
output: 
  html_document:
    toc: true
    fig_width: 7
    fig_height: 7
theme: cosmo
---


# Introduction
In this report, I will take you through a re-analysis methylation data first described by 
[Choufani et al (2019)](https://academic.oup.com/hmg/article/28/3/372/5102486).

In their study, they analysed the DNA methylation patterns of 44 placenta samples conceived
naturally (NAT) and 44 placenta samples conceivced through artificial reproductive technologies (ART) (23 in vitro, 21 in vivo).

The platform used in the study is the Illumina Infinium HumanMethylation450k BeadChip assay.

The methylation data have been deposited to NCBI GEO repository accession number 
[GSE120250](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120250).

# Loading packages and functions
These packackes will help us to perform vital steps such as normalisation, filtering, 
differential analysis, etc, and provide information about the array probe annotaions.

These functions provide shortcuts to help with charts and other analysis. They are 
now available in the separate document "meth_functions.R"

```{r,packages}
source("meth_functions.R")
# Annotation
ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
myann <- data.frame(ann450k[,c("UCSC_RefGene_Name","Regulatory_Feature_Group")])
promoters <- grep("Prom",myann$Regulatory_Feature_Group)
```


# Data Import
Importing Data from GEO.
Raw IDAT files and accompanying Series Matrix publicly available on GEO.
```{r,data}
#Create Dirtectory
dir.create("GSE120250")
ARRAY_SAMPLESHEET="GSE120250/GSE120250_series_matrix.txt.gz"
# only download it if it is not present on the system
if ( !file.exists(ARRAY_SAMPLESHEET ) ) {
    DLFILE=paste(ARRAY_SAMPLESHEET,".gz",sep="")
    download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE120nnn/GSE120250/matrix/GSE120250_series_matrix.txt.gz",
        destfile = DLFILE)
    gunzip(DLFILE)
}



# Reading in GEO series matrix
gse <- getGEO(filename="/mnt/mnorris/ART_methylation/GSE120250/GSE120250_series_matrix.txt.gz")


ARRAY_DATA="GSE120250/GSE120250_RAW.tar"
# only download it if it is not present on the system
if ( !dir.exists("GSE120250/IDAT") ) {
  dir.create("GSE120250/IDAT")
  download.file("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE120250&format=file",
    destfile = ARRAY_DATA)
    untar(exdir = "GSE120250/IDAT", tarfile = "GSE120250/GSE120250_RAW.tar",)
}

baseDir <- "./GSE120250"
targets <- pData(phenoData(gse))


#Include Extra Metadata
extrameta <- read.table("ARTdetails.tsv", sep = "\t", header = TRUE)
extrameta <- extrameta[,c("title","Details.ART")]
targets <- merge(extrameta, targets)
rownames(targets) <- targets$geo_accession

#Reorder dataframe

targets <- targets[order(rownames(targets)),]
mybase <- unique(gsub("_Red.idat.gz" ,"", gsub("_Grn.idat.gz", "" ,list.files("./GSE120250",pattern = "GSM",recursive = TRUE))))
targets$Basename <- paste("GSE120250/", mybase, sep = "")
rgSet <- read.metharray.exp(targets = targets)
```




# Data Normalisation
```{r,norm,fig.cap="Figure 1. Normalisation of bead-array data with SWAN."}
mSet <- preprocessRaw(rgSet)
mSetSw <- SWAN(mSet,verbose=FALSE) 
par(mfrow=c(1,2))
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")
```
## Filter probes

Here we are running parallel analyses, both including and excluding sex chromosomes.

```{r,filterprobes}
# include sex chromosomes
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetSw <- mSetSw[keep,]

# exclude SNP probes
mSetSw <- mapToGenome(mSetSw)
mSetSw_nosnp <- dropLociWithSnps(mSetSw)
mSetSw <- mSetSw_nosnp

# exclude sex chromosomes
keep <- !(featureNames(mSetSw) %in% ann450k$Name[ann450k$chr %in% c("chrX","chrY")])
mSetFlt <- mSetSw[keep,]
```

## Extracting Beta and M-values

```{r,beta_m_vals}
# include sex chromosomes
meth <- getMeth(mSetSw)
unmeth <- getUnmeth(mSetSw)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mSetSw)

# exclude sex chromosomes
meth <- getMeth(mSetFlt)
unmeth <- getUnmeth(mSetFlt)
Mval_flt <- log2((meth + 100)/(unmeth + 100))
beta_flt <- getBeta(mSetFlt)
```


## Removal of Outliers
```{r, outliers}
targets <- subset(targets,`outlier:ch1`=="N")
colnames(Mval) <- sapply(strsplit(colnames(Mval), "_"), "[[", 1)
colnames(Mval_flt) <- sapply(strsplit(colnames(Mval_flt), "_"), "[[", 1)
Mval <- Mval[,which(colnames(Mval) %in% rownames(targets))]
Mval_flt <- Mval_flt[,which(colnames(Mval_flt) %in% rownames(targets))]

```

# MDS analysis

[Multidimensional scaling(https://en.wikipedia.org/wiki/Multidimensional_scaling)  plot is a method used
to identify the major sources of variation in a dataset. In the MDS plots below, I will be plotting the 
first two dimensions (principal components [PCs]), with each sample label coloured either by ART 
classification, sex, ART and sex, and then array chip and then sample plate. 

We wil begin with MDS analysis including the sex chromosomes and then exclude them.

First, let's quantify the contribution of the major principal components. with a scree plot, we can see 
whether most of the variation is captured in the first two PCs or whether it is spread over more PCs.

## Scree Plot
```{r,scree,fig.cap="Figure 2. Screeplot shows contribution of top principal components to the overall variation in the experiment."}

par(mfrow=c(2,1))  
myscree(Mval,main="incl sex chr")
myscree(Mval_flt,main="excl sex chr")

```


## MDS by ART
```{r,mds1,fig.width = 8 ,fig.height = 8,fig.cap="Figure 3. MDS plot coloured by ART classification."}

sample_groups <- factor(targets$characteristics_ch1.2)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$art))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by ART type")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
mydist <- plotMDS(Mval, labels=rownames(targets),col=colours,main="sex chromosomes included")
mydist_flt <- plotMDS(Mval_flt, labels=rownames(targets),col=colours,main="sex chromosomes excluded")

```
## MDS by ART with new metadata
```{r,mds2,fig.width = 8 ,fig.height = 8,fig.cap="Figure 3. MDS plot coloured by ART classification."}

sample_groups <- factor(targets$Details.ART)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$Details.ART))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by ART type")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
plotMDS(mydist, labels=rownames(targets),col=colours,main="sex chromosomes included")
plotMDS(mydist_flt, labels=rownames(targets) ,col=colours,main="sex chromosomes excluded")

```

## MDS by Sex
```{r,mds3,fig.width = 8 ,fig.height = 8, fig.cap="Figure 4. MDS plot coloured by sex."}
sample_groups <- factor(targets$characteristics_ch1.1)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$characteristics_ch1.1))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by sex")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
plotMDS(mydist, labels=rownames(targets),col=colours,main="sex chromosomes included")
plotMDS(mydist_flt, labels=rownames(targets) ,col=colours,main="sex chromosomes excluded")
```

## MDS by ART and sex
```{r,mds4,fig.width = 8 ,fig.height = 8, fig.cap="Figure 5. MDS plot coloured by ART classification and sex."}
targets$ARTgender <- paste(targets$`art treatment:ch1`,targets$`gender:ch1`, sep=("_"))

sample_groups <- factor(targets$ARTgender)

colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$ARTgender))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by ART and sex")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
plotMDS(mydist, labels=rownames(targets), col=colours,main="sex chromosomes included")
plotMDS(mydist_flt, labels=rownames(targets), col=colours,main="sex chromosomes excluded")

```


# Differential analysis

There are a few differential contrasts that would be of interest to us in this study:

* __NAT(40) vs in vivo(15)__: 

* __NAT(40) vs in vitro(18) __: .

* __NAT(40) vs ART(33)__:

* __in vivo(15) vs in vitro(18)__: 

* __NAT(40) vs ICSI(14)__: 

* __NAT(40) vs IUI(9)__:

* __ICSI(14) vs IUI(9)__: 

* __NAT vs IVF__: 

* __NAT(40) vs HYPO(6)__:

The differential analysis is centred around limma to identify differentially methylated probes. 
TopConfects was also run to obtain the probes with the largest confident effect (topconfect).
There are five outputs below:

Potentially Remove in vitro and in vivo comparisons

1. Volcano plot (limma result).

2. Beeswarm plot (top probes by limma p-value).

3. Heatmap (top probes by limma p-value).

4. Beeswarm plot (top probes by topconfect ranking).

5. Heatmap (top probes by topconfect ranking). (TODO)

## Nat vs In vivo

{r,nat_vs_in_vivo_incl,fig.width = 8 ,fig.height = 8, fig.cap="Figure 6: Differential analysis: Naturally conceived vs in vivo including sex chromosomes"}

# Include Sex Chromosomes
samplesheet <- subset(targets,`art treatment:ch1`=="NA" | `art treatment:ch1`=="in vivo")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$`art treatment:ch1`,levels=c("NA","in vivo"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_NA_vs_invivo_inc <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval,name="top_NA_vs_invivo_inc",
    myann=myann ,beta= beta)





```{r,nat_vs_in_vivo_exc, fig.width = 8 , fig.height = 8, fig.cap="Figure 7: Differential analysis: Naturally conceived vs in vivo excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- subset(targets,`art treatment:ch1`=="NA" | `art treatment:ch1`=="in vivo")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$`art treatment:ch1`,levels=c("NA","in vivo"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_NA_vs_invivo_exc <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_NA_vs_invivo_exc",
    myann=myann ,beta= beta)


# Allele
top_NA_vs_invivo_exc$dma$unmeth <- "T"
top_NA_vs_invivo_exc$dma$meth <- "C"
top_NA_vs_invivo_exc$fit$SE <- sqrt(top_NA_vs_invivo_exc$fit$s2.post) * top_NA_vs_invivo_exc$fit$stdev.unscaled


# Extract required columns from dma
top_NA_vs_invivo_metal <-top_NA_vs_invivo_exc$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_NA_vs_invivo_exc_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_NA_vs_invivo_exc$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsin vivo")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_NA_vs_invivo_exc$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsin vivo")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_NA_vs_invivo_metal <- merge(top_NA_vs_invivo_exc_metal, fitCE)
top_NA_vs_invivo_exc_metal <- merge(top_NA_vs_invivo_exc_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$`art treatment:ch1` == "NA",]
cases <- samplesheet[samplesheet$`art treatment:ch1` == "in vivo",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))
top_NA_vs_invivo_metal$N <- Neff


# Outputs for html
head(top_NA_vs_invivo_metal)
head(top_NA_vs_invivo_exc$dma)
#head(top_NA_vs_invivo__exc$dmr)
#top_NA_vs_invivo_exc$comp
#top_NA_vs_invivo_exc$cgi
#names(top_NA_vs_invivo_exc)
head(top_NA_vs_invivo_exc$confects$table$name)

# Output for Meta-analysis
dir.create("misc")
write.table(top_NA_vs_invivo_metal,file="misc/choufani_top_NA_vs_invivo_exc.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```

## Nat vs In vitro

```{r,nat_vs_in_vitro_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 9: Differential analysis: Naturally conceived vs in vitro excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- subset(targets,`art treatment:ch1`=="NA" | `art treatment:ch1`=="in vitro")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$`art treatment:ch1`,levels=c("NA","in vitro"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_NA_vs_invitro_exc <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_NA_vs_invitro_exc",
    myann=myann ,beta= beta)

# Allele
top_NA_vs_invitro_exc$dma$unmeth <- "T"
top_NA_vs_invitro_exc$dma$meth <- "C"
top_NA_vs_invitro_exc$fit$SE <- sqrt(top_NA_vs_invitro_exc$fit$s2.post) * top_NA_vs_invitro_exc$fit$stdev.unscaled


# Extract required columns from dma
top_NA_vs_invitro_metal <-top_NA_vs_invitro_exc$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_NA_vs_invitro_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_NA_vs_invitro_exc$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsin vitro")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_NA_vs_invitro_exc$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsin vitro")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_NA_vs_invitro_metal <- merge(top_NA_vs_invitro_metal, fitCE)
top_NA_vs_invitro_metal <- merge(top_NA_vs_invitro_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$`art treatment:ch1` == "NA",]
cases <- samplesheet[samplesheet$`art treatment:ch1` == "in vitro",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))
top_NA_vs_invivo_metal$N <- Neff


# Outputs for html
head(top_NA_vs_invitro_metal)
head(top_NA_vs_invitro_exc$dma)
#head(top_NA_vs_invitro__exc$dmr)
#top_NA_vs_invitro_exc$comp
#top_NA_vs_invitro_exc$cgi
#names(top_NA_vs_invitro_exc)
head(top_NA_vs_invitro_exc$confects$table$name)

# Output for Meta-analysis
dir.create("misc")
write.table(top_NA_vs_invitro_metal, file="misc/choufani_top_NA_vs_invitro.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```

## Nat vs ART

```{r,ART_vs_Control_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 11: Differential analysis: Naturally conceived vs ART excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- targets
samplesheet$title <- as.character(samplesheet$title)
samplesheet$title <- (sapply(strsplit(samplesheet$title, " "), "[[", 1))
samplesheet$title <- (sapply(strsplit(samplesheet$title, "_"), "[[", 1))
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$title,levels=c("Control","ART"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_ART_vs_Control_exc <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_ART_vs_Natural_exc",
    myann=myann ,beta= beta)

head(top_ART_vs_Control_exc$dma)
head(top_ART_vs_Control_exc$dma)

# Allele
top_ART_vs_Control_exc$dma$unmeth <- "T"
top_ART_vs_Control_exc$dma$meth <- "C"
top_ART_vs_Control_exc$fit$SE <- sqrt(top_ART_vs_Control_exc$fit$s2.post) * top_ART_vs_Control_exc$fit$stdev.unscaled


# Extract required columns from dma
top_ART_vs_Control_metal <-top_ART_vs_Control_exc$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_ART_vs_Control_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_ART_vs_Control_exc$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsART")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_ART_vs_Control_exc$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsART")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_ART_vs_Control_metal <- merge(top_ART_vs_Control_metal, fitCE)
top_ART_vs_Control_metal <- merge(top_ART_vs_Control_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$title == "Control",]
cases <- samplesheet[samplesheet$Details.ART == "ART",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_ART_vs_Control_metal$N <- Neff

head(top_ART_vs_Control_metal)

# Output for Meta-analysis
dir.create("NATvsART")

write.table(top_ART_vs_Control_metal,file="NATvsART/choufani_top_ART_vs_Control.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```

## in vitro vs in vivo


```{r,in_vivo_vs_in_vitro_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 13: Differential analysis: in vitro vs in vivo excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- subset(targets,`art treatment:ch1`=="in vivo" | `art treatment:ch1`=="in vitro")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$`art treatment:ch1`,levels=c("in vivo","in vitro"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_invivo_vs_invitro_exc <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_invivo_vs_invitro_exc",
    myann=myann ,beta= beta)


# Allele
top_invivo_vs_invitro_exc$dma$unmeth <- "T"
top_invivo_vs_invitro_exc$dma$meth <- "C"
top_invivo_vs_invitro_exc$fit$SE <- sqrt(top_invivo_vs_invitro_exc$fit$s2.post) * top_invivo_vs_invitro_exc$fit$stdev.unscaled


# Extract required columns from dma
top_invivo_vs_invitro_metal <-top_invivo_vs_invitro_exc$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_invivo_vs_invitro_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_invivo_vs_invitro_exc$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsin vitro")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_invivo_vs_invitro_exc$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsin vitro")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_invivo_vs_invitro_metal <- merge(top_invivo_vs_invitro_metal, fitCE)
top_invivo_vs_invitro_metal <- merge(top_invivo_vs_invitro_metal, fitSE)

# Number of effective participants
controlss <- samplesheet[samplesheet$`art treatment:ch1` == "in vivo",]
cases <- samplesheet[samplesheet$`art treatment:ch1` == "in vitro",]
ctrl <- nrow(controlss)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))
top_invivo_vs_invitro_metal$N <- Neff


# Outputs for html
head(top_invivo_vs_invitro_metal)
head(top_invivo_vs_invitro_exc$dma)
#head(top_NA_vs_invivo__exc$dmr)
#top_NA_vs_invivo_exc$comp
#top_NA_vs_invivo_exc$cgi
#names(top_NA_vs_invivo_exc)
head(top_invivo_vs_invitro_exc$confects$table$name)

# Output for Meta-analysis
dir.create("misc")
write.table(top_invivo_vs_invitro_metal,file="misc/choufani_top_invivo_vs_invitro.tsv",sep="\t",quote=FALSE, row.names = FALSE) 

```

## Natural vs ICSI
```{r,NA_vs_ICSI, fig.width = 8 ,fig.height = 8, fig.cap="Figure 9: Differential analysis: Naturally conceived vs ICSI excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- subset(targets, `Details.ART`=="NAT" | `Details.ART`=="ICSI")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$Details.ART,levels=c("NAT","ICSI"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_NAT_vs_ICSI <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_NAT_vs_ICSI",
    myann=myann ,beta= beta)

# Allele
top_NAT_vs_ICSI$dma$unmeth <- "T"
top_NAT_vs_ICSI$dma$meth <- "C"
top_NAT_vs_ICSI$fit$SE <- sqrt(top_NAT_vs_ICSI$fit$s2.post) * top_NAT_vs_ICSI$fit$stdev.unscaled


# Extract required columns from dma
top_NAT_vs_ICSI_metal <-top_NAT_vs_ICSI$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_NAT_vs_ICSI_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_NAT_vs_ICSI$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsICSI")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_NAT_vs_ICSI$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsICSI")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_NAT_vs_ICSI_metal <- merge(top_NAT_vs_ICSI_metal, fitCE)
top_NAT_vs_ICSI_metal <- merge(top_NAT_vs_ICSI_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$Details.ART == "NAT",]
cases <- samplesheet[samplesheet$Details.ART == "ICSI",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_NAT_vs_ICSI_metal$N <- Neff





head(top_NAT_vs_ICSI$dma)
head(top_NAT_vs_ICSI$dmr)
head(top_NAT_vs_ICSI_metal)

# Output for Meta-analysis
dir.create("NATvsICSI")
write.table(top_NAT_vs_ICSI_metal, file="NATvsICSI/choufani_top_nat_vs_metal.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```


## Natural vs IUI
```{r,NA_vs_IUI, fig.width = 8 ,fig.height = 8, fig.cap="Figure 9: Differential analysis: Naturally conceived vs IUI excluding sex chromosomes"}
# Exclude Sex Chromosomes
samplesheet <- subset(targets, `Details.ART`=="NAT" | `Details.ART`=="IUI")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$Details.ART,levels=c("NAT","IUI"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_NAT_vs_IUI <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_NAT_vs_IUI",
    myann=myann ,beta= beta)

# Allele
top_NAT_vs_IUI$dma$unmeth <- "T"
top_NAT_vs_IUI$dma$meth <- "C"
top_NAT_vs_IUI$fit$SE <- sqrt(top_NAT_vs_IUI$fit$s2.post) * top_NAT_vs_IUI$fit$stdev.unscaled


# Extract required columns from dma
top_NAT_vs_IUI_metal <-top_NAT_vs_IUI$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_NAT_vs_IUI_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_NAT_vs_IUI$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsIUI")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_NAT_vs_IUI$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsIUI")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_NAT_vs_IUI_metal <- merge(top_NAT_vs_IUI_metal, fitCE)
top_NAT_vs_IUI_metal <- merge(top_NAT_vs_IUI_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$Details.ART == "NAT",]
cases <- samplesheet[samplesheet$Details.ART == "IUI",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_NAT_vs_IUI_metal$N <- Neff





head(top_NAT_vs_IUI$dma)
head(top_NAT_vs_IUI$dmr)
head(top_NAT_vs_IUI_metal)

# Output for Meta-analysis
dir.create("NATvsIUI")
write.table(top_NAT_vs_IUI_metal, file="NATvsIUI/choufani_top_nat_vs_IUI_metal.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```

## ICSI vs IUI
```{r,ICSI_vs_IUI_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 9: Differential analysis: ICSI vs IUI excluding sex chromosomes"}
#Exclude Sex Chromosomes
samplesheet <- subset(targets, `Details.ART`=="ICSI" | `Details.ART`=="IUI")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$Details.ART,levels=c("ICSI","IUI"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_ICSI_vs_IUI <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_ICSI_vs_IUI",
    myann=myann ,beta= beta)

# Allele
top_ICSI_vs_IUI$dma$unmeth <- "T"
top_ICSI_vs_IUI$dma$meth <- "C"
top_ICSI_vs_IUI$fit$SE <- sqrt(top_ICSI_vs_IUI$fit$s2.post) * top_ICSI_vs_IUI$fit$stdev.unscaled


# Extract required columns from dma
top_ICSI_vs_IUI_metal <-top_ICSI_vs_IUI$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_ICSI_vs_IUI_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_ICSI_vs_IUI$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsIUI")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_ICSI_vs_IUI$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsIUI")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_ICSI_vs_IUI_metal <- merge(top_ICSI_vs_IUI_metal, fitCE)
top_ICSI_vs_IUI_metal <- merge(top_ICSI_vs_IUI_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$Details.ART == "ICSI",]
cases <- samplesheet[samplesheet$Details.ART == "IUI",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_ICSI_vs_IUI_metal$N <- Neff





head(top_ICSI_vs_IUI$dma)
head(top_ICSI_vs_IUI$dmr)
head(top_ICSI_vs_IUI_metal)

# Output for Meta-analysis
write.table(top_NAT_vs_IUI_metal, file="NATvsICSI/choufani_top_nat_vs_metal.tsv",sep="\t",quote=FALSE, row.names = FALSE) 

head(top_ICSI_vs_IUI$dma)
head(top_ICSI_vs_IUI$dmr)

# Output for Meta-analysis
write.table(top_ICSI_vs_IUI_metal, file="misc/choufani_top_ICSI_vs_IUI.tsv",sep="\t",quote=FALSE) 
```

## NAT vs HYPO
```{r,Nat_vs_Hypo_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 10: Differential analysis: Natural vs Hypofertile, excluding sex chromosomes"}
#Exclude Sex Chromosomes
samplesheet <- subset(targets, `Details.ART`=="NAT" | `Details.ART`=="Hypofertile")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$Details.ART,levels=c("NAT","Hypofertile"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_na_vs_hypo <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_na_vs_hypo",
    myann=myann ,beta= beta)

# Allele
top_na_vs_hypo$dma$unmeth <- "T"
top_na_vs_hypo$dma$meth <- "C"
top_na_vs_hypo$fit$SE <- sqrt(top_na_vs_hypo$fit$s2.post) * top_na_vs_hypo$fit$stdev.unscaled


# Extract required columns from dma
top_na_vs_hypo_metal <-top_na_vs_hypo$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_na_vs_hypo_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_na_vs_hypo$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsHypofertile")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_na_vs_hypo$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsHypofertile")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_na_vs_hypo_metal <- merge(top_na_vs_hypo_metal, fitCE)
top_na_vs_hypo_metal <- merge(top_na_vs_hypo_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$Details.ART == "NAT",]
cases <- samplesheet[samplesheet$Details.ART == "Hypofertile",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_na_vs_hypo_metal$N <- Neff





head(top_na_vs_hypo$dma)
head(top_na_vs_hypo$dmr)
head(top_na_vs_hypo_metal)

# Output for Meta-analysis
dir.create("NATvsHypo")
write.table(top_na_vs_hypo_metal, file="NATvsHypo/choufani_top_na_vs_hypo.tsv",sep="\t",quote=FALSE, row.names = FALSE)
```
## NAT vs IVF
```{r,Nat_vs_IVF_exc, fig.width = 8 ,fig.height = 8, fig.cap="Figure 10: Differential analysis: Natural vs IVF, excluding sex chromosomes"}
#Exclude Sex Chromosomes
samplesheet <- subset(targets, `Details.ART`=="NAT" | `Details.ART`=="IVF")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "_"), "[[", 1)

groups <- factor(samplesheet$Details.ART,levels=c("NAT","IVF"))
sex <- factor(samplesheet$`gender:ch1`,levels=c("M","F"))

top_na_vs_ivf <- dm_analysis(samplesheet=samplesheet,
    sex=sex,groups=groups,mx=Mval_flt,name="top_na_vs_ivf",
    myann=myann ,beta= beta)

# Allele
top_na_vs_ivf$dma$unmeth <- "T"
top_na_vs_ivf$dma$meth <- "C"
top_na_vs_ivf$fit$SE <- sqrt(top_na_vs_ivf$fit$s2.post) * top_na_vs_ivf$fit$stdev.unscaled


# Extract required columns from dma
top_na_vs_ivf_metal <-top_na_vs_ivf$dma[,c("Row.names", "meth", "unmeth", "AveExpr", "P.Value")]
head(top_na_vs_ivf_metal)

# Convert fit outputs to dataframes
fitCE <- as.data.frame(top_na_vs_ivf$fit$coefficients)
fitCE$Row.names <- row.names(fitCE)
fitCE <- fitCE[,c("Row.names", "groupsIVF")]
names(fitCE)[2]<- "coefficient"


fitSE <- as.data.frame(top_na_vs_ivf$fit$SE)
fitSE$Row.names <- row.names(fitSE)
fitSE <- fitSE[,c("Row.names", "groupsIVF")]
names(fitSE)[2]<- "SE"

# Merge Datasets
top_na_vs_ivf_metal <- merge(top_na_vs_ivf_metal, fitCE)
top_na_vs_ivf_metal <- merge(top_na_vs_ivf_metal, fitSE)

# Number of effective participants
controls <- samplesheet[samplesheet$Details.ART == "NAT",]
cases <- samplesheet[samplesheet$Details.ART == "IVF",]
ctrl <- nrow(controls)
cses <- nrow(cases)
Neff <- 4/((1/cses)+(1/ctrl))

top_na_vs_ivf$N <- Neff





head(top_na_vs_ivf$dma)
head(top_na_vs_ivf$dmr)
head(top_na_vs_ivf_metal)

# Output for Meta-analysis
dir.create("NATvsIVF")
write.table(top_na_vs_ivf_metal, file="NATvsIVF/choufani_top_na_vs_ivf.tsv",sep="\t",quote=FALSE, row.names = FALSE) 
```


# Venn diagrams of the differential methylated probes for each contrast
First, we look at the similarity of DMPs altered between in vivo and in vitro. 

```{r,venn1,fig.cap="Figure 14:"}
#v1 <- list("in vivo up" = top_NA_vs_invivo_exc$dm_up , 
#           "in vitro up" = top_NA_vs_invitro_exc$dm_up ,
#          "in vivo dn" = top_NA_vs_invivo_exc$dm_dn ,
#         "in vitro dn" = top_NA_vs_invitro_exc$dm_dn )
#plot(euler(v1, shape = "ellipse"), quantities = TRUE)
```


# Spearman correlations of each contrast

```{r,spearman}
myranks <- function(x) {
  xx<-x[[1]]
  xx$score <- sign(xx$logFC)/log10(xx$adj.P.Val)
  y <- xx[,"score",drop=FALSE]
  y$rn <- xx$Row.names
  return(y)
}

mycontrasts <- list("top_NA_vs_invivo_exc" = top_NA_vs_invivo_exc,"top_NA_vs_invitro_exc" = top_NA_vs_invitro_exc, "top_ART_vs_Control_exc" = top_ART_vs_Control_exc, "top_invivo_vs_invitro_exc" = top_invivo_vs_invitro_exc, "top_NAT_vs_ICSI" = top_NAT_vs_ICSI, "top_ICSI_vs_IUI" = top_ICSI_vs_IUI)

myrnks <- lapply(X = mycontrasts, FUN = myranks)
df <- join_all(myrnks,by="rn")
rownames(df) <- df$rn
df$rn=NULL
colnames(df) <- names(mycontrasts)
head(df)

mycors <- cor(df,method = "spearman")

my_palette <- colorRampPalette(c("darkred","red", "orange", "yellow","white"))(n = 25)
heatmap.2(mycors,scale="none",margin=c(10, 10),cexRow=0.8,trace="none",cexCol=0.8,
    col=my_palette,main="Spearman correlations")

mycors

```

# Enrichment analysis

We will be using the recently published package mitch to perform enrichment analysis, using
average promoter methylation change as an indicator of gene activity.
Enrichment will be tested with the mitch package.

```{r,genesets}
library("mitch")
# gene sets
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", 
    destfile="ReactomePathways.gmt.zip")
unzip("ReactomePathways.gmt.zip",overwrite = TRUE)
genesets <- gmt_import("ReactomePathways.gmt")
```

## Run 1D enrichment analysis

One dimensional enrichment analysis with REACTOME gene sets using average promoter methylation t-statistic.

```{r,1danalysis}
top_NA_vs_invitro_exc_mitch <- run_mitch_1d(dma= top_NA_vs_invitro_exc$dma, name="top_NA_vs_invitro_mitch")
head(top_NA_vs_invitro_exc_mitch,50)

top_NA_vs_invivo_exc_mitch <- run_mitch_1d(dma= top_NA_vs_invivo_exc$dma, name="top_NA_vs_invivo_mitch")
head(top_NA_vs_invitro_exc_mitch,50)

top_invivo_vs_invitro_exc_mitch <- run_mitch_1d(dma= top_invivo_vs_invitro_exc$dma, name="top_invivo_vs_invitro_mitch")
head(top_invivo_vs_invitro_exc_mitch,50)

top_ART_vs_Control_exc_mitch <- run_mitch_1d(dma= top_ART_vs_Control_exc$dma, name="top_ART_vs_Control_exc_mitch")
head(top_ART_vs_Control_exc_mitch,50)

top_NAT_vs_ICSI_mitch <- run_mitch_1d(dma= top_NAT_vs_ICSI$dma, name="top_top_NAT_vs_ICSI_mitch")
head(top_NAT_vs_ICSI_mitch,50)

top_ICSI_vs_IUI_mitch <- run_mitch_1d(dma= top_ICSI_vs_IUI$dma, name="top_top_ICSI_vs_IUI_mitch")
head(top_ICSI_vs_IUI_mitch,50)
```

## Run multi-dimensional enrichment analysis


```{r,mitch}
#One dimensional enrichment analysis with REACTOME gene sets using average promoter methylation t-statistic.

#xl <- list("NA_vs_invitro"=top_NA_vs_invitro_exc, "top_NA_vs_invivo"=top_NA_vs_invivo_exc$dma, "top_invivo_vs_invitro"=top_invivo_vs_invitro_exc$dma,
#           "top_ART_vs_Control"=top_ART_vs_Control_exc$dma)
#xxl <- lapply(X = xl,run_mitch_rank)  
#xxxl <- lapply(xxl,function(xxl) { xxl$genenames <- rownames(xxl) ; xxl} )
#xxll <- join_all(xxxl,by="genenames")
#rownames(xxll) <- xxll$genenames
#xxll$genenames=NULL
#colnames(xxll) <- names(xl)
#head(xxll)

#capture.output(
#        res <- mitch_calc(xxll,genesets = genesets,priority = "significance")
#        , file = "/dev/null", append = FALSE,
#        type = c("output", "message"), split = FALSE)

#head(res$enrichment_result,20)

#    capture.output(
#        mitch_report(res,outfile=paste("mitch_report.html"))
#        , file = "/dev/null", append = FALSE,
#        type = c("output", "message"), split = FALSE)

#```

## Session Information

```{r,session}

sessionInfo()

```