---
title: "Overview of the Data-set"
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(lme4)
#library(nlme)
#source("http://bioconductor.org/biocLite.R")
#require("GO.db")
#require("limma")
#require("topGO")
#require("biomaRt")
#require("org.Rn.eg.db")
#library(grDevices)
library(Rsamtools)
library(IRanges)
library(GenomicRanges)
library(rtracklayer)
library(geneplotter)
library("ape")
require("multtest")
require("mclust")

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")
```

```{r selection_expressed_isoforms,eval=FALSE}
# An isoforms is considered to be reliably expressed if in all the samples the probability of it belonging to the non-expressed class is below 0.05 (less than 5% chance to below to the background in both replicates) -- soft threshold
#OR
# An isoforms is considered to be expressed if in at least one of the replicate the probability of it belonging to the expressed class is above 0.1 (more than 10% chance to belong to the foreground in at least one replicate) -- hard threshold

SelectExpressed <- function(dat=log2(htseq[,1]+1),frac.bg=0.6,frac.fg=0.1){
  bimdens       <- densityMclust(data=dat,G=2)
  x   <- seq(from=0, to=max(dat),length=100)
  if(is.na(bimdens$parameters$variance$sigmasq[2])){
    Lim.fg           <- qnorm(frac.fg,mean=bimdens$parameters$mean[2],sd=sqrt(bimdens$parameters$variance$sigmasq[1]))
    Lim.bg           <- qnorm(frac.bg,mean=bimdens$parameters$mean[1],sd=sqrt(bimdens$parameters$variance$sigmasq[1]))
    hx2              <- dnorm(x,mean=bimdens$parameters$mean[2],sd=sqrt(bimdens$parameters$variance$sigmasq[1]))
  }
  if(!is.na(bimdens$parameters$variance$sigmasq[2])){
    Lim.fg           <- qnorm(frac.fg,mean=bimdens$parameters$mean[2],sd=sqrt(bimdens$parameters$variance$sigmasq[2]))
    Lim.bg           <- qnorm(frac.bg,mean=bimdens$parameters$mean[1],sd=sqrt(bimdens$parameters$variance$sigmasq[1]))
    hx2              <- dnorm(x,mean=bimdens$parameters$mean[2],sd=sqrt(bimdens$parameters$variance$sigmasq[2]))
  }
  
  hx1 <- dnorm(x,mean=bimdens$parameters$mean[1],sd=sqrt(bimdens$parameters$variance$sigmasq[1]))
  hist(dat,breaks=50,col=rgb(0,0,0,alpha=0.2),freq=FALSE,xlab="",ylab="",main="",ylim=c(0,0.4),las=1,xlim=c(0,20))
  lines(x,hx2 , lwd=2, col="blue")
  lines(x,hx1 , lwd=2, col="red")
  abline(v=Lim.fg,col="blue",lty=2)
  abline(v=Lim.bg,col="red",lty=2)
  return(c(Lim.bg,Lim.fg))
}

#Selection of reliably expressed 3' UTR
#par(mfrow=c(2,4))
#mydat<-data.sum
#for(i in c(1:ncol(mydat))){
#  mydat[,i] <- log2(data.sum[,i]+1)
#}
#lims      <- apply(mydat,2,function(Z)2^SelectExpressed(dat=Z,frac.bg=0.999,frac.fg=0.1))

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

```



# Analysis of the distribution of the raw count {.tabset}
The distributions of the expression of 12'771 3' UTR isoforms detected in cell body and/or distal axons of sympathetic neurons. As shown below, the detected expression levels in axonal compartment is systematically lower than in cell body compartment. As discussed in materials and method section, the amount of mRNA material obtained from axonal compartment is much lower. Furthermore the large amount of mithondria located in axonal compartments arise in a library maily composed of mithocondrial mRNAs. This leads to technical challenges in library preparation which were overcomed by performing two rounds of mRNA linear amplification. We alos find that one of the axonal samples exposed to NT3 was of lower quality (exhibiting systematically lower expression level), which we take into account in our analysis.


```{r analysis_distributions,fig.width=10, fig.height=5,warning=FALSE}
data     <- anno_tot[,grep(".raw.corrected",colnames(anno_tot))]#Raw count corrected for reads arising from anti-sense transcripts which have leaked 
row_info <- anno_tot

#Here focus on the isoforms annotated in Rn5
data.sub              <- data
data.sub              <- data[row_info$is.conservative==TRUE,]
data.log2             <- log2(data.sub)
data.log2[data.sub<=1]<-0
colnames(data.log2)   <- gsub(colnames(data.log2),pattern="\\.raw\\.corrected",repl="")
temp                  <- lapply(colnames(data.log2),function(Z)return(strsplit(Z,split="\\.")[[1]]))         
col_info              <- data.frame(do.call(args=temp,what=rbind))
colnames(col_info)    <- c("treatment","compartment","replicate")
col_info$group        <- factor(paste(col_info$treatment,col_info$compartment,sep="_"),levels=c("NGF_cb","NGF_axon","NT3_cb","NT3_axon"))
mycols_treatment      <- c("#81A4D6","#B072B0")
names(mycols_treatment)<- c("NGF","NT3")



#Look at the distribution of the data
par(mfrow=c(1,2))
boxplot(data.log2,outline=FALSE,las=2)
multidensity(data.log2)

```


As shown below, although one the NT3 axonal sample does exhibit lower quality, the correlation between technical replicates is similar as for the NGF samples.


```{r correlation_between_samples,warning=FALSE,fig.width=10, fig.height=3,warning=FALSE}
par(mfrow=c(1,4))
#CB NGF
plot(data.log2[,"NGF.cb.1"],data.log2[,"NGF.cb.2"],cex=0.2,col=rgb(0,0,0,0.2),pch=19,main="NGF Cell Body",las=1,frame=FALSE,xlab="Cell body 1 [log2(read counts)]",ylab="Cell body 2 [log2(read counts)]")
abline(0,1,col="red")
mtext(side=3,line=0,text=paste("PCC=",round(cor(data.log2[,"NGF.cb.1"],data.log2[,"NGF.cb.2"]),digits=2)),cex=1.0)
#CB NT3
plot(data.log2[,"NT3.cb.1"],data.log2[,"NT3.cb.2"],cex=0.2,col=rgb(0,0,0,0.2),pch=19,main="NT3 Cell Body",las=1,frame=FALSE,xlab="Cell body 1 [log2(read counts)]",ylab="Cell body 2 [log2(read counts)]")
abline(0,1,col="red")
mtext(side=3,line=0,text=paste("PCC=",round(cor(data.log2[,"NT3.cb.1"],data.log2[,"NT3.cb.2"]),digits=2)),cex=1.0)
#Axons NGF
plot(data.log2[,"NGF.axon.1"],data.log2[,"NGF.axon.2"],cex=0.2,col=rgb(0,0,0,0.2),pch=19,main="NGF Axons",las=1,frame=FALSE,xlab="Axons 1 [log2(read counts)]",ylab="Axons 2 [log2(read counts)]")
abline(0,1,col="red")
mtext(side=3,line=0,text=paste("PCC=",round(cor(data.log2[,"NGF.axon.1"],data.log2[,"NGF.axon.2"]),digits=3)),cex=1.0)
#Axons NT3
plot(data.log2[,"NT3.axon.1"],data.log2[,"NT3.axon.2"],cex=0.2,col=rgb(0,0,0,0.2),pch=19,main="NT3 Axons",las=1,frame=FALSE,xlab="Axon 1 [log2(read counts)]",ylab="Axon 2 [log2(read counts)]")
abline(0,1,col="red")
mtext(side=3,line=0,text=paste("PCC=",round(cor(data.log2[,"NT3.axon.1"],data.log2[,"NT3.axon.2"]),digits=3)),cex=1.0)
```

While the above plots are derived from expression quantification based on the Rn5 annotation, below is shown a similar analysis using the total pool of 35,534 3' UTR isoforms which contain 12,771 Rn5 3' UTRs and 22,763 novel 3' UTRs. Using the entire pool of transcripts, the correlation between samples is lower which may be due to increased diversity in transcriptome.

```{r correlation_between_samples_all,warning=FALSE,fig.width=10, fig.height=3,warning=FALSE}

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[,match(coloi,colnames(anno_tot))]
fact                <- 10^8/apply(mydat,2,sum)
for(i in c(1:ncol(mydat))){
  mydat[,i]             <- log2(fact[i]*mydat[,i])
  mydat[mydat[,i]<=0,i] <- 0
}


par(mfrow=c(1,4),mar=c(4,4,2,2))
for(IX in c(3,7,1,5)){
  plot(mydat[,IX],mydat[,IX+1],pch=19,col=rgb(0,0,0,0.2),cex=0.2,frame=F,las=1,xlab=gsub(colnames(mydat)[IX],pattern="\\.raw\\.corrected",repl=""),ylab=gsub(colnames(mydat)[IX+1],pattern="\\.raw\\.corrected",repl=""))
  mtext(side=3,line=0,cex=0.7,text=paste("SRC=",round(cor(mydat[,IX],mydat[,IX+1],method="spearman"),digit=2)))
}
```


# Unsupervised clustering analysis of the samples {.tabset .tabset-fade .tabset-pills}

## All data 
As shown below using hierarchical clustering and principal component analyses of the 8 samples, the first source of gene expression variation is the compartment while the treatment effect is strongly detected in cell body compared to axonal samples. From this analysis, it is also visible that the NT3 axonal samples are clustering together.  
```{r unsupervised_analysis_all,warning=FALSE,fig.width=10, fig.height=5,warning=FALSE}
layout(matrix(c(1,2,3,1,4,5),ncol=3,nrow=2,byrow =TRUE))
#HC
dist1               <- as.matrix(cor(data.log2,metho="spearman",use="na.or.complete"))
mydist              <- as.dist(1-dist1)
hcl                 <- hclust(mydist, method = "complete", members = NULL)
plot(as.phylo(hcl),tip.color=mycols_treatment[match(col_info$treatment,names(mycols_treatment))] ,cex=1.0,label.offset = 0.001,no.margin = FALSE,use.edge.length=TRUE,direction="rightwards",plot=TRUE,font=1)

#SVD on the entire set
#1. SVD init
SVD_eset<- svd(data.log2)
pi      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
d       <- -1*sum(pi*log10(pi))/log10(length(pi))
#2. Remvove first component
datm    <- data.log2-SVD_eset$d[1]*(SVD_eset$u[,1]%*%t(SVD_eset$v[,1]))
SVD_eset<- svd(datm)
pip      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
dp       <- -1*sum(pip*log10(pip))/log10(length(pip))


mypoints       <- c(4,19)
names(mypoints)<- c("cb","axon")

mypchs=rep(mypoints,2)
mycols=c(rep(mycols_treatment[1],2),rep(mycols_treatment[2],2))
mylegs=paste(names(mypchs),names(mycols),sep=":")

par(mar=c(4,4,1,1))
barplot(pip[c(1:which(cumsum(pip)>=0.99)[1])],las=1,ylim=c(0,1),cex.main=0.7,cex.axis=1.0,col="black")
mtext(side = 1, line = 2, text ="principal components", col = "black",cex=0.7, font=1)
mtext(side = 2, line = 2, text ="Fraction of explained variance", col = "black",cex=0.7, font=1)
mtext(side = 3, line = -2, text = paste("Shannon Entropy = ",round(dp,digits=3)), col = "black",cex=0.7, font=1)
i=1
for(j in c(2:4)){
  plot(SVD_eset$v[,i],SVD_eset$v[,j],pch=mypoints[match(col_info$compartment,names(mypoints))],col=mycols_treatment[match(col_info$treatment,names(mycols_treatment))],xlab=paste("PC",i,sep=""),ylab=paste("PC",j,sep=""),cex=1.5)
  abline(h=0,col="grey",lty=2)
  abline(v=0,col="grey",lty=2)
  if(j==2){
    legend("bottomright",ncol=1,pch=mypchs,col=mycols,leg=mylegs,cex=0.8)
  }
}
```


## Cell body samples

Further looking at the samples from cell body compartment only reveals that while treatment effect is captured by PC1, PC2 and PC3 might reveal some technical bias which should not impact the downstream analysis as distributed across different conditions.


```{r unsupervised_analysis_cb,warning=FALSE,fig.width=10, fig.height=2.4,warning=FALSE}
par(mfrow=c(1,3))
ix_oi <- col_info$compartment=="cb"
dat_cb <- data.log2[,ix_oi]
#SVD on the entire set
#1. SVD init
SVD_eset<- svd(dat_cb)
pi      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
d       <- -1*sum(pi*log10(pi))/log10(length(pi))
#2. Remvove first component
datm    <- dat_cb-SVD_eset$d[1]*(SVD_eset$u[,1]%*%t(SVD_eset$v[,1]))
SVD_eset<- svd(datm)
pip      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
dp       <- -1*sum(pip*log10(pip))/log10(length(pip))


mypoints       <- c(4,19)
names(mypoints)<- c("cb","axon")

mypchs=rep(mypoints,2)
mycols=c(rep(mycols_treatment[1],2),rep(mycols_treatment[2],2))
mylegs=paste(names(mypchs),names(mycols),sep=":")

par(mar=c(4,4,1,1))
barplot(pip[c(1:which(cumsum(pip)>=0.99)[1])],las=1,ylim=c(0,1),cex.main=0.7,cex.axis=1.0,col="black")
mtext(side = 1, line = 2, text ="principal components", col = "black",cex=0.7, font=1)
mtext(side = 2, line = 2, text ="Fraction of explained variance", col = "black",cex=0.7, font=1)
mtext(side = 3, line = -2, text = paste("SE = ",round(dp,digits=3)), col = "black",cex=0.7, font=1)
i=1
for(j in c(2:3)){
  plot(SVD_eset$v[,i],SVD_eset$v[,j],col="white",xlab=paste("PC",i,sep=""),ylab=paste("PC",j,sep=""),xlim=c(min(SVD_eset$v[,i])-0.1,max(SVD_eset$v[,i])+0.1),ylim=c(min(SVD_eset$v[,j])-0.1,max(SVD_eset$v[,j])+0.1))
  text(x=SVD_eset$v[,i],y=SVD_eset$v[,j],labels=colnames(dat_cb),col=mycols_treatment[match(col_info$treatment[ix_oi],names(mycols_treatment))],cex=0.7)
  abline(h=0,col="grey",lty=2)
  abline(v=0,col="grey",lty=2)
}
```

## Axonal samples 

Next looking at the samples from axonal compartment only shows that  treatment effect is detected in PC2 and PC3 while PC1 identifies the NT3 sample with lower amount of reads.

```{r unsupervised_analysis_ax,warning=FALSE,fig.width=10, fig.height=2.5,warning=FALSE}
par(mfrow=c(1,3))
ix_oi <- col_info$compartment=="axon"
dat_cb <- data.log2[,ix_oi]
#SVD on the entire set
#1. SVD init
SVD_eset<- svd(dat_cb)
pi      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
d       <- -1*sum(pi*log10(pi))/log10(length(pi))
#2. Remvove first component
datm    <- dat_cb-SVD_eset$d[1]*(SVD_eset$u[,1]%*%t(SVD_eset$v[,1]))
SVD_eset<- svd(datm)
pip      <- SVD_eset$d*SVD_eset$d/sum(SVD_eset$d*SVD_eset$d)
dp       <- -1*sum(pip*log10(pip))/log10(length(pip))


mypoints       <- c(4,19)
names(mypoints)<- c("cb","axon")

mypchs=rep(mypoints,2)
mycols=c(rep(mycols_treatment[1],2),rep(mycols_treatment[2],2))
mylegs=paste(names(mypchs),names(mycols),sep=":")

par(mar=c(4,4,1,1))
barplot(pip[c(1:which(cumsum(pip)>=0.99)[1])],las=1,ylim=c(0,1),cex.main=0.7,cex.axis=1.0,col="black")
mtext(side = 1, line = 2, text ="principal components", col = "black",cex=0.7, font=1)
mtext(side = 2, line = 2, text ="Fraction of explained variance", col = "black",cex=0.7, font=1)
mtext(side = 3, line = -2, text = paste("Shannon Entropy = ",round(dp,digits=3)), col = "black",cex=0.7, font=1)
i=1
for(j in c(2:3)){
  plot(SVD_eset$v[,i],SVD_eset$v[,j],col="white",xlab=paste("PC",i,sep=""),ylab=paste("PC",j,sep=""),xlim=c(min(SVD_eset$v[,i])-0.1,max(SVD_eset$v[,i])+0.1),ylim=c(min(SVD_eset$v[,j])-0.1,max(SVD_eset$v[,j])+0.1))
  text(x=SVD_eset$v[,i],y=SVD_eset$v[,j],labels=colnames(dat_cb),col=mycols_treatment[match(col_info$treatment[ix_oi],names(mycols_treatment))],cex=0.7)
  abline(h=0,col="grey",lty=2)
  abline(v=0,col="grey",lty=2)
}
```

In summary this preliminary unsupervised analysis shows that both treatment and compartment effects are detected at the expression level however that axonal samples exhibit weaker difference between NGF and NT3 condition.

