---
title: "Specificity of the axonal transcriptome in NT3 and NGF conditions"
author: "Raphaelle Luisier, Idiap Research Institute"
date: "February 2023"
output:
  html_document:
    theme: paper
#paper readable
    highlight: monochrome
    code_folding: hide
    toc: true
    toc_depth: 2
#    toc_float: true
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache =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)
#require(gplots)
#library(Rsamtools)
#library(IRanges)
#library(GenomicRanges)
library(rtracklayer)
library(topGO)
#library(geneplotter)


rm(list = ls())
```

# Load the data {.tabset}
3' UTR isoform quantification from 8 samples originating from rat neuronal sympathetic cell cultures where the cell body is treated with low does of NGF and distal axons are either exposed to NGF or NT3. There are a total of 8 samples: cell body versus distal axons in NGF versus NT3 culture condition across 2 technical replicates.

```{r load_data_load_scripts}
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")
load("./data/files_analysis_Nov2017.RData")

```



# Diversity of mRNA pool in cell body versus axons {.tabset}
We next analyse the reliably expressed 3' UTR isoforms and associated transcript ID in each compartment.
```{r analysis_content,warning=FALSE,fig.width=10, fig.height=3.5}
# A. Characterise the expression in each compartment

#Characterise pipeline in terms of how many new annotated: longer, shorter,...
no.isoforms     <- nrow(anno_tot)
is.I0           <- tapply(anno_tot$isoform,INDEX=anno_tot$txID,function(Z)return(sum(Z==0)==0))
no.new.isoform  <- sum(anno_tot$isoform!=0)#13'314 
no.new.longer   <- sum(anno_tot$newL-anno_tot$initL>=30)
no.new.shorter  <- sum(-anno_tot$newL+anno_tot$initL>=30)
focus.iso       <- c(no.new.isoform,no.new.longer,no.new.shorter)

#Novel 3' UTR isoforms
no.txID.modif     <- sum(tapply(abs(anno_tot$newL-anno_tot$initL),INDEX=factor(as.character(anno_tot$txID)),function(Z)return(sum(Z>30)>0)))#10418 txID modified
no.txID.extended  <- sum(unlist(lapply((anno_tot$maxL-anno_tot$initL)[match(unique(anno_tot$txID),anno_tot$txID)],function(Z)return(Z>30))))#8430 txID extended
no.txID.shortened <- sum(unlist(lapply((anno_tot$initL-anno_tot$minL)[match(unique(anno_tot$txID),anno_tot$txID)],function(Z)return(Z>30))))#4664 txID shortened
no.txID.both      <- sum(
  unlist(lapply((anno_tot$initL-anno_tot$minL)[match(unique(anno_tot$txID),anno_tot$txID)],function(Z)return(Z>30)))&
    unlist(lapply((anno_tot$maxL-anno_tot$initL)[match(unique(anno_tot$txID),anno_tot$txID)],function(Z)return(Z>30))))#2111 txID both
focus.txID       <- c(no.txID.modif,no.txID.extended,no.txID.shortened,no.txID.both)


maxdL.pos        <- (anno_tot$maxL-anno_tot$initL)[match(unique(anno_tot$txID),anno_tot$txID)]
maxdL.pos        <- maxdL.pos[maxdL.pos>0]
maxdL.neg        <- abs(anno_tot$minL-anno_tot$initL)[match(unique(anno_tot$txID),anno_tot$txID)]
maxdL.neg        <- maxdL.neg[maxdL.neg>0]


axons                                              <- unique(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso]) 
cb                                                 <- unique(anno_tot$txID[anno_tot$NGF.cb.is.expressed.iso]) 
neurons                                            <- union(axons,cb)
cb.only                                            <- setdiff(cb,axons)
axon.only                                          <- setdiff(axons,cb)
both                                               <- intersect(axons,cb)
detect.txID.ngf                                    <- rep("no",nrow(anno_tot))
detect.txID.ngf[which(anno_tot$txID%in%both)]      <- "both"
detect.txID.ngf[which(anno_tot$txID%in%cb.only)]   <- "cb.only"
detect.txID.ngf[which(anno_tot$txID%in%axon.only)] <- "axon.only"
anno_tot$detect.txID.ngf                           <- detect.txID.ngf

axons                                              <- unique(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso]) 
cb                                                 <- unique(anno_tot$txID[anno_tot$NT3.cb.is.expressed.iso]) 
neurons                                            <- union(axons,cb)
cb.only                                            <- setdiff(cb,axons)
axon.only                                          <- setdiff(axons,cb)
both                                               <- intersect(axons,cb)
detect.txID.nt3                                    <- rep("no",nrow(anno_tot))
detect.txID.nt3[which(anno_tot$txID%in%both)]      <- "both"
detect.txID.nt3[which(anno_tot$txID%in%cb.only)]   <- "cb.only"
detect.txID.nt3[which(anno_tot$txID%in%axon.only)] <- "axon.only"
anno_tot$detect.txID.nt3                           <- detect.txID.nt3


#uniqueID.expression
axons                                             <- anno_tot$uniqueID[anno_tot$NGF.axon.is.expressed.iso]
cb                                                <- anno_tot$uniqueID[anno_tot$NGF.cb.is.expressed.iso]
neur.neur                                         <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%intersect(cb,axons)&anno_tot$txID%in%both] 
neur.cb                                           <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(cb,axons)&anno_tot$txID%in%both] 
neur.ax                                           <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(axons,cb)&anno_tot$txID%in%both]
ax.ax                                             <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(axons,cb)&anno_tot$txID%in%axon.only]
cb.cb                                             <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(cb,axons)&anno_tot$txID%in%cb.only]

detect.iso.ngf                                    <- rep("no",nrow(anno_tot))
detect.iso.ngf[which(anno_tot$uniqueID%in%neur.neur)] <- "neur.neur"
detect.iso.ngf[which(anno_tot$uniqueID%in%neur.cb)]   <- "neur.cb"
detect.iso.ngf[which(anno_tot$uniqueID%in%neur.ax)]   <- "neur.ax"
detect.iso.ngf[which(anno_tot$uniqueID%in%ax.ax)]     <- "ax.ax"
detect.iso.ngf[which(anno_tot$uniqueID%in%cb.cb)]     <- "cb.cb"
anno_tot$detect.iso.ngf                           <- detect.iso.ngf


axons                                             <- anno_tot$uniqueID[anno_tot$NT3.axon.is.expressed.iso]
cb                                                <- anno_tot$uniqueID[anno_tot$NT3.cb.is.expressed.iso]
neur.neur                                         <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%intersect(cb,axons)&anno_tot$txID%in%both] 
neur.cb                                           <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(cb,axons)&anno_tot$txID%in%both] 
neur.ax                                           <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(axons,cb)&anno_tot$txID%in%both]
ax.ax                                             <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(axons,cb)&anno_tot$txID%in%axon.only]
cb.cb                                             <- as.character(anno_tot$uniqueID)[anno_tot$uniqueID%in%setdiff(cb,axons)&anno_tot$txID%in%cb.only]

detect.iso.nt3                                        <- rep("no",nrow(anno_tot))
detect.iso.nt3[which(anno_tot$uniqueID%in%neur.neur)] <- "neur.neur"
detect.iso.nt3[which(anno_tot$uniqueID%in%neur.cb)]   <- "neur.cb"
detect.iso.nt3[which(anno_tot$uniqueID%in%neur.ax)]   <- "neur.ax"
detect.iso.nt3[which(anno_tot$uniqueID%in%ax.ax)]     <- "ax.ax"
detect.iso.nt3[which(anno_tot$uniqueID%in%cb.cb)]     <- "cb.cb"
anno_tot$detect.iso.nt3                              <- detect.iso.nt3

no.txID.ngf <- c(cb=length(unique(anno_tot$txID[anno_tot$NGF.cb.is.expressed.iso])),
                 axons=length(unique(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso]) ))
no.iso.ngf  <- c(cb=length(unique(anno_tot$uniqueID[anno_tot$NGF.cb.is.expressed.iso])),
                 axons=length(unique(anno_tot$uniqueID[anno_tot$NGF.axon.is.expressed.iso]) ))

no.txID.nt3 <- c(cb=length(unique(anno_tot$txID[anno_tot$NT3.cb.is.expressed.iso])),
                 axons=length(unique(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso]) ))
no.iso.nt3  <- c(cb=length(unique(anno_tot$uniqueID[anno_tot$NT3.cb.is.expressed.iso])),
                 axons=length(unique(anno_tot$uniqueID[anno_tot$NT3.axon.is.expressed.iso]) ))

mycols   <- c(rgb(23/255,71/255,120/255),rgb(119/255,192/255,68/255))
par(mfrow=c(1,4))
mp<- barplot(no.iso.ngf,las=1,col=mycols,cex.axis=0.7)
mtext(side=3,line=0,text=no.iso.ngf,cex=0.7,at=mp)
mtext(side=3,line=1,text="NGF",cex=0.7)
mtext(side=2,line=3,text="# 3' UTR isoforms",cex=0.7)

mp<- barplot(no.iso.nt3,las=1,col=mycols,cex.axis=0.7)
mtext(side=3,line=0,text=no.iso.nt3,cex=0.7,at=mp)
mtext(side=3,line=1,text="NT3",cex=0.7)
mtext(side=2,line=3,text="# 3' UTR isoforms",cex=0.7)

mp<- barplot(no.txID.ngf,las=1,col=mycols,cex.axis=0.7)
mtext(side=3,line=0,text=no.txID.ngf,cex=0.7,at=mp)
mtext(side=3,line=1,text="NGF",cex=0.7)
mtext(side=2,line=3,text="# txID",cex=0.7)



mp<- barplot(no.txID.nt3,las=1,col=mycols,cex.axis=0.7)
mtext(side=3,line=0,text=no.txID.nt3,cex=0.7,at=mp)
mtext(side=3,line=1,text="NT3",cex=0.7)
mtext(side=2,line=3,text="# txID",cex=0.7)

```

We now compare the two neurotrophic conditions and find that although the axonal compartments are similar, there are about 2,000 transcripts that are unique to each condition.

```{r analysis_content_bis,warning=FALSE,fig.width=10, fig.height=3.5}
#Compare NGF and NT3 in terms of the number of overlapping transcript ID per compartment
cb.ngf <- unique(anno_tot$txID[anno_tot$NGF.cb.is.expressed.iso])#13174
cb.nt3 <- unique(anno_tot$txID[anno_tot$NT3.cb.is.expressed.iso])#13324
ax.ngf <- unique(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso])#6536
ax.nt3 <- unique(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso])#5693

comp.cb <- list(intersect(cb.ngf,cb.nt3),setdiff(cb.ngf,cb.nt3),setdiff(cb.nt3,cb.ngf))
comp.ax <- list(intersect(ax.ngf,ax.nt3),setdiff(ax.ngf,ax.nt3),setdiff(ax.nt3,ax.ngf))

#Compare NGF and NT3 in terms of the number of overlapping 3' UTR isoforms per compartment
cb.ngf.i <- unique(anno_tot$uniqueID[anno_tot$NGF.cb.is.expressed.iso])#28824
cb.nt3.i <- unique(anno_tot$uniqueID[anno_tot$NT3.cb.is.expressed.iso])#28813
ax.ngf.i <- unique(anno_tot$uniqueID[anno_tot$NGF.axon.is.expressed.iso])#9474
ax.nt3.i <- unique(anno_tot$uniqueID[anno_tot$NT3.axon.is.expressed.iso])#7893

comp.cb.i <- list(intersect(cb.ngf.i,cb.nt3.i),setdiff(cb.ngf.i,cb.nt3.i),setdiff(cb.nt3.i,cb.ngf.i))
comp.ax.i <- list(intersect(ax.ngf.i,ax.nt3.i),setdiff(ax.ngf.i,ax.nt3.i),setdiff(ax.nt3.i,ax.ngf.i))


#Here I perform slightly different analysis in order to compare what is restricted in CB and those targeted to axons.
cb.ngf <- unique(anno_tot$txID[anno_tot$detect.txID.ngf=="cb.only"])#6775
cb.nt3 <- unique(anno_tot$txID[anno_tot$detect.txID.nt3=="cb.only"])#7'798
ax.ngf <- unique(anno_tot$txID[anno_tot$detect.txID.ngf=="both"])#6'399
ax.nt3 <- unique(anno_tot$txID[anno_tot$detect.txID.nt3=="both"])#5'526



comp.cb <- list(intersect(cb.ngf,cb.nt3),setdiff(cb.ngf,cb.nt3),setdiff(cb.nt3,cb.ngf))
comp.ax <- list(intersect(ax.ngf,ax.nt3),setdiff(ax.ngf,ax.nt3),setdiff(ax.nt3,ax.ngf))

col_treatment <- c("grey","#81A4D6","#AE72B0")
par(mfrow=c(1,4))

mp<-barplot(unlist(lapply(comp.cb.i,length)),las=1,ylim=c(0,30000),col=col_treatment,cex.axis=0.7)
mtext(side=1,line=0,text=c("both","NGF only","NT3 only"),at=mp,cex=0.7)
mtext(side=3,line=0,text=unlist(lapply(comp.cb.i,length)),at=mp,cex=0.7)
mtext(side=2,line=2,text="# transcript ID detected in CB",cex=0.7)

mp<-barplot(unlist(lapply(comp.ax.i,length)),las=1,ylim=c(0,6000),col=col_treatment,cex.axis=0.7)
mtext(side=1,line=0,text=c("both","NGF only","NT3 only"),at=mp,cex=0.7)
mtext(side=3,line=0,text=unlist(lapply(comp.ax.i,length)),at=mp,cex=0.7)
mtext(side=2,line=2,text="# 3' UTR isoforms detected in axons",cex=0.7)

mp<-barplot(unlist(lapply(comp.cb,length)),las=1,ylim=c(0,6000),col=col_treatment,cex.axis=0.7)
mtext(side=1,line=0,text=c("both","NGF only","NT3 only"),at=mp,cex=0.7)
mtext(side=3,line=0,text=unlist(lapply(comp.cb,length)),at=mp,cex=0.7)
mtext(side=2,line=2,text="# transcript ID detected in CB",cex=0.7)

mp<-barplot(unlist(lapply(comp.ax,length)),las=1,ylim=c(0,6000),col=col_treatment,cex.axis=0.7)
mtext(side=1,line=0,text=c("both","NGF only","NT3 only"),at=mp,cex=0.7)
mtext(side=3,line=0,text=unlist(lapply(comp.ax,length)),at=mp,cex=0.7)
mtext(side=2,line=2,text="# transcript ID detected in axons",cex=0.7)




```

# Characterise compartment-specific 3' UTR characteristics {.tabset}
As shown below, and as previously reported in Andreassi et al. (Cell Reports, 2021), the transcripts ID that localise to the axons exhibit a significant higher diversity in 3' UTR isoforms and significant longer 3' UTR length. 
```{r analysis_3UTR_length,warning=FALSE,fig.width=10, fig.height=4}

#Extract the txID of all transcripts detected in axons irrespective of their isoform
Axons.txID.ngf <- unique(as.character(anno_tot$txID)[anno_tot$NGF.axon.is.expressed.iso])
Axons.txID.nt3 <- unique(as.character(anno_tot$txID)[anno_tot$NT3.axon.is.expressed.iso])

#Extract the txID of all transcripts detected in CB irrespective of their isoform
CB.txID.ngf    <- unique(as.character(anno_tot$txID)[anno_tot$NGF.cb.is.expressed.iso])
CB.txID.nt3    <- unique(as.character(anno_tot$txID)[anno_tot$NT3.cb.is.expressed.iso])

#Extract the isoform ID of all transcripts detected in axons irrespective of their isoform
Axons.iso.ngf  <- unique(as.character(anno_tot$uniqueID)[anno_tot$NGF.axon.is.expressed.iso])
Axons.iso.nt3  <- unique(as.character(anno_tot$uniqueID)[anno_tot$NT3.axon.is.expressed.iso])

#Extract the isoform ID of all transcripts detected in CB irrespective of their isoform
CB.iso.ngf     <- unique(as.character(anno_tot$uniqueID)[anno_tot$NGF.cb.is.expressed.iso])
CB.iso.nt3     <- unique(as.character(anno_tot$uniqueID)[anno_tot$NT3.cb.is.expressed.iso])


#Extract the list of these and their differences
L1 <- list(
  cb.only.gf=setdiff(CB.txID.ngf,Axons.txID.ngf),
  cb.axons.gf=intersect(CB.txID.ngf,Axons.txID.ngf),
  axon.only.gf=setdiff(Axons.txID.ngf,CB.txID.ngf),
  cb.only.t3=setdiff(CB.txID.nt3,Axons.txID.nt3),
  cb.axons.t3=intersect(CB.txID.nt3,Axons.txID.nt3),
  axon.only.t3=setdiff(Axons.txID.nt3,CB.txID.nt3)
)

L2 <- list(
  cb.only.gf=setdiff(CB.iso.ngf,Axons.iso.ngf),
  cb.axons.gf=intersect(CB.iso.ngf,Axons.iso.ngf),
  axon.only.gf=setdiff(Axons.iso.ngf,CB.iso.ngf),
  cb.only.t3=setdiff(CB.iso.nt3,Axons.iso.nt3),
  cb.axons.t3=intersect(CB.iso.nt3,Axons.iso.nt3),
  axon.only.t3=setdiff(Axons.iso.nt3,CB.iso.nt3)
)


#Extract maximum 3' UTR length of the txID detected in either axons or CV
maxL.ngf        <- tapply(anno_tot$newL[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso])),FUN=max)
maxL.nt3        <- tapply(anno_tot$newL[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso])),FUN=max)
#Extract initi 3' UTR length of the txID detected in either axons or CV
initL.ngf       <- tapply(anno_tot$initL[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso])),FUN=function(X)return(X[1]))                     
initL.nt3       <- tapply(anno_tot$initL[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso])),FUN=function(X)return(X[1]))

no.iso.ngf      <- tapply(anno_tot$initL[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NGF.axon.is.expressed.iso|anno_tot$NGF.cb.is.expressed.iso])),FUN=function(X)return(length(X))) 
no.iso.nt3      <- tapply(anno_tot$initL[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso],INDEX=factor(as.character(anno_tot$txID[anno_tot$NT3.axon.is.expressed.iso|anno_tot$NT3.cb.is.expressed.iso])),FUN=function(X)return(length(X))) 

anno_tot$maxL.ngf  <- maxL.ngf[match(anno_tot$txID,names(maxL.ngf))]
anno_tot$initL.ngf <- initL.ngf[match(anno_tot$txID,names(initL.ngf))]
anno_tot$no.iso.ngf<- no.iso.ngf[match(anno_tot$txID,names(no.iso.ngf))]

anno_tot$maxL.nt3  <- maxL.nt3[match(anno_tot$txID,names(maxL.nt3))]
anno_tot$initL.nt3 <- initL.nt3[match(anno_tot$txID,names(initL.nt3))]
anno_tot$no.iso.nt3<- no.iso.nt3[match(anno_tot$txID,names(no.iso.nt3))]



with.multiple.ngf <- unlist(lapply(L1,function(Z)return(sum(no.iso.ngf[match(Z,names(no.iso.ngf))]>1)/length(Z))))[-c(3:6)]
with.multiple.nt3 <- unlist(lapply(L1,function(Z)return(sum(no.iso.nt3[match(Z,names(no.iso.nt3))]>1)/length(Z))))[-c(1,2,3,6)]


distr.no.iso.ngf  <- rbind(
  as.vector(table(no.iso.ngf[match(L1[[1]],names(no.iso.ngf))]))[c(1:7)]/length(L1[[1]]),
  (as.vector(table(no.iso.ngf[match(L1[[2]],names(no.iso.ngf))]))[c(1:7)])/length(L1[[2]]))

distr.no.iso.nt3  <- rbind(
  as.vector(table(no.iso.nt3[match(L1[[4]],names(no.iso.nt3))]))[c(1:7)]/length(L1[[4]]),
  (as.vector(table(no.iso.nt3[match(L1[[5]],names(no.iso.nt3))]))[c(1:7)])/length(L1[[5]]))

utrL.ngf          <- list(
  Ensembl=initL.ngf[match(unique(c(Axons.txID.ngf,CB.txID.ngf)),names(initL.ngf))],
  cb.only=maxL.ngf[match(L1[[1]],names(maxL.ngf))],
  axons=maxL.ngf[match(L1[[2]],names(maxL.ngf))]
)

utrL.nt3          <- list(
  Ensembl=initL.nt3[match(unique(c(Axons.txID.nt3,CB.txID.nt3)),names(initL.nt3))],
  cb.only=maxL.nt3[match(L1[[4]],names(maxL.nt3))],
  axons=maxL.nt3[match(L1[[5]],names(maxL.nt3))]
)

par(mfrow=c(1,4),mar=c(3,4,3,3))
mycols <- c(rgb(254/255,218/255,0),rgb(150/255,150/255,150/255),"black")
#mp<- barplot(unlist(lapply(L1,length)),las=1,col=mycols,frame=F)
#mtext(side=3,line=0,text=unlist(lapply(L1,length)),at=mp,cex=0.6)
#mtext(side=2,line=3,text="no.txID")

#mp<- barplot(unlist(lapply(L2,length)),las=1,col=mycols,frame=F)
#mtext(side=3,line=0,text=unlist(lapply(L2,length)),at=mp,cex=0.5)
#mtext(side=2,line=3,text="no.3' UTR isoforms")

#Statistical test to investigate whether the fraction of txID with mult iso is different between CB and AX
Is.axons    <- union(as.character(anno_tot$txID)[anno_tot$NGF.axon.is.expressed.iso],
                     as.character(anno_tot$txID)[anno_tot$NGF.cb.is.expressed.iso])%in%L1[[2]]
Is.multiple <- anno_tot$no.iso.ngf[match(union(as.character(anno_tot$txID)[anno_tot$NGF.axon.is.expressed.iso],
                                               as.character(anno_tot$txID)[anno_tot$NGF.cb.is.expressed.iso]),anno_tot$txID)]>1

pval_NGF <- format(fisher.test(Is.axons,Is.multiple)$p.value,scientific=TRUE,digit=1)


Is.axons    <- union(as.character(anno_tot$txID)[anno_tot$NT3.axon.is.expressed.iso],
                     as.character(anno_tot$txID)[anno_tot$NT3.cb.is.expressed.iso])%in%L1[[4]]
Is.multiple <- anno_tot$no.iso.nt3[match(union(as.character(anno_tot$txID)[anno_tot$NT3.axon.is.expressed.iso],
                                               as.character(anno_tot$txID)[anno_tot$NT3.cb.is.expressed.iso]),anno_tot$txID)]>1
pval_NT3 <- format(fisher.test(Is.axons,Is.multiple)$p.value,scientific=TRUE,digit=1)


mp<-barplot(cbind(with.multiple.ngf*100,with.multiple.nt3*100),col=mycols[c(1,2)],las=1,frame=F,beside=T,ylim=c(0,100))
mtext(side=2,line=3,text="fraction of txID with multiple iso",cex=0.7)
mtext(side=1,line=0,at=mp,text=c("NGF","NT3"),cex=0.7)
legend("top",ncol=2,pch=15,col=mycols[c(1,2)],leg=c("CB","axons"),cex=0.7)
mtext(side=3,line=0,text=c("P(Fisher):",pval_NGF,pval_NT3),at=c(0,2,5),cex=0.7)
arrows(x0=1.0,y0=80,x1=3,y1=80,length = 0)
arrows(x0=4,y0=80,x1=6,y1=80,length = 0)

barplot(distr.no.iso.ngf*100,beside=T,col=mycols[c(1,2)],las=1,,frame=F,main="NGF")
mtext(side=2,line=3,text="fraction of txID with multiple iso",cex=0.7)
barplot(distr.no.iso.nt3*100,beside=T,col=mycols[c(1,2)],las=1,,frame=F,main="NT3")
mtext(side=2,line=3,text="fraction of txID with multiple iso",cex=0.7)


#boxplot(utrL.ngf,outline=F,col=c("white",mycols[c(1,2)]),las=1,frame=F,main="ngf")
#mtext(side=3,line=0,text=unlist(lapply(utrL.ngf,length)),at=c(1,2,3),cex=0.5)

#boxplot(utrL.nt3,outline=F,col=c("white",mycols[c(1,2)]),las=1,frame=F,main="nt3")
#mtext(side=3,line=0,text=unlist(lapply(utrL.nt3,length)),at=c(1,2,3),cex=0.5)

boxplot(c(utrL.nt3,utrL.ngf),outline=F,col=c("white",mycols[c(1,2)]),las=1,frame=F,xaxt="n",ylim=c(0,10000))
mtext(side=3,line=2,text=c(unlist(lapply(utrL.ngf,length)),unlist(lapply(utrL.nt3,length))),at=c(1:6),cex=0.5)
mtext(side=3,line=0,at=c(-1,2.5,5.5),text=c("P(WT):",format(wilcox.test(utrL.ngf[[2]],utrL.ngf[[3]], paired=F)$p.value,scientific=TRUE,digit=1),format(wilcox.test(utrL.nt3[[2]],utrL.nt3[[3]], paired=F)$p.value,scientific=TRUE,digit=1)),cex=0.6)
arrows(x0=1.8,y0=10000,x1=3.4,y1=10000,length = 0)
arrows(x0=4.8,y0=10000,x1=6.4,y1=10000,length = 0)
mtext(side=2,line=3,text="3' UTR length [nt]",cex=0.7)
mtext(side=1,line=0,at=c(1:6),text=rep(c("Ens","CB","AX"),2),cex=0.7)
mtext(side=1,line=1,at=c(2,5),text=c("NGF","NT3"),cex=0.7)


```

# Comparison of the NGF and NT3 axonal transcriptome {.tabset}
Here we characterise the differences in the two compartment in terms of 3' UTR isoform diversity, 3' UTR length and biological pathways.

```{r GO_analysis_tx_uniquely_localised,warning=FALSE,fig.width=10, fig.height=4,eval=FALSE}
t2g                     <- read.csv("./data/t2g_biomaRt.csv")
txID2GO                 <- tapply(t2g$ensembl_transcript_id,INDEX=t2g$go_id,FUN=function(x)return(x))
noNodes                 <- 300

mysampleGO.cb           <- lapply(comp.cb,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))
})

mysampleGO.ax           <- lapply(comp.ax,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))
})

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 = ', '))
}



tardir<- "./zenodo/compart/enrichment/"
Name   <- c("cb.both","cb.only.ngf","cb.only.nt3")
myenrich            <- lapply(c(1:length(mysampleGO.cb)),function(Z){
  
  mysampleGO              <- mysampleGO.cb[[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=""))
  
})

Name   <- c("ax.both","ax.only.ngf","ax.only.nt3")
myenrich            <- lapply(c(1:length(mysampleGO.ax)),function(Z){
  
  mysampleGO              <- mysampleGO.ax[[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=""))
  
})

```

```{r,characterise_axonal_pool,warning=FALSE,fig.width=12, fig.height=8}

par(mfrow=c(2,4))

layout(matrix(c(1,2,3,4,5,5,6,6),nrow=2,ncol=4,byrow = TRUE))
mymat.notx <- rbind(c(length(cb.ngf),length(ax.ngf)),
                    c(length(cb.nt3),length(ax.nt3)))
colnames(mymat.notx) <- c("CB","AX")
rownames(mymat.notx) <- c("NGF","NT3")
mp<- barplot(t(mymat.notx),beside=TRUE,border=unlist(lapply(col_treatment[c(2,3)],function(Z)rep(Z,2))),col=c(col_treatment[2],"white",col_treatment[3],"white"),ylim=c(0,8000),ylab="# expressed txID",las=1)
mtext(side=1,line=0,text=c("CB","AX"),at=as.vector(mp),cex=0.7)
mtext(side=3,line=0,text=as.vector(mymat.notx),at=as.vector(mp),cex=0.7)

venn(list(NT3=ax.nt3,NGF=ax.ngf))

mp<-barplot(cbind(with.multiple.ngf*100,with.multiple.nt3*100),border=unlist(lapply(col_treatment[c(2,3)],function(Z)rep(Z,2))),col=c(col_treatment[2],"white",col_treatment[3],"white"),las=1,frame=F,beside=T,ylim=c(0,100))
mtext(side=2,line=3,text="fraction of txID with multiple iso",cex=0.7)
mtext(side=1,line=0,at=mp,text=c("NGF","NT3"),cex=0.7)
mtext(side=3,line=-2,text=c("P(Fisher):",pval_NGF,pval_NT3),at=c(0,2,5),cex=0.7)
arrows(x0=1.0,y0=80,x1=3,y1=80,length = 0)
arrows(x0=4,y0=80,x1=6,y1=80,length = 0)


boxplot(c(utrL.nt3,utrL.ngf[-1]),outline=F,col=c("white",col_treatment[2],"white",col_treatment[3],"white"),las=1,frame=F,xaxt="n",ylim=c(0,10000),border=c("black","black",col_treatment[2],"black",col_treatment[3]))
mtext(side=3,line=2,text=c(unlist(lapply(utrL.ngf,length)),unlist(lapply(utrL.nt3,length))),at=c(1:6),cex=0.5)
mtext(side=3,line=0,at=c(-1,2.5,5.5),text=c("P(WT):",format(wilcox.test(utrL.ngf[[2]],utrL.ngf[[3]], paired=F)$p.value,scientific=TRUE,digit=1),format(wilcox.test(utrL.nt3[[2]],utrL.nt3[[3]], paired=F)$p.value,scientific=TRUE,digit=1)),cex=0.6)
arrows(x0=1.8,y0=10000,x1=3.4,y1=10000,length = 0)
arrows(x0=4.8,y0=10000,x1=6.4,y1=10000,length = 0)
mtext(side=2,line=3,text="3' UTR length [nt]",cex=0.7)
mtext(side=1,line=0,at=c(1:6),text=rep(c("Ens","CB","AX"),2),cex=0.7)
mtext(side=1,line=1,at=c(2,5),text=c("NGF","NT3"),cex=0.7)

par(mar=c(3,15,3,3))
#Compare the GO terms that are specific to Axons NGF versus NT3
enrich_ngf <- read.csv("./zenodo/compart/enrichment/W0_ax.only.ngf_f.csv")
enrich_nt3 <- read.csv("./zenodo/compart/enrichment/W0_ax.only.nt3_f.csv")
terms_oi <- union(as.character(enrich_ngf$GO.ID[enrich_ngf$Significant>=10&enrich_ngf$weight0Fisher<=0.01]),
                  as.character(enrich_nt3$GO.ID[enrich_nt3$Significant>=10&enrich_nt3$weight0Fisher<=0.01]))

terms_oi <- union(as.character(enrich_ngf$GO.ID[enrich_ngf$Significant>=5&enrich_ngf$weight0Fisher<=0.01]),
                  as.character(enrich_nt3$GO.ID[enrich_nt3$Significant>=5&enrich_nt3$weight0Fisher<=0.01]))
enrich_ngf <- read.csv("./zenodo/compart/enrichment/W0_ax.only.ngf.csv")
enrich_nt3 <- read.csv("./zenodo/compart/enrichment/W0_ax.only.nt3.csv")
ix_ngf <- match(terms_oi,enrich_ngf$GO.ID)
ix_nt3 <- match(terms_oi,enrich_nt3$GO.ID)
temp                   <- cbind(as.character(enrich_ngf[ix_ngf,2]),as.character(enrich_ngf[ix_ngf,3]))
temp[!is.na(ix_nt3),1] <- as.character(enrich_nt3[ix_nt3[!is.na(ix_nt3)],2])
temp[!is.na(ix_nt3),2] <- as.character(enrich_nt3[ix_nt3[!is.na(ix_nt3)],3])
temp                   <- data.frame(temp,NGF=ifelse(!is.na(ix_ngf),-log10(enrich_ngf$weight0Fisher[ix_ngf]),0))
temp                   <- data.frame(temp,NT3=ifelse(!is.na(ix_nt3),-log10(enrich_nt3$weight0Fisher[ix_nt3]),0))
vec_ngf         <- -log10(enrich_ngf$weight0Fisher[enrich_ngf$Significant>=5&enrich_ngf$weight0Fisher<=0.01])
names(vec_ngf)  <- as.character(enrich_ngf$Term[enrich_ngf$Significant>=5&enrich_ngf$weight0Fisher<=0.01])
vec_ngf         <- sort(vec_ngf,decreasing=TRUE)
vec_nt3         <- -log10(enrich_nt3$weight0Fisher[enrich_nt3$Significant>=5&enrich_nt3$weight0Fisher<=0.01])
names(vec_nt3)  <- as.character(enrich_nt3$Term[enrich_nt3$Significant>=5&enrich_nt3$weight0Fisher<=0.01])
vec_nt3         <- sort(vec_nt3,decreasing=TRUE)
barplot(rev(vec_ngf[c(1:5)]),col="#81A4D6",las=1,horiz=TRUE,cex.names=1.0)
barplot(rev(vec_nt3[c(1:5)]),col="#AF72B0",las=1,horiz=TRUE,cex.names=1.0)
```

# Dominant factors underlying axonal localisation {.tabset}
As shown below, the fraction of 3' UTR reliably expressed in the axonal compartment depends on the expression in the cell body, supporting the hypothesis that the axonal localisation of mRNA is determined by both specific regulators enabling active transports as well as transcriptional level as highly expressed transcripts are more likely to be detected in axons. This does not mean that the latter are passively transported but rather that these are actively transported to the axons through a machiney which 'randomly' recruit mRNA to actively transport these to axons. 


Notably the transcript length is another key factor that impacts the likelyhood to be detected in the axonal compartment; the longer the transcript the less likely it will be localised. 
```{r analysis_fraction_detected_per_range_expression,warning=FALSE,fig.width=8, fig.height=4}
#Analysis of the dependendence between what is detected and 
colsNGF <- c("#81A4D6","#2D598E","#083872")
colsNT3 <- c("#AE73B1","#79387C","#57055B")
#Calculate per range of expression in the cell body, the fraction of transcript detected in the axonal compartment

avg.cb.ngf        <- apply(log2(1+cbind(anno_tot$NGF.cb.1.raw.corrected,anno_tot$NGF.cb.2.raw.corrected)),1,mean)
avg.cb.nt3        <- apply(log2(1+cbind(anno_tot$NT3.cb.1.raw.corrected,anno_tot$NT3.cb.2.raw.corrected)),1,mean)
avg.ngf.axons     <- apply(log2(1+cbind(anno_tot$NGF.axon.1.raw.corrected,anno_tot$NGF.axon.2.raw.corrected)),1,mean)
avg.nt3.axons     <- apply(log2(1+cbind(anno_tot$NT3.axon.1.raw.corrected,anno_tot$NT3.axon.2.raw.corrected)),1,mean)
par(mfrow=c(1,2))
myLims            <- seq(from=0,by=1.5,to=20)
bins.ngf          <- cut(avg.cb.ngf,myLims,include.lowest = FALSE)
bins.nt3          <- cut(avg.cb.nt3,myLims,include.lowest = FALSE)
frac_detected_ngf <- unlist(tapply(X=anno_tot$detect.txID.ngf,INDEX=bins.ngf,FUN=function(A)return(sum(A=="both")/length(A))))
frac_detected_nt3 <- unlist(tapply(X=anno_tot$detect.txID.nt3,INDEX=bins.nt3,FUN=function(A)return(sum(A=="both")/length(A))))
myx <-unlist(lapply(c(2:length(myLims)),function(IX)return(mean(myLims[c(IX,IX-1)]))))
plot(myx,frac_detected_ngf,las=1,col=colsNGF[1],pch=19,cex=0.7,frame=FALSE,ylab="fraction of detected tX in axons",xlab="expression in cb [log2]",ylim=c(0,1),xlim=c(0,20))
lines(myx,frac_detected_ngf,las=1,col=colsNGF[1])
points(myx,frac_detected_nt3,las=1,col=colsNT3[1],pch=19,cex=0.7)
lines(myx,frac_detected_nt3,las=1,col=colsNT3[1])
legend("bottomright",col=c(colsNGF[1],colsNT3[1]),pch=19,leg=c("NGF","NT3"))


myLims            <- seq(from=2,by=0.25,to=4)
bins.ngf          <- cut(log10(anno_tot$txLength),myLims,include.lowest = FALSE)
bins.nt3          <- cut(log10(anno_tot$txLength),myLims,include.lowest = FALSE)
frac_detected_ngf <- unlist(tapply(X=anno_tot$detect.txID.ngf,INDEX=bins.ngf,FUN=function(A)return(sum(A=="both")/length(A))))
frac_detected_nt3 <- unlist(tapply(X=anno_tot$detect.txID.nt3,INDEX=bins.nt3,FUN=function(A)return(sum(A=="both")/length(A))))
myx <-unlist(lapply(c(2:length(myLims)),function(IX)return(mean(myLims[c(IX,IX-1)]))))
plot(myx,frac_detected_ngf,las=1,col=colsNGF[1],pch=19,cex=0.7,frame=FALSE,ylab="fraction of detected tX in axons",xlab="txlength [log10]",ylim=c(0.4,0.8),xlim=c(2,4))
lines(myx,frac_detected_ngf,las=1,col=colsNGF[1])
points(myx,frac_detected_nt3,las=1,col=colsNT3[1],pch=19,cex=0.7)
lines(myx,frac_detected_nt3,las=1,col=colsNT3[1])
```









