---
title: "SAM Gene Modelling"
author: "Lilly Dipnall"
date: "2023-06-29"
output: html_document
---

#Knit
This chunk sets the directory for when the file knits to a document
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require("knitr")
opts_knit$set(root.dir = "/Users/cogld/Documents/Thesis/epignetics/st3gal3_analysis/st3gal3_repo")
```

```{r}
setwd("/Users/cogld/Documents/Thesis/epignetics/st3gal3_analysis/st3gal3_repo")
```

## Sourcing Packages
```{r}
suppressPackageStartupMessages({
    library("zoo")
    library("gplots")
    library("beeswarm")
    library("topconfects")
    library("FlowSorted.Blood.EPIC")
    library("mitch")
    library("missMethyl")
    library("limma")
    library("minfi")
    library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
    library("IlluminaHumanMethylationEPICmanifest")
    library("RColorBrewer")
    library("matrixStats")
    library("minfiData")
    library("Gviz")
    library("DMRcate")
    library("stringr")
})
```

#Import Data
```{r}
load("~/SyncDrive/RDS27534-Epigenetics/All NICAP EWAS Cleaned/m.08JUN2023.Rdata")
load("~/SyncDrive/RDS27534-Epigenetics/All NICAP EWAS Cleaned/beta.08JUN023.Rdata")

targets <- read.csv("/Users/cogld/Documents/Thesis/epignetics/P3/data/Sample_annotation_demo.csv")
#Convert ID column into string 
targets$id <- as.character(targets$id)
rownames(targets) <- targets$id
targets$adhd <- ifelse(targets$group == 'ADHD', 1, 0)
targets$control <- ifelse(targets$group == 'Control|High Risk', 1, 0)
```

#Subsetting Mval and B DFs for SAM Genes
```{r}
rownames(Mval) <- rownames(beta)
colnames(Mval) <- sub("X", "", colnames(Mval))

#Subsetting Mval to only include W1 IDs
Mval <- as.data.frame(Mval)
Mval1 <- Mval %>% select(-(contains("_4") | contains("_5")))
colnames(Mval1) <- sub("_3", "", colnames(Mval1))
Mval1 <- Mval1 %>% select(-(contains("dup")))

#Subsetting Beta to only include W1 IDs
beta <- as.data.frame(beta)
beta1 <- beta %>% select(-(contains("_4") | contains("_5")))
colnames(beta1) <- sub("_3", "", colnames(beta1))
beta1 <- beta1 %>% select(-(contains("dup")))

#Subsetting Mval and Beta for SAM probes
SAM_Mval <- Mval1[which(rownames(Mval1) %in% SAM_probes_all_list2),]
SAM_beta <- beta1[which(rownames(beta1) %in% SAM_probes_all_list2),]

#Match IDs in targets file to Mval + Beta
SAM_probes_all_list2 <- unique(SAM_probes_all_list)
rownames(targets) <- targets$id

#Match Mval and Beta to Targets file 
Mval2 <- Mval1[,rownames(targets)]
beta2 <- beta1[,rownames(targets)]
beta2$cpg <- rownames(beta2)
SAM_Mval2 <- SAM_Mval[,rownames(targets)] 
SAM_beta2 <- SAM_beta[,rownames(targets)]

dim(targets)
dim(Mval2)
dim(beta2) #162? Come back and check - 163 due to the cpg variable being created
dim(SAM_Mval2)
dim(SAM_beta2)

#compare Mval and Beta probes to SAM probe list
library(arsenal)
summary(comparedf(SAM_Mval, SAM_probes_all))


```



#Regression Models for Diagnostic Status
#1. ADHD diagnosis predicting DNA methylation
```{r}
## ----dmps---------------------------------------------------------------------
# this is the factor of interest
adhd <- ifelse(targets$group == 'ADHD', 1, 0)
control <- ifelse(targets$group == 'Control|High Risk', 1, 0)


group <- factor(adhd)
# this is the individual effect that we need to account for
sex <- factor(targets$ChildGender) 
epi <- as.numeric((targets$Epi))
sample_group <- as.factor((targets$Sample_Group))
array <- as.factor((targets$Array))

# use the above to create a design matrix
#ask about the order and if it matters where each variable sits
#batch correction - run a PCA
design1b <- model.matrix(~epi + array + sex + group)
 
# fit the linear model 
fit1b <- lmFit(Mval2, design1b)

# Fitting to an estimated Baysian distribution
fit1b <- eBayes(fit1b)

# look at the numbers of DM CpGs at FDR < 0.05
summary(decideTests(fit1b))


## ----annotatedmps-------------------------------------------------------------

# get the table of results for the first contrast 
ann450kSub1b <- anno[match(rownames(Mval2),anno$Name),
                      c(1:4,12:19,24:ncol(anno))]
DMPs1b <- topTable(fit1b, num=Inf, coef=11)
head(DMPs1b)

# ----write, eval=FALSE--------------------------------------------------------
write.table(DMPs1b, file="DMPs1b.csv", sep=",", row.names=FALSE)

# ----figure9, fig.width=10, fig.height=10, fig.cap="Plotting the top few differentially methylated CpGs is a good way to check whether the results make sense."----

# plot the top 4 most significantly differentially methylated CpGs 
par(mfrow=c(2,2))
sapply(rownames(DMPs1b)[1:4], function(cpg){
  plotCpg(beta2, cpg=cpg, pheno=targets$adhd, ylab = "Beta values")
})

## ----cpgannotate--------------------------------------------------------------

myAnnotation1b <- cpg.annotate(object = Mval2, datatype = "array", what = "M", 
                             analysis.type = "differential", design = design1b,
                             coef = 11, arraytype = "EPIC")
str(myAnnotation1b)

## ----ST3GAL3 Subset--------------------------------------------------------------
#Run one that is just the gene bodies and doesn't include the promoter and enhancers
SAM_DMPs1b <- DMPs1b[which(rownames(DMPs1b) %in% SAM_probes_all_list2),]
SAM_DMPs1b$adj.P.Val <- p.adjust(SAM_DMPs1b$P.Value, method = "fdr")
SAM_DMPs1b

write.table(SAM_DMPs1b, file = "SAM_DMPs1b.csv", sep = ",", quote = FALSE)

# plot the top 4 most significantly differentially methylated CpGs 
par(mfrow=c(2,2))
sapply(rownames(SAM_DMPs1b)[1:4], function(cpg){
  plotCpg(beta, cpg=cpg, pheno=targets$adhd, ylab = "Beta values")
})

```
#Symptom Severity
```{r}
targets_adhd <- read.csv("/Users/cogld/Documents/Thesis/epignetics/P3/data/Sample_annotation_adhd.csv")
targets_adhd$id <- as.character(targets_adhd$id)
rownames(targets_adhd) <- targets_adhd$id
targets_adhd$adhd <- ifelse(targets_adhd$group == 'ADHD', 1, 0)


#Match IDs in targets_adhd file to Mval + Beta

#Match Mval and Beta to Targets file 
Mval_adhd <- Mval1[,rownames(targets_adhd)]
beta_adhd <- beta1[,rownames(targets_adhd)]
beta_adhd$cpg <- rownames(beta2)
SAM_Mval_adhd <- SAM_Mval[,rownames(targets_adhd)] 
SAM_beta_adhd <- SAM_beta[,rownames(targets_adhd)]

dim(targets_adhd)
dim(Mval_adhd)
dim(beta_adhd) #162? Come back and check - 163 due to the cpg variable being created
dim(SAM_Mval_adhd)
dim(SAM_beta_adhd)
```
#2. ADHD Symptom Severity predicting DNA methylation
```{r}
## ----dmps---------------------------------------------------------------------

group <- factor(targets_adhd$group)
symp <- targets_adhd$p3_cpadhdi_raw
# this is the individual effect that we need to account for
sex <- factor(targets_adhd$ChildGender) 
epi <- as.numeric((targets_adhd$Epi))
sample_group <- as.factor((targets_adhd$Sample_Group))
array <- as.factor((targets_adhd$Array))

# use the above to create a design matrix
#ask about the order and if it matters where each variable sits
#batch correction - run a PCA
design2 <- model.matrix(~epi + array + sex + symp)
 
# fit the linear model 
fit2 <- lmFit(Mval_adhd, design2)

# Fitting to an estimated Baysian distribution
fit2 <- eBayes(fit2)

# look at the numbers of DM CpGs at FDR < 0.05
summary(decideTests(fit2))


## ----annotatedmps-------------------------------------------------------------

# get the table of results for the first contrast 
ann450kSub2 <- anno[match(rownames(Mval_adhd),anno$Name),
                      c(1:4,12:19,24:ncol(anno))]
DMPs2 <- topTable(fit2, num=Inf, coef=11)
head(DMPs2)

# ----write, eval=FALSE--------------------------------------------------------
write.table(DMPs2, file="DMPs2.csv", sep=",", row.names=FALSE)

# ----figure9, fig.width=10, fig.height=10, fig.cap="Plotting the top few differentially methylated CpGs is a good way to check whether the results make sense."----

# plot the top 4 most significantly differentially methylated CpGs 
par(mfrow=c(2,2))
sapply(rownames(DMPs2)[1:4], function(cpg){
  plotCpg(beta_adhd, cpg=cpg, pheno=targets$adhd, ylab = "Beta values")
})

## ----cpgannotate--------------------------------------------------------------

myAnnotation2 <- cpg.annotate(object = Mval_adhd, datatype = "array", what = "M", 
                             analysis.type = "differential", design = design2,
                             coef = 11, arraytype = "EPIC")
str(myAnnotation2)

## ----ST3GAL3 Subset--------------------------------------------------------------
#Run one that is just the gene bodies and doesn't include the promoter and enhancers
SAM_DMPs2 <- DMPs2[which(rownames(DMPs2) %in% SAM_probes_all_list2),]
SAM_DMPs2$adj.P.Val <- p.adjust(SAM_DMPs2$P.Value, method = "fdr")
SAM_DMPs2

write.table(SAM_DMPs2, file = "SAM_DMP2.csv", sep = ",", quote = FALSE)

# plot the top 4 most significantly differentially methylated CpGs 
par(mfrow=c(2,2))
sapply(rownames(SAM_DMPs1b)[1:4], function(cpg){
  plotCpg(beta, cpg=cpg, pheno=targets_adhd$adhd, ylab = "Beta values")
})

```


