Upload R file

RNA half-life estimates/half_life_calculation20231108.R

+library(MASS)
+library(modelr)
+library(reshape2)
+library(ggplot2)
+library(cowplot)
+
+exampledf <- data.frame(time=c(0, 5, 10, 20, 40), measure=c(0.4, 0.31, 0.22, 0.15, 0.08))
+
+estimate_half_life_modelr_catch <- function(gname, adf, initial_k){
+  # using confm from the MASS package to estimate confidence interval
+  # using rsquared from the modelr package to estimate goodness of fit
+  # Initial guess values for speed of decay and experimental lag
+  initialestimates<-c(initial_k)
+  names(initialestimates) <- c("k")
+  #print(adf$measure)
+  y <- adf$measure/adf$measure[1]
+  x <- adf$time
+  xydf <- data.frame(x=x, y=y)
+  if(class(try((m <- nls(formula= y ~ exp(-(x*k)), 
+                         start=initialestimates))[[2]],
+               silent = TRUE)) != "try-error"){
+    k_estimated <- coef(summary(m))[1]
+    half_life <- log(2)/k_estimated
+    if(class(try((confm <- confint(m, level=0.95))[[2]],
+                 silent = TRUE)) != "try-error"){
+      confmin <- confm[1]
+      h_l_max <- log(2)/confmin
+      confmax <- confm[2]
+      h_l_min <- log(2)/confmax
+      conf_int <- max(c(half_life-h_l_min, h_l_max-half_life))
+    }
+    else{conf_int <- NA}
+    #print(xydf)
+    rsq <- rsquare(model=m, data=xydf)
+    fres <-list(gname=gname, k_est = k_estimated, half_life=half_life, conf_int=conf_int,
+                rsquared=rsq)
+  }
+  else{
+    fres <- list(gname=gname, k_est=NA, half_life=NA, conf_int=NA, rsquared=NA)
+  }
+  return(fres)
+}
+
+test <- estimate_half_life_modelr_catch("mytest", exampledf, 0.1)
+print(test)
+# $gname
+# [1] "mytest"
+# 
+# $k_est
+# [1] 0.0496424
+# 
+# $half_life
+# [1] 13.9628
+# 
+# $conf_int
+# [1] 3.341715
+# 
+# $rsquared
+# [1] 0.9815205
+
+#A simple plot to see how the data are distributed:
+times_fit <- seq(from=0, to=40, by=1)
+fitdf <- data.frame(timesf=times_fit, mydata=NA)
+fitdf$mydata <- exp(-test[["k_est"]]*fitdf$timesf)
+fitdf.m <- melt(fitdf, id.vars = c("timesf"))
+
+dftoplot <- exampledf
+dftoplot$measure <- dftoplot$measure/dftoplot$measure[1]
+
+ggplot(data=dftoplot, aes(x=time, y=measure))+
+  geom_jitter(width=0.3, shape=21, fill=NA, stroke=2)+
+  geom_line(data=fitdf.m, aes(x=timesf, y=value), linewidth=2, alpha=0.5)+
+  scale_fill_viridis()+
+  theme_cowplot()+
+  theme(panel.grid = element_line(color = "gray",
+                                  linewidth = 0.75,
+                                  linetype = 2))+
+  xlab("time (min.)")+
+  ylab("fraction of RNA")+
+  ylim(0, 1)+
+  theme(legend.position="none")