BCB 718 Modeling: HW 3

Homework 3 – BCB 718 – Spring 2021

Construct 5 different functions (.m files) to calculate 5 different features of a time series (10 points). Calculate and report all five features for the time series measurements in Figure 5 of PNAS 99:2399. Use these 15 time points shown for 10 uM IP3

---

---

DATA INITIALIZATION

knitr::opts_chunk$set(echo = TRUE)

There were 25 warnings (use warnings() to see them)

library(ggplot2);library(gridExtra);library(kableExtra);library(gt);library(ggpubr);library(ggseg);library(pROC)

Type 'citation("pROC")' for a citation.

Attaching package: ‘pROC’

The following objects are masked from ‘package:stats’:

    cov, smooth, var

gghisto <- list(
  theme(axis.text.x = element_text(face="bold", color="cornflowerblue", size=14,angle=17),
          axis.text.y = element_text(face="bold", color="royalblue4", 
          size=16, angle=25),
          axis.title=element_text(size=17,face="italic"),
          plot.title = element_text(size=20,face="bold.italic")))
ggbar <- list(
  theme(axis.text.x = element_text(face="bold", color="gray18", size=10,angle=20),
          axis.text.y = element_text(face="bold", color="royalblue4", 
          size=12, angle=25),
          axis.title=element_text(size=12,face="italic"),
          plot.title = element_text(size=16,face="bold.italic")))

x = c(0, 0.07, 0.14, 0.21, 0.28, 0.35, 0.42, 0.49, 0.56, 0.63, 0.7, 0.77, 0.84, 0.91, 0.98) 
y = c(0.02, 0.20, 0.34, 0.23, 0.14, 0.10, 0.08, 0.07, 0.06, 0.05, 0.05, 0.04, 0.04, 0.03, 0.03)
df = as.data.frame(cbind(x,y))

xmax <- df$x[df$y==max(df$y)]
x1 <- df$x[df$x < xmax][which.min(abs(df$y[df$x < xmax]-max(df$y)/2))]
x2 <- df$x[df$x > xmax][which.min(abs(df$y[df$x > xmax]-max(df$y)/2))]

---

---

BASE PLOT WITH 10uM IP3 OVER 0.98 SECONDS

p <- ggplot(df, aes(x=x, y=y)) + 
  geom_point(size=7,shape=21,stroke=1.4,color="lightcoral",fill="lightcoral",alpha=.6)+
  geom_line(size=2,color="lightcoral",alpha=.3) +
  scale_y_continuous(limits = c(0, .4), labels=function(y) paste0(y*100,"%"))+
  scale_x_continuous(limits = c(0, 1), breaks=seq(from=0,to=1,by=0.1))+
  ylab("Open Probability") + xlab("Time (s)") +
  gghisto

p+
  ggtitle("10uM IP3 Dose Response \nover 0.98 seconds")

---

---

1. PROBABILITY AT 0.98 SECONDS

y[x==max(x)]

[1] 0.03

p +
  geom_point(aes(x=max(x),y=y[x==max(x)]), size= 4,colour="red")+
  annotate("text", x=.9,y=.12, label= "Probability at \n0.98 sec = 3%")+
  ggtitle("10uM IP3 - Probability of Receptor\nBeing Open At 0.98 Seconds")

---

---

2. TIME FROM PEAK PROBABILITY

df2 <- df[x>=xmax,]
max(df2$x)-min(df2$x)

[1] 0.84

p +
  geom_point(data=df2, size= 4,colour="red")+
  geom_line(data=df2, size=.6,color="red",alpha=.5)+
  annotate("text", x=.4,y=.29, label= "Time from Peak = 0.84 Sec")+
  ggtitle("10uM IP3 - Time remainder \nfrom Maximum Probability")

---

---

3. LATENCY TO MAXIMUM PROBABILITY

df3 <- df[x<=xmax,]
max(df3$x)-min(df3$x)

[1] 0.14

p +
  geom_point(data=df3, size= 4,colour="red")+
  geom_line(data=df3, size=.6,color="red",alpha=.5)+
  geom_segment(aes(x = 0.14, y = 0, xend = 0.14, yend = 0.34),
  color="lightcoral",linetype="dashed")+
  annotate("text", x=.3,y=.32, label= "Time to Maximum \n= 0.14 Sec")+
  ggtitle("10uM IP3 - Latency to Maximum \nProbability of Open Receptors")

---

---

4. MAXIMUM PROBABILITY

max(y)

[1] 0.34

p +
  geom_point(x=x[y==max(y)],y=max(y), size= 4,colour="red")+
  annotate("text", x=x[y==max(y)]+.19,y=max(y), label= "Peak Probability = 34%")+
  geom_segment(aes(x = 0.14, y = 0, xend = 0.14, yend = 0.34),
  color="lightcoral",linetype="dashed")+
  ggtitle("10uM IP3 - Peak Probability of \nOpen Receptors over 0.98 Sec")

---

---

5. FULL WIDTH AT HALF MAXIMUM

FWHM <- x2-x1
p +
  geom_point(aes(x=x1,y=y[x==x1]), size= 4,colour="red")+
  geom_point(aes(x=x2,y=y[x==x2]), size= 4,colour="red")+
  geom_segment(aes(x = x1, y = y[x==x1], xend = x2, yend = y[x==x2]),
  color="red",linetype="dashed")+
  annotate("text", x=.39,y=.23, label= "Full Width \nHalf Maximum\n= 0.21 sec")+
  ggtitle("10uM IP3 - Full Width at \nHalf Maximum over 0.98 sec")

---

---

5. AREA UNDER THE CURVE

myroc <- roc(df$y,df$x)

'response' has more than two levels. Consider setting 'levels' explicitly or using 'multiclass.roc' insteadSetting levels: control = 0.02, case = 0.03
Setting direction: controls < cases

plot(myroc, auc.polygon=TRUE)