---
title: "Differential expression analysis of combined wildtype and APOE-TR datasets"
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)
```


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

data2 <- read_xlsx("~/Google Drive (ahaghani@usc.edu)/Amin documents/Finch lab projects/My Phd project/Research projects/GSM RNAseq/GSM rnaseq individual sample.xlsx")
groups2 <- read_xlsx("~/Google Drive (ahaghani@usc.edu)/Amin documents/Finch lab projects/My Phd project/Research projects/GSM RNAseq/groups.xlsx", sheet = 2)

total <- merge(data2, data[,c(5,8:39)], by='Gene Symbol', all=TRUE)
total <- total[complete.cases(total),]

groups$drug <- "vehicle"
groups$batch <- 2
groups2$batch <- 1
groups <- rbind(groups2, groups)


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

```{r group description}
table(groups$Sex, groups$ApoE, groups$Treatment, groups$batch)
```


```{r convert to edgeR}
library(limma)
library(edgeR)
x <- DGEList(counts = geneCount[,which(groups$drug=="vehicle")], samples = groups[which(groups$drug=="vehicle"),], group = groups$group[which(groups$drug=="vehicle")])
class(x)
dim(x)
genes <- map
x$genes <- genes
dim(genes)
```

# Data pre-processing

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

```{r zeroes}
table(rowSums(cpm!=0)==48)
```

__Analysis of the common genes__

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

keep.exprs <- rowSums(cpm!=0)>=48 
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(x$samples$group)
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 design}
airPol <- factor(x$samples$Treatment)

female <- factor(x$samples$Sex)
levels(female) <- c(1,0)
female <- factor(female, levels = c(0,1), labels = c("Male","Female"))

apoe <- factor(x$samples$ApoE)
batch <- factor(x$samples$batch)

  
design <- model.matrix(~female+apoe+airPol)
colnames(design) <- c("(Intercept)","Female", "E4", "Wt","nPM")
design
```


```{r voom combat batch correction, fig.height=4, fig.width=8}
library(sva)
par(mfrow=c(1,2))
batch = unclass(x$samples$batch)
v <- voom(x, design, plot=TRUE)
v.bc = ComBat(dat=v$E, batch=batch, mod=design[,5])
mds<-plotMDS(v.bc,labels=x$samples$group, ndim=3)

vfit <- lmFit(v.bc, design) #this is least square method

efit <- eBayes(vfit)
plotSA(efit)

norm.data <- cbind(map[rownames(v.bc),],v.bc)
#write.csv(norm.data, "./Results/wt+apoe normalized data.csv")
```


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

vennDiagram(dt[,c(3,5)], circle.col=c("turquoise", "salmon", "red"), include=c("up", "down"), counts.col=c("red", "blue"))

apoe4 <- topTreat(tfit, coef=3, n=Inf)
apoe4 <- cbind(map[rownames(apoe4),], apoe4)

nPM <- topTreat(tfit, coef=5, n=Inf)
nPM <- cbind(map[rownames(nPM),], nPM)

```


