---
title: "Differential expression analysis of female APOE3 cerebral cortex"
Author: Amin Haghani
Date: April 2020
output: html_document
---

Libraries
```{r, message=FALSE}
#if(require(data.table)) install.packages("data.table")
library(data.table)
library(RColorBrewer)
library(ComplexHeatmap)
library(Glimma)
library(readxl)
library(edgeR)
```


Load the files
```{r}
data <- read_xlsx("RNAseq TR-APOE exposure 2012.xlsx")
groups <- read_xlsx("groups.xlsx")
groups$group <- paste(groups$Sex, groups$ApoE, groups$Treatment, sep = "_")


map <- data[,1:7]
geneCount <- data[,-1:-7]
rownames(geneCount) <- make.names(map$`Gene Symbol`, unique = TRUE)

```


```{r convert to edgeR}
library(limma)
library(edgeR)

x <- DGEList(counts = geneCount[, which(groups$Sex=="Female"&groups$ApoE=="E3")], 
             samples = groups[ which(groups$Sex=="Female"&groups$ApoE=="E3"),], 
             group = groups$group[which(groups$Sex=="Female"&groups$ApoE=="E3")])
class(x)
genes <- map
x$genes <- genes
dim(x)
```

# Data pre-processing

```{r cpm and log cpm}
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE)
```


```{r filter, drop the outlier}
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE)


keep.exprs <- rowSums(cpm!=0)>=8 
x <- x[keep.exprs, keep.lib.sizes=FALSE]
dim(x)
```

```{r filterplot1, fig.height=4, fig.width=8}

samplenames <- x$samples$group
nsamples <- ncol(x)
col.group <- factor(groups$group[which(groups$Sex=="Female"&groups$ApoE=="E3")])
par(mfrow=c(1,2))
plot(density(lcpm[,1]), col=col.group, lwd=2, ylim=c(0,0.4), las=2, 
     main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=0, lty=3)
for (i in 2:nsamples){
 den <- density(lcpm[,i])
 lines(den$x, den$y, col=col.group[i], lwd=2)
}
#legend("topright", samplenames, text.col=col, bty="n")
x <- calcNormFactors(x, method = "TMM")
lcpm <- cpm(x, log=TRUE)
plot(density(lcpm[,1]), col=col.group, lwd=2, ylim=c(0,0.4), las=2, 
     main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=0, lty=3)
for (i in 2:nsamples){
   den <- density(lcpm[,i])
   lines(den$x, den$y, col=col.group[i], lwd=2)
}
#legend("topright", samplenames, text.col=col, bty="n")
```

```{r normalize}
x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors
```

```{r normalizemodifieddata}
x2 <- x
x2$samples$norm.factors <- 1
```

```{r plotmodifieddata, fig.height=4, fig.width=8}
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col.group, main="", cex.axis=0.75, names=col.group)
title(main="A. Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2, method = "TMM")  
x2$samples$norm.factors
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col.group, main="", cex.axis=0.75, names=col.group)
title(main="B. Normalised data",ylab="Log-cpm")
```


```{r MDS1, fig.height=4, fig.width=8}
lcpm <- cpm(x, log=TRUE)
par(mfrow=c(1,2))

levels(col.group) <-  brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
plotMDS(lcpm, labels=x$samples$Groups, col=col.group)
title(main="A. Sample groups")
```

```{r design}
airPol <- factor(x$samples$Treatment)
  
design <- model.matrix(~airPol)
colnames(design) <- c("(Intercept)","nPM")
design

```


```{r voom, fig.height=4, fig.width=8}
par(mfrow=c(1,2))
v <- voom(x, design, plot=TRUE)
v
vfit <- lmFit(v, design) 
efit <- eBayes(vfit)
plotSA(efit)
```


```{r decidetests prior to adjustment}
summary(decideTests(efit, p.value = 0.05, lfc = 0, adjust.method = "fdr"))
summary(decideTests(efit, p.value = 0.005, lfc = 0, adjust.method = "none"))
```


```{r treat and look at the sig genes}
tfit <- treat(efit, lfc=0) 
dt <- decideTests(tfit, adjust.method = "none", p.value = 0.005)
summary(dt)
```

```{r MDplot, fig.keep='none'}
plotMD(tfit, column=2, status=dt[,2], main=colnames(tfit)[2], 
       xlim=c(-8,13))
```


```{r treat and look at the sig genes}
tfit <- treat(efit, lfc=0) 
nPM <- topTreat(tfit, coef=2, n=Inf)
head(nPM)
write.csv(nPM, "./Results/Female, E3, nPM/nPM DEG.csv")
```


