---
title: "Differential Gene Expression Analysis"
author: "Raphaëlle Luisier"
date: "Idiap Research Institute -- February 2023"
output:
  github_document:
    toc: true
    toc_depth: 3
    dev: jpeg
---

```{r setup, include=FALSE}
working_dir = "~/Documents/Scripts/AxonLoc/"
setwd(working_dir)

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(fig.width=16, fig.height=10, figure.align = "center") 
library(ggplot2) # devtools::install_github('hadley/ggplot2')
library(RColorBrewer)
library(wesanderson)
library(knitr)
library(mclust)
#library(limma)
library(gplots)
#library(lme4)
#library(nlme)
warning = FALSE
library("GO.db")
require("limma")
library("topGO")
#require("biomaRt")
#require("org.Rn.eg.db")
#library(grDevices)
library(Rsamtools)
library(IRanges)
library(GenomicRanges)
library(rtracklayer)
library(geneplotter)
require("multtest")
require("mclust")
library("edgeR")

rm(list = ls())
## 

```

```{r load_data,echo=FALSE}

load("./data/myfiles.RData")#"anno_tot","myUTR"
#myUTR        <- import.gff("./data/myUTR.sorted.gtf",format="gtf")
#anno_tot     <- read.delim("./data/anno_tot.tab",header=TRUE,sep="\t")
names(myUTR) <- as.character(myUTR$ID)
source("./scripts/nested_fun_differential_PAS.R")

```

```{r perform_differential_GE,echo=FALSE}
#A. Select only those reliably expressed in CB
coloi                <- c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected","NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")
data.con.cb          <- anno_tot[anno_tot$NT3.cb.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso,match(coloi,colnames(anno_tot))]
row_info             <- anno_tot[anno_tot$NT3.cb.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso,]

#B. Sum-up all reads from same tX
data.sum             <- apply(data.con.cb,2,function(Z)return(tapply(Z,INDEX=factor(as.character(row_info$txID)),FUN=sum)))
row_info             <- row_info[match(rownames(data.sum),row_info$txID),]

#C. Remove redundant tX for same GS selecting only those which have longest 3' UTR
temp              <- as.vector(row_info$maxL)
names(temp)       <- as.character(row_info$txID)
txOI              <- tapply(temp,INDEX=factor(as.character(row_info$geneSymbol)),function(Z)return(names(Z)[which(Z==max(Z))[1]]))
mysel             <- which(row_info$txID%in%txOI)
data.sum          <- data.sum[mysel,]
row_info          <- row_info[mysel,]

#D. edgeR normalisation (all together)
count.matrix          <- data.sum
group                 <- factor(c("NGF","NGF","NT3","NT3"),levels=c("NGF","NT3"))
my.dgelist            <- DGEList(counts = count.matrix,group = group)


#E.  Calculate normalization factors; stored in my.dgelist$samples
my.dgelist  <- calcNormFactors(my.dgelist)

#F. Estimate the common dispersion; stored in my.dgelist$common.disp
my.dgelist  <- estimateCommonDisp(my.dgelist)

#G. Estimate tagwise dispersions
my.dgelist  <- estimateTagwiseDisp(my.dgelist)

#H. The square root of the common dispersion gives the coefficient of variation of biological variation (BCV).
bcv_samples<-round(sqrt(my.dgelist$common.disp),digits=2)


#I. Estimate mean expression
TD               <- group
design           <- model.matrix(~0+TD)
colnames(design) <- levels(TD)
#The GLM likelihood ratio test is based on the idea of fitting negative binomial GLMs with the Cox-Reid dispersion estimates.
#Given raw counts, a fixed value for the dispersion parameter and a design matrix, the function glmFit() fits the negative binomial GLM for each tag and produces an object of class DGEGLM with some new components
fit              <- glmFit(my.dgelist,design)
et               <- exactTest(my.dgelist,pair=c("NGF","NT3"))#Will compare NT3 versus NGF
detags           <- rownames(topTags(et)$table)
de               <- decideTestsDGE(et, p=0.1, adjust="BH")
detags           <- rownames(my.dgelist)[as.logical(de)]
res              <- et$table
res              <- data.frame(res,cpm(my.dgelist)[, order(my.dgelist$samples$group)])
res$GS           <- as.character(row_info$geneSymbol)[match(as.character(rownames(res)),as.character(row_info$txID))]

#Create output of differentially expressed genes
sel1           <- abs(res$logFC)>log2(1.5)
sel1.1         <- res$logFC>log2(1.5)
sel1.2         <- res$logFC<(log2(1/1.5))
sel2           <- res$PValue<=0.01
sel3           <- res$is.diff
sel3           <- rep(TRUE,length(sel2))
selt           <- sel1&sel2&sel3
selt.1         <- sel1.1&sel2&sel3
selt.2         <- sel1.2&sel2&sel3
```

We find 320 up-regulated genes in NT3 compared to NGF and 385 up-regulated genes in NGF compared to NT3.


```{r continue,echo=FALSE}
#sum(selt.1)#320 NT3>NGF
#sum(selt.2)#385 NGF>NT3

Col            <- rep(rgb(0.2,0.2,0.2,0.4),nrow(res))
Col[sel2]      <- "black"
Col[selt]      <- "red"
#res            <- res[selt,]
colnames(res)[1]<- "logFC[NT3-NGF]"
sub <- cpm(fit$fitted.values)
sub <- sub[match(rownames(res),rownames(sub)),]
colnames(sub)<- paste(colnames(sub),"[fitted CPM]",sep="")
init <- cpm(fit$counts)
init <- init[match(rownames(res),rownames(init)),]
colnames(init)<- paste(colnames(init),"[raw CPM]",sep="")
res <- data.frame(row_info$geneSymbol[match(rownames(res),row_info$txID)],res,row_info[match(rownames(res),row_info$txID),],init,sub)


#J. GO enrichment analysis
t2g                     <- read.csv("./data/t2g_biomaRt.csv")
txID2GO                 <- tapply(t2g$ensembl_transcript_id,INDEX=t2g$go_id,FUN=function(x)return(x))
noNodes                 <- 300

GetGOI <- function(sampleGO=temp.sampleGO,mygoID=as.character(temp$GO.ID)[1]){
  entrez.oi       <- intersect(genesInTerm(sampleGO, mygoID)[[1]],sigGenes(sampleGO))
  gs.oi           <- unique(as.character(t2g$external_gene_name[match(entrez.oi,t2g$ensembl_transcript_id)]))
  return(paste(gs.oi, collapse = ', '))
}

mysampleGO.diff           <- lapply(list(res$txID[selt.1],res$txID[selt.2]),function(Z){
  geneNames              <- unique(t2g$ensembl_transcript_id)
  geneList               <- factor(as.integer(geneNames %in% Z))
  names(geneList)        <- geneNames
  return(new("topGOdata",description = "Simple session", ontology = "BP",allGenes = geneList, geneSel = Z,nodeSize = 10,annot = annFUN.GO2genes,GO2gene=txID2GO))
})



tardir <- "./data/enrichment/"
dir.create(tardir)
Name   <- c("up.DEG.NT3","up.DEG.NGF")
myenrich            <- lapply(c(1:length(mysampleGO.diff)),function(Z){
  mysampleGO              <- mysampleGO.diff[[Z]]
  resultFisher            <- runTest(mysampleGO, algorithm = "classic", statistic = "fisher")
  resultFisher.weight01   <- runTest(mysampleGO, algorithm = "weight01", statistic = "fisher")
  temp                    <- GenTable(mysampleGO, classicFisher = resultFisher,weight0Fisher=resultFisher.weight01,orderBy = "classicFisher", ranksOf = "classicFisher", topNodes = noNodes)
  temp                    <- temp[temp$Significant>5,]
  temp                    <- data.frame(temp,GOI=unlist(lapply(as.character(temp$GO.ID),function(Z)return(GetGOI(sampleGO=mysampleGO,mygoID=Z)))))
  write.csv(temp,file=paste(tardir,"Fisher_",Name[Z],".csv",sep=""))
  
  temp                    <-GenTable(mysampleGO, classicFisher = resultFisher,weight0Fisher=resultFisher.weight01,orderBy = "weight0Fisher", ranksOf = "weight0Fisher", topNodes = noNodes)
  temp                    <- temp[temp$Significant>5,]
  temp                    <- data.frame(temp,GOI=unlist(lapply(as.character(temp$GO.ID),function(Z)return(GetGOI(sampleGO=mysampleGO,mygoID=Z)))))
  write.csv(temp,file=paste(tardir,"W0_",Name[Z],".csv",sep=""))
})


mysampleGO.diff           <- lapply(list(as.character(res.filt.nt3$txID),as.character(res.filt.ngf$txID)),function(Z){
  geneNames              <- unique(t2g$ensembl_transcript_id)
  geneList               <- factor(as.integer(geneNames %in% Z))
  names(geneList)        <- geneNames
  return(new("topGOdata",description = "Simple session", ontology = "BP",allGenes = geneList, geneSel = Z,nodeSize = 10,annot = annFUN.GO2genes,GO2gene=txID2GO))
})


Name   <- c("up.DEG.NT3.f","up.DEG.NGF.f")
myenrich            <- lapply(c(1:length(mysampleGO.diff)),function(Z){
  mysampleGO              <- mysampleGO.diff[[Z]]
  resultFisher            <- runTest(mysampleGO, algorithm = "classic", statistic = "fisher")
  resultFisher.weight01   <- runTest(mysampleGO, algorithm = "weight01", statistic = "fisher")
  temp                    <- GenTable(mysampleGO, classicFisher = resultFisher,weight0Fisher=resultFisher.weight01,orderBy = "classicFisher", ranksOf = "classicFisher", topNodes = noNodes)
  temp                    <- temp[temp$Significant>5,]
  temp                    <- data.frame(temp,GOI=unlist(lapply(as.character(temp$GO.ID),function(Z)return(GetGOI(sampleGO=mysampleGO,mygoID=Z)))))
  write.csv(temp,file=paste(tardir,"Fisher_",Name[Z],".csv",sep=""))
  
  temp                    <-GenTable(mysampleGO, classicFisher = resultFisher,weight0Fisher=resultFisher.weight01,orderBy = "weight0Fisher", ranksOf = "weight0Fisher", topNodes = noNodes)
  temp                    <- temp[temp$Significant>5,]
  temp                    <- data.frame(temp,GOI=unlist(lapply(as.character(temp$GO.ID),function(Z)return(GetGOI(sampleGO=mysampleGO,mygoID=Z)))))
  write.csv(temp,file=paste(tardir,"W0_",Name[Z],".csv",sep=""))
})




#Are these longer for those up-regulated in NGF (selt.2) compared to NT3 (selt.1) ?
myL<-lapply(list(res$txID[selt.1],res$txID[selt.2]),function(Z)return(row_info$maxL[match(Z,row_info$txID)]))
wilcox.test(myL[[1]],myL[[2]],paired = FALSE)
	Wilcoxon rank sum test with continuity correction
data:  myL[[1]] and myL[[2]]
W = 49815, p-value = 1.202e-05
alternative hypothesis: true location shift is not equal to 0



par(mfrow=c(2,2),mar=c(3,3,3,3))
boxplot(myL,outline=F,las=1,main="",frame=F,ylab="3' UTR length",xaxt="n",col=rev(c("#81A4D6","#AF71AF")))
mtext(side=1,at=c(1,2),line=0,text=c("NT3","NGF"),cex=0.7)
Col <- rep("grey",length(res$PValue))
Col[res$PValue<0.01]<- "black"
Col[res$logFC.NT3.NGF.<(-log2(1.5))&res$PValue<0.01]<-"#81A4D6"
Col[res$logFC.NT3.NGF.>(log2(1.5))&res$PValue<0.01]<-"#AF71AF"
plot(x=res$logFC,y=-log10(res$PValue),col=Col,pch=19,cex=0.05,xlab="",ylab="",xlim=c(-4,4),frame=F,las=1)
mtext(side=2,line=2,text="-log10(P-Value)",cex=0.7)
mtext(side=1,line=2,text="log2FC",cex=0.7)
abline(h=-log10(0.01),col="grey",lwd=1,lty=2)
abline(v=log2(1.5),col="grey",lwd=1,lty=2)
abline(v=-log2(1.5),col="grey",lwd=1,lty=2)
text(-3,15,paste("n(NGF>NT3)=",sum(selt.2),sep=""),cex=0.5)
text(3,15,paste("n(NGF<NT3)=",sum(selt.1),sep=""),cex=0.5)
par(mfrow=c(2,2),mar=c(5,4,2.5,1),cex=0.7)
plotMDS(my.dgelist, method="bcv",xlab="",ylab="")
mtext(side=1,line=2,text="BCV distance 1")
mtext(side=2,line=2,text="BCV distance 2")
plotBCV(my.dgelist,xlab="",ylab="")
mtext(side=3,line=0.5,text=paste("BVC=",bcv_samples,sep=""))
mtext(side=1,line=2,text="Average log CPM",cex=0.7)
mtext(side=2,line=2,text="Biological Coefficient of Variation",cex=0.7)
plotSmear(et, de.tags=detags,xlab="",ylab="")
mtext(side=1,line=2,text="Average logCPM",cex=0.7)
mtext(side=2,line=2,text="logFC",cex=0.7)
abline(h = c(-1, 1), col = "blue")



#H. Plot example genes

coloi                <- c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected","NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")
data.con.cb          <- anno_tot[anno_tot$NT3.cb.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso,match(coloi,colnames(anno_tot))]
row_info             <- anno_tot[anno_tot$NT3.cb.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso,]
data.sum             <- apply(data.con.cb,2,function(Z)return(tapply(Z,INDEX=factor(as.character(row_info$txID)),FUN=sum)))
data.sum.log2        <- log2(1+data.sum)
plotGenesPoints   <- function(mytx="ENSRNOT00000022880",raw=data.sum.log2){
  
  grp   <- factor(c("NGF","NGF","NT3","NT3"),levels=c("NGF","NT3")) 
  mycols <- rep(c("#81A4D6","#AF72B0"),2)
  myX   <- c(1,3)
  
  
  all  <-tapply(raw[match(mytx,rownames(raw)),],INDEX=grp,FUN=function(X)return(X))
  mym  <-tapply(raw[match(mytx,rownames(raw)),],INDEX=grp,FUN=function(X)return(mean(X)))
  
  plot(xlim=c(0,max(myX)+0.5),ylim=c(floor(min(unlist(all))),ceiling(max(unlist(all)))),x=rep(myX[1],length(all[[1]])),y=all[[1]],col=mycols[1],pch=19,frame=F,cex=0.8,ylab="",xlab="",las=1,xaxt="n")
  #plot(xlim=c(0,max(myX)+0.5),ylim=c(8,16),x=rep(myX[1],length(all[[1]])),y=all[[1]],col=mycols[1],pch=19,frame=F,cex=0.8,ylab="",xlab="",las=1,xaxt="n")
  
  arrows(x0=myX[1]-0.5,y0=mym[1], x1=myX[1]+0.5, angle=90, code=3, length=0,col=mycols[1])
  for(i in c(2:length(myX))){
    points(x=rep(myX[i],length(all[[i]])),y=all[[i]],col=mycols[i],pch=19,cex=0.8)
    arrows(x0=myX[i]-0.5,y0=mym[i], x1=myX[i]+0.5, angle=90, code=3, length=0,col=mycols[i])
  }
  mtext(side=2,line=3,text="read count [log2]",cex=0.7)
  mtext(side=1,line=0,text=c("NGF","NT3"),at=c(1.5,3.5),cex=0.7)  
  mtext(side=3,line=0,text=mytx,cex=0.7)
  mtext(side=3,line=1,text=as.character(anno_tot$geneSymbol)[match(mytx,anno_tot$txID)],cex=0.7)
}

error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
  if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}



diff.tx.NGF <- c("ENSRNOT00000025362","ENSRNOT00000022880","ENSRNOT00000025372","ENSRNOT00000005866","ENSRNOT00000066932","ENSRNOT00000018718","ENSRNOT00000007698","ENSRNOT00000026009","ENSRNOT00000051757","ENSRNOT00000016265")

diff.tx.NT3 <- c("ENSRNOT00000022943","ENSRNOT00000022943","ENSRNOT00000004320","ENSRNOT00000022635","ENSRNOT00000024900","ENSRNOT00000007246","ENSRNOT00000030344","ENSRNOT00000072040","ENSRNOT00000068102","ENSRNOT00000022276","ENSRNOT00000025221","ENSRNOT00000044483","ENSRNOT00000018697","ENSRNOT00000003633","ENSRNOT00000076589","ENSRNOT00000050549","ENSRNOT00000027047","ENSRNOT00000017453","ENSRNOT00000043085")

par(mfrow=c(3,4),mar=c(3,3,4,2))
lapply(diff.tx.NGF,function(myTX)plotGenesPoints(mytx=myTX,raw=data.sum.log2))
par(mfrow=c(3,4),mar=c(2,4,4,2))
lapply(diff.tx.NT3,function(myTX)plotGenesPoints(mytx=myTX,raw=data.sum.log2))


```