calc_dCT.R

#Load libraries
library(data.table)
library(dplyr)
library(tidyr)
library(tibble)
library(stringr)

#Load files (spliced and retained are arg[1:2] output is arg[3])
file_input = commandArgs(trailingOnly = T)
normR = str_detect(file_input[3], "normR")
#Function to determine which is spliced vs retained
if(str_detect(file_input[1], "ss")==TRUE & str_detect(file_input[2], "iu")==TRUE){
  spliced = fread(file_input[1])
  retained = fread(file_input[2])
}

if(str_detect(file_input[1], "iu")==TRUE & str_detect(file_input[2], "ss")==TRUE){
  retained = fread(file_input[1])
  spliced = fread(file_input[2])
}

#Split Name column into the ID data
spliced = spliced %>% separate(Name, c("gene", "isoform", "time", "rep"), " ")
retained = retained %>% separate(Name, c("gene", "isoform", "time", "rep"), " ")

#Get average of spliced and retained for each timepoint so can work out which is the most stable and then do deltaCT
#mean_spliced = spliced %>% group_by(time) %>% summarize(meanCT_spliced = mean(Ct)) %>% mutate(time=as.numeric(str_remove(time, "t"))) %>% arrange(time)
#mean_retained = retained %>% group_by(time) %>% summarize(meanCT_retained = mean(Ct)) %>% mutate(time=as.numeric(str_remove(time, "t"))) %>% arrange(time)

#spliced_first_time = mean_spliced[1,2] %>% unlist() %>% unname()
#spliced_last_time = mean_spliced[nrow(mean_spliced),2] %>% unlist() %>% unname()
#dif_spliced = spliced_last_time/spliced_first_time

#retained_first_time = mean_retained[1,2] %>% unlist() %>% unname()
#retained_last_time = mean_retained[nrow(mean_retained),2] %>% unlist() %>% unname()
#dif_retained = retained_last_time/retained_first_time

#if(dif_spliced>dif_retained){
#  most_stable = spliced
#}else if(dif_retained>dif_spliced){
#  most_stable = retained
#}

if(normR){
  most_stable=retained
}else{
  most_stable=spliced
}

colnames(most_stable)[5]="most_stable_CT"
most_stable = most_stable %>% select(-c(isoform, gene, SEM))

#Do delta CT against most stable
spliced = left_join(spliced, most_stable, by=c("time", "rep")) %>% mutate(dCT = Ct-most_stable_CT)
retained = left_join(retained, most_stable, by=c("time", "rep")) %>% mutate(dCT = Ct-most_stable_CT)

#Output
output = rbind(spliced, retained)

#Write out
fwrite(output, file_input[3])