#Title: Estimate fitting decay curve
#Auther: Naoto Imamachi
#ver: 1.0.0
#Date: 2015-10-08
###Estimate_normalization_factor_function###
BridgeRHalfLifeCalcRaw <- function(filename = "BridgeR_4_Normalized_expression_data.txt",
group,
hour,
InforColumn = 4,
CutoffRelExp = 0.1,
CutoffDataPoint = 3,
OutputFile = "BridgeR_5_HalfLife_calculation.txt"){
###Prepare_file_infor###
time_points <- length(hour)
group_number <- length(group)
input_file <- fread(filename, header=T)
output_file <- OutputFile
###print_header###
cat("",file=output_file)
hour_label <- NULL
for(a in 1:group_number){
if(!is.null(hour_label)){
cat("\t", file=output_file, append=T)
}
hour_label <- NULL
for(x in hour){
hour_label <- append(hour_label, paste("T", x, "_", a, sep=""))
}
infor_st <- 1 + (a - 1)*(time_points + InforColumn)
infor_ed <- (InforColumn)*a + (a - 1)*time_points
infor <- colnames(input_file)[infor_st:infor_ed]
cat(infor,hour_label, sep="\t", file=output_file, append=T)
cat("\t", sep="", file=output_file, append=T)
cat("Model","R2","half_life", sep="\t", file=output_file, append=T)
}
cat("\n", sep="", file=output_file, append=T)
###calc_RNA_half_lives###
gene_number <- length(input_file[[1]])
for(x in 1:gene_number){
data <- as.vector(as.matrix(input_file[x,]))
for(a in 1:group_number){
if(a != 1){
cat("\t", sep="", file=output_file, append=T)
}
infor_st <- 1 + (a - 1)*(time_points + InforColumn)
infor_ed <- (InforColumn)*a + (a - 1)*time_points
exp_st <- infor_ed + 1
exp_ed <- infor_ed + time_points
gene_infor <- data[infor_st:infor_ed]
cat(gene_infor, sep="\t", file=output_file, append=T)
cat("\t", file=output_file, append=T)
exp <- as.numeric(data[exp_st:exp_ed])
cat(exp, sep="\t", file=output_file, append=T)
cat("\t", file=output_file, append=T)
time_point_exp <- data.frame(hour,exp)
time_point_exp <- time_point_exp[time_point_exp$exp >= CutoffRelExp, ]
data_point <- length(time_point_exp$exp)
if(!is.null(time_point_exp)){
if(data_point >= CutoffDataPoint){
model <- lm(log(time_point_exp$exp) ~ time_point_exp$hour - 1)
model_summary <- summary(model)
coef <- -model_summary$coefficients[1]
#coef_error <- model_summary$coefficients[2]
#coef_p <- model_summary$coefficients[4]
r_squared <- model_summary$r.squared
#adj_r_squared <- model_summary$adj.r.squared
#residual_standard_err <- model_summary$sigma
half_life <- log(2)/coef
if(coef < 0 || half_life >= 24){
half_life <- 24
}
cat("Raw",r_squared,half_life, sep="\t", file=output_file, append=T)
}else{
cat("Notest","NA","NA", sep="\t", file=output_file, append=T)
}
}else{
cat("Notest","NA","NA", sep="\t", file=output_file, append=T)
}
}
cat("\n", file=output_file, append=T)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.