---
title: "Autism Spectrum Disorder Methylation analysis - Guthrie card analysis"
author: "Namitha Mohandras, Mark Ziemann"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    fig_width: 7
    fig_height: 7
theme: cosmo
---

Source: https://github.com/markziemann/asd_meth

## Introduction

Here we are analysing Guthrie card methylation in 9 twin pairs.

```{r, packages}

baseDir=getwd()
dataDir=paste(baseDir,"/ASD_EPIC_DATA",sep="")
library("missMethyl")
library("limma")
library("minfi")
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
library("ruv")
library("RColorBrewer")
library("matrixStats")
library("gplots")
library("WGCNA")
library("FlowSorted.Blood.450k")
library("reshape2")
library("ggplot2")
library("missMethyl")
library("DMRcate")
library("FlowSorted.Blood.EPIC")
library("mitch")
library("kableExtra")

```

## Load data

Load the annotation data and the Epic methylation data.

This analysis is to be conducted on Ubuntu with R4.

```{r,load_data1}

ann = getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b2.hg19)
targets_gen = read.metharray.sheet(dataDir, pattern = "ASD_guthrie_summarysheet.csv")
#targets$ID = paste(targets$Sample_Group,targets_gen$Sample_Name,sep=".")
rgSet = read.metharray.exp(targets = targets_gen)
sampleNames(rgSet) = targets_gen$SampleID

```

## Testing MZ status

```{r,snptest}

snpBetas = getSnpBeta(rgSet)
d = dist(t(snpBetas))
hr = hclust(d, method = "complete", members=NULL)
plot(hr)

```

## Quality control

```{r,qc1}

detP = detectionP(rgSet)
qcReport(rgSet, sampNames = targets_gen$Sample_Name,
  pdf="qc-report_ASD_guthrie_nov27.pdf")
cols=brewer.pal(4,"Set1")

barplot(apply(detP,2,mean),
  col=as.numeric(factor(targets_gen$Sample_Name)),
  las=2,cex.names= 0.8, cex.axis=0.75,
  main="Mean detection p-values of probe signals",
  ylab="Mean detection p-value")

barplot(apply(detP,2,mean),
  col=as.numeric(factor(targets_gen$Sample_Name)),
  las=2,cex.names= 0.8, cex.axis=0.75,ylim=c(0,0.010),
  main="Mean detection p-values of probe signals",
  ylab="Mean detection p-value")

```

## Preprocessing

```{r,preprocess1}

mset.raw = preprocessRaw(rgSet)

```

## Data exploration

Using Multi-dimensional scaling (MDS) plots before filtering.

```{r,mds1}

mdsPlot(mset.raw, sampGroups = targets_gen$Sample_Name, 
  sampNames=targets_gen$Social_interaction_on_ADOS,legendPos="bottom")

mdsPlot(mset.raw, sampGroups = targets_gen$Sex, 
  sampNames=targets_gen$SampleID,legendPos="bottom")

```


## Normalisation

SQN is used as it gave a more stable result.

```{r,norm1}

#SQN method
mSetSq = preprocessQuantile(rgSet)

densityPlot(rgSet, sampGroups = targets_gen$Social_interaction_on_ADOS,
  main="Raw", legend = FALSE)

densityPlot(getBeta(mSetSq), 
  sampGroups = targets_gen$Social_interaction_on_ADOS,main="SQN", 
  legend = FALSE)

```

## Cell type composition analysis

Using old method.

```{r,celltypes1}

rgSet$Slide <- as.numeric(rgSet$Slide)
rgSet$Sex <- as.character(rgSet$Sex)
#rgSet$Sample_Name <- as.character(rgSet$Sample_Name)
#sampleNames(rgSet) = targets_gen$SampleID

cellCounts_new <- estimateCellCounts(rgSet, compositeCellType = "Blood", 
  processMethod = "auto", probeSelect = "auto", 
  cellTypes = c("CD8T","CD4T", "NK","Bcell","Mono","Gran"),
  referencePlatform = c("IlluminaHumanMethylation450k"),
  returnAll = FALSE, meanPlot = TRUE)

#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")

cellCounts_long <- melt(cellCounts_new, id = "celltype")

cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1,  targets_gen$SampleID),"Diagnosis"]

colnames(cellCounts_long) <- c("diagnosis","celltype","value")

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) + 
  geom_boxplot() 

wx <- lapply(1:ncol(cellCounts_new) , function(i) {
  wilcox.test(a[,i],c[,i], paired=FALSE)
} )
names(wx) <- colnames(cellCounts_new)
wx

```

Try estimatecellcounts2.

```{r,celltypes2}

cells <- estimateCellCounts2(rgSet, referencePlatform= "IlluminaHumanMethylationEPIC",
  returnAll = TRUE)

mset <- cells$normalizedData

cellCounts_new <- cells[[1]]
#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")

cellCounts_long <- melt(cellCounts_new, id = "celltype")

cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1,  targets_gen$SampleID),"Diagnosis"]

colnames(cellCounts_long) <- c("diagnosis","celltype","value")

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
  geom_boxplot()

pdf("cells_guthrie.pdf")

par(mar=c(3,8,2,2))

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
  geom_boxplot()

dev.off()

wx <- lapply(1:ncol(cellCounts_new) , function(i) {
  wilcox.test(a[,i],c[,i], paired=FALSE)
} )
names(wx) <- colnames(cellCounts_new)
wx

#densityByProbeType(mSetSq, main = "Raw")
#densityByProbeType(mset[,1], main = "estimateCellCounts2")
#mset_reduced <- mset[which(rownames(mset) %in% names(keep[keep==TRUE])),]
#meth <- getMeth(mset_reduced)
#unmeth <- getUnmeth(mset_reduced)
#Mval <- log2((meth + 100)/(unmeth + 100))
#beta <- getBeta(mset_reduced)
#dim(Mval)

```

## Filtering

None of the p-values are less than 0.05.

Now filter for high detection p-value and overlap with SNP.

Need to get a copy of the Xreact probes from Namitha. 
In the mean time I have downloaded from github - link below.

```{r,filt1}

detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mset <- mset[keep,]
gmset <- mapToGenome(mset)

#remove SNPs
gmset_flt = dropLociWithSnps(gmset, snps = c("CpG", "SBE"))

#Removing cross-reactive probes
XURL="https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/EPIC/13059_2016_1066_MOESM1_ESM.csv"
Xreact <- read.csv(XURL)

#Xreact = read.csv(file="/group/canc2/puumba/Data/InfiniumData/NamithaData/Rprojects/Autism/Analysis_Sept11/EPIC_850k_crossreactiveProbes.csv", stringsAsFactors=FALSE)
#Xreact = read.csv(file="~/48639-non-specific-probes-Illumina450k.csv", stringsAsFactors=FALSE)
noXreact <-  !(featureNames(gmset) %in% Xreact$X)

gmset <- gmset[noXreact,]

#Removing probes on X and Y chromosomes
autosomes <- !(featureNames(gmset) %in% ann$Name[ann$chr %in% c("chrX","chrY")])
gmset_flt <- gmset[autosomes,]

#Relative log expression (RLE plot)
mvals = getM(gmset_flt)
medSq = apply(mvals, 1, median)
YSq = mvals - medSq

boxplot(YSq,outline=FALSE,ylim=c(-1.5,1.5), 
  ylab="Relative Log Methylation Value", 
  cols=as.character(factor(targets_gen$Social_interaction_on_ADOS,)) )

```


## MDS plots generation after filtering

```{r,mds2}

pal = brewer.pal(8, "Dark2")
mds1Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(1,2))
mds2Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(1,3))
mds3Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(2,3))
mds4Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(3,4))

plotMDS(mds1Sq, xlab="Dimension 1", ylab="Dimension 2",
  col=pal[as.factor(targets_gen$Diagnosis)],
  dim=c(1,2), labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds2Sq, xlab="Dimension 1", ylab="Dimension 3",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(1,3),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds3Sq, xlab="Dimension 2", ylab="Dimension 3",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(2,3),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds4Sq, xlab="Dimension 3", ylab="Dimension 4",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(3,4),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

```

## Principal Component Analysis (PCA)

```{r,pca1}

fit <- prcomp(t(mvals),center = TRUE, scale = TRUE,retx=TRUE)
loadings = fit$x
plot(fit,type="lines")
nGenes = nrow(mvals)
nSamples = ncol(mvals)
datTraits = targets_gen[,7:15]
moduleTraitCor = cor(loadings[,1:6], datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
par(cex=0.75, mar = c(6, 8.5, 3, 3))
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", 
  signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)

labeledHeatmap(Matrix = t(moduleTraitCor), 
  xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
  yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
  textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
  cex.lab.y = 0.6, zlim = c(-1,1), 
  main = paste("PCA-trait relationships: Top principal components"))

pdf("pca_guthrie.pdf")

par(mar=c(3,8,2,2))

labeledHeatmap(Matrix = t(moduleTraitCor), 
  xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
  yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
  textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
  cex.lab.y = 0.6, zlim = c(-1,1), 
  main = paste("PCA-trait relationships: Top principal components"))

dev.off()

#blah()

```

## Regression analysis - within twin pair on ADOS

ADOS is the outcome variable (continuous)

```{r,reg1}

mDat = mvals
Twin_Group <- factor(targets_gen$Family_ID)
ADOS <- targets_gen$Social_interaction_on_ADOS
names(ADOS) <- targets_gen$SampleID
design_tw <- model.matrix(~ Twin_Group + ADOS )

fit_tw <- lmFit(mDat, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
top <- topTable(fit2_tw,coef="ADOS",num=Inf, sort.by = "P")
head(top)
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)

ann_sub = 
  ann[,c("chr","pos","strand","Name","Islands_Name",
    "Relation_to_Island","UCSC_RefGene_Name","UCSC_RefGene_Group")]

output <-merge(ann_sub,top,by.x="Name",by.y="row.names")

write.csv(output[order(output$P.Value),],
  file="ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv",row.names=FALSE)

top1000 <- topTable(fit2_tw,coef="ADOS",num=1000, sort.by = "P")
output_1000 <- merge(ann_sub,top1000,by.x="Name",by.y="row.names")

write.csv(output_1000[order(output_1000$P.Value),],
  file="ASD_guthrie_top_dmps_onADOS_Nov27_withintw_limma_top1k.csv",row.names=FALSE)

saveRDS(design_tw, "gu_design.rds")
saveRDS(mvals, "gu_mvals.rds")

output <- output[order(output$P.Value),]

output <- subset(output,P.Value<1e-4)

output %>% kbl() %>% kable_paper("hover", full_width = F)

```

## Converting to beta from M values

```{r,betavals1}

bDat = ilogit2(mDat)
bDat_new = getBeta(gmset_flt)
#View(bDat)
write.csv(bDat,file="ASD_guthrie_beta_onADOS_withintw_Nov27.csv",row.names=TRUE)

```

## RUV

Extract illumina negative control

```{r,ruv1}

#Extract illumina negative control
INCs <- getINCs(rgSet)
head(INCs)
#add negative control data to M-values
Mc <- rbind(mvals,INCs)
ctl <- rownames(Mc) %in% rownames(INCs)
table(ctl)

#old RUV syntax is updated
#rfit1 <- RUVfit(data=Mc, design=design_tw, coef=2, ctl=ctl) # Stage 1 analysis

# this is the updatedcommand - not sure if correct
rfit1 <- RUVfit(Y = Mc, X = design_tw[,"ADOS"] , ctl = ctl) # Stage 1 analysis

rfit2 <- RUVadj(Y = Mc, fit = rfit1)

top1 <- topRUV(rfit2, num=Inf)

head(top1)

#ctl <- rownames(mValsSq) %in% rownames(top1[top1$p.ebayes > 0.5,])
#bottom 80% used as negative controls
ctl <- rownames(mvals) %in% rownames(top1[ceiling(0.2*nrow(top1)):nrow(top1),])
table(ctl)

# Perform RUV adjustment and fit
rfit1 <- RUVfit(Y = mvals, X=design_tw[,"ADOS"], ctl=ctl) # Stage 2 analysis 
rfit2 <- RUVadj(Y = mvals, fit=rfit1)
top_ruv <- topRUV(rfit2, number = Inf)
top_ruv_1000 <- topRUV(rfit2, number = 1000)
head(top_ruv_1000)
#dat_subset_ruv = as.matrix(subset(mValsSq, rownames(mValsSq) %in% rownames(top_ruv_1000)))
#heatmap.2(dat_subset_ruv, cexRow=0.8, cexCol=0.8, col=colorRampPalette(c("Blue","Yellow")), density.info="none", trace="none", scale="row", main="diff meth probes due to epilepsy")
output_ruv = merge(ann,top_ruv,by.x="Name",by.y="row.names")
output_ruv_1000 = merge(ann,top_ruv_1000,by.x="Name",by.y="row.names")
#output_full_epicanno = merge(ann_epic,top,by.x="Name",by.y="row.names")

write.csv(output_ruv[order(output_ruv$F.p),],
  file="ASD_guthrie_topruv80%_dmps_onADOS_withintw_Nov29.csv",
  row.names=FALSE)

write.csv(output_ruv_1000[order(output_ruv_1000$F.p),],
  file="ASD_guthrie_topruv80%_1kdmps_onADOS_withintw_Nov29.csv",
  row.names=FALSE)

output_ruv <- data.frame(output_ruv[order(output_ruv$F.p),])

output_ruv <- output_ruv[order(output_ruv$p_X1),]

output_ruv <- subset(output_ruv,p_X1<1e-4)

output_ruv %>% kbl() %>% kable_paper("hover", full_width = F)

```

## Visualise the effect of RUV adjustment

```{r,ruv2}

#Madj <- getAdjusted(mValsSw, rfit2)
# cant get this part to work
mvals <- mvals[which(rownames(mvals) %in% colnames(rfit2$Y) ),]
#Madj <- missMethyl::getAdj(Y = mValsSq, fit = rfit2)

plotMDS(mvals, labels=targets_gen$Sample_Name,
  col=as.integer(factor(targets_gen$Social_interaction_on_ADOS)),
  main="Unadjusted")

#legend(legend=c("1","0"),pch=7,cex=0.7,col=1:2)

# i had to comment the below stuff out because Madj couldn't be generated.
#plotMDS(Madj, labels=targets_gen$Sample_Name, 
#  col=as.integer(factor(targets_gen$Social_interaction_on_ADOS)),
#  main="Adjusted: RUV-inverse")
#legend(legend=c("1","0"),pch=7,cex=0.7,col=1:2)

#Trying different k values
# mark fixed a major bug here
rfit_k1 <- RUVfit(Y=mvals, X=design_tw[,"ADOS"], ctl=ctl, 
  method = "ruv4", k=1) # Stage 2 analysis 

rfit_k1 <- RUVadj(Y=mvals, fit=rfit_k1)

rfit_k2 <- RUVfit(Y=mvals, X=design_tw[,"ADOS"], ctl=ctl,
  method = "ruv4", k=2) # Stage 2 analysis 

rfit_k2 <- RUVadj(Y=mvals, fit=rfit_k2)

rfit_k8 <- RUVfit(Y=mvals, X=design_tw[,"ADOS"], ctl=ctl,
  method = "ruv4", k=8) # Stage 2 analysis 

rfit_k8 <- RUVadj(Y=mvals, fit=rfit_k8)

#get adjusted values - doesn't work as at Apr 2022.
#Madj1 <- missMethyl::getAdj(Y=mValsSq, fit=rfit_k1)
#Madj2 <- missMethyl::getAdj(Y=mValsSq, fit=rfit_k2)
#Madj3 <- missMethyl::getAdj(Y=mValsSq, fit=rfit_k8)

#Visualise plots
par(mfrow=c(1,1))
plotMDS(mvals, labels=targets_gen$Sample_Name, 
  col=as.integer(factor(targets_gen$Diagnosis)),
  dim=c(14,15),
  main="Adjusted")

#legend("topleft",legend=c("Cancer","Normal"),pch=16,cex=1,col=1:2)

#plotMDS(Madj, labels=targets_gen$Sample_Name, col=as.integer(factor(targets_gen$Diagnosis..L.0..M.1..H.2.)),dim=c(14,15),
#        main="Adjusted: RUV-inverse")
#legend("topleft",legend=c("Cancer","Normal"),pch=16,cex=1,col=1:2)

#plotMDS(Madj1, labels=targets_gen$Sample_Name, col=as.integer(factor(targets_gen$Diagnosis..L.0..M.1..H.2.)),dim=c(14,15),
#        main="Adjusted: RUV-4, k=1")
#legend("bottomleft",legend=c("Cancer","Normal"),pch=16,cex=1,col=1:2)

#plotMDS(Madj2, labels=targets_gen$Sample_Name, col=as.integer(factor(targets_gen$Diagnosis..L.0..M.1..H.2.)),dim=c(14,15),
#        main="Adjusted: RUV-4, k=2")
#legend("bottomright",legend=c("Cancer","Normal"),pch=16,cex=1,col=1:2)

#plotMDS(Madj3, labels=targets_gen$Sample_Name, col=as.integer(factor(targets_gen$Diagnosis..L.0..M.1..H.2.)),dim=c(14,15),
#        main="Adjusted: RUV-4, k=8")
#legend

```

## Plotting effect sizes of top DMPs

```{r,effect1}

res <- read.csv("ASD_guthrie_topruv80%_dmps_onADOS_withintw_Nov29.csv", header = TRUE)
TOPPROBES <- head(res$Name,9)

effect_plot <- function(PROBE) {
  GENE <- res[which(res$Name==PROBE),"UCSC_RefGene_Name"]
  GENE <- paste0(unique(unlist(strsplit(GENE,";"))),collapse=";")
  HEADER=paste(GENE,PROBE)
  BETAS <- ilogit2(mDat[rownames(mDat) == PROBE,])
  plot(BETAS, ADOS, main = HEADER, ylab = "ADOS", xlab = "Beta Value")
  regl5 <- lm(ADOS ~ BETAS)
  SLOPE=signif(regl5$coefficients[2],3)
  abline(regl5)
  mycor <- cor.test(ADOS, BETAS, method = "pearson")
  mycor_stat <- signif(mycor$estimate,3)
  mycor_p <- signif(mycor$p.value,3)
  SUBHEAD <- paste("R=",mycor_stat," p=",mycor_p," slope=",SLOPE, sep="")
  mtext(SUBHEAD,cex=0.6)

}

par(mfrow=c(3,3))

sapply(TOPPROBES,effect_plot)

pdf("effect_guthrie.pdf")
par(mfrow=c(3,3))
sapply(TOPPROBES,effect_plot)
dev.off()

par(mfrow=c(1,1))

```


## Gene ontology of top DMPs

Gene Ontology analysis (gometh): top 1000 probes

```{r,gometh1}

res <- read.csv("ASD_guthrie_topruv80%_dmps_onADOS_withintw_Nov29.csv", header = TRUE)
sigCpGs_1k = res$Name[1:1000]

sigCpGs_1k = as.character(sigCpGs_1k)
all = res$Name
length(all)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)
write.csv(gometh_kegg, file = "GOmeth_kegg.csv", row.names = TRUE)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)
write.csv(gometh_go, file = "GOmeth_go.csv", row.names = TRUE)

```

Now specifically analyse hypermethylated probes

```{r,gometh2}

res2 <- subset(res,b_X1>0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

```

Now specifically analyse hypomethylated probes

```{r,gometh3}

res2 <- subset(res,b_X1<0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

```

## Gene ontology of top DMPs LIMMA data

Gene Ontology analysis (gometh): top 1000 probes (limma data).

```{r,gometh4}

res <- read.csv("ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv", header = TRUE)
sigCpGs_1k = res$Name[1:1000]

sigCpGs_1k = as.character(sigCpGs_1k)
all = res$Name
length(all)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

```

Now specifically analyse hypermethylated probes

```{r,gometh5}

res2 <- subset(res,logFC>0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

```

Now specifically analyse hypomethylated probes

```{r,gometh6}

res2 <- subset(res,logFC<0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

```

## DMRCate - differentially methylated region analysis

```{r,dmrcate1}

#design matrix in regression 
design_tw <- model.matrix(~ADOS+Twin_Group)
design_tw_dmrc <- model.matrix(~ADOS)
#myannotation <- cpg.annotate("array", mDat, analysis.type="differential", design=design_tw, coef=2, fdr=1)
myannotation <- cpg.annotate("array", mDat, what = "M", 
  arraytype = "EPIC", analysis.type="differential", design_tw, coef=2)
dmrcoutput <- dmrcate(myannotation, lambda=1000, C=2, pcutoff=0.001)
results.ranges <- data.frame(extractRanges(dmrcoutput, genome = "hg19") )
head(results.ranges, 20)
write.csv(results.ranges, file = "dmrcoutput_guthrie.csv", row.names = TRUE)

```

## Session information

```{r,sessioninfo}

sessionInfo()

```
