4_Axonal_Localisation_Scoring.Rmd

---
title: "Modelling of axonal locasation"
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}
#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)
library("multtest")
library("mclust")
library("plotly")
library("roll")
library("fitdistrplus")
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,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")
source("./scripts/nested_fun_nt3_ngf.R")
source("./scripts/nested_fun_APA_CB.R")
#myUTR        <- import.gff("./data/myUTR.sorted.gtf",format="gtf")
load("./zenodo/files_analysis_Nov2017.RData")#anno_tot: 35534:52

source("./scripts/nested_fun_differential_PAS.R")
source("./scripts/nested_fun_nt3_ngf.R")
diff.transport.neurotrophins<-read.csv("./zenodo/transport/diff_transport.csv")

colNGF <- c("#81A4D6",rgb(129/255,164/255,214/255,0.25))
colNT3 <- c("#AE72B0",rgb(174/255,114/255,176/255,0.25))
load("./zenodo/transport/regulation/welch.RData")
load("./zenodo/transport/regulation/fisher.RData")
```


We first need to select the reliably expressed genes and transcripts i.e. those which are reliably expressed in the cell body

```{r selection_reliably_expressed_isoforms,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=10, fig.height=3}
par(mfrow=c(2,4))
coloi <- c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected","NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected","NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected","NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")
mydat                <- anno_tot[which(anno_tot$is.conservative),match(coloi,colnames(anno_tot))]
mydatall             <- anno_tot[,match(coloi,colnames(anno_tot))]
for(i in c(1:ncol(mydat))){
  mydat[,i] <- log2(mydat[,i]+1)
}
lims      <- apply(mydat,2,function(Z)2^SelectExpressed(dat=Z,frac.bg=0.95,frac.fg=0.1))

tempsel   <- apply(mydatall,2,function(Z)return(Z>=mean(lims[1,c(1,2)])))
soft.sel  <- cbind(tempsel[,1]&tempsel[,2],tempsel[,3]&tempsel[,4],tempsel[,5]&tempsel[,6],tempsel[,7]&tempsel[,8])
tempsel   <- do.call(what=cbind,lapply(c(1:8),function(Z)return(mydatall[,Z]>=lims[2,Z])))
hard.sel  <- cbind(tempsel[,1]|tempsel[,2],tempsel[,3]|tempsel[,4],tempsel[,5]|tempsel[,6],tempsel[,7]|tempsel[,8])
final.sel <- cbind(soft.sel[,1]|hard.sel[,1],soft.sel[,2]|hard.sel[,2],soft.sel[,3]|hard.sel[,3],soft.sel[,4]|hard.sel[,4])

anno_tot$NGF.axon.is.expressed.iso          <- final.sel[,1]
anno_tot$NGF.cb.is.expressed.iso            <- final.sel[,2]
anno_tot$NT3.axon.is.expressed.iso          <- final.sel[,3]
anno_tot$NT3.cb.is.expressed.iso            <- final.sel[,4]

#Select only those reliably expressed in CB of both NGF and NT3
coloi                <- grep(colnames(anno_tot),pattern=".raw.corrected")
seloi                <- anno_tot$NT3.cb.is.expressed.iso&anno_tot$NGF.cb.is.expressed.iso#If we want to study differential transport, we need to make sure that things are expressed at least in both NGF and NT3 CB
#seloi <- c(1:nrow(anno_tot))
data.con.cb          <- anno_tot[seloi,coloi]
row_info             <- anno_tot[seloi,]
anno_tot             <- anno_tot[seloi,]
mydatall             <- log2(1+data.con.cb)
myUTR                <- myUTR[seloi,]

```
 
```{r compute_relevant_values_for_downstream_analysis}
#B. Compute average expression in each compartment
avg.ngf.axons <- apply(mydatall[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected")],1,mean)
avg.nt3.axons <- apply(mydatall[,c("NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)
avg.ngf.cb    <- apply(mydatall[,c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.nt3.cb    <- apply(mydatall[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")],1,mean)
avg.cb        <- apply(mydatall[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected","NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.ax        <- apply(mydatall[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected","NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)

avg.ngf.axons.raw <- apply(data.con.cb[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected")],1,mean)
avg.nt3.axons.raw <- apply(data.con.cb[,c("NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)
avg.ngf.cb.raw    <- apply(data.con.cb[,c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.nt3.cb.raw    <- apply(data.con.cb[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")],1,mean)
```
 
# Differential transport using classic abundance ratios
```{r difference_abundance_ratios,eval=FALSE}
#A. Compute relative log ratio
RUD.ngf  <- avg.ngf.axons-avg.ngf.cb
RUD.nt3  <- avg.nt3.axons-avg.nt3.cb
diff.RUD <- RUD.ngf-RUD.nt3


#B. Compute relative usage (PUD has larger dynamic range than SUD; with SUD lowly expressed genes in CB have higher chances of getting high score so better PUD which detects changes across the entire spectrum of CB expression level)
PUD.ngf <- avg.ngf.axons/(avg.ngf.cb+avg.ngf.axons)
PUD.nt3 <- avg.nt3.axons/(avg.nt3.cb+avg.nt3.axons)
PUD.nt3[(avg.nt3.cb+avg.nt3.axons)==0] <-0
PUD.ngf[(avg.ngf.cb+avg.ngf.axons)==0] <-0
diff.PUD <- PUD.ngf-PUD.nt3


SUD.ngf <- avg.ngf.axons.raw/(avg.ngf.cb.raw+avg.ngf.axons.raw)
SUD.nt3 <- avg.nt3.axons.raw/(avg.nt3.cb.raw+avg.nt3.axons.raw)
SUD.nt3[(avg.nt3.cb.raw+avg.nt3.axons.raw)==0] <-0
SUD.ngf[(avg.ngf.cb.raw+avg.ngf.axons.raw)==0] <-0
diff.SUD <- SUD.ngf-SUD.nt3


LUD.ngf <- avg.ngf.axons.raw/avg.ngf.cb.raw
LUD.nt3 <- avg.nt3.axons.raw/avg.nt3.cb.raw
LUD.nt3[avg.nt3.cb.raw==0] <-0
LUD.ngf[avg.ngf.cb.raw==0] <-0
diff.LUD <- LUD.ngf/SUD.nt3



#C. Compute Fisher count P-value --> P-value on differential transport between NGF and NT3 while previous analysis is about transport efficiency
mytempdat    <- data.con.cb
sum.dat      <- t(apply(mytempdat,1,function(Z)return(tapply(Z,INDEX=c("NGF.axon","NGF.axon","NGF.cb","NGF.cb","NT3.axon","NT3.axon","NT3.cb","NT3.cb"),FUN=sum))))
p.val.diff   <- apply(sum.dat,1,function(Z)return(fisher.test(round(cbind(Z[c(1,2)],Z[c(3,4)])))$p.value))
FDR.diff     <-  p.adjust(p.val.diff,method="fdr")

#D. Prepare output
PUD                 <- data.frame(PUD.ngf=PUD.ngf,PUD.nt3=PUD.nt3)
RUD                 <- data.frame(RUD.ngf=RUD.ngf,RUD.nt3=RUD.nt3)
SUD                 <- data.frame(SUD.ngf=SUD.ngf,SUD.nt3=SUD.nt3)
dPUD.rel            <- PUD$PUD.ngf*(avg.ngf.cb/max(avg.ngf.cb))-PUD$PUD.nt3*(avg.nt3.cb/max(avg.nt3.cb))
dSUD.rel            <- SUD$SUD.ngf*(avg.ngf.cb/max(avg.ngf.cb))-SUD$SUD.nt3*(avg.nt3.cb/max(avg.nt3.cb))

myOut               <- data.frame(anno_tot[,c("txID","geneSymbol","uniqueID")],FDR=FDR.diff,dRUD=diff.RUD,dPUD=diff.PUD,dSUD=diff.SUD,dPUD.rel=dPUD.rel,dSUD.rel=dSUD.rel,PUD,SUD,
                                  PUD.ngf.rel=PUD$PUD.ngf*(avg.ngf.cb/max(avg.ngf.cb)),
                                  SUD.ngf.rel=SUD$SUD.ngf*(avg.ngf.cb/max(avg.ngf.cb)),
                                  PUD.nt3.rel=PUD$PUD.nt3*(avg.nt3.cb/max(avg.nt3.cb)),
                                  SUD.nt3.rel=SUD$SUD.nt3*(avg.nt3.cb/max(avg.nt3.cb)),
                                  dRUD=diff.RUD,RUD,data.con.cb)


#G. Extract differentially transported txID based on RUD and PUD
#Test for overlap with Repetitive regions -- test whether this is required
myRpM                           <- import.gff("./zenodo/rmsk_rn5_ucsc.sorted.gff",format="gff")
myRpM                           <- reduce(myRpM,min.gapwidth=20)
POS                             <- as.character(strand(myUTR))=="+"
myUTR.focus                     <- myUTR
start(myUTR.focus)[POS]         <- end(myUTR)[POS]-499
end(myUTR.focus)[!POS]          <- start(myUTR)[!POS]+499
gOver                           <- findOverlaps(query=myUTR.focus,subject=myRpM,ignore.strand=FALSE)
idxU                            <- queryHits(gOver)
idxG                            <- subjectHits(gOver)
new.start                       <- apply(cbind(start(myUTR.focus)[idxU],start(myRpM)[idxG]),1,max)
new.end                         <- apply(cbind(end(myUTR.focus)[idxU],end(myRpM)[idxG]),1,min)
my.overlapping.utr              <- myUTR.focus[idxU,]
start(my.overlapping.utr)       <- new.start
end(my.overlapping.utr)         <- new.end
mywidth                         <- tapply(width(my.overlapping.utr),INDEX=factor(as.character(my.overlapping.utr$ID)),FUN=sum)
oi.rpm                          <- names(mywidth)[mywidth>=300]

#save(list="myOut",file="./zenodo/transport/diff_transport_ratios/Differential_abundance_ratios.RData")

sel.A.up  <- myOut$dRUD>1
sel.A.do  <- myOut$dRUD<(-1)
sel.B.up  <- myOut$dPUD.rel>0.10
sel.B.do  <- myOut$dPUD.rel<(-0.10)
sel.C     <- myOut$FDR<0.01
sel.D     <- !(myOut$uniqueID%in%oi.rpm)

sel.more.NGF <- sel.A.up&sel.B.up&sel.D&sel.C&anno_tot$NGF.axon.is.expressed.iso&data.con.cb[,1]>20&data.con.cb[,2]>20
sel.more.NT3 <- sel.A.do&sel.B.do&sel.D&sel.C&anno_tot$NT3.axon.is.expressed.iso&(data.con.cb[,5]>20|data.con.cb[,6]>20)#More laxist on NT3 since it is not as covered as the NGF

txID.more.NGF <- unique(as.character(anno_tot$txID)[sel.more.NGF])#643
txID.more.NT3 <- unique(as.character(anno_tot$txID)[sel.more.NT3])#322

#write.csv(myOut[sel.more.NGF,],"./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
#write.csv(myOut[sel.more.NT3,],"./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")
```

```{r load_abundance_ratios_analysis}
more_NGF <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
more_NT3 <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")
load("./zenodo/transport/diff_transport_ratios/Differential_abundance_ratios.RData")
```

 
# Development of the axonal localisation score {.tabset}

## Selection of linear model that best fit the data {.tabset}
We next want to test whether the read count in the axonal compartment depends on the read count in the cell body as well as on transcript length (log10). We do this by fitting a linear model of iteratively increasing order for NGF and NT3 and then extract AIC.

Below is the result for the NGF condition:

```{r, lm_NGF,echo=TRUE}
# C. Modelling axonal transport in function of 1) read counts in CB and 2) tx length
# With the linear model --using both size and expression in CB -- the aim is to identify those which are not behaving in the same way as most other species of the same size 
# Linear model to measure the extend of deviation from the linear model

dat.ngf           <- data.frame(txL=log10(anno_tot$txLength[match(names(avg.ngf.axons),anno_tot$uniqueID)]),
                                cb=avg.ngf.cb,
                                axons=avg.ngf.axons)
dat.nt3           <- data.frame(txL=log10(anno_tot$txLength[match(names(avg.nt3.axons),anno_tot$uniqueID)]),
                                cb=avg.nt3.cb,
                                axons=avg.nt3.axons)

dat.ngf           <- data.frame(txL=log10(anno_tot$txLength),
                                cb=avg.ngf.cb,
                                axons=avg.ngf.axons)
dat.nt3           <- data.frame(txL=log10(anno_tot$txLength),
                                cb=avg.nt3.cb,
                                axons=avg.nt3.axons)

polyfit <- function(i) x <- AIC(lm(axons~poly(cb,i)+txL,dat=dat.ngf))

fit1.1 <- lm(axons~cb,dat=dat.ngf)
fit1.2 <- lm(axons~cb+txL,dat=dat.ngf)
fit2.1 <- lm(axons~poly(cb,2),dat=dat.ngf)
fit2.2 <- lm(axons~poly(cb,2)+txL,dat=dat.ngf)
fit3.1 <- lm(axons~poly(cb,3),dat=dat.ngf)
fit3.2 <- lm(axons~poly(cb,3)+txL,dat=dat.ngf)
fit4.1 <- lm(axons~poly(cb,4),dat=dat.ngf)
fit4.2 <- lm(axons~poly(cb,4)+txL,dat=dat.ngf)#selected one
fit5.1 <- lm(axons~poly(cb,5),dat=dat.ngf)
fit5.2 <- lm(axons~poly(cb,5)+txL,dat=dat.ngf)
print(AIC(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2))
print(anova(fit4.2,fit5.2,test="Chisq"))
print(anova(fit4.2,fit4.1,test="Chisq"))
```

Below is the result for the NT3 condition.
```{r, lm_NT3,echo=TRUE}

fit1.1 <- lm(axons~cb,dat=dat.nt3)
fit1.2 <- lm(axons~cb+txL,dat=dat.nt3)
fit2.1 <- lm(axons~poly(cb,2),dat=dat.nt3)
fit2.2 <- lm(axons~poly(cb,2)+txL,dat=dat.nt3)
fit3.1 <- lm(axons~poly(cb,3),dat=dat.nt3)
fit3.2 <- lm(axons~poly(cb,3)+txL,dat=dat.nt3)
fit4.1 <- lm(axons~poly(cb,4),dat=dat.nt3)
fit4.2 <- lm(axons~poly(cb,4)+txL,dat=dat.nt3)#selected one
fit5.1 <- lm(axons~poly(cb,5),dat=dat.nt3)
fit5.2 <- lm(axons~poly(cb,5)+txL,dat=dat.nt3)
print(AIC(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2))
print(anova(fit4.2,fit5.2,test="Chisq"))
print(anova(fit4.2,fit4.1,test="Chisq"))


myfits <- list(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2)
#par(mfrow=c(5,2))
#lapply(myfits,function(Z){
#  pdt     <- predict(Z)
#  plot(avg.ngf.cb,avg.ngf.axons,pch=19,col=rgb(0,0,0,0.1),cex=0.2,main=,cex.main=0.7)
#  points(avg.ngf.cb,pdt,pch=19,col=rgb(1,0,0,0.1),cex=0.2)
#})
```

The selected model is polynom of degree 4 with transcript length. Here we check that the predictions (red dots) are close to the observed read counts for the selected model.

```{r relation_CB_axons_model,echo=TRUE,warning=FALSE,fig.width=6, fig.height=3}
par(mfrow=c(1,1))
Z=myfits[[8]]
pdt     <- predict(Z)
plot(avg.ngf.cb,avg.ngf.axons,pch=19,col=rgb(0,0,0,0.1),cex=0.2,cex.main=0.7,las=1,frame=FALSE,main="NGF")
points(avg.ngf.cb,pdt,pch=19,col=rgb(1,0,0,0.1),cex=0.2)
```

Interestingly, transcripts that are highly expressed in cell bodies are not necessarily enriched in axons or dendrites, indicating that at least a subset of RNAs do not reach the peripheral compartments by 'passive' transport but are sorted and delivered with an selective mechanism.

## Compute the localisation score  {.tabset}
We next use use a Bayesian approach to get a localisation score with positive values indicating over-transport and negative values indicating under-transport.  We first group 3' UTR isoforms of similar expression level in the CB and of similar transcript length together. Assuming that isoforms within a group should exhibit similar expression in the cell body, of mean and standard-deviation derived from the linear model learned in the previous step. For each 3' UTR isoform, we calculate the likelihood of that the observed axonal read count are higher or lower than the predicted read counts given the prior distribution.

```{r transport_locasation,echo=TRUE,eval=FALSE}
#Extract the probability of each data-point being generated to identify those which are under- or over-transported

#A. Fit polynomial of degree 4 to each data-set
fit4.2.ngf      <- lm(axons~poly(cb,4)+txL,dat=dat.ngf)#selected one
fit4.2.nt3      <- lm(axons~poly(cb,4)+txL,dat=dat.nt3)#selected one

#B. Model prediction 
pdt.ngf         <- predict(fit4.2.ngf)
pdt.nt3         <- predict(fit4.2.nt3)
plot(pdt.ngf,pdt.nt3,pch=19,col=rgb(0,0,0,0.2),cex=0.2)


RUD.ngf         <- avg.ngf.axons-avg.ngf.cb
RUD.nt3         <- avg.nt3.axons-avg.nt3.cb
out.ngf         <- data.frame(dat.ngf,prediction=pdt.ngf,RUD=RUD.ngf,anno_tot[match(names(avg.ngf.axons),anno_tot$uniqueID),],avg.ngf.cb,avg.ngf.axons)
out.nt3         <- data.frame(dat.nt3,prediction=pdt.nt3,RUD=RUD.nt3,anno_tot[match(names(avg.nt3.axons),anno_tot$uniqueID),],avg.nt3.cb,avg.nt3.axons)

#D. Create intervals of expression according to expression in CB and transcript length 
CreateBins <- function(myout){
  bins.cb                        <- cut(myout$cb,breaks=seq(from=3,to=20,by=1),iclude.lowest=T)
  bins.txL                       <- cut(myout$txL,breaks=seq(from=2,to=4.5,by=0.25),iclude.lowest=T)
  bins.both                      <- paste(bins.cb,bins.txL,sep=".")
  temp                           <- data.frame(table(bins.both))
  tokeep                         <- temp$bins.both[temp$Freq>=5] 
  #tokeep                         <- temp$bins.both[temp$Freq>=1] 
  bins.both[!bins.both%in%tokeep]<-NA
  return(factor(as.character(bins.both)))
}
bins.both.ngf <- CreateBins(myout=out.ngf)
bins.both.nt3 <- CreateBins(myout=out.nt3)

#E.Get the probability of the observed read count in the axons larger or smaller than the observed read counts in axons given their observed read counts in the CB and their transcript length

GetProbs <- function(mybins=bins.both.ngf,myfit=fit4.2.ngf,out=out.ngf){
  
  myOut <- do.call(lapply(c(1:length(levels(mybins))),function(IX){
    print(IX)
    myprediction         <- compute.prediction(IX,myfit=myfit,bins=mybins)
    subdat               <- out[which(mybins==levels(mybins)[IX]),]
    if(nrow(subdat)>=10){
      #Fit prior normal distribution of the read counts in the axons with prior mean and sd from prediction given the intervals
      fit.norm    <- fitdist(subdat$axons, "norm",start=list(mean=mean(myprediction),sd=sd(myprediction)))
      #Sample from the prior
      my.rdm.norm <- rnorm(n=10^4, mean=fit.norm$estimate["mean"], sd=fit.norm$estimate["sd"])
      #Compute the posterior prob. i.e. empirical P-value
      temp        <- do.call(what=rbind,lapply(subdat$axons,function(Z)return(c(sum(my.rdm.norm<=Z)/10^4,sum(my.rdm.norm>=Z)/10^4))))
      colnames(temp)<- c("prob(N).smaller","prob(N).bigger")
      temp          <-data.frame(temp,subdat)
      return(temp)
    }
    else{
      temp <- matrix(NA,ncol=2,nrow=nrow(subdat))
      colnames(temp)<- c("prob(N).smaller","prob(N).bigger")
      temp          <-data.frame(temp,subdat)
      return(temp)
    }
    }),what=rbind)
  
  
  p.value.norm <- apply(cbind(-log10(myOut$prob.N..smaller+min(myOut$prob.N..smaller[myOut$prob.N..smaller>0],na.rm=TRUE)),
                              -log10(myOut$prob.N..bigger+min(myOut$prob.N..bigger[myOut$prob.N..bigger>0],na.rm=TRUE))),1,function(Z){
                                M<- max(Z)
                                if(is.na(M))return(NA)
                                if(which(Z==M)[1]==1)return(-M)
                                else
                                  return(M)
                              })

  myOut <- data.frame(p.value.norm,myOut)
  return(myOut)
}


myOut.ngf <- GetProbs(mybins=bins.both.ngf,myfit=fit4.2.ngf,out=out.ngf)
myOut.nt3 <- GetProbs(mybins=bins.both.nt3,myfit=fit4.2.nt3,out=out.nt3)
ix        <- match(rownames(myOut.nt3),rownames(myOut.ngf))
ix        <- ix[!is.na(ix)]
myOut.ngf <- myOut.ngf[ix,]

ix        <- match(rownames(myOut.ngf),rownames(myOut.nt3))
ix        <- ix[!is.na(ix)]
myOut.nt3 <- myOut.nt3[ix,]



transport <- data.frame(uniqueIDs=rownames(myOut.ngf),
                                           GS=as.character(anno_tot$geneSymbol[match(rownames(myOut.ngf),anno_tot$uniqueID)]),
                                           txID=as.character(myOut.ngf$txID),
                                           diff.transport=myOut.ngf$p.value.norm-myOut.nt3$p.value.norm,
                                           diff.RUD=myOut.ngf$RUD-myOut.nt3$RUD,
                                           transport.ngf=myOut.ngf$p.value.norm,
                                           transport.nt3=myOut.nt3$p.value.norm,
                                           RUD.ngf=myOut.ngf$RUD,
                                           RUD.nt3=myOut.nt3$RUD,
                                           cb.ngf=myOut.ngf$cb,
                                           cb.nt3=myOut.nt3$cb,
                                           axons.ngf=myOut.ngf$axons,
                                           axons.nt3=myOut.nt3$axons,
                                           anno_tot$NGF.cb.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NGF.cb.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NGF.axon.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NGF.axon.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NT3.cb.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NT3.cb.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NT3.axon.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
                                           anno_tot$NT3.axon.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)])
transport$txLength         <- anno_tot$txLength[match(rownames(myOut.ngf),anno_tot$uniqueID)]


transport             <- transport[!is.na(transport$transport.ngf),]
mygr                 <- myUTR[match(transport$uniqueIDs,myUTR$ID),]
transport$coord      <- paste(seqnames(mygr),":",start(mygr),"-",end(mygr),sep="")
tempax <- do.call(what=rbind,args=lapply(c(1:nrow(transport)),function(ix){
  sel=transport$txID==transport$txID[ix]&transport$wL>transport$wL[ix]
  return(c(ifelse(sum(sel)>0,max(transport$transport.ngf[sel]),NA),ifelse(sum(sel)>0,max(transport$transport.nt3[sel]),NA)))
}))

transport$tranport.ngf.longer=tempax[,1]
transport$tranport.nt3.longer=tempax[,2]
transport$diff_transport_longer=tempax[,1]-tempax[,2]
#write.csv(transport,""./zenodo/transport/diff_localisationscore/transport.csv")
```

Let's now check that the axonal localisation score does not depend on the read count in cell body nor on the transcript length. 

```{r stats_localisation_scores,echo=TRUE,fig.width=10, fig.height=10}

transport <- read.csv("./zenodo/transport/diff_localisationscore/transport.csv")
#transport <- read.csv("./data/transport.csv")


PlotWithAVG <- function(myxval=transport$cb.ngf,myyval=transport$transport.ngf,myxlab="avg.ngf,cb [log2]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20)){
  BY=0.5
  plot(myxval,myyval,pch=19,col=rgb(0,0,0,0.1),cex=0.2,xlab="",ylab="",las=1,main="",frame=F,xlim=myxlim,ylim=myylim,cex.lab=0.7,cex.axis=0.7)
  grid()
  myY<-tapply(myyval,cut(myxval,breaks=seq(from=min(myxval),to=max(myxval),by=BY),include.lowest=T),FUN=function(W)return(mean(W,na.rm=TRUE)))
  myX<- seq(from=min(myxval),to=max(myxval),by=BY)[-1]
  lines(myX,myY,col="red")
  mtext(side=1,line=2,text=myxlab,cex=0.7)
  mtext(side=2,line=3,text=myylab,cex=0.7)
}

layout(matrix(c(1,1,2,2,3,4,5,6,7,8,9,10),ncol=4,nrow=3,byrow = TRUE))

hist(transport$transport.ngf,col=rgb(129/255,164/255,214/255,0.5),breaks=20,las=1,cex.lab=0.7,cex.axis=0.7,xlim=c(-4,4),main="",xlab="localisation score [NGF]")
hist(transport$transport.nt3,col=rgb(175/255,113/255,175/255,1.0),breaks=20,las=1,cex.lab=0.7,cex.axis=0.7,xlim=c(-4,4),main="",xlab="localisation score [NT3]")


PlotWithAVG(myxval=transport$cb.ngf,myyval=transport$RUD.ngf,myxlab="avg.ngf.cb [log2]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(0,20))
PlotWithAVG(myxval=transport$cb.ngf,myyval=transport$transport.ngf,myxlab="avg.ngf,cb [log2]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20))


PlotWithAVG(myxval=transport$cb.nt3,myyval=transport$RUD.nt3,myxlab="avg.cb [log2] in NT3",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(0,20))
PlotWithAVG(myxval=transport$cb.nt3,myyval=transport$transport.nt3,myxlab="avg.cb [log2] in NT3",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20))




PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$RUD.ngf,myxlab="tx length [log10]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(2.5,4))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$transport.ngf,myxlab="tx length [log10]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(2.5,4))

PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$RUD.nt3,myxlab="tx length [log10]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(2.5,4))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$transport.nt3,myxlab="tx length [log10]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(2.5,4))



```


# Differential localisation analysis  {.tabset}

We can now identify those 3' UTR isoforms which are over-transported in NGF and NT3.

```{r identify_over_under_transported,echo=TRUE,eval=FALSE}

#A. More transported in NGF axons (avoid selecting candidate of axonal remodeling i.e. transport longer should as well be in the same direction)
#Should be more transported in NGF samples
sel_diff <- transport$diff.transport>1.0
#Should be transport or over-transported in NGF
sel_NGF <- transport$transport.ngf>=(1.0)
#Should be under-transported in NT3
sel_NT3 <- transport$transport.nt3<=(0.0)
#If longer, then should be in the same range of differential transport otherwise risk to be axonal remodelling candidate
sel_rem <- is.na(transport$tranport.ngf.longer)#|transport$diff_transport_longer>1.0
#selection for coverage
sel_expression <- transport$anno_tot.NGF.axon.1.raw.corrected.match.rownames.myOut.ngf...>=40&transport$anno_tot.NGF.axon.2.raw.corrected.match.rownames.myOut.ngf...>=40
more_transported_NGF <- data.frame(coord=transport$coord,transport)[sel_diff&sel_NGF&sel_NT3&sel_rem,]#482


#B. More transported in NT3 axons
#Should be more transported in NT3 samples
sel_diff <- transport$diff.transport<=(-1.0)
#Should be transport or over-transported in NT3
sel_NGF <- transport$transport.nt3>=(1.0)
#Should be under-transported in NT3
sel_NT3 <- transport$transport.ngf<=(0.0)
#If longer, then should be in the same range of differential transport otherwise risk to be axonal remodelling candidate
sel_rem <- is.na(transport$tranport.nt3.longer)#|transport$diff_transport_longer>1.0
#selection for coverage
sel_expression <- transport$anno_tot.NT3.axon.1.raw.corrected.match.rownames.myOut.ngf...>=40&transport$anno_tot.NT3.axon.2.raw.corrected.match.rownames.myOut.ngf...>=40
more_transported_NT3 <- data.frame(coord=transport$coord,transport)[sel_diff&sel_NGF&sel_NT3&sel_rem,]


```


```{r compare_diff_transport_methods,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=9, fig.height=6}
#LS
more_transported_NGF <- read.csv("./zenodo/transport/diff_localisationscore/more_transport_NGF.csv")
more_transported_NT3 <- read.csv("./zenodo/transport/diff_localisationscore/more_transport_NT3.csv")

#Relative Abundance Ratios
more_NGF <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
more_NT3 <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")


par(mfrow=c(2,2),mar=c(4,4,4,4))
colsNGF <- c("#81A4D6","#2D598E","#083872")
colsNT3 <- c("#AE73B1","#79387C","#57055B")
mycols <- rep(rgb(0,0,0,0.1),nrow(transport))
mycols[transport$uniqueIDs%in%more_transported_NT3$uniqueIDs]<- colsNT3[1]
mycols[transport$uniqueIDs%in%more_transported_NGF$uniqueIDs]<- colsNGF[1]
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],transport$diff.transport,cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on Localisation Score")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff transport efficiency",cex=0.7)
grid()
text(x=10,y=5,col=colsNGF[1],lab=paste("#over in NGF=",nrow(more_transported_NGF)),cex=0.7)
text(x=10,y=-5,col=colsNT3[1],lab=paste("#over in NT3=",nrow(more_transported_NT3)),cex=0.7)

plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],myOut$dRUD[match(transport$uniqueIDs,myOut$uniqueID)],cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on Localisation Score")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff abundance ratios",cex=0.7)
grid()


mycols <- rep(rgb(0,0,0,0.1),nrow(transport))
mycols[transport$uniqueIDs%in%more_NT3$uniqueID]<- colsNT3[1]
mycols[transport$uniqueIDs%in%more_NGF$uniqueID]<- colsNGF[1]
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],transport$diff.transport,cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on abundance ratios")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff transport efficiency",cex=0.7)
grid()
text(x=10,y=5,col=colsNGF[1],lab=paste("#over in NGF=",nrow(more_NGF)),cex=0.7)
text(x=10,y=-5,col=colsNT3[1],lab=paste("#over in NT3=",nrow(more_NT3)),cex=0.7)

plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],myOut$dRUD[match(transport$uniqueIDs,myOut$uniqueID)],cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on abundance ratiose")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff abundance ratios",cex=0.7)
grid()

```

# Pathway enrichment analysis
```{r GO_analysis,eval=FALSE}

#
#GO on transported BPs compared to random -- on continuous values
#
#
require("topGO")
t2g                <- read.csv("./data/t2g_biomaRt.csv")
x                  <- org.Rn.egGO
mapped_genes       <- mappedkeys(x)
GO2geneID          <- as.list(x[mapped_genes])
geneID2GO          <- as.list(org.Rn.egGO2EG)
goterms            <- Term(GOTERM)
geneNames          <- names(GO2geneID)
myGS               <- anno_tot$geneSymbol
eid                <- t2g$entrezgene[match(myGS,t2g$external_gene_name)]
eid                <- eid[!is.na(eid)]
geneList           <- factor(as.numeric(geneNames%in%as.character(eid)))
names(geneList)    <- geneNames
mysampleGO         <- new("topGOdata",description = "Simple session", ontology = "BP",allGenes = geneList, geneSel = as.character(eid),nodeSize = 10,annot = annFUN.GO2genes,GO2gene=geneID2GO)
myterms            <- data.frame(goID=names(goterms),Term=unlist(goterms))
size.terms         <- unlist(lapply(names(goterms),function(goID){
  if(length(genesInTerm(mysampleGO, goID))==0){return(0)}
  else{
    return(length(genesInTerm(mysampleGO, goID)[[1]]))
  }
}))
goterms <- goterms[size.terms>20]



myMat.rdm           <- list(replicate(1000,myOut.ngf$p.val.max[sample(c(1:nrow(myOut.ngf)))]),
                            replicate(1000,myOut.ngf$p.value.norm[sample(c(1:nrow(myOut.ngf)))]),
                            replicate(1000,myOut.ngf$p.value.bn[sample(c(1:nrow(myOut.ngf)))]),
                            replicate(1000,myOut.ngf$error[sample(c(1:nrow(myOut.ngf)))]),
                            replicate(1000,myOut.nt3$p.val.max[sample(c(1:nrow(myOut.nt3)))]),
                            replicate(1000,myOut.nt3$p.value.norm[sample(c(1:nrow(myOut.nt3)))]),
                            replicate(1000,myOut.nt3$p.value.bn[sample(c(1:nrow(myOut.nt3)))]),
                            replicate(1000,myOut.nt3$error[sample(c(1:nrow(myOut.nt3)))]))
                            
                            
myVec               <- list(myOut.ngf$p.val.max,
                            myOut.ngf$p.value.norm,
                            myOut.ngf$p.value.bn,
                            myOut.ngf$error,
                            myOut.nt3$p.val.max,
                            myOut.nt3$p.value.norm,
                            myOut.nt3$p.value.bn,
                            myOut.nt3$error)
names(myVec[[1]])<-  names(myVec[[2]])<- names(myVec[[3]])<- names(myVec[[4]])<-  myOut.ngf$geneSymbol 
names(myVec[[5]])<-  names(myVec[[6]])<- names(myVec[[7]])<- names(myVec[[8]])<-  myOut.nt3$geneSymbol 


temp                <- do.call(lapply(c(1:length(goterms)),function(Z){
  print(Z)
  return(unlist(lapply(c(1:length(myVec)),function(IX)return(GetEnrichGO(goID=names(goterms)[Z],myvec=myVec[[IX]],myMat=myMat.rdm[[IX]])))))
}),what=rbind)

enrich.GO.transport <- data.frame(names(goterms),goterms,temp)
colnames(enrich.GO.transport)[c(3:10)]<- c("p.val.max.ngf","p.value.norm.ngf","p.value.bn.ngf","error.ngf",
                                           "p.val.max.nt3","p.value.norm.nt3","p.value.bn.nt3","error.nt3")

enrich.GO.transport                <- enrich.GO.transport[,c(1,2,3,7,4,8,5,9,6,10)]
enrich.GO.transport$diff.pval.max  <- enrich.GO.transport$p.val.max.ngf-enrich.GO.transport$p.val.max.nt3
enrich.GO.transport$diff.pval.norm <- enrich.GO.transport$p.val.norm.ngf-enrich.GO.transport$p.val.norm.nt3
enrich.GO.transport$diff.pval.bn   <- enrich.GO.transport$p.val.bn.ngf-enrich.GO.transport$p.val.bn.nt3
enrich.GO.transport$diff.error     <- enrich.GO.transport$error.ngf-enrich.GO.transport$error.nt3
#write.csv(enrich.GO.transport,file="./zenodo/transport/enrich/enrich_transport.csv")

#Difference in transported GO between NGF and NT3
temp                <- do.call(lapply(c(1:length(goterms)),function(Z){
  print(Z)
  return(GetEnrichDiff(goID=names(goterms)[Z]))
}),what=rbind)

enrich.GO.diff.transport          <- data.frame(names(goterms),goterms,temp)
colnames(enrich.GO.diff.transport)<- c("GO.ID","Term","my.diff.norm","zscore.norm","my.diff.error","no.terms")
#write.csv(enrich.GO.diff.transport,file="./zenodo/transport/enrich/enrich_diff_transport.csv")
```


```{r bp_over_transported,,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=9, fig.height=6}
enrich.GO.transport<- read.csv("./zenodo/transport/enrich/enrich_transport.csv")
enrich.GO.diff.transport <- read.csv("./zenodo/transport/enrich/enrich_diff_transport.csv")

#To check whether the p-value makes sense, we check the biological pathways which are best transported, or more transported than expected. And it seems to work pretty weel :-)
myid                    <-"GO:0030574"
myVec.nt3               <- myOut.nt3$p.value.norm
names(myVec.nt3)        <- myOut.nt3$geneSymbol 
myVec.ngf               <- myOut.ngf$p.value.norm
names(myVec.ngf)        <- myOut.ngf$geneSymbol 

PlotEnrichGO(goID=myid,myvec=myVec[[1]],name="-log10(P-value) transport")
PlotEnrichGO(goID=myid,myvec=myVec[[5]],name="-log10(P-value) transport")



                            
                            
```

# Regulation of transport
Here we integrate with CLIP-sequencing data to look into the potential of 126 RBPs to regulate axonal localisation.

```{r load_CLIP_data}
myranges         <- list(c(100,-100),c(3000,2950),c(2750,2700),c(2500,2450),c(2250,2200),c(2000,1950),c(1750,1700),c(1500,1450),c(1400,1350),c(1300,1250),c(1200,1150),c(1100,1050),c(1000,950),c(900,850),c(800,750),c(700,650),c(600,550),c(500,450),c(450,400),c(400,350),c(350,300),c(300,250),c(250,200),c(200,150),c(150,100),c(100,50),c(50,0),c(0,-50),c(-50,-100),c(-100,-150))
mynames_rages   <- unlist(lapply(myranges,function(Z)return(paste("[",Z[[1]],":",Z[[2]],"]",sep=""))))

load("./zenodo/clips.RData")#"myClips","myranges","mynames_rages"
#List of 30 objects == for 30 regions along 3' UTR region, we collect the number of cross-link events per 3' UTR isoforms. Matrix [37513 X 335] with 37,513 isoforms and 335 clip data
#

no_rbp_bound            <- matrix(nrow=nrow(myClips[[1]]),ncol=length(myranges),0)
rownames(no_rbp_bound)  <- rownames(myClips[[1]])
colnames(no_rbp_bound)  <- mynames_rages

for(IX in c(1:length(mynames_rages))){
  clipdat          <- myClips[[IX]]
  no_rbp_bound[,IX]  <- apply(clipdat>0,1,sum)   
  print(IX)
}

```



## RBPome underlying cell body restriction
In order to identify negative regulators of axonal localisation, we focus on the region [-250;-50] preceding the 3' end and perform Fisher count enrichment test to look for RBPs exhibiting significant enrichment in cross-link regions in those 3' UTR isoforms exhibiting restricted axonal localisation. 



```{r prepare_data_for_fisher,eval=TRUE}
#diff.transport.neurotrophins<-read.csv("~/Desktop/DataAnalsyisRiccio/AxonalTransport/diff_transport.csv")
source("./scripts/nested_fun_nt3_ngf.R")

q_NGF           <- quantile(diff.transport.neurotrophins$transport.ngf,seq(from=0,by=0.1,to=1.0))
q_NT3           <- quantile(diff.transport.neurotrophins$transport.nt3,seq(from=0,by=0.1,to=1.0))
myVal           <- list(ifelse(diff.transport.neurotrophins$transport.ngf>q_NGF[10],1,0),
                        ifelse(diff.transport.neurotrophins$transport.ngf<q_NGF[2],1,0),
                        ifelse(diff.transport.neurotrophins$transport.nt3>q_NT3[10],1,0),
                        ifelse(diff.transport.neurotrophins$transport.nt3<q_NT3[2],1,0),
                        ifelse(diff.transport.neurotrophins$uniqueIDs%in%more_transported_NGF$uniqueIDs,1,0),
                        ifelse(diff.transport.neurotrophins$uniqueIDs%in%more_transported_NT3$uniqueIDs,1,0))



#clip_1   <- read.csv("/Users/raphaelle/Desktop/DataAnalsyisRiccio/AxonalTransport/[50:0]eCLIP_hepg2.csv")
#no_reads <- clip_1[match(as.character(diff.transport.neurotrophins$uniqueIDs),as.character(clip_1$X)),-1]

colNGF <- c("#81A4D6",rgb(129/255,164/255,214/255,0.25))
colNT3 <- c("#AE72B0",rgb(174/255,114/255,176/255,0.25))

ROI       <- c("[250:200]","[200:150]","[150:100]","[100:50]")
clipdatOI <- myClips[[match(ROI[1],mynames_rages)]]
clipdatOI[is.na(clipdatOI)]<-0
for(reg in ROI[-1]){
  clipdatOI <- clipdatOI+myClips[[match(reg,mynames_rages)]]
}
clipdatOI <- clipdatOI[match(diff.transport.neurotrophins$uniqueIDs,rownames(clipdatOI)),]

diff.transport.neurotrophins$wL <- anno_tot$newL[match(diff.transport.neurotrophins$uniqueIDs,anno_tot$uniqueID)]

is_over_NGF    = factor(ifelse(myVal[[1]]==1,"over_NGF","none"))
is_under_NGF   = factor(ifelse(myVal[[2]]==1,"under_NGF","none"))
is_over_NT3    = factor(ifelse(myVal[[3]]==1,"over_NT3","none"))
is_under_NT3   = factor(ifelse(myVal[[4]]==1,"under_NT3","none"))
is_more_NGF    = factor(ifelse(myVal[[5]]==1,"over_NGF","none"))
is_more_NT3    = factor(ifelse(myVal[[6]]==1,"over_NT3","none"))
diff_NGF       = unlist(mapply(A=as.character(is_over_NGF),B=as.character(is_under_NGF),function(A,B){
  if(sum(c(A,B)!="none")>0){return(c(A,B)[which(c(A,B)!="none")])}
  else{
    return("none")
  }
}))

diff_NT3       = unlist(mapply(A=as.character(is_over_NT3),B=as.character(is_under_NT3),function(A,B){
  if(sum(c(A,B)!="none")>0){return(c(A,B)[which(c(A,B)!="none")])}
  else{
    return("none")
  }
}))

diff_neuro       = unlist(mapply(A=as.character(is_more_NGF),B=as.character(is_more_NT3),function(A,B){
  if(sum(c(A,B)!="none")>0){return(c(A,B)[which(c(A,B)!="none")])}
  else{
    return("none")
  }
}))

```

```{r compute_fisher,eval=FALSE}
myFracsROI       <- matrix(0,ncol=7,nrow=ncol(clipdatOI))
myPvalsFisherROI <- matrix(0,ncol=9,nrow=ncol(clipdatOI))

colnames(myPvalsFisherROI)<-c("fisher_NGF_over","fisher_NGF_under","fisher_NGF_diff",
                           "fisher_NT3_over","fisher_NT3_under","fisher_NT3_diff",
                           "fisher_more_NGF","fisher_more_NT3","fisher_diff")

colnames(myFracsROI)<-c("frac_bg","frac_NGF_over","frac_NGF_under",
                     "frac_NT3_over","frac_NT3_under",
                     "frac_more_NGF","frac_more_NT3")

rownames(myPvalsFisherROI)<-rownames(myFracsROI)<-colnames(clipdatOI)

for(IX in  c(1:ncol(clipdatOI))){

  is_bound              <- factor(ifelse(clipdatOI[,IX]>0,"bound","unbound"))
  
  frac_bg       <- sum(clipdatOI[,IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[,IX]))
  frac_NGF_over <- sum(clipdatOI[is_over_NGF!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_over_NGF!="none",IX]))
  frac_NGF_under <- sum(clipdatOI[is_under_NGF!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_under_NGF!="none",IX]))
  frac_NT3_over <- sum(clipdatOI[is_over_NT3!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_over_NT3!="none",IX]))
  frac_NT3_under <- sum(clipdatOI[is_under_NT3!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_under_NT3!="none",IX]))
  frac_more_NGF <- sum(clipdatOI[is_more_NGF!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_more_NGF!="none",IX]))
  frac_more_NT3 <- sum(clipdatOI[is_more_NT3!="none",IX]>0,na.rm=TRUE)/sum(!is.na(clipdatOI[is_more_NT3!="none",IX]))
  
  myFracsROI[IX,] <- c(frac_bg,
                       frac_NGF_over,frac_NGF_under,
                       frac_NT3_over,frac_NT3_under,
                       frac_more_NGF,frac_more_NT3)

  
  fisher_over_NGF           <- -log10(fisher.test(x=is_bound,y=is_over_NGF)$p.value)*sign(frac_NGF_over-frac_bg)
  fisher_under_NGF          <- -log10(fisher.test(x=is_bound,y=is_under_NGF)$p.value)*sign(frac_NGF_under-frac_bg)
  fisher_underover_NGF      <- -log10(fisher.test(x=is_bound[diff_NGF!="none"],y=factor(diff_NGF[diff_NGF!="none"]))$p.value)*sign(frac_NGF_over-frac_NGF_under)
  
  fisher_over_NT3           <- -log10(fisher.test(x=is_bound,y=is_over_NT3)$p.value)*sign(frac_NT3_over-frac_bg)
  fisher_under_NT3          <- -log10(fisher.test(x=is_bound,y=is_under_NT3)$p.value)*sign(frac_NT3_under-frac_bg)
  fisher_underover_NT3      <- -log10(fisher.test(x=is_bound[diff_NT3!="none"],y=factor(diff_NT3[diff_NT3!="none"]))$p.value)*sign(frac_NT3_over-frac_NT3_under)
  
  fisher_more_NGF           <- -log10(fisher.test(x=is_bound,y=is_more_NGF)$p.value)*sign(frac_more_NGF-frac_bg)
  fisher_more_NT3           <- -log10(fisher.test(x=is_bound,y=is_more_NT3)$p.value)*sign(frac_more_NT3-frac_bg)
  fisher_diff               <- -log10(fisher.test(x=is_bound[which(diff_neuro!="none")],y=factor(as.character(diff_neuro[which(diff_neuro!="none")])))$p.value)*sign(frac_more_NGF-frac_more_NT3)
  
  myPvalsFisherROI[IX,] <- c(fisher_over_NGF,fisher_under_NGF,fisher_underover_NGF,
                             fisher_over_NT3,fisher_under_NT3,fisher_underover_NT3,
                             fisher_more_NGF,fisher_more_NT3,fisher_diff)
  print(IX)
}


save(list=c("myPvalsFisherROI","myFracsROI"),file="./zenodo/transport/regulation/fisher.RData")
```

```{r identify_regulators_stability,eval=TRUE}
IX=match("U2AF2_1_K652",colnames(clipdatOI))
myFracsROI <- cbind(myFracsROI,myFracsROI[,match("frac_more_NGF",colnames(myFracsROI))]-myFracsROI[,match("frac_more_NT3",colnames(myFracsROI))])

mylims <- apply(myPvalsFisherROI,2,function(Z)return(boxplot(Z,plot=FALSE)$stats))

#Negative regulators of mRNA stability--those are enriched in CB restricted and depleted in over transported
sel_NGF_fisher <- rownames(myPvalsFisherROI)[which(myPvalsFisherROI[,2]>mylims[5,2])]
sel_NT3_fisher <- rownames(myPvalsFisherROI)[which(myPvalsFisherROI[,5]>mylims[5,5])]

myRBP                       <- read.csv("./zenodo/RBPs_clip_description_full_updated.csv")
myregulators_stability      <- myRBP[myRBP$Clip%in%intersect(sel_NGF_fisher,sel_NT3_fisher),]
#write.csv(myregulators_stability,file="./zenodo/transport/regulation/regulators_stabiliyt.csv")

```

Here we perform Wech's t-test to identify regulators of transport globally.

```{r welchs_ttest, eval=FALSE}
names_files     <- c("transport_NGF","transport_NT3","NT3_vs_NGF_transport")

myTtest_transport_reg <- list()
for(ix in c(1:length(mynames_rages))){
  clipdat            <- myClips[[ix]]
  clipdat            <- clipdat[match(as.character(diff.transport.neurotrophins$uniqueIDs),rownames(clipdat)),]
  
  
  myTtest_transport_reg[[ix]]           <- matrix(0,ncol=3,nrow=ncol(clipdat))
  colnames(myTtest_transport_reg[[ix]]) <- c("NGF_ttest","NT3_ttest","diff_transport_ttest")
  rownames(myTtest_transport_reg[[ix]]) <- colnames(clipdat)
  for(colix in c(1:ncol(clipdat))){
    no.sites <- clipdat[,colix]
    
    myX               <- diff.transport.neurotrophins$transport.ngf[no.sites==0]
    myY               <- diff.transport.neurotrophins$transport.ngf[no.sites>0] 
    myTtest_transport_reg[[ix]][colix,1] <- -log10(t.test(x=myX,y=myY,var.equal = FALSE)$p.value)*sign(mean(myY,na.rm=TRUE)-mean(myX,na.rm=TRUE))
    
    
    myX <- diff.transport.neurotrophins$transport.nt3[no.sites==0]
    myY <- diff.transport.neurotrophins$transport.nt3[no.sites>0] 
    myTtest_transport_reg[[ix]][colix,2] <- -log10(t.test(x=myX,y=myY,var.equal = FALSE)$p.value)*sign(mean(myY,na.rm=TRUE)-mean(myX,na.rm=TRUE))
    
    
    myX <- diff.transport.neurotrophins$diff.transport[no.sites==0]
    myY <- diff.transport.neurotrophins$diff.transport[no.sites>0] 
    myTtest_transport_reg[[ix]][colix,3] <- -log10(t.test(x=myX,y=myY,var.equal = FALSE)$p.value)*sign(mean(myY,na.rm=TRUE)-mean(myX,na.rm=TRUE))
    print(colix)
  }
  #write.csv(myTtest_transport_reg[[ix]],paste("~/Desktop/DataAnalsyisRiccio/AxonalTransport/binary_motif_enrichment/",mynames_rages[ix],"_",names_files[i],"ttest.csv",sep=""))
  print(ix)
}
#transport efficient -- regulators overall
ttest_transport_NGF              <- do.call(lapply(myTtest_transport_reg,function(TTest)return(TTest[,1])),what=cbind)
colnames(ttest_transport_NGF)    <- mynames_rages
rownames(ttest_transport_NGF)    <- rownames(myTtest_transport_reg[[1]])

ttest_transport_NT3               <- do.call(lapply(myTtest_transport_reg,function(TTest)return(TTest[,2])),what=cbind)
colnames(ttest_transport_NT3)    <- mynames_rages
rownames(ttest_transport_NT3)    <- rownames(myTtest_transport_reg[[1]])

ttest_difftransport             <- do.call(lapply(myTtest_transport_reg,function(TTest)return(TTest[,3])),what=cbind)
colnames(ttest_difftransport)    <- mynames_rages
rownames(ttest_difftransport)    <- rownames(myTtest_transport_reg[[1]])
 
#save(list=c("ttest_difftransport","ttest_transport_NT3","ttest_transport_NGF"),file="./zenodo/transport/regulation/welch.RData")
```


```{r analysis_regulatory_regions,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=8}

plotLineWithSD<-function(x=(-unlist(lapply(myranges,mean))[-1]),mymat=fisher_proxi_NGF_reg,col1="#81A4D6",col2=rgb(129/255,164/255,214/255,0.3),YLIM=c(-2,2),YLAB="Ip Fisher Enrichment",MAIN="Proximal"){
  mean_force=apply(mymat,2,function(W)return(median(W,na.rm=TRUE)))
  sd=apply(mymat,2,function(W)return(sd(W,na.rm=TRUE)))
  psd<-mean_force+sd
  nsd<-mean_force-sd
  plot(x, mean_force, ty="l", col=col1, ylab="", lty=1,lwd=3,las=1,frame=FALSE,ylim=YLIM,xlab="")
  lines(x, psd,col=col2)
  lines(x, nsd,col=col2)
  polygon(x=c(x, rev(x)), y=c(psd, rev(nsd)), col=col2,border=col2)
  lines(x, mean_force, col=col1,lwd=3)
  abline(h=0,col="black")
  mtext(side=3,line=0,text=MAIN,cex=0.7)
  mtext(side=1,line=2,text="distance from 3' end",cex=0.7)
  mtext(side=2,line=2,text=YLAB,cex=0.7)
}

par(mfrow=c(2,3))
myrangesL         <- list(c(100,-100),c(3000,2950),c(2750,2700),c(2500,2450),c(2250,2200),c(2000,1950),c(1750,1700),c(1500,1450),c(1400,1350),c(1300,1250),c(1200,1150),c(1100,1050),c(1000,950),c(900,850),c(800,750),c(700,650),c(600,550),c(500,450),c(450,400),c(400,350),c(350,300),c(300,250),c(250,200),c(200,150),c(150,100),c(100,50),c(50,0),c(0,-50),c(-50,-100),c(-100,-150))
mynames_ragesL   <- unlist(lapply(myrangesL,function(Z)return(paste("[",Z[[1]],":",Z[[2]],"]",sep=""))))

plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=ttest_transport_NGF[,-1],col1=colNGF[1],col2=colNGF[2],YLIM=c(-5,50),YLAB="P-value [-log10 -- Welch]",MAIN="Global regulators transport NGF")
abline(v=0,col="darkred")

plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=SelectPositive(ttest_transport_NGF[,-1]),col1=colNGF[1],col2=colNGF[2],YLIM=c(0,60),YLAB="P-value [-log10 -- Welch]",MAIN="Positive regulators transport NGF")
abline(v=0,col="darkred")
plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=SelectNegative(ttest_transport_NGF[,-1]),col1=colNGF[1],col2=colNGF[2],YLIM=c(0,60),YLAB="P-value [-log10 -- Welch]",MAIN="Negative regulators transport NGF")
abline(v=0,col="darkred")

plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=ttest_transport_NT3[,-1],col1=colNT3[1],col2=colNT3[2],YLIM=c(-5,50),YLAB="P-value [-log10 -- Welch]",MAIN="Global regulators transport NT3")
abline(v=0,col="darkred")

plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=SelectPositive(ttest_transport_NT3[,-1]),col1=colNT3[1],col2=colNT3[2],YLIM=c(-5,50),YLAB="P-value [-log10 -- Welch]",MAIN="Positive regulators transport NT3")
abline(v=0,col="darkred")

plotLineWithSD(x=(-unlist(lapply(myrangesL,mean))[-1]),mymat=SelectNegative(ttest_transport_NT3[,-1]),col1=colNT3[1],col2=colNT3[2],YLIM=c(-5,50),YLAB="P-value [-log10 -- Welch]",MAIN="Negative regulators transport NT3")
abline(v=0,col="darkred")


```


The preferential localisation of the negative regulators of axonal localisation is shown below:
```{r analysis_negative_regulators,eval=FALSE,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=5, fig.height=6}
myregulators_stability      <- read.csv("./zenodo/transport/regulation/regulators_stabiliyt.csv")
par(mfrow=c(3,3))
PlotHeatmapRegulators(mots=myregulators_stability$Clip,mymats=list(ttest_transport_NGF[,-1],ttest_transport_NT3[,-1]),scaling=FALSE)
```



The preferential localisation of the positive regulators of axonal localisation is shown below:
```{r analysis_positive_regulators,eval=FALSE,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=5, fig.height=6}
myregulators_transport <- read.csv("./zenodo/transport/regulation/regulators_transport.csv")
par(mfrow=c(3,3))
PlotHeatmapRegulators(mots=myregulators_transport$Clip,mymats=list(ttest_transport_NGF[,-1],ttest_transport_NT3[,-1]),scaling=FALSE)
```




```{r condensate,eval=FALSE,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=5, fig.height=3}


mytempregs <- list(positive_reg <- rownames(myPvalsFisherROI)[which(myPvalsFisherROI[,1]>mylims[4,1])],
                   negative_reg <- rownames(myPvalsFisherROI)[which(myPvalsFisherROI[,2]>mylims[4,2])])


require(scales)
BxPlotComps <- function(myList=lapply(mytempregs,function(Z)return(myRBP$transport.ngf[which(myRBP$Clips%in%Z)])),
                        myname="transport NGF",
                        myCol=c(colsNGF[1],colsNGF[2]),
                        YLIM=c(-2,2)){
  boxplot(myList,xlab="",ylab="",xaxt="n",col=myCol,frame=FALSE,outline=FALSE,las=1,ylim=YLIM)
  mtext(side=1,line=0,at=c(1,2),text=c("(+)reg","(-)reg"),cex=0.7)
  mtext(side=3,line=0,text=paste("P=",scientific(t.test(myList[[1]],myList[[2]],var.equal = FALSE)$p.value,digit=2),sep=""),cex=0.7)
  mtext(side=2,line=2,text=myname,cex=0.7)
}


par(mar=c(2,3,2,1))
mytempregs <- list(positive_reg <- unique(myregulators_transport$RBPs),
                   negative_reg <- unique(myregulators_stability$RBPs))
layout(matrix(c(1,2,3,1,4,5),ncol=3,nrow=2,byrow = TRUE))
par(mar=c(2,3,2,1),mfrow=c(1,1))
BxPlotComps(myList=lapply(mytempregs,function(Z)return(myRBP$Condensate_Zscore[match(Z,myRBP$RBPs)])),
            myname="condensate",
            myCol=c("white","grey"),
            YLIM=c(1.5,3.5))

```


```{r TDP43_iClip_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="TDP43_iClip",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```


```{r IGF2BP1_2_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="IGF2BP1_2_K652",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```


```{r FMR1_1_K652_K652_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="FMR1_1_K652",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```


```{r KHSRP_2_K652_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="KHSRP_2_K652",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```

```{r HuR_iClip_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="HuR_iClip",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```

```{r FUBP3_1_HepG2_ex,echo=TRUE,warning=FALSE,eval=TRUE,fig.width=12, fig.height=6}
PlotTransportAnalysisGlobal(mot="FUBP3_1_HepG2",by_te=0.1,by_cb=0.1,limfrac=0.6,myMAX=0.6)
```