paperFigures.R

library("rootSolve")
library("data.table")
## use library("SingleSampleRNAdynamics") if the package is installed
source("SingleSampleRNAdynamics/R/SingleSampleRNAdynamics.R")
library("ggplot2")
library("EnsDb.Mmusculus.v79")
library("org.Mm.eg.db")
library("gridExtra")


## The figures of the paper were generated using the following sessionInfo()
##
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS

## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so

## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_CH.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=de_CH.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=de_CH.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C       

## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     

## other attached packages:
##  [1] gridExtra_2.3              org.Mm.eg.db_3.12.0       
##  [3] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.14.0          
##  [5] AnnotationFilter_1.14.0    GenomicFeatures_1.42.3    
##  [7] AnnotationDbi_1.52.0       Biobase_2.50.0            
##  [9] GenomicRanges_1.42.0       GenomeInfoDb_1.26.6       
## [11] IRanges_2.24.1             S4Vectors_0.28.1          
## [13] BiocGenerics_0.36.0        data.table_1.14.0         
## [15] rootSolve_1.8.2.1         

## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6                  lattice_0.20-41            
##  [3] prettyunits_1.1.1           Rsamtools_2.6.0            
##  [5] Biostrings_2.58.0           assertthat_0.2.1           
##  [7] utf8_1.2.1                  BiocFileCache_1.14.0       
##  [9] R6_2.5.0                    RSQLite_2.2.5              
## [11] httr_1.4.2                  pillar_1.5.1               
## [13] zlibbioc_1.36.0             rlang_0.4.10               
## [15] progress_1.2.2              lazyeval_0.2.2             
## [17] curl_4.3                    rstudioapi_0.13            
## [19] blob_1.2.1                  Matrix_1.3-2               
## [21] BiocParallel_1.24.1         stringr_1.4.0              
## [23] ProtGenerics_1.22.0         RCurl_1.98-1.3             
## [25] bit_4.0.4                   biomaRt_2.46.3             
## [27] DelayedArray_0.16.3         compiler_4.0.5             
## [29] rtracklayer_1.50.0          pkgconfig_2.0.3            
## [31] askpass_1.1                 openssl_1.4.3              
## [33] tidyselect_1.1.0            SummarizedExperiment_1.20.0
## [35] tibble_3.1.0                GenomeInfoDbData_1.2.4     
## [37] matrixStats_0.58.0          XML_3.99-0.6               
## [39] fansi_0.4.2                 crayon_1.4.1               
## [41] dplyr_1.0.5                 dbplyr_2.1.1               
## [43] GenomicAlignments_1.26.0    bitops_1.0-6               
## [45] rappdirs_0.3.3              grid_4.0.5                 
## [47] gtable_0.3.0                lifecycle_1.0.0            
## [49] DBI_1.1.1                   magrittr_2.0.1             
## [51] stringi_1.5.3               cachem_1.0.4               
## [53] XVector_0.30.0              xml2_1.3.2                 
## [55] ellipsis_0.3.1              generics_0.1.0             
## [57] vctrs_0.3.7                 tools_4.0.5                
## [59] bit64_4.0.5                 glue_1.4.2                 
## [61] purrr_0.3.4                 hms_1.0.0                  
## [63] MatrixGenerics_1.2.1        fastmap_1.1.0              
## [65] memoise_2.0.0              

with.inspect.comp <- TRUE # set to FALSE to avoid the lengthy rate generation from INSPEcT

figdir  <- "mfigs"
figpath <- function(fname){
    return(paste(figdir,fname,sep="/"))
}

######### SIMULATED DATA (Fig.2 ) ###################3



nsim <- 50000  # number of data points

tp <- t(exp(seq(-5,5))) # time points for trajectories

ntp  <- length(tp)

## generating random rates
ra <- exp(rnorm(nsim,0,2)) # production rate
rb <- exp(rnorm(nsim,0,2)) # degradation rate
rc <- exp(rnorm(nsim,0,2)) # processing rate

pre.p <- ra/rc*exp(-rc %*% tp) # pre-existing (or unlabeled) premature RNA
pre.m <- ra/(rb-rc)*exp(-rc %*% tp) - ra*rc/(rb*(rb-rc))*exp(-rb %*% tp) # pre-existing mature RNA

lab.p <- ra/rc*(1-exp(-rc %*% tp)) # labeled premature RNA
lab.m <- ra/rb*(1+rc/(rb-rc)*exp(-rb %*% tp))-ra/(rb-rc)*exp(-rc %*% tp)# labeled mature RNA

pre.frac <- pre.p/(pre.p+pre.m) ## unlabeled observables 
lab.frac <- lab.p/(lab.p+lab.m) ## labeled observables


ersd = 0.00 ## simulated noise level (not used because we also have real data)

pre.frac <- exp(log(pre.frac)+rnorm(nsim,sd=ersd))
lab.frac <- exp(log(lab.frac)+rnorm(nsim,sd=ersd))

pre.frac[pre.frac==0] <-NA
lab.frac[lab.frac==0] <- NA
sc <- 5 # using a single time point to infer rates.
t <- tp[sc]
frac <-cbind(pre.frac[,sc],lab.frac[,sc])

## estimating rates from observables
my.ratesn  <-  apply(frac,1,solve.rates,t=t)
rates <- cbind(sapply(my.ratesn,"[[",1),sapply(my.ratesn,"[[",2))

## classifying number of solutions
to.show <- sapply(my.ratesn,function(x)length(x)==2)
get.table <- function(lrb,lrc,rates,ambig,frac){
    diffb = sqrt((log(rb)-rates[,1])^2)
    diffc = sqrt((log(rc)-rates[,2])^2)
    res.data <- data.table(rb=rb,rc=rc,lrb=log(rb), lrc=log(rc),diffb=diffb,diffc=diffc,preex.ratio=frac[,1],label.ratio=frac[,2])
    res.data[,col:=ifelse(!ambig,"A",ifelse(diffc>0.1,"B","C"))]
    res.data[,k:= lrc-lrb]
    return(res.data)
}
resn  <- get.table(log(rb),log(rc),rates,to.show,frac)

p1 <- ggplot(data=resn,aes(x=preex.ratio,y=label.ratio,color=col))  
p1 <- p1 +geom_point(aes(stroke=0))
p1 <- p1 +geom_line(data=data.frame(x=seq(0,1,0.05),col="green"),aes(x=x,y=1/(2-x),colour=as.factor(col)),size=1.5,show.legend=FALSE,guide=FALSE) + theme_classic() + scale_color_manual(name ="rate estimate", labels=c("ambiguous","wrong","correct","b=1/(2-a)"), values=c(rgb(0,0,1,0.5),rgb(1,0,0,0.5),rgb(0.5,0.5,0.5,0.5),"green"))+labs(x="unlabeled ratio (a)",y = "labeled ratio (b)")
p1 <- p1 + scale_x_continuous(limits=c(0,1)) + scale_y_continuous(limits=c(0,1)) 
p1  <- p1 + theme(text = element_text(size = 20), legend.position=c(0.7,0.3))
pp <- p1 + geom_abline(slope=1,color="red")+geom_abline(intercept=1,slope=0,color="blue")
pp

### plotting trajectories Fig 2  bottom center
trajs <- data.table(lra=log(ra),lrb=log(rb),lrc=log(rc),k=log(rc)-log(rb),preex.ratio=pre.frac,label.ratio=lab.frac)
p1  <- ggplot(trajs,aes(x=preex.ratio.V2,y=label.ratio.V2,z=k))+theme_classic()
for (i in seq(-6,6)){
    for (j in seq_along(tp)){
        cx  <- paste0("preex.ratio.V",j)
        cy  <- paste0("label.ratio.V",j)
        topl <- trajs[abs(2*k-i)<0.01]
        topl$col <- i
        p1  <- p1+ geom_point(topl,mapping=aes_string(x=cx,y=cy,col="col"),shape=20,size=1)
    }
}
p1  <- p1 + scale_color_gradient2(low="darkolivegreen1",high="darkorange1",mid="black",midpoint=0,name= "k [log]")
p1 <- p1 +labs(x="unlabeled ratio (a)",y = "labeled ratio (b)")
p1  <-  p1 +geom_line(data=data.frame(x=seq(0,1,0.05),col=1.0,k=100),aes(x=x,y=1/(2-x)),color="green",show.legend=FALSE)
p1  <- p1 +  geom_abline(slope=1,color="red")+geom_abline(intercept=1,slope=0,color="blue")
kk  <- seq(-4,4,4)
pt  <- p1 + annotate("text",x = 1/(1+exp(kk/2)),y =1.05,label=paste("log(k)=",kk/2))
pt <- pt + theme(text = element_text(size = 20), legend.position=c(0.8,0.3))
pt  




#### noise analysis ##################

#all.sdrel  <-c(0,1,5,10,50)
all.sd  <- c(0,0.01,0.05,0.1,0.5)
sim.corr1S <- matrix(NA,nrow=length(all.sd),ncol=3)
init.corr <- sim.corr1S
nsol <- rep(NA,length(all.sd))
for (i in rev(seq(length(all.sd)))){
ersd = log(1+all.sd[i]) ## simulated noise level

### simulated introns and exons
pre.int = pre.p 
pre.ex = pre.p + pre.m
lab.int = lab.p
lab.ex = lab.p + lab.m

## adding gaussian noise in log space
npre.int <- exp(log(pre.int)+rnorm(nsim,sd=ersd))
npre.ex <- exp(log(pre.ex)+rnorm(nsim,sd=ersd))
nlab.int <- exp(log(lab.int)+rnorm(nsim,sd=ersd))
nlab.ex <- exp(log(lab.ex)+rnorm(nsim,sd=ersd))



sc <- 5 # using a single time point to infer rates.
t <- tp[sc]




## building required dataframe
sim.data  <-  data.frame(sim.prod=log(ra), sim.deg=log(rb),sim.proc=log(rc),unlabeled.intron=npre.int[,sc], unlabeled.exon=npre.ex[,sc],labeled.intron=nlab.int[,sc],labeled.exon=nlab.ex[,sc])

## estimating rates from observables
res.data <-singleSampleRateEstimation(sim.data)

## computing correlations
keep <- which(! res.data$twosols)
init.corr[i,1] <- cor(c(npre.ex[keep,sc],npre.int[keep,sc],nlab.ex[keep,sc],nlab.int[keep,sc]),c(pre.ex[keep,sc],pre.int[keep,sc],lab.ex[keep,sc],lab.int[keep,sc]))

sim.corr1S[i,1] <- res.data[twosols==FALSE,cor(x=sim.prod,y=prod.rate,use="pairwise.complete")]
sim.corr1S[i,2] <- res.data[twosols==FALSE,cor(x=sim.deg,y=deg.rate,use="pairwise.complete")]
sim.corr1S[i,3] <- res.data[twosols==FALSE,cor(x=sim.proc,y=proc.rate,use="pairwise.complete",)]

nsol[i] <- res.data[!is.na(deg.rate),.N]/sum(!is.na(rates[,1]))
}

corr.sim <- data.frame(sd.val = all.sd,init.cor=init.corr[,1],frac.sol=nsol,prod.cor=sim.corr1S[,1],deg.cor=sim.corr1S[,2],proc.cor=sim.corr1S[,3])

## Fig 2 - bottom right 
pn <- ggplot(data=corr.sim[1:5,], aes(x=sd.val)) + geom_line(aes(y=deg.cor,color="degradation rate"),size=1) + geom_line(aes(y=proc.cor,color="processing rate"),size=1)+ geom_line(aes(y=prod.cor,color="synthesis rates"),size=1)+labs(x="noise relative standard deviation",y="correlation with true value",color="") +theme_classic()+theme(legend.position=c(0.5,0.8),text=element_text(size=20)) + scale_x_continuous(labels = scales::percent_format()) + scale_y_continuous(labels = scales::percent_format())
mfrac.sol <- corr.sim$frac.sol[1]
pn  <- pn + geom_line(aes(y=frac.sol/mfrac.sol,color=" fraction of estimated rates"),size=1,linetype="dashed")
pn <- pn + scale_color_manual(values = c("black", "red","blue", "green"))


## fig.2  - top row 

ps <- ggplot(data=res.data[solvable=="C"],aes(x=sim.prod,y=prod.rate-log(t))) + geom_point(alpha=0.1)+theme_classic() +geom_abline(color="red")+ theme(text=element_text(size=20)) + scale_y_continuous(limits=c(-7.5,7.5))+ scale_x_continuous(limits=c(-7.5,7.5)) + labs(x="true value [log]",y = "estimated value [log]", title="synthesis rate")

pd <- ggplot(data=res.data[solvable=="C"],aes(x=sim.deg,y=deg.rate-log(t))) + geom_point(alpha=0.1)+theme_classic()+geom_abline(color="red")+ theme(text=element_text(size=20)) + scale_y_continuous(limits=c(-7.5,7.5))+ scale_x_continuous(limits=c(-7.5,7.5)) + labs(x="true value [log]",y = "estimated value [log]", title="degradation rate")

pb <- ggplot(data=res.data[solvable=="C"],aes(x=sim.proc,y=proc.rate-log(t))) + geom_point(alpha=0.1)+theme_classic()+geom_abline(color="red")+ theme(text=element_text(size=20)) + scale_y_continuous(limits=c(-7.5,7.5))+ scale_x_continuous(limits=c(-7.5,7.5)) +labs(x="true value [log]",y = "estimated value [log]", title="processing rate")

grid.arrange(ps,pd,pb,pp,pt,pn,nrow=2,heights=c(1,2))

dev.copy2pdf(file = figpath("simul_dat.pdf"),onefile=T)



##### REAL DATA ####################3

##loading  and preparing data ####

data.WT10 <- fread("./data/transcripts_tpm.csv")



## getting transcript annotations
edb  <- EnsDb.Mmusculus.v79
biot <- select(edb, keys=data.WT10$rn, columns=c("TXID", "TXBIOTYPE"),keytype="TXID")
idx <- match(biot$TXID,data.WT10$rn)
data.WT10$biotype <- NA
data.WT10$biotype[idx] <- biot$TXBIOTYPE

rownames(data.WT10) <- data.WT10$rn
data.WT10[,rn:=NULL]

## getting gene symbols


gene.info <- select(edb, keys=rownames(data.WT10), columns=c("TXID", "TXBIOTYPE", "SYMBOL"),keytype="TXID")
data.WT10$txid = rownames(data.WT10)
data.WT10[,symbol:=gene.info$SYMBOL[match(txid,gene.info$TXID)]]




###### Figure 3  ##############

## Fig 3 left - plotting the raw data
p1 <- ggplot(data.WT10,aes(x=log(P_1_WT_10_8d.intron/P_1_WT_10_8d.exon),y = log(L_1_WT_10_8d.intron/L_1_WT_10_8d.exon)))
##p1 <- p1+ geom_point(alpha=0.01,pch=20)+scale_x_continuous(limits = c(-7, 3)) + scale_y_continuous(limits = c(-7, 3))
p1 <- p1 + geom_rect(xmin=-10,xmax=-7,ymin=0,ymax=10,fill="#CCFFFF",alpha=0.1)
p1 <- p1 + geom_rect(xmin=-10,xmax=-7,ymin=-10,ymax=0,fill="#FFCCCC",alpha=0.1)
p1 <- p1+ geom_point(aes(alpha=log(1+P_1_WT_10_8d.exon + P_1_WT_10_8d.intron)),pch=20)+scale_x_continuous(limits = c(-7, 3)) + scale_y_continuous(limits = c(-7, 3)) + scale_alpha(range=c(0,0.1),guide=F)

p1 <- p1 + geom_density_2d(color="wheat")

p1 <- p1 + geom_line(data=data.frame(x=seq(0,1,0.05)),aes(x=log(x),y=log(1/(2-x))),colour="green") + geom_abline(slope=1,colour=rgb(1,0,0,0.5)) + geom_abline(slope=0,colour=rgb(0,0,1,0.5))+theme_classic()+theme(text = element_text(size = 20))+ ylab("labeled ratio (b) [log]")+ xlab("unlabeled ratio (a) [log]") + xlab(bquote("observed unlabeled ratio "*r[u]~ "[log]")) + ylab(bquote("observed labeled ratio "*r[l]~ "[log]"))

p1

nrep  <- 3

pre.frac <- data.WT10[,seq(nrep+1,2*nrep),with=F]/data.WT10[,seq(1:nrep),with=F]
lab.frac  <- data.WT10[,seq(3*nrep+1,4*nrep),with=F]/data.WT10[,seq(2*nrep+1,3*nrep),with=F]
data.WT10[,paste("pre.frac",seq(nrep),sep=".") := pre.frac]
data.WT10[,paste("lab.frac",seq(nrep),sep=".") := lab.frac]

### Fig 3 left - inset
p2 <- ggplot(data.WT10[biotype=="protein_coding" & is.finite(lab.frac.1)],aes(x=lab.frac.1>1,y=log(P_1_WT_10_8d.exon),fill=lab.frac.1>1))+geom_boxplot(notch=T)+theme_classic()
p2 <- p2 + xlab("")+ylab("expression [log TPM]")+scale_x_discrete(labels=c(expression(r[l] <= 1),expression(r[l] > 1))) +theme(legend.position = "none",text = element_text(size = 16))+scale_fill_manual(values=c("#FFCCCC", "#CCFFFF"))
p3  <- p1+ annotation_custom(ggplotGrob(p2),xmin=-1.5,xmax=3,ymin=-8,ymax=-3)
p3
#dev.copy2pdf(file=figpath("data_phase.pdf"),onefile=T)


## classifying transcripts:
## C: solvable with a single solution
## A: potentially ambiguous solutions
## B: bad and discarded
## L: limit cases, k may be infered but not the rates (rates are too fast)

data.WT10[,solvable1:=ifelse(lab.frac.1>1/(2-pre.frac.1),"C","A")]
data.WT10[,solvable2:=ifelse(lab.frac.2>1/(2-pre.frac.2),"C","A")]
data.WT10[,solvable3:=ifelse(lab.frac.3>1/(2-pre.frac.3),"C","A")]


EPS = exp(-30)
data.WT10[pre.frac.1 < EPS,solvable1:="L"]
data.WT10[pre.frac.2 < EPS,solvable2:="L"]
data.WT10[pre.frac.3 < EPS,solvable3:="L"]

data.WT10[lab.frac.1>1 | pre.frac.1 >= lab.frac.1,solvable1:="B"]
data.WT10[lab.frac.2>1 | pre.frac.2 >= lab.frac.2,solvable2:="B"]
data.WT10[lab.frac.3>1 | pre.frac.3 >= lab.frac.3,solvable3:="B"]


pl <- length(setdiff(data.WT10[biotype=="protein_coding" &  P_1_WT_10_8d.exon>1,symbol],data.WT10[biotype=="protein_coding" &  P_1_WT_10_8d.exon>1 & solvable1 %in% c("A","C","L"),symbol]))/length(unique(data.WT10[biotype=="protein_coding" &  P_1_WT_10_8d.exon>1,symbol]))
print(paste("proportion of genes with bad ratio",pl))

table(data.WT10$solvable1)
table(data.WT10$solvable2)
table(data.WT10$solvable3)


tt <- table(data.WT10[biotype=="protein_coding"& L_1_WT_10_8d.exon>10 & L_1_WT_10_8d.intron>10,solvable1])
(tt[2]+tt[4])/sum(tt)

tt <- table(data.WT10[ L_1_WT_10_8d.exon>10 & L_1_WT_10_8d.intron>10,solvable1])
(tt[2]+tt[4])/sum(tt)

cl <- NULL ## change this to cl  <- init.cluster(10) if you want to work in parallel

## computing rates ####

refactor <- function(mat){return(list(mat[1,],mat[2,]))}
refactor2 <- function(mat){return(list(mat[,1],mat[,2]))}
refactor4 <- function(mat){return(list(mat[1,],mat[2,],mat[3,],mat[4,]))}
refactor3 <- function(mat){return(list(mat[1,],mat[2,],mat[3,]))}

## computing processing and degradation rates for solvable cases
data.WT10[solvable1=="C",c("deg.rate1","proc.rate1"):=refactor(mapply(solve.rates2,pre.frac.1,lab.frac.1))]
data.WT10[solvable2=="C",c("deg.rate2","proc.rate2"):=refactor(mapply(solve.rates2,pre.frac.2,lab.frac.2))]
data.WT10[solvable3=="C",c("deg.rate3","proc.rate3"):=refactor(mapply(solve.rates2,pre.frac.3,lab.frac.3))]

## computing processing and rates for potentially ambiguous cases
data.WT10[solvable1=="A",c("deg.rate1","proc.rate1","deg.rate1b","proc.rate1b"):=refactor4(mapply(solve.rates3,pre.frac.1,lab.frac.1))]
data.WT10[solvable2=="A",c("deg.rate2","proc.rate2","deg.rate2b","proc.rate2b"):=refactor4(mapply(solve.rates3,pre.frac.2,lab.frac.2))]
data.WT10[solvable3=="A",c("deg.rate3","proc.rate3","deg.rate3b","proc.rate3b"):=refactor4(mapply(solve.rates3,pre.frac.3,lab.frac.3))]

## finding out ambiguous cases
data.WT10[,twosols.1:= !is.na(deg.rate1b)]
data.WT10[,twosols.2:= !is.na(deg.rate2b)]
data.WT10[,twosols.3:= !is.na(deg.rate3b)]

## optimization procedure for unsolvable cases
data.WT10[solvable1=="A" & !twosols.1, c("deg.rate1","proc.rate1"):=refactor2(eval.rates.vec(pre.frac.1,lab.frac.1,cl=cl))]
data.WT10[solvable2=="A" & !twosols.2, c("deg.rate2","proc.rate2"):=refactor2(eval.rates.vec(pre.frac.2,lab.frac.2,cl=cl))]
data.WT10[solvable3=="A" & !twosols.3, c("deg.rate3","proc.rate3"):=refactor2(eval.rates.vec(pre.frac.3,lab.frac.3,cl=cl))]


## computing production rates
data.WT10[,prod.rate1:= get.production.rate(log(P_1_WT_10_8d.exon-P_1_WT_10_8d.intron),deg.rate1,proc.rate1)]
data.WT10[,prod.rate2:= get.production.rate(log(P_2_WT_10_8d.exon-P_2_WT_10_8d.intron),deg.rate2,proc.rate2)]
data.WT10[,prod.rate3:= get.production.rate(log(P_3_WT_10_8d.exon-P_3_WT_10_8d.intron),deg.rate3,proc.rate3)]

## also for ambiguous cases
data.WT10[twosols.1==T,prod.rate1b:= get.production.rate(log(P_1_WT_10_8d.exon-P_1_WT_10_8d.intron),deg.rate1b,proc.rate1b)]
data.WT10[twosols.2==T,prod.rate2b:= get.production.rate(log(P_2_WT_10_8d.exon-P_2_WT_10_8d.intron),deg.rate2b,proc.rate2b)]
data.WT10[twosols.3==T,prod.rate3b:= get.production.rate(log(P_3_WT_10_8d.exon-P_3_WT_10_8d.intron),deg.rate3b,proc.rate3b)]

## a simple heuristics to guess which of the two solutions is correct, take the one with the smaller the production rate
swap.ambiguous <- TRUE
if(swap.ambiguous){
    swap.elmt <- function(a1,a2,idx){
        tmp <- a1[idx]
        a1[idx] <- a2[idx]
        a2[idx] <- tmp
        return(list(a1,a2))
    }
    swap <- function(a1,b1,c1,a2,b2,c2){
        idx <- which(a1>a2)
        la <- swap.elmt(a1,a2,idx)
        lb <- swap.elmt(b1,b2,idx)
        lc <- swap.elmt(c1,c2,idx)
        return(list(la[[1]],lb[[1]],lc[[1]],la[[2]],lb[[2]],lc[[2]]))
    }
    data.WT10[twosols.1==T,c("prod.rate1","deg.rate1","proc.rate1","prod.rate1b","deg.rate1b","proc.rate1b"):=swap(prod.rate1,deg.rate1, proc.rate1,prod.rate1b,deg.rate1b, proc.rate1b)]
    data.WT10[twosols.2==T,c("prod.rate2","deg.rate2","proc.rate2","prod.rate2b","deg.rate2b","proc.rate2b"):=swap(prod.rate2,deg.rate2, proc.rate2,prod.rate2b,deg.rate2b, proc.rate2b)]
    data.WT10[twosols.3==T,c("prod.rate3","deg.rate3","proc.rate3","prod.rate3b","deg.rate3b","proc.rate3b"):=swap(prod.rate3,deg.rate3, proc.rate3,prod.rate3b,deg.rate3b, proc.rate3b)]
}
k.name <- paste("k",seq(3),sep=".")
k.name  <- c(k.name,paste0(k.name,"b"))
data.WT10[, c("k.1","k.2","k.3","k.1b","k.2b","k.3b") := list(proc.rate1-deg.rate1,proc.rate2-deg.rate2,proc.rate3-deg.rate3,proc.rate1b-deg.rate1b,proc.rate2b-deg.rate2b,proc.rate3b-deg.rate3b)]

## gathering some stats to compare across replicates
data.WT10[, sd.k:= apply(cbind(proc.rate1-deg.rate1,proc.rate2-deg.rate2,proc.rate3-deg.rate3),1,sd)]
data.WT10[, sd.proc:= apply(cbind(proc.rate1,proc.rate2,proc.rate3),1,sd)]
data.WT10[, sd.deg:= apply(cbind(deg.rate1,deg.rate2,deg.rate3),1,sd)]
data.WT10[, sd.prod:= apply(cbind(prod.rate1,prod.rate2,prod.rate3),1,sd)]
data.WT10[, avg.deg:= apply(cbind(deg.rate1,deg.rate2,deg.rate3),1,mean)]
data.WT10[, avg.k:= apply(cbind(proc.rate1-deg.rate1,proc.rate2-deg.rate2,proc.rate3-deg.rate3),1,mean)]
data.WT10[, avg.proc:= apply(cbind(proc.rate1,proc.rate2,proc.rate3),1,mean)]
data.WT10[, avg.prod:= apply(cbind(prod.rate1,prod.rate2,prod.rate3),1,mean)]

## Fig 3, right
p1 <- ggplot(data.WT10[biotype=="protein_coding"],aes(x=prod.rate1,y=proc.rate1))+ scale_x_continuous(limits = c(-2, 8)) + scale_y_continuous(limits = c(-5, 5))
p1 <- p1 +  geom_point(pch=19,alpha = 0.05,cex=1) +scale_color_gradient2(low="black",mid="blue",high="red",midpoint=5,na.value = "grey")
p1 <- p1 + theme_classic() +theme(text = element_text(size=20))+coord_fixed(ratio = 1)
p1 <- p1 + labs(x="synthesis rate [log]", y="processing rate [log]",color="intron size [log]")
p1 <- p1 + geom_density_2d(color="wheat", bins=8)
p1

grid.arrange(p3,p1,nrow=1,widths=c(1,1))
dev.copy2pdf(file=figpath("real_data_rates.pdf"),onefile=T)



### Suppl. Figure 

data.WT10[,prod.rate1b:= get.production.rate.labeled(log(L_1_WT_10_8d.exon-L_1_WT_10_8d.intron),deg.rate1,proc.rate1)]

## data.WT10[,prod.rate1b:= get.norm.factor(log(L_1_WT_10_8d.exon),prod.rate1,deg.rate1,proc.rate1)+prod.rate1]

data.WT10[solvable1=="A",exact1:= abs(func(k.1,pre.frac.1,lab.frac.1)) < 0.0001]
data.WT10[solvable1=="C",exact1:=TRUE]

rs <- data.WT10[exact1==T & is.finite(prod.rate1)& is.finite(prod.rate1b),cor(prod.rate1,prod.rate1b,use="p")]

p1 <- ggplot(data.WT10[exact1==T],aes(x=prod.rate1,y=prod.rate1b,color=deg.rate1))+geom_point(alpha = 0.08) +scale_x_continuous(limits = c(-6, 7)) + scale_y_continuous(limits = c(-6, 7))+theme_classic() + scale_color_gradient2(low="green",mid="blue",high="red",midpoint=-4,limits=c(-12,NA)) + geom_abline()+theme(text = element_text(size = 20),legend.position=c(0.2,0.8))+ labs(y="synthesis rate estimated from labeled RNA [log]",x="synthesis rate estimated from unlabeled RNA [log]",color = "degradation rate [log]") + annotate("text",x=4,y=-4,label=paste("R =",format(100*rs,digits=0),"%"),size=8)

p1

dev.copy2pdf(file=figpath("synthesis_rates2.pdf"),onefile=T)

### Figure 4  ################33


## projecting rates on the abundance and reactivity axes for 
data.WT10[biotype=="protein_coding" &  !is.na(deg.rate1),concentration:= (prod.rate1-deg.rate1)/2]
data.WT10[biotype=="protein_coding" & !is.na(deg.rate1),reactivity:= (prod.rate1+deg.rate1)/2]

print("correlations:")
data.WT10[biotype=="protein_coding" & !is.na(deg.rate1),cor(concentration,reactivity)]
data.WT10[biotype=="protein_coding" &  !is.na(deg.rate1),cor(deg.rate1,prod.rate1)]


go.cat <- c("transcription", "monosaccharide metabolism")
go.fnames <- list()
go.fnames[["transcription"]] <- "transcription_go.csv"
go.fnames[["monosaccharide metabolism"]] <- "monosaccharide_met_go.csv"

go.dir <- "./data/"
evid <- c("EXP","IMP","IDA")

cat.pos <- data.table(cat=go.cat)
data.WT10[,go_func:=as.character(NA)]

for (gcat in go.cat){
#    data.WT10[,go_func:=as.character(NA)]
    gene.list <- fread(paste0(go.dir,go.fnames[[gcat]]))
    keep <- which(gene.list[["GO EVIDENCE CODE"]] %in% evid)
    data.WT10[symbol %in% gene.list$SYMBOL[keep],go_func:=ifelse(is.na(go_func),gcat,"pleiotropic")]

    cat.pos[cat==gcat,reactivity:= data.WT10[biotype=="protein_coding" & go_func==gcat,median(reactivity,na.rm=T)]]
    ss <- wilcox.test(data.WT10[biotype=="protein_coding"  & go_func==gcat,reactivity],data.WT10[biotype=="protein_coding" & is.na(go_func),reactivity])
    cat.pos[cat==gcat,reactivity.pval:=ss$p.value]
    cat.pos[cat==gcat,concentration:= data.WT10[biotype=="protein_coding" & go_func==gcat,median(concentration,na.rm=T)]]
    ss <- wilcox.test(data.WT10[biotype=="protein_coding" & go_func==gcat,concentration],data.WT10[biotype=="protein_coding"  & is.na(go_func),concentration])
    cat.pos[cat==gcat,concentration.pval:=ss$p.value]
    cat.pos[cat==gcat,n:=data.WT10[biotype=="protein_coding" & go_func==gcat,.N]]
}
data.WT10[go_func=="pleiotropic",go_func:=NA]

## plotting Fig. 4 left
ax2 <- data.table(from.x=c(0,0),from.y=c(0,0),to.x=c(3,4),to.y=c(3,-4),y.off=c(1,-1),label=c("responsiveness","abundance"))
p1 <- ggplot(data.WT10[biotype=="protein_coding"],aes(x=prod.rate1,y=deg.rate1)) + scale_x_continuous(limits = c(-3, 10)) + scale_y_continuous(limits = c(-5.5, 4.5))
p1 <- p1 + geom_abline(slope=1,intercept=seq(-8,0,2),alpha=0.1,linetype="solid",color="brown")+geom_abline(slope=-1,intercept=seq(-4,6,2),alpha=0.1,linetype="dashed",color="brown")
p1 <- p1 + geom_point(alpha = 0.02,color="black")
p1 <- p1 + geom_point(data=data.WT10[biotype=="protein_coding"& !is.na(go_func)],aes(color=go_func),pch=19,alpha = 1,cex=1)
p1 <- p1 + labs(x="synthesis rate [log]", y="degradation rate[log]",color="GO categories")
p1 <- p1 + theme_classic() +theme(legend.position=c(0.8,0.2),text = element_text(size=20),legend.text = element_text(size = 12),legend.title = element_text(size = 16))+coord_fixed(ratio = 1)
p1 <- p1 + geom_segment(aes(x=from.x,y=from.y,xend=to.x,yend=to.y),color="black",data=ax2,arrow = arrow(length = unit(0.2, "cm")),lineend = "round") + geom_text(aes(x=to.x,y=to.y+0.5*sign(to.y),label=label),data=ax2,size=6)
p1
med = data.WT10[biotype=="protein_coding",list(concentration=mean(concentration,na.rm=T),reactivity=mean(reactivity,na.rm=T),prod.rate=mean(prod.rate1,na.rm=T),deg.rate=mean(deg.rate1,na.rm=T) ) ,by=go_func]
p1 <- p1 + geom_point(data=med,aes(x=concentration,y=-concentration,color=go_func),pch=15,cex=3)
p1 <- p1 + geom_point(data=med,aes(x=reactivity,y=reactivity,color=go_func),pch=15,cex=3)
p1

## plotting fig 4 right
p2 <- ggplot(data.WT10[biotype=="protein_coding"],aes(x=concentration,y=reactivity)) + scale_x_continuous(limits = c(-0.5, 4.5)) + scale_y_continuous(limits = c(-5, 5))
p2 <- p2 + geom_point(data=data.WT10[biotype=="protein_coding"& !is.na(go_func)],aes(color=go_func),pch=19,alpha = 0.5,cex=1)
p2 <- p2 + labs(x="steady-state abundance [log]", y="responsiveness [log]")
p2 <- p2 + theme_classic() +theme(legend.position = "none",text = element_text(size=20),legend.text = element_text(size = 12),legend.title = element_text(size = 16))+coord_fixed(ratio = 1,clip="off",expand=F)
p2 <- p2 + geom_point(data=med,aes(x=concentration,y=reactivity,color=go_func),pch=17,cex=3)
p2 <- p2 + geom_hline(aes(yintercept=reactivity,color=go_func),data=med,linetype="dotted")+ geom_vline(aes(xintercept=concentration,color=go_func),data=med,linetype="dotted")
p2 <- p2 + geom_point(data=med,aes(y=reactivity,x=-Inf,color=go_func),pch=15,cex=3)+geom_point(data=med,aes(y=-Inf,x=concentration,color=go_func),pch=15,cex=3)
p2

grid.arrange(p1,p2,nrow=1,widths=c(2.3,1))

dev.copy2pdf(file=figpath("deg_synt_rates.pdf"),onefile=T)

## checking statistical significance

wilcox.test(data.WT10[biotype=="protein_coding" & go_func=="monosaccharide metabolism",reactivity],data.WT10[biotype=="protein_coding" & is.na(go_func),reactivity])
wilcox.test(data.WT10[biotype=="protein_coding" & go_func=="monosaccharide metabolism",concentration],data.WT10[biotype=="protein_coding" & is.na(go_func),concentration])
wilcox.test(data.WT10[biotype=="protein_coding" & go_func=="transcription",reactivity],data.WT10[biotype=="protein_coding" & is.na(go_func),reactivity])
wilcox.test(data.WT10[biotype=="protein_coding" & go_func=="transcription",concentration],data.WT10[biotype=="protein_coding" & is.na(go_func),concentration])


data.WT10[biotype=="protein_coding",cor(prod.rate1,proc.rate1,use="p")]
data.WT10[biotype=="protein_coding",cor(prod.rate1,deg.rate1,use="p")]


############ Figure 5 ##########################

## comparing with Herzog et al.

data.dir <- "./data/"
am.fname <- "rates_herzog.txt"
adata <- fread(paste(data.dir,am.fname,sep="/"),fill=T)
colnames(adata)[which(colnames(adata)=="Half-life (h)")] <- "Half.life"
res <- select(org.Mm.eg.db, keys=adata$Name , columns=c("SYMBOL","ENSEMBLTRANS"), keytype="SYMBOL")
ma.idx <- match(rownames(data.WT10),res$ENSEMBLTRANS)
data.WT10[,gene.name:= res$SYMBOL[ma.idx]]
data.WT10[,sum1:=P_1_WT_10_8d.exon+L_1_WT_10_8d.exon]
data.WT10[,sum2:=P_2_WT_10_8d.exon+L_2_WT_10_8d.exon]
data.WT10[,sum3:=P_3_WT_10_8d.exon+L_3_WT_10_8d.exon]
data.WT10b = merge(data.WT10[biotype=="protein_coding"],adata,by.x="gene.name",by.y="Name",all.y=T)



aa1 <- data.WT10b[solvable1 %in% c("A","C"),list(fin.deg=sum(deg.rate1*sum1,na.rm=T)/sum(sum1,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),msum=sum(sum1,na.rm=T),solvable=solvable1[which.max(sum1)],frac.solv=sum(sum1[solvable1=="A"])/sum(sum1),lab.frac=sum(lab.frac.1*sum1)/sum(sum1),pre.frac=sum(pre.frac.1*sum1)/sum(sum1),k=sum(k.1*sum1)/sum(sum1)),by=gene.name]

aa2 <- data.WT10b[solvable2 %in% c("A","C"),list(fin.deg=sum(deg.rate2*sum2,na.rm=T)/sum(sum2,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),msum=sum(sum2,na.rm=T),solvable=solvable2[which.max(sum2)],frac.solv=sum(sum2[solvable2=="A"])/sum(sum2),lab.frac=sum(lab.frac.2*sum2)/sum(sum2),pre.frac=sum(pre.frac.2*sum2)/sum(sum2),k=sum(k.2*sum2)/sum(sum2)),by=gene.name]

aa3 <- data.WT10b[solvable3 %in% c("A","C"),list(fin.deg=sum(deg.rate3*sum3,na.rm=T)/sum(sum3,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),msum=sum(sum3,na.rm=T),solvable=solvable3[which.max(sum3)],frac.solv=sum(sum3[solvable3=="A"])/sum(sum3),lab.frac=sum(lab.frac.3*sum3)/sum(sum3),pre.frac=sum(pre.frac.3*sum3)/sum(sum3),k=sum(k.3*sum3)/sum(sum3)),by=gene.name]

data.WT10b[,asum:=(sum1+sum2+sum3)/3]
data.WT10[,rn := rownames(data.WT10)]

# Fig 5, left

thresh <- 200
ss <- summary(lm(fin.deg+log(6) ~ slam.deg,data=aa1[msum>thresh],weights=1-lab.frac))

pv <- ggplot(aa1[msum>thresh],aes(y=fin.deg+log(6),x=slam.deg))+ geom_point(aes(alpha=1-lab.frac))
pv <- pv + geom_abline(slope=ss$coefficient[2,1],intercept=ss$coefficient[1,1],colour="red")+ theme_classic()+annotate("text",x=-0.75,y=-4,label=paste("R =",format(100*sqrt(ss$r.squared),digits=0),"%"),size=8)+ labs(x="slam-seq estimate", y="single sample estimate ", title="degradation rate [log]")+ guides(alpha=FALSE)+theme(text = element_text(size=20))
pv



#### Fig. 5, right 

cor.func <- function(th,aa){
     ss <- summary(lm(fin.deg+log(6) ~ slam.deg,data=aa[msum>th & is.finite(fin.deg) & lab.frac<1],weights=1-lab.frac))
     return(c(sqrt(ss$r.squared)*sign(coef(ss)[2,1]),ss$df[2], coef(ss)[2,1]))
    ##return (c(cor(aa[msum>th,fin.deg],aa[msum>th,slam.deg],method="spearman", use= "pairwise.complete.obs"),sum(aa[msum>th,is.finite(fin.deg+slam.deg) ])))  
}
cor.data <- data.table(thresh=c(0,100,200,500,1000))
cort <- matrix(NA,5,3)
dft  <- matrix(NA,5,3)
##laa  <- list(aa1,aa2,aa3,aa4,aa5)
laa  <- list(aa1,aa2,aa3)
for (j in seq(length(laa))){
    var =paste0(c("cor","df","slope"),j)
    fres  <- sapply(seq(5),function(i)cor.func(cor.data$thresh[i],laa[[j]]))
    for (i in seq(3)){
        cor.data[,(var[i]) := fres[i,]]
        }
}
cor.data[,quant1:=df1/df1[1]]
cor.data[,quant2:=df2/df2[1]]
cor.data[,quant3:=df3/df3[1]]

pc <- ggplot(cor.data,aes(x=thresh/2))+geom_line(aes(y=cor1,color="replicate 1"))+geom_line(aes(y=cor2,color="replicate 2"))+geom_line(aes(y=cor3,color="replicate 3"))+scale_y_continuous(lim=c(0,0.7))+theme_classic()+theme(text = element_text(size=20),legend.position=c(0.8,0.5),axis.text=element_text(angle=0,size=12)) +geom_point(data=cor.data[thresh==200,],mapping=aes(y=cor1, color="replicate 1"),size=2)
pc <- pc+ scale_x_continuous(breaks=cor.data[,thresh/2],label=cor.data[,paste0(thresh/2,"\n(",format(100-100*(quant1+quant2+quant3)/3,digits=2),"%)")],name="expression threshold [TPM (and quantiles)]") + labs(x="expression threshold [TPM (and quantiles)]",y="correlation",color="")
pc

grid.arrange(pv,pc,nrow=1)

dev.copy2pdf(file=figpath("slamseq_comp.pdf"),onefile=T)


############# Figure 6  and sup Fig D8 ##############################333

### generate inspect rates

if (with.inspect.comp){
    library(INSPEcT)
    sset = c(1:nrow(data.WT10))


    data.WT10[,p.exon.var:= apply(cbind(P_1_WT_10_8d.exon,P_2_WT_10_8d.exon,P_3_WT_10_8d.exon),1,var)]
 lfit <- loess(log(p.exon.var)~log(P_1_WT_10_8d.exon),data = data.WT10[is.finite(log(p.exon.var)+log(P_1_WT_10_8d.exon))],span=0.1)
 p.smooth.exon.var=predict(lfit,newdata=data.WT10)
 data.WT10[,p.smooth.exon.var:=exp(p.smooth.exon.var)]
 
 data.WT10[,p.intron.var:= apply(cbind(P_1_WT_10_8d.intron,P_2_WT_10_8d.intron,P_3_WT_10_8d.intron),1,var)]
 lfit <- loess(log(p.intron.var)~log(P_1_WT_10_8d.intron),data = data.WT10[is.finite(log(p.intron.var)+log(P_1_WT_10_8d.intron))],span=0.1)
 p.smooth.intron.var=predict(lfit,newdata=data.WT10)
 data.WT10[,p.smooth.intron.var:=exp(p.smooth.intron.var)]
 


 data.WT10[,l.exon.var:= apply(cbind(L_1_WT_10_8d.exon,L_2_WT_10_8d.exon,L_3_WT_10_8d.exon),1,var)]
 lfit <- loess(log(l.exon.var)~log(L_1_WT_10_8d.exon),data = data.WT10[is.finite(log(l.exon.var)+log(L_1_WT_10_8d.exon))],span=0.1)
 l.smooth.exon.var=predict(lfit,newdata=data.WT10)
 data.WT10[,l.smooth.exon.var:=exp(l.smooth.exon.var)]

 data.WT10[,l.intron.var:= apply(cbind(L_1_WT_10_8d.intron,L_2_WT_10_8d.intron,L_3_WT_10_8d.intron),1,var)]
 lfit <- loess(log(l.intron.var)~log(L_1_WT_10_8d.intron),data = data.WT10[is.finite(log(l.intron.var)+log(L_1_WT_10_8d.intron))],span=0.1)
 l.smooth.intron.var=predict(lfit,newdata=data.WT10)
 data.WT10[,l.smooth.intron.var:=exp(l.smooth.intron.var)]



#plot(data.WT10[is.finite(log(p.exon.var)+log(P_1_WT_10_8d.exon)),log(P_1_WT_10_8d.exon)],logvar)
ggplot(data.WT10,aes(x=log(P_1_WT_10_8d.exon),y=0.5*log(p.exon.var)))+geom_point(alpha=0.05)+theme_classic() + geom_smooth(span=0.1)+ geom_abline(aes(slope=1,intercept=0),color="red")
    
 ggplot(data.WT10,aes(x=log(P_1_WT_10_8d.intron),y=0.5*log(p.intron.var)))+geom_point(alpha=0.05)+theme_classic() + geom_smooth(span=0.1) + geom_abline(aes(slope=1,intercept=0),color="red")
 ggplot(data.WT10,aes(x=log(L_1_WT_10_8d.exon),y=0.5*log(l.exon.var)))+geom_point(alpha=0.05)+theme_classic() + geom_smooth(span=0.1)+ geom_abline(aes(slope=1,intercept=0),color="red")
 ggplot(data.WT10,aes(x=log(L_1_WT_10_8d.intron),y=log(l.intron.var)))+geom_point(alpha=0.05)+theme_classic() + geom_smooth(span=0.1)+ geom_abline(aes(slope=1,intercept=0),color="red")





## gfit <- gam(log(p.exon.var)~log(P_1_WT_10_8d.exon),data = data.WT10[is.finite(log(p.exon.var)+log(P_1_WT_10_8d.exon))])
## sm <- gam(gfit,data=data.WT10)

 generate.inspect.rates <- function(fname,idata.WT10,sset,benchmark=FALSE){
     insp.rates <- list()
     for (i in seq(1,3)){
         p.ex <- paste("P",i,"WT_10_8d.exon",sep="_")
         p.intr <- paste("P",i,"WT_10_8d.intron",sep="_")
         l.ex <- paste("L",i,"WT_10_8d.exon",sep="_")
         l.intr <- paste("L",i,"WT_10_8d.intron",sep="_")
         
         exonExpMat <- matrix(idata.WT10[[p.ex]][sset],ncol=1)
         intronExpMat  <- matrix(idata.WT10[[p.intr]][sset],ncol=1)
         exonVarMat  <- matrix(idata.WT10[sset,p.smooth.exon.var] ,ncol=1,nrow=length(sset))
         intronVarMat  <- matrix(idata.WT10[sset,p.smooth.intron.var] ,ncol=1,nrow=length(sset))
         rownames(exonExpMat) <- idata.WT10$rn[sset]
         rownames(intronExpMat) <- idata.WT10$rn[sset] 
         rownames(exonVarMat) <- idata.WT10$rn[sset] 
         rownames(intronVarMat) <- idata.WT10$rn[sset] 

         exonExpNac <- matrix(idata.WT10[[l.ex]][sset],ncol=1)
         intronExpNac  <- matrix(idata.WT10[[l.intr]][sset],ncol=1)
         exonVarNac  <- matrix(idata.WT10[sset,l.smooth.exon.var] ,ncol=1,nrow=length(sset))
         intronVarNac  <- matrix(idata.WT10[sset,l.smooth.intron.var] ,ncol=1,nrow=length(sset))
         
         rownames(exonExpNac) <- idata.WT10$rn[sset]
         rownames(intronExpNac) <- idata.WT10$rn[sset]
         rownames(exonVarNac) <- idata.WT10$rn[sset]
         rownames(intronVarNac) <- idata.WT10$rn[sset]
         
         
         matureExp  <- list(exonsExpressions=exonExpMat, intronsExpressions=intronExpMat, exonsVariance=exonVarMat,intronsVariance=intronVarMat)
         nacentExp  <- list(exonsExpressions=exonExpNac, intronsExpressions=intronExpNac, exonsVariance=exonVarNac,intronsVariance=intronVarNac)
         atoc <- system.time({
         inspObj<-newINSPEcT(tpts=c(1/6),labeling_time=1/6,nascentExpressions=nacentExp,matureExpressions=matureExp,degDuringPulse=FALSE, preexisting=TRUE)
         insp.rates[[i]] <- data.frame(rn = ratesFirstGuess(inspObj,'name'))
         insp.rates[[i]][[paste("isynthesis",i,sep=".")]] <- ratesFirstGuess(inspObj,'synthesis')
         insp.rates[[i]][[paste("iprocessing",i,sep=".")]] <- ratesFirstGuess(inspObj,'processing')           
         insp.rates[[i]][[paste("idegradation",i,sep=".")]] <- ratesFirstGuess(inspObj,'degradation')
         })
         if (benchmark){
           return(c(n.inspect=length(sset),elapsed.inspect=atoc[1]+atoc[2],n.out.inspect=nrow(insp.rates[[i]])))
        }else{
           write.csv(insp.rates[[i]],file=paste0(fname,i,".csv"))
        }
        
     }
 }

    bench.ssre <- function(data.WT10,sset){
        TPM.data <-data.frame(unlabeled.intron=data.WT10[sset,P_1_WT_10_8d.intron],unlabeled.exon=data.WT10[sset,P_1_WT_10_8d.exon],labeled.intron=data.WT10[sset,L_1_WT_10_8d.intron],labeled.exon=data.WT10[sset,L_1_WT_10_8d.exon])
        #tic()
        atoc <- system.time(ssre <- singleSampleRateEstimation(TPM.data))
        #atoc <- toc()
        return(c(n.ssre=length(sset),elapsed.ssre=atoc[1]+atoc[2],n.out.ssre=ssre[is.finite(deg.rate),.N]))
    }
    
benchmarking  <- FALSE ## set to TRUE to compare cpu time
inspfname <- "insp_rates_loess_p" # base filename for inspect rates
nbtrans <- nrow(data.WT10)

 if (benchmarking){
    sample.sizes  <- c(10,100,1000,10000,nbtrans)
    bench <- c()
    for (samp.size in sample.sizes[1:4]){
        sset <- sample(nbtrans,samp.size)
        b.insp <- generate.inspect.rates(inspfname,data.WT10,sset,benchmark=T)
        b.insp
        b.ssre <- bench.ssre(data.WT10,sset)
        b.ssre
        bench <- rbind(bench,c(b.insp,b.ssre))
    }
    time.ratio <- (bench[,2]*bench[,6])/(bench[,5]*bench[,3])    
    bench <- cbind(bench,time.ratio)
    bench
 }else{
     sset  <- seq(nbtrans)
     generate.inspect.rates(inspfname,data.WT10,sset,benchmark=F)
 }

 
 tmp1  <- data.WT10
 for (i in seq(3)){
     insp.rates <- fread(file=paste0(inspfname,i,".csv"))
     names(insp.rates)[1] <- "txid"
     insp.rates[,2:4]  <-log(insp.rates[,2:4])-log(6) ## - log(6) to move from hour^-1 to 10min ^-1 units
     tmp1 <- merge(tmp1,insp.rates,by=c("txid"),all=T)
 }
 comp.insp <- tmp1
 
 al.val  <- 0.05
 thresh  <-  0
 cond  <- which(comp.insp[,biotype=="protein_coding"])
 irate.name  <- c("isynthesis","iprocessing","idegradation")
 rate.name  <- c("prod.rate","proc.rate","deg.rate")
 irate.label <- c("synthesis","processing", "degradation")
 mlims  <- list(c(-7,8),c(-4,3),c(-6,3))
 imlims  <- list(c(0,8),c(-1,5),c(-2,2))
 replid <- c(2,3,1)
 cpi <- list()
 cp <- list()
 ccp  <- list()
 for (i in seq(3)){ # for all three rates    
     ## plot replicates for inspect
     icol1 <- paste(irate.name[i],replid[1],sep=".")
     icol2 <- paste(irate.name[i],replid[2],sep=".")
     icolx <- sym(icol1)
     icoly <- sym(icol2)
     ilabx <- paste("repl.",replid[1], "INSPEcT","rate [log]")
     ilaby <- paste("repl.",replid[2], "INSPEcT","rate [log]")
     rr  <- cor(comp.insp[cond,get(icol1)],comp.insp[cond,get(icol2)], method="spearman", use= "pairwise.complete.obs")
     cpi[[i]] <-  ggplot(comp.insp[cond],aes(x=!!icolx,y=!!icoly))+geom_point(alpha=al.val)+theme_classic()+geom_abline(slope=1,color="red")+scale_x_continuous(limits=mlims[[i]])+scale_y_continuous(limits=mlims[[i]])+annotate("text",x=Inf,y=-Inf,vjust=0,hjust=1,label=paste("R =",format(100*rr,digits=0),"%"),size=8)+xlab(ilabx)+ylab(ilaby)+theme(axis.title.x = element_text(size = 16),axis.title.y = element_text(size = 16), plot.title=element_text(hjust=0.5,size=26))+labs(title=irate.label[i])
  
    ## plot replicates for SSRE
     col1 <- paste0(rate.name[i],replid[1])
     col2 <- paste0(rate.name[i],replid[2])
     colx <- sym(col1)
     coly <- sym(col2)
     labx <- paste("repl.",replid[1], "SSRE","rate [log]")
     laby <- paste("repl.",replid[2], "SSRE","rate [log]")
     rr  <- cor(comp.insp[cond,get(col1)],comp.insp[cond,get(col2)], method="spearman", use= "pairwise.complete.obs")
     cp[[i]] <-  ggplot(comp.insp[cond],aes(x=!!colx,y=!!coly))+geom_point(alpha=al.val)+theme_classic()+geom_abline(slope=1,color="red")+scale_x_continuous(limits=mlims[[i]])+scale_y_continuous(limits=mlims[[i]])+annotate("text",x=Inf,y=-Inf,vjust=0,hjust=1,label=paste("R =",format(100*rr,digits=0),"%"),size=8)+xlab(labx)+ylab(laby)+theme(axis.title.x = element_text(size = 16),axis.title.y = element_text(size = 16))
     
    ## plot inspect vs SSRE estimates
     col1 <- paste0(rate.name[i],replid[3])
     icol2 <- paste(irate.name[i],replid[3],sep=".")
     colx <- sym(col1)
     icoly <- sym(icol2)
     labx <- paste("repl.",replid[3], "SSRE",irate.label[i],"rate [log]")
     ilaby <- paste("repl.",replid[3], "INSPEcT",irate.label[i],"rate [log]")
     labx <- "SSRE"
     ilaby <- "INSPEcT"
     rr  <- cor(comp.insp[cond,get(col1)],comp.insp[cond,get(icol2)], method="spearman", use= "pairwise.complete.obs")
     ccp[[i]] <- ggplot(comp.insp[cond],aes(x=!!colx,y=!!icoly))+geom_point(alpha=al.val)+theme_classic()+geom_abline(slope=1,color="red")+scale_x_continuous(limits=mlims[[i]])+scale_y_continuous(limits=imlims[[i]])+annotate("text",x=Inf,y=-Inf,vjust=0,hjust=1,label=paste("R =",format(100*rr,digits=0),"%"),size=8)+xlab(labx)+ylab(ilaby)+theme(axis.title.x = element_text(size = 16),axis.title.y = element_text(size = 16), plot.title=element_text(hjust=0.5,size=26) )+labs(title=irate.label[i])
}


 spi <- list()
 icol1 <- paste(irate.name[1],replid[3],sep=".")
 icol2 <- paste(irate.name[2],replid[3],sep=".")
 icolx <- sym(icol1)
 icoly <- sym(icol2)
 ilabx <- paste("INSPEcT",irate.label[1],"rate [log]")
 ilaby <- paste("INSPEcT",irate.label[2],"rate [log]")
 rr  <- cor(comp.insp[cond,get(icol1)],comp.insp[cond,get(icol2)], method="spearman", use= "pairwise.complete.obs")
 spi[[1]] <- ggplot(comp.insp[cond],aes(x=!!icolx,y=!!icoly))+geom_point(alpha=al.val)+theme_classic()+scale_x_continuous(limits=imlims[[1]])+scale_y_continuous(limits=imlims[[2]])+annotate("text",x=Inf,y= -Inf,vjust=0,hjust=1,label=paste("R =",format(100*rr,digits=0),"%"),size=8)+xlab(ilabx)+ylab(ilaby)+theme(axis.title.x = element_text(size = 16),axis.title.y = element_text(size = 16))
 
 icol2 <- paste(irate.name[3],replid[3],sep=".")
 icoly <- sym(icol2)
 ilabx <- paste("repl.",replid[3], "INSPEcT",irate.label[1],"rate [log]")
 ilaby <- paste("repl.",replid[3], "INSPEcT",irate.label[3],"rate [log]")
 rr  <- cor(comp.insp[cond,get(icol1)],comp.insp[cond,get(icol2)], method="spearman", use= "pairwise.complete.obs")
 spi[[2]] <- ggplot(comp.insp[cond],aes(x=!!icolx,y=!!icoly))+geom_point(alpha=al.val)+theme_classic()+scale_x_continuous(limits=imlims[[1]])+scale_y_continuous(limits=imlims[[3]])+annotate("text",x=Inf,y=-Inf,vjust=0,hjust=1,label=paste("R =",format(100*rr,digits=0),"%"),size=8)+xlab(ilabx)+ylab(ilaby)+theme(axis.title.x = element_text(size = 16),axis.title.y = element_text(size = 16))
 

## compare rates with herzog et al
 comp1 <- merge(comp.insp[biotype=="protein_coding"],adata,by.x="gene.name",by.y="Name",all.y=T)
 comp1[,sum1:=P_1_WT_10_8d.exon+L_1_WT_10_8d.exon]
 comp1[,sum2:=P_2_WT_10_8d.exon+L_2_WT_10_8d.exon]
 comp1[,sum3:=P_3_WT_10_8d.exon+L_3_WT_10_8d.exon]
 
 aa1i <- comp1[,list(fin.deg=sum(idegradation.1*sum1,na.rm=T)/sum(sum1,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),fin.syn=sum(isynthesis.1*sum1,na.rm=T)/sum(sum1,na.rm=T),msum=sum(sum1,na.rm=T),solvable=solvable1[which.max(sum1)],frac.solv=sum(sum1[solvable1=="A"])/sum(sum1),lab.frac=0*sum(lab.frac.1*sum1)/sum(sum1),pre.frac=sum(pre.frac.1*sum1)/sum(sum1),k=sum(k.1*sum1)/sum(sum1)),by=gene.name]
 aa1i[fin.deg==0,fin.deg:=NA]


 aa2i <- comp1[,list(fin.deg=sum(idegradation.2*sum2,na.rm=T)/sum(sum2,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),fin.syn=sum(isynthesis.2*sum2,na.rm=T)/sum(sum2,na.rm=T),msum=sum(sum2,na.rm=T),solvable=solvable2[which.max(sum2)],frac.solv=sum(sum2[solvable2=="A"])/sum(sum2),lab.frac=0*sum(lab.frac.2*sum2)/sum(sum2),pre.frac=sum(pre.frac.2*sum2)/sum(sum2),k=sum(k.2*sum2)/sum(sum2)),by=gene.name]
 aa2i[fin.deg==0,fin.deg:=NA]

 aa3i <- comp1[,list(fin.deg=sum(idegradation.3*sum3,na.rm=T)/sum(sum3,na.rm=T),slam.deg=mean(log(log(2))-log(Half.life),na.rm=T),fin.syn=sum(isynthesis.3*sum3,na.rm=T)/sum(sum3,na.rm=T),msum=sum(sum3,na.rm=T),solvable=solvable3[which.max(sum3)],frac.solv=sum(sum3[solvable3=="A"])/sum(sum3),lab.frac=0*sum(lab.frac.3*sum3)/sum(sum3),pre.frac=sum(pre.frac.3*sum3)/sum(sum3),k=sum(k.3*sum3)/sum(sum3)),by=gene.name]
 aa3i[fin.deg==0,fin.deg:=NA]
 

 icor.data <- data.table(thresh=c(0,100,200,500,1000))
 cort <- matrix(NA,5,3)
 dft  <- matrix(NA,5,3)
laa  <- list(aa1i,aa2i,aa3i)
 for (j in seq(length(laa))){
     var =paste0(c("cor","df","slope"),j)
     fres  <- sapply(seq(5),function(i)cor.func(cor.data$thresh[i],laa[[j]]))
     for (i in seq(3)){
         icor.data[,(var[i]) := fres[i,]]
     }
 }
 icor.data[,quant1:=df1/df1[1]]
 icor.data[,quant2:=df2/df2[1]]
 icor.data[,quant3:=df3/df3[1]]
 
 colsbp <-  c(paste0("cor",seq(3)),paste0("slope",seq(3)))
 bpdata <- data.table(x=seq(3),cor=icor.data[1,c(get(colsbp[1]),get(colsbp[2]),get(colsbp[3]))], slope=icor.data[1,c(get(colsbp[4]),get(colsbp[5]),get(colsbp[6]))], method="INSPEcT")
 bpdata <- rbind(bpdata,data.table(x=seq(3),cor=cor.data[1,c(get(colsbp[1]),get(colsbp[2]),get(colsbp[3]))],slope=cor.data[1,c(get(colsbp[4]),get(colsbp[5]),get(colsbp[6]))], method="SSRE"))
 bpdata[method=="SSRE",si:=1]
 bpdata[method=="INSPEcT",si:=-1]
 pb  <- ggplot(bpdata,aes(x=x,y=cor,fill=method))+geom_bar(stat="identity",position="dodge")+theme_classic()+theme(text = element_text(size=20),legend.position="top",legend.title = element_blank())+xlab("replicate")+ylab("correlation and slope") + scale_fill_manual(values=c("blue", "turquoise"))+geom_point(aes(x=x+si*0.25,y=slope,color=method),size=5)
 
 grid.arrange(cpi[[1]],cpi[[2]],cpi[[3]],cp[[1]],cp[[2]],cp[[3]],nrow=2)
 dev.copy2pdf(file=figpath("replcor.pdf"),onefile=T)
 grid.arrange(ccp[[1]],ccp[[2]],ccp[[3]], pb, spi[[1]],spi[[2]],nrow=2)
 dev.copy2pdf(file=figpath("inspcomp.pdf"),onefile=T)
}