draw_stability_plots.R

#Load libraries
library(dplyr)
library(ggplot2)
library(data.table)
library(stringr)

#Take file input on command line to allow integration into pipelines (Rscript draw_graphs.R 'file.csv' 'outputs.png')
file_input = commandArgs(trailingOnly = T)

#Read as a data.table
file = fread(file_input[1])

#Force time to be numeric 
file = file %>% mutate(time=as.numeric(str_remove_all(time,"t")))

#Summarize and get mean +/- SEM
file_summary = file %>% group_by(isoform, time) %>% summarize(mean=mean(`2^-ddCT`), sd = sd(`2^-ddCT`), n=n()) %>% mutate(sem=sd/sqrt(n))

#Get gene name
gene = file[1,1] %>% unlist() %>% unname() %>% toupper()

#Replace iu with retained and ss with spliced
file_summary = file_summary %>% mutate(Isoform=ifelse(isoform=="iu", "Retained", ifelse(isoform=="ss", "Spliced", NA)))

#Plot with ggplot2
file_summary %>% 
  mutate(Gene=gene) %>%
  ggplot(aes(x=time, y=mean, group=Isoform, color=Isoform)) + 
  geom_pointrange(aes(ymin=mean-sem, ymax=mean+sem)) + geom_line(size=1.5) + 
  facet_wrap(~Gene) + 
  labs(x = "Time (hours)",
       y = expression(paste(2^-Delta,""^Delta,""^"CT"))) + 
  theme_bw(base_size=8) + 
  theme(axis.text=element_text(size=8),
        axis.title=element_text(size=8),
        legend.position = "top",
        legend.text = element_text(size=8),
        strip.background = element_rect(fill="white"),
        strip.text=element_text(size=8)) +
  scale_x_continuous(breaks=c(0,3,6,9,12)) +
  scale_color_manual(values=c("#999999", "orange")) -> saveme

#Save as GENENAME_stability_plot.png
png(file_input[2], width=3, height=2, units="in", res=1200)
saveme
dev.off()