Gordon2012.Rmd

---
title: "Gordon 2012 Reanalysis"
author: "The methylation in ART meta-analysis group"
date:  "`r Sys.Date()`"
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 [Gordon et al (2012)](https://genome.cshlp.org/content/early/2012/07/05/gr.136598.111).

In their study, they analysed the DNA methylation patterns of CBMCs (18MZ, 9DZ), HUVECS (14MZ, 10DZ) and Placenta (8MZ, 7DZ)

The platform used in the study is the Illumina Infinium HumanMethylation27k BeadChip assay.
The authors used packages from the Bioconductor Project.

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

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

```{r,packages}
suppressPackageStartupMessages({
    library("R.utils")
    library("missMethyl")
    library("limma")
    library("minfi")
    library("RColorBrewer")
    library("eulerr")
    library("plyr")
    library("gplots")
    library("reshape2")
    library("beeswarm")
    library("GEOquery")
    library("Biobase")
    library("readxl")
  })
```

## Loading Functions
These functions provide shortcuts to help with charts and other analysis. They will eventually be shoved into another Rscript or package but can stay here for now.
```{r,load_functions}
# scree plot shows the amount of variation in a dataset that is accounted
# for by the first N principal components
myscree <- function(mx,n=10,main="") {
  pc<-princomp(Mval)$sdev
  pcp <- pc/sum(pc)*100
  pcp <- pcp[1:10]
  barplot(pcp,cex.names = 1,las=2,ylim=c(0,60),
      ylab="percent (%) variance explained", main=main)
  text((0.5:length(pcp)*1.2),pcp,label=signif(pcp,3),pos=3,cex=0.8)
}

# Here is a function to make a volcano plot
make_volcano <- function(dm,name,mx) {
    sig <- subset(dm,adj.P.Val<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,logFC>0))
    N_DN=nrow(subset(sig,logFC<0))
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn")
    plot(dm$logFC,-log10(dm$P.Val),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log FC", ylab="-log10 pval")
    mtext(HEADER)
    grid()
    points(sig$logFC,-log10(sig$P.Val),cex=0.5,pch=19,col="red")
}

# Here is a function to make heatmaps 
make_heatmap <- function(dm,name,mx,n) {
  topgenes <-  rownames(head(dm[order(dm$P.Value),],n))
  ss <- mx[which(rownames(mx) %in% topgenes),]
  my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 25)
  heatmap.2(ss,scale="row",margin=c(10, 10),cexRow=0.4,trace="none",cexCol=0.4,
      col=my_palette, main=name)
}

# make beeswarm charts
# dm = a limma differential meth object
# name = character name of the limma dm object
# mx = matrix of normalised data
# groups = a vector of factors corresponding to the cols in mx
# n = the number of top significant genes to plot (default = 15) 
make_beeswarms <- function(dm,name,mx,groups,n=15) {
    par(mar=c(3,3,1,1))
    NCOLS=5
    NROWS=floor(n/NCOLS)
    if (n %% NCOLS > 0) { NROWS <- NROWS + 1 }
    par(mfrow=c(NROWS, NCOLS))
    topgenes <-  rownames(head(dm[order(dm$P.Value),],n))
    ss <- mx[which(rownames(mx) %in% topgenes),]
    n <- 1:n
    g1name=levels(groups)[1]
    g2name=levels(groups)[2]
    g1dat <- ss[n,which(groups == g1name)]
    g2dat <- ss[n,which(groups == g2name)]
    g1l <-lapply(split(g1dat, row.names(g1dat)), unlist)
    g2l <-lapply(split(g2dat, row.names(g2dat)), unlist)

    for (i in n) {
      mydat <- list(g1l[[i]],g2l[[i]])
        beeswarm(mydat,ylim=c(-8,8),cex=0.2, pch=19,
        las=2, cex.lab=0.6, main=names( g1l )[i] , 
        ylab="",labels = c(g1name,g2name))
      grid()
    }
}


make_dm_plots <- function(dm,name,mx,groups=groups) {
    make_volcano(dm,name,mx)
    make_beeswarms(dm ,name , mx , groups , n= 15)
    make_heatmap(dm , name , mx ,n = 50)
} 
```


## Data Import

```{r,download}

ARRAY_DATA="pets_27k" 
# only download it if it is not present on the system


# Series Matrix
SERIES_MATRIX="GSE36642_SERIES_MATRIX.txt.gz"
download.file("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE36nnn/GSE36642/matrix/GSE36642_series_matrix.txt.gz",destfile=SERIES_MATRIX)

# Reading in GEO series matrix
gse <- getGEO(filename="/mnt/mnorris/ART_methylation/GSE36642_SERIES_MATRIX.txt.gz")
targets <- pData(phenoData(gse))

untar(exdir = "IDAT", tarfile = "GSE36642_RAW.tar",)
SampleTable <- read_excel("~/SampleTable Infinium 27k(2).xlsx")
SampleTableSentrix <- paste(SampleTable$Sentrix_ID,SampleTable$Sentrix_Position, sep = "_")
View(Sample)
baseDir <- "."

IDATS <- list.files("/mnt/mnorris/pets_27k/",recursive = TRUE, pattern="idat") 
IDATS <- IDATS[grep("Error",IDATS,invert = TRUE)] 
IDATS
basename(IDATS)
SentrixID <- sapply(strsplit(basename(IDATS),"_"),"[[",1) 
SentrixPOS <- sapply(strsplit(basename(IDATS),"_"),"[[",2) 
Sentrix <- paste(SentrixID,SentrixPOS,sep = "_")

intersect(SampleTableSentrix,SentrixID)

rgSet <- read.metharray.exp(targets = targets)
rgSet

head(targets)
```

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


## Fliter 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 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)
dim(Mval)

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

```

## MDS Analysis

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Session Information

```{r,session}

sessionInfo()

```