diff --git a/GCAT_users_manual.odt b/GCAT_users_manual.odt
index ee3a033..3d67d97 100755
Binary files a/GCAT_users_manual.odt and b/GCAT_users_manual.odt differ
diff --git a/GCAT_users_manual_draft-6.docx b/GCAT_users_manual_draft-6.docx
deleted file mode 100644
index c3e6d7b..0000000
Binary files a/GCAT_users_manual_draft-6.docx and /dev/null differ
diff --git a/R/.Rhistory b/R/.Rhistory
deleted file mode 100644
index e69de29..0000000
diff --git a/R/GCAT/.RData b/R/GCAT/.RData
deleted file mode 100644
index e5b917d..0000000
Binary files a/R/GCAT/.RData and /dev/null differ
diff --git a/R/GCAT/.Rbuildignore b/R/GCAT/.Rbuildignore
index 91114bf..13e458c 100644
--- a/R/GCAT/.Rbuildignore
+++ b/R/GCAT/.Rbuildignore
@@ -1,2 +1,2 @@
-^.*\.Rproj$
-^\.Rproj\.user$
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/R/GCAT/.Rhistory b/R/GCAT/.Rhistory
deleted file mode 100644
index 8feb3e6..0000000
--- a/R/GCAT/.Rhistory
+++ /dev/null
@@ -1,159 +0,0 @@
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4_old/trunk/R/GCAT/R/addingParams.R')
-library(GCAT)
-source('~/Documents/GCAT4_old/trunk/R/GCAT/R/addingParams.R')
-t
-warnings()
-library(GCAT)
-source('~/Documents/GCAT4_old/trunk/R/GCAT/R/addingParams.R')
-warnings()
-library(GCAT)
-source('~/Documents/GCAT4_old/trunk/R/GCAT/R/addingParams.R')
-warnings()
-library(GCAT)
-source('~/Documents/GCAT4_old/trunk/R/GCAT/R/addingParams.R')
-library(GCAT)
-library(GCAT)
-library("GCAT")
-library("codetools")
-checkUsagePackage("GCAT")
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-checkUsagePackage("GCAT")
-library("codetoolds")
-library("codetoolss")
-library("codetools")
-checkUsagePackage("GCAT")
-library(GCAT)
-library("codetools")
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-checkUsagePackage("GCAT")
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/fit.model.R')
-library(GCAT)
-checkUsagePackage("GCAT")
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-checkUsagePackage("GCAT")
-source('~/Documents/GCAT4/trunk/R/GCAT/R/class.model.R')
-library(GCAT)
-checkUsagePackage("GCAT")
-library(GCAT)
-checkUsagePackage("GCAT")
-library(GCAT)
-checkUsagePackage("GCAT")
-?? checkusagePackage
-checkUsagePackage("GCAT", all + T)
-checkUsagePackage("GCAT", all = T)
-checkUsagePackage("GCAT")
-checkUsagePackage("GCAT")
-source('~/Documents/GCAT4/trunk/R/GCAT/R/plot.fit.R')
-library(GCAT)
-checkUsagePackage("GCAT")
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-source('~/Documents/GCAT4/trunk/R/GCAT/R/plot.fit.R')
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/plot.fit.R')
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/table2well.R')
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library("codetools")
-library("GCAT")
-checkUsagePackage(GCAT)
-checkUsagePackage("GCAT)
-checkUsagePackage("GCAT")
-checkUsagePackage("GCAT")
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R')
-t
-checkUsagePack(GCAT)
-checkUsagePackage(GCAT)
-checkUsagePackage("GCAT")
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/normalize.and.transform.R')
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-source('~/Documents/GCAT4/trunk/R/GCAT/R/plot.fit.R')
-library(GCAT)
-library(GCAT)
-library(GCAT)
-?? plot
-? plot
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library(GCAT)
-library("codetools")
-checkPackageUsage("GCAT")
-checkUsagePackage("GCAT")
-checkUsagePackage("GCAT", all = T)
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/build_options b/R/GCAT/.Rproj.user/7CB73FCA/build_options
deleted file mode 100644
index 5b5e33b..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/build_options
+++ /dev/null
@@ -1,4 +0,0 @@
-auto_roxygenize_for_build_and_reload="1"
-auto_roxygenize_for_build_package="1"
-auto_roxygenize_for_check="1"
-makefile_args=""
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/persistent-state b/R/GCAT/.Rproj.user/7CB73FCA/persistent-state
deleted file mode 100644
index f8c5903..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/persistent-state
+++ /dev/null
@@ -1,9 +0,0 @@
-build-last-errors="[]"
-build-last-errors-base-dir="~/Documents/GCAT4/trunk/R/GCAT/"
-build-last-outputs="[{\"output\":\"==> roxygenize('.', roclets=c('rd'))\\n\\n\",\"type\":0},{\"output\":\"* checking for changes ... \",\"type\":1},{\"output\":\"DONE\\n\",\"type\":1},{\"output\":\"\\n\",\"type\":1},{\"output\":\"==> R CMD INSTALL --no-multiarch --with-keep.source GCAT\\n\\n\",\"type\":0},{\"output\":\"* installing to library ‘/home/minh/R/x86_64-pc-linux-gnu-library/3.0’\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"* installing *source* package ‘GCAT’ ...\\n\",\"type\":1},{\"output\":\"** R\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"** inst\\n\",\"type\":1},{\"output\":\"** preparing package for lazy loading\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"\\n\",\"type\":1},{\"output\":\"Reading input files...\",\"type\":1},{\"output\":\"Error in try(gcat.load.data(file.name = file.name, input.data = input.data, : \\n\",\"type\":1},{\"output\":\" could not find function \\\"gcat.load.data\\\"\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"Creating a generic function for 'plot' from package 'graphics' in package 'GCAT'\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"No man pages found in package 'GCAT' \\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"** help\\n\",\"type\":1},{\"output\":\"*** installing help indices\\n\",\"type\":1},{\"output\":\"** building package indices\\n\",\"type\":1},{\"output\":\"** testing if installed package can be loaded\\n\",\"type\":1},{\"output\":\"\",\"type\":1},{\"output\":\"* DONE (GCAT)\\n\",\"type\":1},{\"output\":\"\",\"type\":1}]"
-compile_pdf_state="{\"errors\":[],\"output\":\"\",\"running\":false,\"tab_visible\":false,\"target_file\":\"\"}"
-console_procs="[]"
-files.monitored-path=""
-find-in-files-state="{\"handle\":\"\",\"input\":\"\",\"path\":\"\",\"regex\":false,\"results\":{\"file\":[],\"line\":[],\"lineValue\":[],\"matchOff\":[],\"matchOn\":[]},\"running\":false}"
-imageDirtyState="1"
-saveActionState="-1"
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/1B6124D6 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/1B6124D6
deleted file mode 100644
index dc615d1..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/1B6124D6
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# Estimate the growth curve slope at each timepoint of a well #\n# #\n########################################################################\n#\n# uses the functions and (see well.class.R) \n# adds estimated slopes as a new column to the \"screen.data\" slot\n\ncalculate.slopes = function(input.well, silent = T){\n # Get the growth curve data (excluding removed points, but not excluding points marked as tanking)\n\tgrowth.data = data.from(input.well, remove = T, remove.tanking = F)\n\tx = growth.data[,1]\n\ty = growth.data[,2]\n\t\n\t# Get a list of timepoint indices \n\tindices = as.numeric(rownames(growth.data))\n\t\n # Default slope is NA, values will be filled in as they are calculated\n slopes = rep(NA, length(input.well))\n\n\tfor (i in 1:length(x)){ \n\t\tif (i == 1)\n\t\t\tslopes[indices[i]] = NA\n # Calculate the slope of the line drawn from each point to the point preceding it. \n\t\telse\n\t\t\tslopes[indices[i]] = (y[i] - y[i-1])/(x[i] - x[i-1])\n\t\t}\n # Add a Slope column to the \"screen.data\" slot \n\tinput.well@screen.data$Slope = slopes\n\tif (!silent)\n\t\tcat(\"slopes filled for\", input.well@position[1], well.name(input.well), \"\\n\")\n\treturn(input.well)\n }\n \n \n########################################################################\n# #\n# Use slope estimates to check growth curves for tanking and OD jumps #\n# #\n########################################################################\n#\n# uses the functions and (see well.class.R) \n# Arguments: \n# ----- stringency parameters ----\n# remove.jumps - should the program remove OD jumps? default F (just report them) - \n# should be set to T if data contains distinct jumps in OD that need to be eliminated \n# otherwise, this might be too stringent and will result in loss of data. \n# check.start - which timepoint should checking for jumps and tanking start at? this is included because early timepoints can be unstable.\n# fall.cutoff - what downward slope should constitute a fall in OD? \n# jump.cutoffs - multipliers to determine whether a curve jumps up or down (see methods 1 and 2, below)\n# tank.limit - how many timepoints in a row can have falling slopes until the curve is marked as tanking?\n# tank.cutoff - what proportion of the maximum OD can the curve go below until it is considered tanking?\n\n# ---- input/output ----\n# silent - output to R console?\n# draw - plot growth curve with curve checking details? \n#\n# Fills the \"curve.par\" slot in the well with the starting index of tanking (NA if none is found)\n# \n\n#check.slopes = function(input.well, check.start = 8, fall.cutoff = -.0025, remove.jumps = F,\n#\t\t\tjump.multipliers = -c(15, 500, 10), tank.cutoff = 1.0, tank.limit = 3, silent = T, draw = T){\n\n#changed default values to parameters to account for settling\ncheck.slopes = function(input.well, check.start = 22, fall.cutoff = -.0025, remove.jumps = F,\n\t\t\tjump.multipliers = -c(15, 500, 10), tank.cutoff = 1.0, tank.limit = 6, silent = T, draw = T){\n\n if (!silent)\n cat(\"\\nNow checking slopes for\", plate.name(input.well), well.name(input.well))\n\n\t# Get estimated slopes and untransformed points from the well \n slopes = input.well@screen.data$Slope\n\tx = data.from(input.well, remove=F, remove.tanking=F)[,1]\n\ty = data.from(input.well, remove=F, remove.tanking=F)[,2] \n\t\n\t# Get a list of indices with valid slope estimates\n\tindices = (1:length(input.well))[!is.na(slopes)]\n\n # Do not report tanking or check for jumps if there are fewer valid points than the number needed to detect it\n\tif (length(indices) < tank.limit){\n\t\tinput.well@curve.par$tanking.start = NA\n\t\tif (!silent)\n\t\t\tcat(\"...does not have enough points to check.\\n\")\n\t\treturn(input.well)\n\t}\n\t\n #######################################################################################\n # Create indicator variables and recalculate cutoffs based on timepoint density. #\n #######################################################################################\n #\n\n # Create , a vector of indicator variables for each timepoint with a valid slope estimate. \n\tslope.indicator = rep(NA, length(slopes))\n\t\n\t# Calculate the mean time difference between two timepoints (this typically doesn't vary too much)\n time.diff = mean(diff(x[-1]))\n \n\t# Use the mean time difference to recalculate what should constitute a fall in OD using (which should be a proportion)\n\t# Honestly I don't remember why the fifth root thing is in here...this is probably going to be revised later. \n fall.cutoff = fall.cutoff * time.diff ^ (1/5) / 0.9506785\n \t\n\t# Recalculate stringency parameters for jump detection based on spread of timepoints\n\tjump.cutoffs = jump.multipliers* fall.cutoff \n\n # Recalculate tanking limit based on spread of timepoints\n\ttank.limit = round(tank.limit / time.diff ^ (1/5) * 0.9506785)\n\n\t# Cycle through the indices of input.wells with valid slope estimate\n\tcounter = 0\n\n for(j in 1:length(indices)){\n #######################################################################################\n # Method #1 for finding OD jumps: compare the slope estimate of each point to the #\n # ones for the closest surrounding points. #\n #######################################################################################\n # Get indices of the two closest surrounding timepoints with valid slope estimates. \n\t\tif (j == 1)\n\t\t\tprev.i = indices[2]\n\t\telse\n\t\t\tprev.i = indices[j-1]\n\t\tif (j == length(indices))\n\t\t\tnext.i = indices[length(indices) - 1]\n\t\telse\n\t\t\tnext.i = indices[j+1]\n\t\ti = indices[j]\n\t\t\n\t\t# How the program determines a jump up:\n\t\t# If slope estimate of current timepoint is larger than times the highest surrounding slope estimate plus \n\t\t# Add a \"2\" to the indicator variable\n\t\t# How the program determines a fall:\n\t\t# If slope estimate of current timepoint is more negative than \n\t\t# Add a \"5\" to the indicator variable\n\t\t# How the program determines a jump down:\n\t\t# If slope estimate is lower than AND is smaller than times the lowest surrounding slope estimate minus \n\t\t# Add a \"6\" to the slope indicator variable\n\t\t#\n\t\t# If none of these are true, add a \"0\" to the indicator variable\n\t\t\n\t\tif (slopes[i] > jump.cutoffs[2] * max(c(slopes[next.i],slopes[prev.i]),0) + jump.cutoffs[1])\n\t\t\tslope.indicator[i] = 2\n\t\telse if (slopes[i] < fall.cutoff){\n\t\t\tslope.indicator[i] = 5\n\t\t\tif (slopes[i] < jump.cutoffs[2] * min(c(slopes[next.i],slopes[prev.i], 0)) - jump.cutoffs[1])\n\t\t\t slope.indicator[i] = 6\t\n\t\t\t}\n\t\telse\n\t\t\tslope.indicator[i] = 0\n\n\n #######################################################################################\n # Method #2 for finding OD jumps: see if each point lies close to a line drawn #\n # between the closest surrounding points. #\n #######################################################################################\n #\n # Use variable to track the location of each point. If two subsequent points lie farther \n # away than the cutoff from their respectively drawn lines AND are on different sides, then count that as a jump. \n \n\t\tif (j > 1 & j < length(indices)){\n\t\t\n\t\t # Make equation (y=mx+b) for line drawn between two surrounding points\n\t\t\tm = (y[next.i] - y[prev.i])/(x[next.i] - x[prev.i])\n\t\t\tb = y[prev.i] - m * x[prev.i]\t\t\n\t\t\t\n # Estimate y from that line. Points will be judged by how much their true y value deviate from this estimate.\n \n # calculate b for perpendicular line from observed point to drawn line (slope is -1/m)\n b2 = y[i] + x[i]/m\n\n # solve equation for intersection to determine the shortest Euclidean distance between the point and line.\n # assign a sign to the distance based on the vertical distance. \n est.x = (b2 - b) / (m + 1/m)\n est.y = est.x * m + b\n #est.y = m * x[i] + b\n #est.x = x[i]\n \n if(m != 0)\n point.distance = sqrt((y[i]-est.y)^2 + (x[i]-est.x)^2) * sign(y[i]-est.y)\n else # horizontal case\n point.distance = y[i] - b\n \n #print(paste(i, point.distance, slopes[i], jump.cutoffs[2] * max(c(slopes[next.i],slopes[prev.i]),0) + jump.cutoffs[1], point.distance > jump.cutoffs[3])) \n \n\t\t\tcolor = \"gray30\"\n\t\t\t# If the true point exceeds that estimate by more than , update to positive.\n\t\t\t# if the counter weas previously negative, mark this as a jump up. \n\t\t\tif (point.distance > jump.cutoffs[3]){\n\t\t\t\tif (counter == -1){\n\t\t\t\t\tslope.indicator[i] = 2\n\t\t\t\t\tcolor = \"red\"\n\t\t\t\t\t}\n\t\t\t\tcounter = 1\n\t\t\t\t}\n\t\t\t# If the true point is under that estimate by more than , update to negative.\n\t\t\t# if the counter was previously positive, mark this as a jump down. \n\t\t\telse if (point.distance < -jump.cutoffs[3]){\n\t\t\t\tif (counter == 1){\n\t\t\t\t\tslope.indicator[i] = 6\n\t\t\t\t\tcolor = \"red\"\n\t\t\t\t\t}\n\t\t\t\tcounter = -1\n\t\t\t\t}\n\t\t\t# If the true point lies within of that estimate, update to zero. \n\t\t\telse \n\t\t\t\tcounter = 0\n \t\n\t\t\tif(draw)\n # Graphic representation: draw each line used in Method #2 as a dotted line, \n # and highlight in red if a jump was detected\n\t\t\t\tlines(x[c(prev.i, next.i)], y[c(prev.i, next.i)], lty = 2, col = color)\n\t\t\t}\n\t\t}\n\n #######################################################################################\n # Check for tanking by looking for unbroken series of points with falling slopes. #\n #######################################################################################\n #\n # Cycle through , adding to until the end of the curve or until reaches the \n\n\ttank = 0\n\ti = 1\n\t\n\twhile(i < length(slope.indicator) & tank < tank.limit){\n\t# If a fall was not detected, reset to 0. \n\t\tif (is.na(slope.indicator[i]))\n\t\t\ttank = 0 \n # If a fall was detected at a point index greater than , add 1 to .\n \telse if (slope.indicator[i] >= 5 & i > check.start)\n\t\t\ttank = tank + 1\n\t\telse\n\t\t\ttank = 0\n\n\t\ti = i + 1\n\t\t}\n\t\t\n\t# If the above loop was terminated because reached , update the \"curve.par\" \n # slot to denote the first point at which tanking started (should be the last index checked minus )\n # also truncate so that it does not include the timepoints after tanking. \n \n\tif (tank == tank.limit){\n\t\tinput.well@curve.par$tanking.start = i - tank.limit\n\t\tslope.indicator = slope.indicator[1:i]\n\t\tif (!silent)\n\t\t cat(\"...tanks at timepoint\", i - tank.limit, \".\\n\")\n\t\t}\n\telse{\n\t\tinput.well@curve.par$tanking.start = NA\n\t\tif (!silent)\n\t\t\tcat(\"...does not tank.\\n\")\n }\n \n #######################################################################################\n # Method #2 of checking for tanking: see if OD falls below cutoff # \n # (as a proportion of max OD) without recovery (according to ) #\n #######################################################################################\n #\n i = check.start\n tanking.start = NA\n while(i < length(y) & is.na(tanking.start)){\n # If the next ODs beyond i do not reach the cutoff, then mark i as tanking. \n if (all(y[i+(1:tank.limit)] < max(y,na.rm=T)*tank.cutoff, na.rm=T)) \n tanking.start = i\n i = i+1\n }\n \n # Graphic representation: draw the indicators used in Method #1 using the pch symbols in \n # slope index = 2: an upward-pointing traingle for an upward jump \n # slope index = 5: a diamond for a fall\n # slope index = 6: a downward-pointing triangle for downward jump \n\tif (draw){\n points(data.from(input.well, remove = F, remove.tanking = F)[which(slope.indicator != 0),], \n \tcol = 2, bg = 2, cex = 1.3, pch = slope.indicator[which(slope.indicator != 0)])\n }\n \n #######################################################################################\n # Decide what to do about any remaining jumps in OD #\n #######################################################################################\n \n\tjump.up = which(slope.indicator == 2)\n\tjump.down = which(slope.indicator == 6)\n\t# is a variable which keeps track of all the jumps, whether up or down. \n\tjump.all = sort(c(match(jump.down, indices), match(jump.up, indices)))\n # commented out; jump not working\n # if (length(jump.all) > 0)\n # add.info = paste(\"Jump(s) detected at timepoint(s)\",paste(indices[jump.all],collapse=\" \"))\n # else\n add.info = \"\"\n \n # If is true, use the following automated process to try and remove OD jumps by selectively removing points from analysis.\n # if not, just return the well with the above slot filled. \n \n if (!remove.jumps)\n input.well@add.info = add.info\n else{\n # Cycle through first few jumps (before ). is a logical that controls this loop. \n \tremove.initial.jump = T\n \twhile (length(jump.all) > 0 & jump.all[1] < check.start & remove.initial.jump){\n \t\n # If any other jumps are also before ...\n \t\tif (any(jump.all[-1] < check.start)){ \n \t\t # ...and the next jump is in a different direction, stop the loop and don't remove the first one. \n \t\t\tif(slope.indicator[indices[min(jump.all[-1])]] != slope.indicator[indices[jump.all[1]]])\n \t\t\t\tremove.initial.jump = F\n \t\t\t\t# ...or if the next two jumps are different, stop the loop and don't remove the first one. \n \t\t\telse if(length(jump.all[-1]) > 1 &\n \t\t\t\tslope.indicator[indices[jump.all[2]]] != slope.indicator[indices[jump.all[3]]])\n \t\t\t\tremove.initial.jump = F\t\n \t\t\t# ...otherwise, remove the jump and keep looping. \t\n \t\t\telse \n \t\t\t\tremove.initial.jump = T\n \t\t\t}\n # If no jumps other than the first one are before , remove it and keep looping. \n \t\telse \n \t\t\tremove.initial.jump = T\n \t\t\t\n \t # If the initial jump is to be removed, remove all points before the jump from analysis. \n # also delete the initial jump from \t\t\n \t\tif (remove.initial.jump){\n \t\t\tinput.well = remove.points(input.well, 1:(indices[jump.all[1]] - 1))\n \t\t\tinput.well@add.info = paste(add.info, \"and removed.\")\n \t\t\tjump.all = jump.all[-1]\n \t\t\t}\n \t\t}\t\n # If greater than 3 jumps remain, discard the curve as uninterpretable\t\n \tif (length(jump.all) >= 4){\n \t\tinput.well = remove.points(input.well, 1:length(input.well))\n \t\tinput.well@add.info = paste(add.info, \" - data was discarded.\")\n \t\t}\t\n \telse{\n \t# If there are 3 jumps, remove all points after the last one from analysis and delete the last jump from \n \t\tif(length(jump.all) == 3){\n \t\t\tinput.well = remove.points(input.well, indices[jump.all[3]]:length(input.well))\n \t\t\tinput.well@add.info = paste(add.info, \"and removed.\")\n \t\t\tjump.all = jump.all[-3] \n \t\t\t}\n \t\t\t\n\t\t# If there are now 2 jumps...\n \t\tif(length(jump.all) == 2){\n \t\t # ...and they are different (one up, one down), remove the points in between them from analysis.\n \t\t\tif (diff(slope.indicator[indices[jump.all]]) != 0 ){\n \t\t \t\tinput.well = remove.points(input.well, indices[jump.all[1]:(jump.all[2] - 1)])\n \t\t\t\tinput.well@add.info = paste(add.info, \"and removed.\")\n \t\t\t\t}\n \t\t\t\t# ...and they are in the same direction, remove all the points after the first one from analysis.\n \t\t\telse{\n \t\t\t\tinput.well = remove.points(input.well, indices[jump.all[1]]:length(input.well))\n\t \t\t\tinput.well@add.info = paste(add.info, \"and removed.\")\n \t\t\t\t}\t\t\t\t\n \t\t\t}\n\t\t# If there is only one jump, remove all points after it from analysis. \n \t\telse if (length(jump.all == 1)){\n \t\t\tinput.well = remove.points(input.well, indices[jump.all[1]]:length(input.well))\n \t\t\tinput.well@add.info = paste(add.info, \"and removed.\")\n \t\t\tjump.all = jump.all[-1] \n \t\t\t}\n \t\t}\n \t}\n \tif(!silent)\n \t cat(\"\\t\", input.well@add.info)\n\treturn(input.well)\n\t}\n\t\n########################################################################\n# #\n# Check wells for growth, remove from analysis if OD is too low # \n# #\n########################################################################\n#\n# The well will be tagged with no.growth = T in the slot \"curve.par\" if raw OD values (except for )\n# do not increase beyond above the specified time of inoculation for that well () \n\ncheck.growth = function(input.well, growth.cutoff, start.index = 2){\n\n # Get raw ODs (not including ) and slope estimates from the well \n # as well as OD at inoculation timepoint \n\traw.ODs = raw.data(input.well)[,2]\n start.OD = raw.ODs[start.index]\n \n raw.ODs[input.well@screen.data$Remove] = NA\n\tslope.estimates = slopes(input.well, remove.tanking = T, na.rm = T)\n\n # If fewer than 3 points remain in the analysis with valid slope estimates, discard the well. \n\tif (length(slope.estimates) < 3 | all(is.na(slope.estimates)))\n\t\tinput.well@curve.par$no.growth= T\t\t\n\telse{\n\t# If there are no points at all in the raw ODs \n\t\tif(all(is.na(raw.ODs)))\n\t\t\tinput.well@curve.par$no.growth = T\n\t\telse if(max(raw.ODs, na.rm=T) - start.OD < growth.cutoff) # See if OD increases by at least above\n\t\t\tinput.well@curve.par$no.growth = T\n\t\telse\n\t\t\tinput.well@curve.par$no.growth = F\n\t\t}\n if(all(raw.data(input.well)[,2] - raw.data(input.well)[start.index,2] < growth.cutoff)) \n input.well@add.info = \"\" # This is simply to reduce the amount of unnecessary info in the output. \n # If the well is below growth cutoff anyway, don't bother reporting other errors. \n\treturn(input.well)\n\t}\n\n\n",
- "created" : 1425413252161.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2162891425",
- "id" : "1B6124D6",
- "lastKnownWriteTime" : 1424208623,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/slope.analysis.R",
- "project_path" : "R/slope.analysis.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/26C12555 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/26C12555
deleted file mode 100644
index 1eedd3e..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/26C12555
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# Functions to calculate various things about wells based on fit model #\n# #\n########################################################################\n\n# S3 generic\nlag <- function(fitted.well, ...)\n{\n UseMethod(\"lag\")\n}\n\n#\n# Common arguments:\n# fitted.well - should be a well containing the results of , most functions will return NA if well has not been fit yet.\n# unlog - should the value be returned on the linear scale as opposed to the log-transformed scale?\n# constant.added - for returning values on the linear scale, what was the constant added before the log transform?\n# digits - passed to the function, default is no rounding (infinity digits)\n\nunlog = function(x, constant.added) {\n ########################################################################\n # Transform values back to OD scale #\n ########################################################################\nexp(x) - constant.added\n}\n\nwell.eval = function(fitted.well, Time = NULL){\n ########################################################################\n # Evaluate estimated OD at any timepoints using the fitted model #\n ########################################################################\n\n # If no timepoints are provided, use the ones collected in the experiment itself.\n\tif(!is.numeric(Time))\n\t\tTime = data.from(fitted.well)$Time\n\n # Use of equation is deprecated. Use nls and loess models stored in the well object instead\n # Attempt to use with the fitted equation and parameters to get estimates for OD at the given timepoints.\n\t#output = try(eval(fitted.well@equation, fitted.well@fit.par), silent = T)\n \n # Predict log.OD value(s) using nls model if present. If no nls model, try using loess.\n if (length(fitted.well@nls)>0) {\n output = try(predict(fitted.well@nls,list(Time=Time)),silent=T)\n } else if (length(fitted.well@loess)>0) {\n output = try(predict(fitted.well@loess,Time),silent=T)\n } else {\n output = NA\n }\n\n # Return values. If OD evaluation failed for any reason, return NULL.\n if (is.numeric(output)){\n return(output)\n } else {\n return(NULL)\n\t}\n}\n\nmodel.residuals = function(fitted.well, unlog = F){\n ########################################################################\n # Evaluate model residuals using the measured vs. fitted log.OD values #\n ########################################################################\n\tmeasured.OD = data.from(fitted.well)[,2]\n\n\t# Use with no Time argument to get fitted OD values at measured timepoints.\n\tpredicted.OD = well.eval(fitted.well)\n\n\t# If all values are valid, return the differences\n\tif (!is.numeric(predicted.OD))\n\t\treturn(NA)\n\telse\n return(measured.OD - predicted.OD)\n\t}\n\ndev.from.mean = function(fitted.well){\n ########################################################################\n # Evaluate deviations of log.OD values from the mean #\n ########################################################################\n measured.ODs = data.from(fitted.well,remove=T,na.rm=T)[,2]\n \n # Get the mean values of these measured ODs.\n mean.ODs = mean(measured.ODs)\n \n if (!is.numeric(mean.ODs))\n return (NA)\n else\n return (measured.ODs - mean.ODs)\n}\n\nrss = function(fitted.well){\n #######################################################################\n # Get the residual sum of square. #\n #######################################################################\n if (length(fitted.well@rss) == 0)\n return (NA)\n else\n return (fitted.well@rss)\n}\n\nmodel.good.fit = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate a metric for fit accuracy using squared residuals #\n ########################################################################\n\n # Sum of squared residuals\n\tRSS = rss(fitted.well)\n \n # Total sum of squared\n tot = sum(dev.from.mean(fitted.well)^2)\n \n # Coefficient of determination\n return (1 - RSS/tot)\n\t}\n\nparameter.text = function(fitted.well){\n ########################################################################\n # Output a string with values of fitted parameters #\n ########################################################################\n \n # Get a list of fitted parameters\n fit.par = fitted.well@fit.par\n \n # Giving the parameter text descriptive names.\n if (length(fitted.well@fit.par) != 0){\n names(fit.par)[1] = \"A\" \n names(fit.par)[2] = \"b\" \n names(fit.par)[3] = \"lambda\" \n names(fit.par)[4] = \"max.spec.growth.rate\" \n \n if (fitted.well@model.name == \"richards sigmoid\"){ \n names(fit.par)[5] = \"shape.par\" \n } \n \n if (fitted.well@model.name == \"richards sigmoid with linear par.\"){ \n names(fit.par)[5] = \"shape.param\" \n names(fit.par)[6] = \"linear term\"\n } \n \n if (fitted.well@model.name == \"logistic sigmoid with linear par.\")\n names(fit.par)[5] = \"linear.term\"\n \n # if loess, just show smoothing param\n if(fitted.well@model.name == \"local polynomial regression fit.\")\n fit.par = fitted.well@fit.par[\"smoothing parameter\"]\n }\n \n # Return nothing if the list is empty. Otherwise, concatenate the terms in the list with the parameter names.\n\tif(!is.list(fit.par))\n\t\treturn()\n else{\n \toutput = \"\"\n \ti = 1\n \twhile(i <= length(fit.par)){\n \t\toutput = paste(output, names(fit.par)[i], \"=\", round(as.numeric(fit.par[i]),3), \"; \", sep = \"\")\n \t\ti = i + 1\n if (i %% 6 == 0)\n output = paste(output, \"\\n\")\n \t\t}\n \toutput\n \t}\n\t}\n\nmax.spec.growth.rate = function(fitted.well, digits = Inf, ...){\n ########################################################################\n # Calculate maximum specific growth rate #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n round(fitted.well@fit.par$u,digits)\n}\n\n\nplateau = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate plateau log.OD from fitted parameters #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n plat = fitted.well@fit.par$A + fitted.well@fit.par$b\n \n\tif (!is.numeric(plat)) {\n\t plat = NA\n\t} else {\n plat = round(plat, digits)\n\t}\n\treturn(plat)\n}\n\nbaseline = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate baseline log.OD from fitted parameters #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n\n base = fitted.well@fit.par$b\n\n # If A (plateau OD) is invalid, return NA.\n\tif (!is.numeric(fitted.well@fit.par$A))\n\t\tbase = NA\n # If b (baseline OD) is invalid but plateau OD was valid, return zero.\n else if (!is.numeric(base))\n\t\tbase = 0\n\telse{\n\t\t base = round(base, digits)\n\t\t}\n\treturn(base)\n\t}\n\ninoc.log.OD = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate log.OD at inoculation from fitted parameters #\n ########################################################################\n\n # Evaluated the fitted model at the inoculation timepoint (should be zero from using from table2wells.R)\n\tif (is.null(well.eval(fitted.well)))\n\t\treturn(NA)\n else{\n inoc.time = fitted.well@screen.data$Time[fitted.well@start.index]\n inoc.log.OD = well.eval(fitted.well, inoc.time)\n if (is.na(inoc.log.OD)) inoc.log.OD = fitted.well@fit.par$b # need this in a special case: loess fits with start.index = 1 \n return(round(inoc.log.OD, digits))\n }\n\t}\n\nmax.log.OD = function(fitted.well, digits = Inf, ...){\n ########################################################################\n # Calculate max log.OD from model fit #\n ########################################################################\n\n # Evaluated the fitted model at the final timepoint (just the last valid timepoint in the experiment)\n\tif (is.null(well.eval(fitted.well)))\n\t\treturn(NA)\n else{\n \treturn(round(max(well.eval(fitted.well),na.rm=T), digits))\n }\n}\n\n\nprojected.growth = function(fitted.well,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n\tplateau(fitted.well,digits) - inoc.log.OD(fitted.well,digits)\n}\n\nprojected.growth.OD = function(fitted.well,constant.added,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n value = unlog(plateau(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added)\n round(value,digits)\n}\n\n\nachieved.growth = function(fitted.well,digits=Inf) {\n ########################################################################\n # Calculate achieved growth: max.log.OD minus the inoculated log.OD #\n ########################################################################\n max.log.OD(fitted.well,digits) - inoc.log.OD(fitted.well,digits)\n}\n\nachieved.growth.OD = function(fitted.well,constant.added,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n value = unlog(max.log.OD(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added)\n round(value,digits)\n}\n\nreach.plateau = function(fitted.well, cutoff = 0.75){\n ########################################################################\n # Did the curve come close to the plateau OD during the experiment? #\n ########################################################################\n\n plat = plateau(fitted.well)\n inoc = inoc.log.OD(fitted.well)\n final = max.log.OD(fitted.well)\n\n\tif (!is.na(final)){\n # If the plateau is the same as the OD at inoculation, return TRUE\n if ((plat - inoc) == 0)\n return(T)\n # If the difference between the final OD and inoculation OD is at least a certain proportion\n # of the difference between the plateau and inoculated ODs, return TRUE.\n else\n return((final - inoc) / (plat - inoc) > cutoff)\n\t\t}\n\telse\n\t\treturn(T)\n\t\t# If no final OD was calculated (if curve was not fit properly) just return T.\n\t}\n\n\nlag.time = function(fitted.well, digits = Inf, ...){\n ########################################################################\n # Calculate the lag time from the fitted OD #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n fitted.well@fit.par$lam\n}\n\n# new params for GCAT 4.0\n\namplitude = function(fitted.well){\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n return(fitted.well@fit.par$A)\n}\n\nshape.par = function(fitted.well){\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.par$v), NA, fitted.well@fit.par$v)\n}\n\nmax.spec.growth.rate.SE = function(fitted.well, ...){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$u), NA, fitted.well@fit.std.err$u)\n}\n\nlag.time.SE = function(fitted.well, ...){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$lam), NA, fitted.well@fit.std.err$lam)\n}\n\nshape.par.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$v), NA, fitted.well@fit.std.err$v)\n}\n\namplitude.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$A), NA, fitted.well@fit.std.err$A)\n}\n\nbaseline.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$b), NA, fitted.well@fit.std.err$b)\n}\n\n# used to calulate the inflection.time value\ninflection.time = function(well){\n if (length(well@loess) == 0 && length(well@nls) == 0) return(NA) # can' compute inflection time in the absence of a fit\n data = data.from(well)\n Time = data[,1]\n t = seq(from = min(Time), to = max(Time), by = (max(Time)-min(Time))/1000)\n y = well.eval(well,t)\n if (is.null(y)) return(NA)\n delta.t = diff(t)\n dydt = diff(y)/delta.t\n infl.index = which.max(dydt)\n t[infl.index]\n}\n",
- "created" : 1425413281936.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2757150562",
- "id" : "26C12555",
- "lastKnownWriteTime" : 1428437531,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/fitted.calculations.R",
- "project_path" : "R/fitted.calculations.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2C8FAF10 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2C8FAF10
deleted file mode 100644
index 2aeaf02..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2C8FAF10
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\nrequire(pheatmap)\nrequire(gplots)\n\n########################################################################\n# #\n# Graphic output functions for fitted well objects. The functions are #\n# fairly complicated and intertwined and may need revision. #\n# #\n########################################################################\n\n# S3 Generic.\nplot <- function(input.well, ...) {\n UseMethod(\"plot\")\n}\n\n########################################################################\n# Basic function plots time vs. OD from a well object #\n########################################################################\n#' plot.data\n#'\n#' Basic function plots time vs. OD from a well object\n#'\n#' @param input.well The well object that need to be plottedd\n#' @param unlog should data be plotted on a linear (vs. logarithmic) scale?\n#' @param view.raw.data should the raw data be plotted? (\n#' @param number.points should points be labeled with numeric indices?\n#' @param scale determines the font scale for the entire graph. all cex values are calculated from this.\n#' @param draw.symbols - should be called on the well and markings drawn on the graph?\n#' @param ... additional arguments passed to plot()\nplot.data = function(input.well, view.raw.data = F, unlog = F, scale = 1, \n main = paste(plate.name(input.well), well.name(input.well)), number.points = T, \n draw.symbols = F, constant.added, ylim, ...){\n \n # Get data as well as a vector showing which points were removed. \n\tinput.data = data.from(input.well, remove = F, remove.tanking = F, raw.data=view.raw.data)\n\tremoved.points = !(rownames(input.data) %in% rownames(data.from(input.well, remove = T, remove.tanking = T)))\n point.colors = as.character(factor(removed.points,levels=c(F,T),labels=c(\"black\",\"gray80\")))\n\n # Draw the axes and all text labels first.\n\tpar(mar = c(5, 4, 4, 5)+0.1)\n\tplot(input.data, main = main, xlab = \"Time(hours)\", ylab = \"log(OD - blank + const)\",\n mex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, cex.lab = 1.2*scale, type =\"n\",...)\n \n\t# Draw a second vertical axis, showing unlogged OD scale\n\t# - Determine the range of the labels: from min.OD to max.OD\n if (class(try(ylim,silent=T)) == \"try-error\") {\n OD = unlog(input.data[,2],constant.added)\n baseline.OD = unlog(baseline(input.well),constant.added)\n min.OD = min(min(OD,na.rm=T),baseline.OD,na.rm=T)\n plateau.OD = unlog(plateau(input.well),constant.added)\n max.OD = max(max(OD,na.rm=T),plateau.OD,na.rm=T)\n } else {\n min.OD = unlog(ylim[1],constant.added)\n max.OD = unlog(ylim[2],constant.added)\n }\n\t# - Compute labels and their positions\n\tOD.labels = seq(from = min.OD, to = max.OD, length.out = 5)\n\tOD.labels = signif(OD.labels,2)\n\tOD.at = log(OD.labels+constant.added)\n # - Draw the axis\n\taxis(side=4, at=OD.at, labels=OD.labels, cex.axis = 1.2*scale, cex.lab = 1.2*scale)\n\tmtext(4, text = \"OD - blank\", line = 3, cex=1.2)\n \n\t# If is true, then label each point with the index of its timepoint and plot removed points in grey, others in black. \n if (number.points)\n\t\ttext(input.data$Time, input.data[,2], rownames(input.data), col = point.colors, cex = 0.5*scale)\n\t# Otherwise plot all points, using a different plotting character for removed points. \t\n else\n\t\tpoints(input.data$Time, input.data[,2], pch = 1 + removed.points*15)\n\n # If is set to T, then draw all the markings that makes to determine curve parameters. \n\tif (draw.symbols & !view.raw.data)\n\t\tcheck.slopes(input.well, draw = T)\n\treturn()\n\t}\n\t\n########################################################################\n# Plots the fitted model curve from a well object if it exists #\n########################################################################\n#\n# time: specify which points (in units of time) to plot fitted OD values for. if not specifies, plot all timepoints in range of well. \n\nplot.model = function(input.well, col = 1, scale = 1, lty = 1, time = NULL, unlog = F, constant.added=1, ...){\n\t \n #input.data = data.from(input.well)\n\t#growth = input.data[,2]\n\n # If no list of timepoints is specified, get a list of 360 timepoints (should be smooth enough) from the well's range. \n\tif (is.null(time)){\n\t\ttime.fin = max(data.from(input.well, raw.data = T, remove = F, remove.tanking = F)$Time)\n\t\ttime = seq(0, time.fin, length.out = 360)\n\t}\n \n # Evaluate the predicted OD at the specified timepoints based on the fitted model. \n\tpredicted.OD = well.eval(input.well, time) \n # If any values were returned, plot the lines on the current graph. Otherwise, just return without doing anything.\n\tif (is.numeric(predicted.OD))\n\t\tlines(time, predicted.OD, col = col, lty = lty, lw = 2 * scale) \n else\n return()\t\n \n}\n\t\n########################################################################\n# Put various parameters and info in text form on the graphs #\n########################################################################\n#\t\t\t\ndraw.text = function(input.well, scale = 0.5, xlim = 0, ylim = 0,...){\n\n\t#input.data = data.from(input.well, remove = F, remove.tanking = F)\n\t#fit = input.well@fit.par\n \n # - fit information (fit status, model name if available, jump detection output, fit parameters if available) from well \n # color = red if no fit, blue if fit, green if skipped\n # - empty or inoculated well.\n # color = green if empty, blue if inoculated, red if inoculated but has no growth or empty but has growth. \n \n \tcol2 = \"blue\" \n\ttext2 = paste(input.well@fit.info, input.well@model.name, input.well@add.info, \"\\n\", parameter.text(input.well))\n\n\tif (length(input.well@fit.par) == 0) # no fit\n\t\tcol2 = \"red\"\n\n\tif (is.empty(input.well)){\n\t\t text1 = \"empty well\"\n\t\tif(!lacks.growth(input.well) | length(input.well@fit.par) == 0) # growth curve fit for an empty well\n\t\t\tcol1 = \"red\"\n\t\telse\n\t\t\tcol1 = \"forestgreen\"\n\t\tif (length(input.well@model.name) == 0) # well was skipped \n\t\t\tcol1 = col2 = \"forestgreen\"\n\t\t}\n\telse{\n\t\ttext1 = \"inoculated well\"\n\t\tif(lacks.growth(input.well) | length(input.well@fit.par) == 0) # failure to fit an inoculated well\n\t\t\tcol1 = \"red\"\n\t\telse\n\t\t\tcol1 = \"forestgreen\"\n\t\t}\n\t\t\n\t# - goodness of fit metric. \n\t# color = red if below 2, yellow if between 2 and 2.72, and green if above 2.72. \n\t\n\tif(!is.na(model.good.fit(input.well))){\n\t\t#if (model.good.fit(input.well, unlog = F) > 2.72)\n\t\t#\tcol1.5 = \"forestgreen\"\n\t\t#else if (model.good.fit(input.well, unlog = F)> 2.0)\n\t\t#\tcol1.5 = \"gold2\"\n\t\t#else\n\t\t#\tcol1.5 = \"red\"\n col1.5 = \"forestgreen\"\n\t\ttext1.5 = paste(\"R squared:\", round(model.good.fit(input.well),3))\n\t\t}\n\telse\n\t col1.5 = text1.5 = NULL\n\t \n # Print all text at the top of the graph with approprate positions and scaling \t\n \ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0.025 * diff(ylim), \n \t\ttext1.5, cex = 1.5*scale, col = col1.5) \n \t\t\n\ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0 * diff(ylim), \n\t\t\ttext1, cex = 1.5*scale, col = col1)\n\n\ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0.03 * diff(ylim), \n\t\t\ttext2, pos = 1, cex = 1.5*scale, col = col2)\t\n\t}\t\n\n########################################################################\n# Draw lines on graph denoting calculated parameters #\n########################################################################\n#\t\n# - should curve parameters be labeled? \n\t\t\ndraw.calc.par = function(input.well, scale = 0.5, unlog = F, constant.added, show.num = T){\n\n # Don't do anything if well was not fit. \n\tif (is.null(well.eval(input.well)))\n\t\treturn()\n \n # Collect values for various curve parameters. \n\tbaseline = baseline(input.well)\n\tinoc.log.OD = inoc.log.OD(input.well)\n\tmax.log.OD = max.log.OD(input.well) \n\tplateau = plateau(input.well) \n\tinflection.time = input.well@inflection.time # was a param in model\n\tfin.time = (inflection.time+max(data.from(input.well)[,1]))/2\n\n\t# = timepoint at greatest growth \n\t# = OD measurement at , minus the constant added before taking the log (if reversing the transformation)\n\t# = slope (on log scale) at (specific growth)\n # had to add the unlog code. was calculated differently before NWD 7/21/14\n\tmax.slope = max.spec.growth.rate(input.well)\n\tmax.y = well.eval(input.well, inflection.time)\n\tlag.x = lag.time(input.well) \n lag.y = baseline\n\t\n\t# ---- Specific growth rate ---- #\n\tlines(c(lag.x, inflection.time), c(lag.y, max.y), lty = 2, col = \"red\")\n\n\n # Blue dotted line at time of maximum growth, with text label for specific growth rate. \n\tabline(v = inflection.time, lty = 2, lw = (scale^2)*2, col = \"blue\")\n if(show.num) text(inflection.time, max.y, round(max.slope,3), col = \"blue\", cex = 1.5*scale, pos = 2)\n\n # inoculation OD and baseline of the fitted model\n abline(h = inoc.log.OD, lw = scale*2, lty = 3)\n abline(h = baseline, col = \"red\", lw = (scale^2)*2, lty = 2)\n if(show.num) {\n text(fin.time, inoc.log.OD, paste(round(inoc.log.OD,3),\"\\n\",sep=\"\") , col = \"black\", cex = 1.5*scale, pos = 2)\n text(fin.time, baseline, paste(\"\\n\\n\", round(baseline,3), sep=\"\") , col = \"red\", cex = 1.5*scale, pos = 2)\n }\n\n # ---- Lag time ---- #\n # Do not draw a horizontal line to lag time if it is 0 or negative. \n # Otherwise draw a red line from the starting point to the lag time, and label with the lag time \n \tif (lag.time(input.well) == 0){\n\t\tif(show.num) text(0, inoc.log.OD, \"\\n\\n0.000\", col = \"red\", cex = 1.5*scale, pos = 4)\n\t\t}\n\telse{\n\t\tlines(c(0, lag.x), c(baseline, baseline), col = \"red\", lw = (scale^2)*2, lty = 2)\n\t\tif(show.num) text(lag.x, lag.y, paste(\"\\n\\n\", round(lag.time(input.well),3)), col = \"red\", cex = 1.5*scale, pos = 2)\n\t}\n\n # ---- Total growth ---- #\n\n # Draw horizontal lines for the max.log.OD in black, the plateau in green and the initial OD in black.\n abline(h = max.log.OD, lty = 3, lw = scale*2)\n abline(h = plateau, lty = 2, lw = (scale^2)*2, col = \"forestgreen\")\n\t\n\t# Draw a vertical line from the initial OD to the final OD in black, and then to the plateau in gray. \n\tlines(c(fin.time, fin.time), c(inoc.log.OD, max.log.OD), lw = (scale^2)*2, lty = 3)\n\tlines(c(fin.time, fin.time), c(max.log.OD, plateau), lw = (scale^2)*2, lty = 3, col = \"grey\")\n\n # Text: plateau and initial ODs (on left), difference between initial and final OD on right\n\tif(show.num){\n text(fin.time, plateau, paste(round(plateau,3),\"\\n\",sep=\"\") , col = \"forestgreen\", cex = 1.5*scale, pos = 2)\n text(fin.time, max.log.OD, paste(\"\\n\\n\\n\",round(max.log.OD,3),sep=\"\") , col = \"black\", cex = 1.5*scale, pos = 2)\n text(fin.time, .5*(max.log.OD-inoc.log.OD)+inoc.log.OD, round(max.log.OD - inoc.log.OD,3), cex = 1.5*scale, pos = 4)\n # difference between final and plateau OD (if large enough) \n if (!reach.plateau(input.well))\n\t\t text(fin.time, .5*(plateau-max.log.OD)+max.log.OD, paste(\"(\", round(plateau - max.log.OD,3), \")\", sep = \"\"), col = \"grey\", cex = 1.5*scale, pos = 2) \n }\n }\n\t\n########################################################################\n# Draw residuals from the nonlinear fit with option for lowess line #\n########################################################################\n#\t\t\nplot.residuals = function(input.well, xlim = NULL, lowess = T, ...){\n well = input.well\n\tdata = data.from(well, remove = F, remove.tanking = F)\n\n\tif (is.null(xlim))\n\t\txlim = c(min(data$Time, 0)-1, max(data$Time))\n\t\t\n\tplot(data.from(well)[,1], model.residuals(well), main = paste(plate.name(well), well.name(well), \"\\n[Residuals]\"),\n\t\txlab = \"Time(hours)\", ylab = paste(\"Residual\", names(data)[2]), xlim = xlim)\n\t\n\tabline(0,0, lty = 2)\n\t\n\tif (lowess)\n\t\tlines(lowess(data.from(well)[,1], model.residuals(well)), lw = 2, col = \"red\")\n\t}\n\n##############################################################################\n# This function is used to create a heatmap using: \n# specific growth, total growth, and lag time\n# for each well on a plate.\n#\n# @params\n# fitted.well.array: matrix containing well array object data\n# attribute: the data type we should use to create a heatmap\n# @returns\n# path of heatmap pdf file\n##############################################################################\ncreate.heatmap = function(fitted.well.array, attribute, unlog=NULL){\n attr.name <- deparse(substitute(attribute))\n pdf.name <- \"\"\n if(class(fitted.well.array) == \"matrix\"){\n #We may want to sub() out periods from plate.ID if it causes problems\n plate.ID = unique(unlist(aapply(fitted.well.array,plate.name)))[1]\n if(is.null(unlog)) {\n spec.growth = unlist(aapply(fitted.well.array, attribute))\n }\n # currently only total growth needs to be unlogged if unlog == T\n else {\n spec.growth = unlist(aapply(fitted.well.array, attribute))\n }\n num.dig = 3 #how many digits should be put on pdf?\n max = round(max(spec.growth, na.rm=T), digits=num.dig)\n min = round(min(spec.growth, na.rm=T), digits=num.dig)\n avg = round(mean(spec.growth, na.rm=T), digits=num.dig)\n heat.text = paste(toupper(sub(\"\\\\.\", \" \", attr.name)), \":\\n\", plate.ID, \"\\n\",\n paste(\"Max:\",max ,\"Min:\" ,min ,\"Avg:\", avg, sep=\"\"))\n \n attr.name <- sub(\"\\\\.\", \"_\", attr.name) #do not want periods in file path\n letters <- attr(fitted.well.array, \"dimnames\")[[1]]\n for(i in 1:length(letters)) letters[i] = paste(\" \", letters[i], \" \")\n nums <- attr(fitted.well.array, \"dimnames\")[[2]]\n for(i in 1:length(nums)) nums[i] = paste(\" \", nums[i], \" \")\n heat <- matrix(spec.growth, nrow=dim(fitted.well.array)[1], ncol=dim(fitted.well.array)[2], dimnames=list(letters,nums))\n pdf.name <- paste(getwd(), \"/\", plate.ID, \"_\", attr.name, \".pdf\", sep=\"\")\n pdf(pdf.name)\n #heatmap(heat, Rowv=NA, Colv=NA, revC=T, scale=\"none\", na.rm=T, main=plate.ID, col=rainbow(100), margins=c(6,6))\n #mtext(paste(\"Max:\", round(max(spec.growth, na.rm=T), digits=4),\"Min:\", round(min(spec.growth, na.rm=T), digits=4), \"Avg:\", round(mean(spec.growth, na.rm=T), digits=4)), side=1, line=3)\n pheatmap(heat, color=colorpanel(100, \"red\", \"orange\", \"yellow\"),\n border_color=\"black\", cell_width=2, cell_height=3,\n cluster_rows=F, cluster_cols=F, scale='none', main=heat.text, fontsize=16)\n dev.off()\n }\n else {\n return(\"Error\") \n }\n return(pdf.name)\n}\n\n########################################################################\n# Draw grids of 96 points as a visual representation of fit status, #\n# and other info for an array of fitted well objects, plate by plate #\n########################################################################\n#\n\t\nplate.overview = function(fitted.well.array, scale = 1, plate.ncol = 12, plate.nrow = 8){\n \n \n # Start with a list of the unique plate names in the fitted well array \n # and an appropriately-sized grid of coordinates to plot wells on.\n\tplates = unique(unlist(aapply(fitted.well.array,plate.name)))\n \n\tgrid = data.frame(x = rep(rep(1:plate.ncol, each = plate.nrow), length(plates)),\n\t\t\t y = rep( rep(-(1:plate.nrow), times = plate.ncol), length(plates)))\n\n\n # Gather information on each well to display on each of the coordinates in :\n # - was it marked as empty in the plate layout?\n # - did the program find it to contain no growth (\"dead\")? \n # - was the fitting procedure successful? \n # - did the curve tank? if so, at what timepoint? if not, or if the curve was marked as dead anyway, do not display the value. \n # - does the \"additional info\" slot indicate that any points were removed or the whole well discarded?\n \n\tempty = unlist(aapply(fitted.well.array, is.empty))\n\tdead = unlist(aapply(fitted.well.array, lacks.growth))\n\tfit = unlist(aapply(fitted.well.array, contains.fit))\n \n\ttanking = unlist(aapply(fitted.well.array, tanking.start))\n\ttanking[is.na(tanking) | tanking == 1 | dead] = \"\"\n\n\terrors = unlist(aapply(fitted.well.array, function(well){\n\t\tif (length(well@add.info) == 0)\n\t\t\t\"\"\n\t\telse if (grepl(\"removed\", well@add.info))\n\t\t\t\"-\"\n else if (grepl(\"detected\", well@add.info))\n \"+\"\n\t\telse if (grepl(\"discarded\", well@add.info))\n\t\t\t\"!\"\n\t\telse\n\t\t\t\"\"\n\t\t}))\n\t\n # Color and plotting character vectors (length = the number of wells in the array)\n # Default = 1 (open point, black)\n\tcolors = char = rep(1, length(tanking))\n\n # Desired colors\n colors[empty & dead] = \"green3\" # Empty well with no growth.\n colors[!empty & fit] = \"blue\" # Inoculated well with successfully fitted growth curve.\n \n # Undesired colors \n colors[empty & !dead] = \"darkolivegreen4\" # Inoculated well with some growth. \n colors[!empty & !fit] = \"red\" # Inoculated well with no successfully fit (either no growth or unsuccessful fit).\n \n char[!dead & fit] = 19 # Filled points for non-empty wells with successful fits \n char[!dead & !fit] = 4 # an X for non-empty wells with failed fits. \n \n char[errors == \"!\"] = 8 # Asterisk for discarded wells. \n char[errors == \"-\" & dead ] = 5 # Open diamond for empty wells (after removing points).\n char[errors == \"-\" & !dead & fit] = 23 # Filled diamond for non-empty wells with removed points and successful fits. \n char[errors == \"-\" & !dead & !fit] = 8 # Asterisk for wells with removed points and failed fits.\n \n \n\tfor (plate in 1:length(plates)){\n \n\t\tindices = (plate - 1) * plate.nrow*plate.ncol + 1:(plate.nrow*plate.ncol)\n\n # Plot the grid using colors and plotting characters determined above. \n\t\tplot(grid[indices,], col = colors[indices], bg = colors[indices], pch = char[indices], \n\t\t\t main = plates[plate], mex = scale, cex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, \n\t\t\t xaxt = \"n\", yaxt = \"n\", xlim = c(-0.5,plate.ncol + 0.5), ylim = c(-(plate.nrow + 1.5), 0.5), xlab = \"\", ylab = \"\")\n \n # Symbol legends\n \n legend.xpos = (c(-1,2.75,6.5,6.86,10.25)+0.5)*(plate.ncol+1)/13 - 0.5\n legend.ypos = -(plate.nrow + 0.5)\n \n legend(x=legend.xpos[1], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Empty, no growth\",\"Empty with growth\"),\n pch = c(1,19),\n pt.bg = c(\"green3\",\"darkolivegreen4\"),\n col = c(\"green3\",\"darkolivegreen4\")\n )\n legend(x=legend.xpos[2], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Inoculated with growth\", \"Inoculated, no growth\"),\n pch = c(19,1),\n pt.bg = c(\"blue\",\"red\"),\n col = c(\"blue\",\"red\")\n )\n legend(x=legend.xpos[3], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Well tanks at specified index\", \"Some points removed\"),\n pch = c(21,23),\n pt.bg = c(\"grey\",\"grey\"),\n col = c(\"black\",\"black\")\n ) \n \n text(x=legend.xpos[4], y=legend.ypos - 0.29,\"#\",cex=0.5*scale) \n \n legend(x=legend.xpos[5], y=legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Model fitting failed\", \"Well discarded\"),\n pch = c(4,8),\n pt.bg = c(\"black\",\"black\"),\n col = c(\"black\",\"black\")\n )\n \n # Add tanking indices if any were found. \n \ttext(grid[indices,] + 0.30, cex = 0.75*scale, \n labels = tanking[indices], col = colors[indices])\n\n # Label rows and columns\n\t\ttext(-1, -1:-plate.nrow, pos = 4, LETTERS[1:plate.nrow], cex = scale) \n\t\ttext( 1:plate.ncol, 0 , 1:plate.ncol, cex = scale) \t\n\t\t}\n \n\t}\n\n########################################################################\n# Draw each well in an array of fitted well objects in succession. #\n# Include options for adding notations, text info and fit parameters. #\n########################################################################\n#\n\nview.fit = function(fitted.data, indices = 1:length(fitted.data), \n unlog = F, constant.added, xlim = NULL, ylim = NULL, display.legend = T, \n\t\t show.text = T, show.calc = T, draw.guess = NULL, draw.symbols = F, number.points = T, \n\t\t\tuser.advance = T, show.residuals = F, scale = 1,...){\n\n if(!is.array(fitted.data))\n fitted.data = list(fitted.data)\n\n # Determine the boundaries for the axes (if user did not specify them)\n if(is.null(ylim)){\n min.y = min(unlist(aapply(fitted.data, function(well){\n \tif (unlog) well@use.log = F\n \tmin.y = min(data.from(well, remove = F, remove.tanking = F)[,2], na.rm = T)\n min(min.y, well@fit.par$b)\n \t})))\n max.y = max(unlist(aapply(fitted.data, function(well){\n if (unlog) well@use.log = F\n \tmax.y = max(data.from(well, remove = F, remove.tanking = F)[,2], na.rm = T)\n max(max.y, well@fit.par$b + well@fit.par$A)\n \t})))\n ylim = c(min.y, min.y + (max.y-min.y)*1.15) - unlog*constant.added\n }\n if(is.null(xlim)){\n min.x = min(unlist(aapply(fitted.data, function(well){\n min(data.from(well, remove = F, remove.tanking = F)[,1], na.rm = T)\n \t})))\n max.x = max(unlist(aapply(fitted.data, function(well){\n \tmax(data.from(well, remove = F, remove.tanking = F)[,1], na.rm = T)\n \t})))\n xlim = c(min.x - 0.05 * (max.x-min.x), max.x)\n }\n \n # Display a figure legend\n if(display.legend){\n well.fit.legend(xlim=xlim,ylim=ylim,scale=scale,constant.added=constant.added)\n if(user.advance){\n prompt = readline(\" to continue or Q to quit >>\")\n if (toupper(prompt) == \"Q\") break\n }\n }\n # Start to cycle through the wells \n\twell.number = 1\n\twhile (well.number <= length(fitted.data)) {\t\t\n\t\t# Only show wells specified by (default all wells)\n if (well.number %in% indices){ \n # plot the well\n fitted.well = fitted.data[[well.number]]\n plot(x=fitted.well, constant.added = constant.added, xlim = xlim, ylim = ylim,\n unlog = unlog, well.number = well.number, scale = scale, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, ...)\n \n if(user.advance)\n cat(\"\\n[\", well.number, \"] \", plate.name(fitted.well), \" \", well.name(fitted.well), \".\", sep = \"\")\n \n if (show.residuals & is.numeric(model.residuals(fitted.well))){\n if(user.advance)\n if (toupper(readline(\" for residuals >>\")) == \"Q\") break\n plot.residuals(fitted.well)\n }\n \n # Allow user to advance the currently shown well if specified. \n\t\t\tif (user.advance){\n \n\t\t\t\tprompt = readline(\" to continue, or type # of next well or Q to quit >>\")\n\t\t\t\tif (toupper(prompt) == \"Q\") break\n\n user.input = suppressWarnings(try(as.numeric(prompt),silent=T))\n \n # Go onto the next well unless input is a number. \n\t\t\t\tif (is.numeric(user.input) & !is.na(user.input) & length(user.input) > 0)\n\t\t\t\t\twell.number = user.input - 1\n\t\t\t }\n\t\t\t}\n # Advance the loop\n well.number = well.number + 1\n\t\t}\t\t\n\t}\t\n\n\nwell.fit.legend = function(xlim, ylim, scale = 1, constant.added){\n par(mar = c(5, 4, 4, 5)+0.1)\n plot(0,0, main = \"[Index] \\n; \",\n xlim = xlim, ylim = ylim, xlab = \"Time\", ylab = \"log(OD - blank + const)\", \n mex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, cex.lab = 1.2*scale, type = \"n\")\n \n # Draw a second vertical axis, showing unlogged OD scale\n min.OD = unlog(ylim[1],constant.added)\n max.OD = unlog(ylim[2],constant.added)\n OD.labels = seq(from = min.OD, to = max.OD, length.out = 5)\n OD.labels = round(OD.labels,1)\n OD.at = log(OD.labels+constant.added)\n axis(side=4, at=OD.at, labels=OD.labels, cex.axis = 1.2*scale, cex.lab = 1.2*scale)\n mtext(4, text = \"OD - blank\", line = 3, cex=1.2)\n \n # Sample max. slope line\n abline(v=min(xlim)+0.5*max(xlim), col=\"blue\", lty=2)\n text(mean(xlim),min(ylim)+0.4*diff(ylim),labels=\"Maximum specific\\ngrowth rate\",col=\"blue\",pos=2,cex=0.75*scale)\n \n # Sample plateau line\n abline(h=min(ylim)+0.8*diff(ylim),col=\"forestgreen\",lty=2)\n text(min(xlim)+0.9*diff(xlim),ylim+0.8*diff(ylim),labels=\"Growth plateau\",col=\"forestgreen\",pos=3,cex=0.75*scale)\n\n # Sample max.log.OD line\n abline(h=min(ylim)+0.7*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.9*diff(xlim),ylim+0.7*diff(ylim),labels=\"max.log.OD\",col=\"black\",pos=1,cex=0.75*scale)\n \n # Sample inoc.log.OD\n abline(h=min(ylim)+0.1*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.1*diff(xlim),min(ylim)+0.1*diff(ylim),labels=\"Fitted growth\\nat inoculation\",col=\"black\",pos=3,cex=0.75*scale)\n \n # Sample baseline\n abline(h=min(ylim)+0.05*diff(ylim),col=\"red\",lty=2)\n text(min(xlim)+0.1*diff(xlim),min(ylim)+0.05*diff(ylim),labels=\"Baseline\",col=\"red\",pos=1,cex=0.75*scale)\n\n # Sample lag time\n lines(min(xlim)+c(0.1,0.25,0.50)*max(xlim),min(ylim)+c(0.05,0.05,0.4)*diff(ylim),col=\"red\",lty=2)\n text(min(xlim)+0.25*max(xlim),min(ylim)+0.05*diff(ylim),labels=\"Lag time\",col=\"red\",pos=1,cex=0.75*scale)\n \n # Sample achieved growth\n lines(min(xlim)+c(0.75,0.75)*max(xlim),min(ylim)+c(0.1,0.7)*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.75*max(xlim),min(ylim)+0.3*diff(ylim),labels=\"Achieved growth\",col=\"black\",cex=0.75*scale)\n \n # Sample plateau - achieved growth\n lines(min(xlim)+c(0.75,0.75)*max(xlim),min(ylim)+c(0.7,0.8)*diff(ylim),col=\"grey\",lty=3)\n text(min(xlim)+0.75*max(xlim),min(ylim)+0.75*diff(ylim),labels=\"Projected minus achieved growth\",col=\"grey\",cex=0.75*scale)\n \n # Symbol legend\n legend(x=\"right\", title = \"Timepoint Symbols\", legend = c(\"Normal point\", \"Ignored point\"),\n cex = 0.75*scale, pt.cex = c(0.6,0.6)*scale, pch = c(35,35), col=c(\"black\",\"gray80\"),\n x.intersp=1, xjust = 1, y.intersp=1.5)\n}\n\npdf.by.plate = function(fitted.data, out.prefix = \"\", upload.timestamp = NULL, \n out.dir = getwd(), unlog = F, constant.added, silent = T, overview.jpgs = T, plate.ncol = 12, plate.nrow = 8,...){\n \n # Prepare timestamp for addition to output file names. \n filename.timestamp = strftime(upload.timestamp, format=\"_%Y-%m-%d_%H.%M.%S\")\n \n # Start file list with the overview pdf\n file.list.out = paste(out.dir,\"/\",out.prefix, \"_overview\", filename.timestamp, \".pdf\",sep=\"\")\n \n # Call to draw a graphic representation of each plate in this file. \n pdf(file.list.out, title = paste(out.prefix, \"plate overview\"))\n plate.overview.out = try(plate.overview(fitted.data),silent=T)\n if(class(plate.overview.out) == \"try-error\")\n stop(\"Error in : \", plate.overview.out)\n \n # Close devices\n while(dev.cur() != 1)\n dev.off() \n \n # Cycle through each plate \n for(i in 1:dim(fitted.data)[3]){\n \n # Get plate ID and position in data array.\n plate.ID = dimnames(fitted.data)[[3]][i]\n plate.indices = (i-1) * plate.nrow*plate.ncol + 1:(plate.nrow*plate.ncol)\n if(overview.jpgs){\n # most be > 1 to partition value breaks for heatmap\n well.matrix <- aapply(fitted.data[,,i], max.spec.growth.rate) \n num.wells <- length(well.matrix[!sapply(well.matrix, is.na)])\n if(num.wells > 1){\n #Heatmap block##########################################################\n #alongside the jpgs file create 3 heatmaps for each plate. NWD\n spec.heat.file = create.heatmap(fitted.data[,,i], max.spec.growth.rate)\n if(spec.heat.file == \"Error\")\n stop(\"Error in for specific growth\")\n lag.heat.file = create.heatmap(fitted.data[,,i], lag.time)\n if(lag.heat.file == \"Error\")\n stop(\"Error in for lag time\")\n total.heat.file = create.heatmap(fitted.data[,,i], achieved.growth)\n if(total.heat.file == \"Error\")\n stop(\"Error in for total growth\")\n # Add name of file if successfully written to file list output. Including heatmap files NWD\n file.list.out = c(file.list.out, spec.heat.file, lag.heat.file, total.heat.file)\n ########################################################################\n }\n jpg.name = paste(out.dir, \"/\", plate.ID, \"_overview\", \".jpg\", sep=\"\")\n jpeg(jpg.name, quality = 90, width = 600, height = 480)\n plate.overview.out = try(plate.overview(fitted.data[,,i]),silent = T)\n if(class(plate.overview.out) == \"try-error\")\n stop(\"Error in : \", plate.overview.out)\n }\n else\n jpg.name = c()\n \n # Open a separate PDF for each plate.\n if(!silent) cat(\"\\nprinting PDF for\", plate.ID)\n pdf.name = paste(out.dir, \"/\", plate.ID, \"_plots\", filename.timestamp, \".pdf\", sep=\"\")\n pdf(pdf.name, title = paste(\"R Graphics output for plate\", plate.ID))\n \n # Call to draw each well on the plate to the pdf. \n view.fit.out = try(view.fit(fitted.data, indices = plate.indices, unlog=unlog, constant.added=constant.added, user.advance=F,...),silent=T) \n \n if(class(view.fit.out) == \"try-error\")\n stop(\"Error in : \", view.fit.out)\n\n # Close all devices\n while(dev.cur() != 1)\n dev.off() \n \n if(!silent) cat(\"...done!\\n\\twritten to\", pdf.name, \"\\n\") \n file.list.out = c(file.list.out, jpg.name , pdf.name)\n }\n return(file.list.out)\n } \n",
- "created" : 1425413258441.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "1576748553",
- "id" : "2C8FAF10",
- "lastKnownWriteTime" : 1428438949,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/plot.fit.R",
- "project_path" : "R/plot.fit.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2FF709FD b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2FF709FD
deleted file mode 100644
index 577cdb9..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/2FF709FD
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n# Notes by Jason\n# 9/07/11\n\n\n########################################################################\n# #\n# Function for loading data from tabular format into an object array #\n# #\n########################################################################\n#' Load tabular data\n#'\n#' This function handles loading data from tabular format (.csv, tab-delimited text or R data frame object)\n#' and returns an array of well objects, each filled with raw Time vs. OD data. \n#' It takes single-plate or multiple-plate format data. For single-plate data, \n#'it calls on the function \\code{gcat.reorganize.single.plate.data} to rearrange the table before creating the output object. \n#'\n#' @param file.name Complete path and file name of a comma-separated values (.csv) file containing growth curve data \n#' in the multiple-plate (long) format. \n#' @param input.data A list of tables representing input files read with \\code{read.table}. Used to save time in cases\n#' of running multiple analyses on the same dataset. If used, the function will ignore \\code{file.name} entirely.\n#' @param load.type .csv by default. \n#' @param plate.laout Specifies the layout of the given plate.\n#' @param single.plate.ID specifies a plate name for a single-plate read. If NULL, this is derived from the file name. \n#' @param blank.value Blank OD measurement for uninoculated wells. By default(NULL), the value of the first OD\n#' measurement in each well is used.\n#' @param add.constant A value for r in the log(OD + r) transformation.\n#' @param plate.nrow The number of rows in the input files.\n#' @param plate.ncol The number of columns in the input files.\n#' @param input.skip.lines specifies a plate name for a single-plate read. If NULL, this is derived from the file name. \n#' @param multi.column.headers The headers of the column when analyzing multiple plates.\n#' @param single.column.headers The headers of the column when analyzing a single plate.\n#' @param layout.sheet.headers The headers of the layout file.\n#' @param silent Surpress all messages.\n#' @param verbose Display all messages when analyzing each well.\n#' \n#' @return A list of well objects.\ngcat.load.data = function(file.name = NULL, load.type = \"csv\", input.data = NULL, single.plate.ID = NULL, \n plate.layout = NULL,plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0,\n multi.column.headers = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\"), single.column.headers = c(\"\",\"A1\"), \n layout.sheet.headers = c(\"Strain\", \"Media Definition\"),\n blank.value = NULL, start.index = 2, single.plate = F, silent = T){\n\n ########################################################################\n # Read from .csv, tab-delimited text file or data frame object #\n ########################################################################\n \n\tif(is.null(input.data)){\n\t\t# Either read from .csv.\n\t\tinput.data = read.csv(file.name, stringsAsFactors=F, skip = input.skip.lines, fileEncoding='UTF-8')\n\n\t\t# Determine single plate name if not specified. \n if (is.null(single.plate.ID)){\n # Split the file name by \".\" and discard the last member (file extension). \n single.plate.ID = strsplit(basename(file.name),\"\\\\.\")[[1]]\n single.plate.ID = paste(single.plate.ID[-length(single.plate.ID)],collapse=\".\")\n }\n\t\t}\n\n # Call to arrange data from a single plate format file\n\tif(single.plate)\n\t\tinput.data = gcat.reorganize.single.plate.data(input.data = input.data, single.column.headers,\n blank.value = blank.value, single.plate.ID = single.plate.ID, plate.nrow = plate.nrow, plate.ncol = plate.ncol, silent=silent)\n\n ########################################################################\n # Search for and standardize column headers in input data #\n ########################################################################\n \n # Go through the specified column headers, determining what their positions are in the \n # input data frame and if any are missing.\n \n # Get names of column headers in input data\n input.colnames = colnames(input.data)\n \n # Create a list denoting the column numbers in input data that match each of the specified column names, \n # and a separate list for any missing columns. \n \n column.matches = c()\n missing.list = NULL\n \n\tfor(i in 1:length(multi.column.headers)){\n\t\tif (multi.column.headers[i] %in% input.colnames)\n\t\t\tcolumn.matches[i] = min(which(input.colnames == multi.column.headers[i]))\n\t\t# Take the first column in input file that matches a specified column header.\n\t\telse{\n\t\t\tmissing.list = c(missing.list, i)\n\t\t}\n\t}\n\n # If any columns are missing, stop and report error with missing column names\n\tif (is.vector(missing.list)){\n\t\tmessage = \"The following columns:\"\n\t\tfor (i in missing.list)\n\t\t\tmessage = paste(message, paste(\" \", multi.column.headers[i]), sep = \"\\n\") \n\t\tstop(message, \"\\n were not found in the data file.\")\n\t\t}\n\n # Reorder and rename the columns, using the list of matching column numbers from above.\n\tinput.data = input.data[,column.matches]\n\tnames(input.data)[1:4] = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\") \n\n # Use 'substring' to split the alphanumeric \"Well\" field into row (letters) and column (numbers)\n\tinput.data$Well.row = substring(input.data$Well, 0,1)\n\tinput.data$Well.col = as.numeric(substring(input.data$Well, 2))\n\n\n ########################################################################\n # Create an array of well objects with the Time and OD data #\n ########################################################################\n #\n # Use the by function to split up the data frame into shorter segments by well (row, column and plate)\n \n\twell.array = by(data = input.data[,c(\"OD\", \"Time\")], \n INDICES = list(input.data$Well.row,input.data$Well.col,input.data$Plate.ID), \n FUN = function(x){data.frame(Time=x$Time, OD=x$OD,stringsAsFactors = F)}, simplify = F)\n\n \n # Then apply the function (found in well.class) to each member to create a well object\n well.array = aapply(well.array,function(x){well(x$Time,x$OD)})\n\n # Differentiate any duplicate plate names in the array's dimnames \n\tnew.dimnames = dimnames(well.array)\n for (i in length(new.dimnames[[3]]):1){\n\t\tif (any(new.dimnames[[3]][-i] == new.dimnames[[3]][i]))\n\t\t\tnew.dimnames[[3]][i] = paste(\"another_\", new.dimnames[[3]][i], sep = \"\")\t\n\t\t}\n\tdimnames(well.array) = new.dimnames\n\t\n # Copy the plate/row/column names found in the dimnames into the array objects themselves (use \"position\" slot)\n\tfor(plate in unique(dimnames(well.array)[[3]])){\n\t\tfor (col in unique(dimnames(well.array)[[2]])){\n\t\t\tfor(row in unique(dimnames(well.array)[[1]])){\n\t\t\t\twell.array[[row,col,plate]]@position = c(plate=plate,row=row,col=col)\n\t\t\t\t}\n\t\t\t}\n\t\t}\n\t\t\n ########################################################################\n # Add plate layout information to well array #\n ########################################################################\n \n # Use the object to add media and strain information to the \"well.info\" slot of each well\n # Also set the value of the parameter in slot \"curve.par\" to T for empty wells. \n ########################################################################\n # Add plate layout information to well array #\n ########################################################################\n \n # Use the object to add media and strain information to the \"well.info\" slot of each well\n # Also set the value of the parameter in slot \"curve.par\" to T for empty wells. \n \n\n # If is not provided, do not add strain information, and assume all wells are inoculated. \n if(is.null(plate.layout)){ \n plate.layout = data.frame(Row=rep(PLATE.LETTERS[1:plate.nrow],plate.ncol),Column=rep(1:plate.ncol,each=plate.nrow),rep(\"Unknown Strain\",96),rep(\"Unknown Media\",96))\n colnames(plate.layout) = c(\"Row\", \"Column\", layout.sheet.headers)\n } \n else\n if(!silent) cat(\"\\n\\t\\tusing plate layout to fill well info.\")\n \n\tfor(plate in unique(dimnames(well.array)[[3]])){\n\t\tfor (col in unique(dimnames(well.array)[[2]])){\n\t\t\tfor(row in unique(dimnames(well.array)[[1]])){\n well = well.array[[row,col,plate]]\n # For each well on each plate, find the corresponding row in . \n # If refers to specific plates, then use those to find the correct row. \n # Otherwise, generalize across all plates. \n if (\"Plate.ID\" %in% names(plate.layout)) \n layout.row.number = which(plate.layout$Column==well@position[\"col\"] & \n plate.layout$Row==well@position[\"row\"] & \n plate.layout$Plate.ID==well@position[\"plate\"] ) \n else \n layout.row.number = which(plate.layout$Column==well@position[\"col\"] & \n plate.layout$Row==well@position[\"row\"])\n \n # Error if either no rows or more than one row matches the well\n if (length(layout.row.number) != 1)\n stop(\"incorrectly formatted plate layout! check names of columns, rows, and plates (if applicable).\")\n \n # Add any additional columns to the well's \"well.info\" slot\n well.info = plate.layout[layout.row.number,!(names(plate.layout) %in% c(\"Row\",\"Column\",\"Plate.ID\",layout.sheet.headers))]\n \n # Fix the column name issue if there is only one additional entry. \n if(length(well.info) == 1){\n well.info = data.frame(well.info,stringsAsFactors=F) \n names(well.info) = names(plate.layout)[!(names(plate.layout) %in% c(\"Row\",\"Column\",\"Plate.ID\",layout.sheet.headers))] \n } \n well@well.info = as.list(well.info)\n \n well@well.info$Strain = plate.layout[layout.row.number, layout.sheet.headers[1]]\n well@well.info$Media = plate.layout[layout.row.number, layout.sheet.headers[2]]\n \n # Set parameter in slot \"curve.par\" accordingly \n well@curve.par$empty.well = (plate.layout$Strain[layout.row.number] == \"Empty\") \n well.array[[row,col,plate]] = well\n\t\t\t\t}\n\t\t\t}\n\t\t}\n \n # Set start index value in each well\n well.array = aapply(well.array, function(x,start.index) { x@start.index = start.index; x }, start.index)\n \n ########################################################################\n # Return values to R #\n ######################################################################## \n # \n\t# Console output if desired, return the completed well array.\n\tif (!silent)\n\t\tcat(\"\\n\\t\", dim(well.array)[[3]], \"plates added to array from\", file.name)\n\treturn(well.array)\n\t}\n\n\n\n########################################################################\n# #\n# Reorganize data from single-plate input format before reading #\n# #\n########################################################################\n#\n# This function reorganizes the data frame from a single-plate format file. \n# input.data - data frame read straight from a single-plate format data file. \n# single.plate.ID - specifies a plate name for a single-plate read, since none is given in the single-plate format\n# The plate will be named Plate_1 unless otherwise specified. \n\ngcat.reorganize.single.plate.data = function(input.data, blank.value = NULL, single.column.headers, single.plate.ID = \"Plate_1\", \n plate.nrow = 8, plate.ncol = 12, silent=T){\n \n ########################################################################\n # Standardize the formatting and search for specified column names #\n ######################################################################## \n # \n # Locate the first and last rows from the table and return errors if not defined \n # Note: this only works if the time column is the first column\n \n\theader.row = min(which(input.data[,1] == single.column.headers[1])) \n if (length(header.row) != 1 | is.infinite(header.row))\n stop(\"could not locate header row in input file!\")\n \n # The last row: where column 2 starts to be blank, or the total number of rows, whichever is smaller \n extra.rows.start = min(which(input.data[-(1:header.row),2] == \"\"), which(is.na(input.data[-(1:header.row),2])), nrow(input.data[-(1:header.row),]))\n if (length(extra.rows.start) != 1 & is.infinite(extra.rows.start)) \n stop(\"could not locate last row in input file!\")\n\n # Use header row to rename the columns, then cut off extra rows (including the ones above header)\n\tnames(input.data) = as.character(unlist(input.data[header.row,]))\n input.data = input.data[(header.row+1):(header.row+extra.rows.start-1),]\n \n # Time column: allow for multiple matches to the name (since it's usually blank) but assume it's the first one\n\tTime.column = which(names(input.data) == single.column.headers[1])\n\tif (length(Time.column) != 1){\n if(!silent) cat(\"No unique time column in input.data file! Using the first one encountered.\")\n\t\tTime.column = min(Time.column)\n\t\t}\n # First well column (default A1): only allow for one match.\t\n\tWell.column.start = which(names(input.data) == single.column.headers[2])\n\tif (length(Well.column.start) != 1)\n\t\tstop(\"No unique start point for well columns in input.data file!\")\n\n # If the time column was found, rename it \"Time\" and reformat it into a numeric value\n # Adjust the blank measurement timestamp to -1 seconds if there is one\n\n names(input.data)[Time.column] = \"Time\"\n \n # Note: Some single plate screens have timepoints ending with \"s\" for seconds. \n # This line removes the \"s\" while maintaining general compatibility. \n\tinput.data$Time = unlist(strsplit(input.data$Time, \"s\"))\n\n # If is NULL (default - takes the first OD as the blank reading), then the first timepoint can labeled something non-numeric. \n # In that case, rename it to match the first real timepoint minus one. \n # when user input blank value, Blank timepoint i.e. input.data$Time[1] == Blank, labeled as \"Blank\" from data input file\n # It also should rename it to match the first real timepoint minus one. \n if(is.null(blank.value) || is.na(as.numeric(input.data$Time[1])))\n input.data$Time[1] = as.numeric(input.data$Time[2]) - 1\n \n ########################################################################\n # Start to fill the reformatted data frame #\n ######################################################################## \n \n # If all columns are present, make a list of all the wells.\n\twell.list = paste(rep(PLATE.LETTERS[1:plate.nrow], each = plate.ncol), rep(formatC(1:plate.ncol, digits = log(plate.ncol, 10), flag = \"0\"), plate.nrow), sep = \"\")\n\n #\tDuplicate the well names times the number of time measurements in each well\t\n Well = rep(well.list, each = length(input.data[,Time.column]))\t\n\t\t\n # Duplicate times the length of the entire output \n Plate.ID = rep(single.plate.ID, length(Well))\n\n # Duplicate the time column times the number of wells and convert to numeric\n\tTime = as.numeric(rep(input.data[,Time.column], times = length(Well.column.start:ncol(input.data))))\n\n # Append OD measurements from each well together and convert to numeric\n\tOD = c()\n\tfor (i in Well.column.start:ncol(input.data)){\n\t\tOD = as.numeric(c(OD, input.data[,i]))\n\t\tOD = unlist(OD)\n\t\t}\n\n # Fill and return the data frame containing the above four vectors.\n\toutput = data.frame(Plate.ID, Well, OD, Time)\t\n\t\n # Include any extra columns that were not Time or OD measurements?\n\tfor(i in (1:length(names(input.data)))[-c(Time.column,Well.column.start:ncol(input.data))]){\n\t\tnew.column = data.frame(rep(input.data[,i], length(Well.column.start:ncol(input.data))))\n\t\tnames(new.column) = names(input.data)[i]\n\t\toutput = cbind(output, new.column)\n\t\t}\t\n\treturn(output)\n}\n\n\n########################################################################\n# #\n# Function to combine two well array datasets by plate #\n# #\n########################################################################\n# ----------------------------------------------------------\n# This function can append together arrays created using \n# Arguments: any number of array objects as long as they are all output straight from \n\ngcat.append.arrays = function(...){\n\n # Transfer arrays to a list\n\targs.arrays = list(...)\n\tfirst.array = args.arrays[[1]]\n\tfirst.dims = dim(first.array)\n plate.nrow = args.arrays[[4]]\n plate.ncol = args.arrays[[3]]\n input.arrays = list(args.arrays[[1]], args.arrays[[2]])\n\tfor (i in 2:length(input.arrays)){\n\t\tnext.array = input.arrays[[i]]\n\t\tnext.dims = dim(next.array)\n\t\n # Check to make sure the arrays have proper dimensions for 96-well plate data\n if (!(all(c(first.dims[1], next.dims[1]) == plate.nrow) & all(c(first.dims[2], next.dims[2]) == plate.ncol)))\n\t\t\tstop(\"Data should have dimensions (\",plate.nrow,\",\",plate.ncol,\",n)!\")\n\t\t\n # If dimensions are alright, append dimensions and dimension names\t\n new.dims = c(plate.nrow,plate.ncol,first.dims[3]+next.dims[3])\n\t\tnew.names = dimnames(first.array)\n\t\tnew.names[[3]] = c(new.names[[3]], dimnames(next.array)[[3]]) \n\n # Differentiate duplicate names\n\t\tfor (i in length(new.names[[3]]):1){\n\t\t\tif (any(new.names[[3]][-i] == new.names[[3]][i]))\n\t\t\t\tnew.names[[3]][i] = paste(\"another_\", new.names[[3]][i], sep = \"\")\n\t\t\t}\n\t\n\t\t# Create a new array\n new.array = c(first.array, next.array)\n\t\tdim(new.array) = new.dims\n\t\tdimnames(new.array) = new.names\n\t\t\n # Update plate name in well objects\n\t\tfor (i in 1:length(unlist(new.array)))\n\t\t\tnew.array[[i]]@position[1] = new.names[[3]][floor((i-1)/96)+1]\n\t\t\n\t\t# Loop until complete\n\t\tfirst.array = new.array\n\t\tfirst.dims = dim(first.array)\n\t\t}\n\treturn(new.array)\n\t}\n\n\n########################################################################\n# #\n# Convert timestamps to hours from start and sort timepoints #\n# #\n########################################################################\n#\n# This function acts on a single well and modifies the raw data stored in slot \"screen.data\"\n# \n# input.well - an object of class well\n# time.format - specifies the time format. allowed values are \"%S\", for seconds, \"%d\", for days, or anything complying with ISO C / POSIX standards; see \n# note: reading directly from excel to R results in timestamps being converted to days.\n# start.index - which timepoint should be used as the starting time at inoculation?\n\ngcat.start.times = function(input.well, time.input, start.index = 2) { \n \n if(start.index > length(input.well))\n stop(\"Start index is greater than length of well!\")\n \n # If using a numeric time format, simply multiply times by the appropriate conversion factor\n\t# Conversion factor should be supplied to convert timestamps to hours. For example, \n # should be equal to 1/3600 if \"time\" is in seconds, 24 if \"time\" is in days, etc.\n\n time.format = time.input # Set default constant from rails user input\n \n if (is.numeric(time.format))\n\t\tinput.well@screen.data$Time = (input.well@screen.data$Time - min(input.well@screen.data$Time)) * time.format\n else{\n # Otherwise, convert timestamps from ISO C / POSIX to numeric values representing seconds (since the dawn of time?) and subtract out the initial value. \n\t\trtimes = input.well@screen.data$Time\n ptimes = strptime(as.character(rtimes),time.format)\n\t\tctimes = as.POSIXct(ptimes)\n\t\tint.times = unclass(ctimes)\n # Time will be in seconds, convert to hours by dividing by 3600\n\t\tinput.well@screen.data$Time = (int.times - min(int.times))/3600\n\t\t}\n\t# Sort raw data by timestamp and return the input.well\n\tinput.well@screen.data = input.well@screen.data[order(input.well@screen.data$Time),]\n \n input.well@screen.data$Time = input.well@screen.data$Time - input.well@screen.data$Time[start.index]\n \n if(all(is.na(input.well@screen.data$Time)))\n stop(\"Error in .\") \n rownames(input.well@screen.data) = NULL\n\treturn(input.well)\n }\t\n",
- "created" : 1425413240737.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2796681086",
- "id" : "2FF709FD",
- "lastKnownWriteTime" : 1428438949,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/table2well.R",
- "project_path" : "R/table2well.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/3DC2D3DA b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/3DC2D3DA
deleted file mode 100644
index be33fea..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/3DC2D3DA
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n# GCAT version 5.00\n# Notes by Jason\n# 08/18/2011\n\n# Initialization\n\n\nPLATE.LETTERS = paste(rep(c(\"\", LETTERS), each = 26), rep(LETTERS, 26), sep=\"\")\nglobal.version.number = packageDescription(pkg=\"GCAT\")$Version \n\n########################################################################\n# #\n# Top-level functions for analysis of screening data from .csv files. #\n# #\n########################################################################\n# This functions is called directly by the user interface.\n# They in turn call the main function (below) multiple times for each data file provided in . \n\n# Arguments:\n# file.list - a list of full paths to .csv files. all files must be in the same format (see )\n# single.plate - are the file in the single plate (wide) format vs. the multi-plate (long) format?\n# layout.file - (optional) provide full path to a layout file with strain and media definitions (applies to all files in list)\n\n# out.dir - name a directory to output the table of curve parameters to (defaults to working directory) \n# graphic.dir - name a directory to output the images of the fitted curves to (defaults to subdirectory \"pics\" of above)\n\n\n# add.constant- should be a numeric constant that will be added to each curve before the log transform (defaults to 1) \n# blank.value - user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well. \n# start.index - which timepoint should be used as the first one after inoculation (defaults to the 2th one)\n# growth.cutoff - minimum threshold for curve growth. \n# points.to.remove - a list of numbers referring to troublesome points that should be removed across all wells.\n# remove.jumps - should the slope checking function be on the lookout for large jumps in OD?\n \n# silent - should messages be returned to the console?\n# verbose - should sub-functions return messages to console? (when I say verbose, I mean it!)\n\n# Returns:\n# if = F (default), avector of full paths to all the files generated by the function.\n# otherwise, the fitted array of well objects. \n\n# Use this function to analyze any set of .csv files using the same plate layout info. \n\n#' Analyze screening growth data from the given .csv files.\n#' \n#' Top-level GCAT function \n#' \n#' @param file.list A list of full paths to .csv files. all files must be in the same format (see )\n#' @param single.plate The file in the single plate (wide) format vs. the multi-plate (long) format?\n#' @param layout.file Full path to a layout file with strain and media definitions (applies to all files in list)\n#' @param out.dir A directory to output the table of curve parameters to (defaults to working directory) \n#' @param graphic.dir A directory to output the images of the fitted curves to (defaults to subdirectory \"pics\" of above)\n#' @param use.linear.param Whether to use linear parameters or not?\n#' @param use.loess Whether to use LOESS model or not?\n#' @param smooth.param Smoothing parameter for LOESS model.\n#' @param add.constant A numeric constant that will be added to each curve before the log transform (defaults to 1) \n#' @param blank.value User can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well. \n#' @param start.index Which timepoint should be used as the first one after inoculation (defaults to the 2th one)\n#' @param growth.cutoff Minimum threshold for curve growth. \n#' @param points.to.remove A list of numbers referring to troublesome points that should be removed across all wells.\n#' @param remove.jumps Should the slope checking function be on the lookout for large jumps in OD?\n#' @param time.input The time setting in which the current system is running?\n#' @param plate.nrow The number of rows in a plate.\n#' @param plate.ncol The number of columns in a plate.\n#' @param input.skip.lines If specified, this number of lines shall be skipped from the top when reading the input file with read.csv \n#' @param multi.column.headers The headers of the result tabular data when analyzing multiple plates at once.\n#' @param single.column.headers The headers of the result tebaular data when analyzaing a single plate.\n#' @param layout.sheet.headers The headers of the layout file?\n#' @param silent Shoulde messages be returned to the console?\n#' @param verbose Should sub-functions return messages to console? (when I say verbose, I mean it!)\n#' @param overview.jpgs Should GCAT enable an overview image?\n#' \n#' @return A list of the output files.\ngcat.analysis.main = function(file.list, single.plate, layout.file = NULL, \n out.dir = getwd(), graphic.dir = paste(out.dir, \"/pics\", sep = \"\"), \n add.constant = 0.1, blank.value = NULL, start.index = 2, growth.cutoff = 0.05,\n use.linear.param = F, use.loess = F, smooth.param=0.1,\n points.to.remove = 0, remove.jumps = F, time.input = NA,\n plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0,\n multi.column.headers = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\"), single.column.headers = c(\"\",\"A1\"), \n layout.sheet.headers = c(\"Strain\", \"Media Definition\"),\n silent = T, verbose = F, return.fit = F, overview.jpgs = T){\n\n # MB: Prototyping system unwanted argument guarding. Proper function \n # will be added in the future.\n # Not the best solution.\n if (is.na(time.input)) {\n if (single.plate)\n time.input = 1/3600\n else\n exception(\"Error: \", \"time.input is NA.\")\n }\n \n if (add.constant < 0)\n exception(\"Error: \", \"The constant r should not be negative.\")\n # End prototyping temporary solution.\n \n upload.timestamp = strftime(Sys.time(), format=\"%Y-%m-%d %H:%M:%S\") # Get a timestamp for the time of upload. \n fitted.well.array.master = list()\n source.file.list = c()\n \n dim(fitted.well.array.master) = c(plate.nrow,plate.ncol,0)\n dimnames(fitted.well.array.master) = list(PLATE.LETTERS[1:plate.nrow], 1:plate.ncol, c())\n \n for(file.name in file.list){\n \n # Call on the file with single plate options\n \tfitted.well.array = try(gcat.fit.main(file.name = file.name, load.type = \"csv\", \n single.plate = single.plate, layout.file = layout.file, start.index = start.index, \n time.input = time.input, add.constant = add.constant, blank.value = blank.value, \n growth.cutoff = growth.cutoff, points.to.remove = points.to.remove, remove.jumps = remove.jumps,\n \t use.linear.param=use.linear.param, use.loess=use.loess, smooth.param=smooth.param,\n plate.nrow = plate.nrow, plate.ncol = plate.ncol, multi.column.headers = multi.column.headers, \n single.column.headers = single.column.headers, layout.sheet.headers = layout.sheet.headers,\n input.skip.lines = input.skip.lines, silent = silent, verbose = verbose), silent = T)\n \n # Return error message if the function fails.\n if(class(fitted.well.array) == \"try-error\")\n return(as.character(fitted.well.array)) \n }\n \n # Add fitted well array onto existing fitted wells\n fitted.well.array.master = gcat.append.arrays(fitted.well.array.master, fitted.well.array, plate.ncol, plate.nrow)\n \n # Remove the \"processed_\" tag from file names and add to the list of source files.\n source.file.list = c(source.file.list, basename(paste(strsplit(file.name, \"processed_\")[[1]],collapse=\"/\"))) \n\n out.files = try(gcat.output.main(fitted.well.array.master, out.prefix = \"output\", \n source.file.list = source.file.list, upload.timestamp = upload.timestamp, \n growth.cutoff = growth.cutoff, add.constant = add.constant, blank.value = blank.value, start.index = start.index, \n points.to.remove = points.to.remove, remove.jumps = remove.jumps, \n out.dir = out.dir, graphic.dir = graphic.dir, overview.jpgs=overview.jpgs,\n use.linear.param=use.linear.param, use.loess=use.loess, plate.ncol = plate.ncol, plate.nrow = plate.nrow,\n silent = silent), silent = T)\n \n # Return file list or error message otherwise return \"successful analysis\" message (?)\n\n # file.list = c(\"Data was successfully analyzed.\", file.list) # <--- yet to be implemented. causes errors downstream right now \n if(class(out.files) == \"try-error\") return(as.character(out.files)) \n \n if(return.fit) return(fitted.well.array.master)\n else return(out.files) \n}\n\n########################################################################\n# #\n# Main function for analysis of screening data from input tables. #\n# #\n########################################################################\n# This is the main function that handles all the analyses for files in both single and multiple plate formats. \n# It is called by the top level function \n#\n# It then calls the following functions on each member of the array: \n# - curve normalization and standardization: , , , , \n# - curve shape analysis before model fitting: , , \n# - to fit a nonlinear model to the growth data: \n# Finally, it returns the fitted array of well objects.\n\n#' Main analysis function for GCAT \n#' \n#' This is the main function that handles all the analyses for data files in both single and multiple plate formats. \n#' It is called by the top level function \\code{gcat.analysis.main} along with \\code{gcat.output.main}.\n#' \n#' @param file.name Complete path and file name of a comma-separated values (.csv) file containing growth curve data \n#' in the multiple-plate (long) format. \n#' @param input.data A list of tables representing input files read with \\code{read.table}. Used to save time in cases\n#' of running multiple analyses on the same dataset. If used, the function will ignore \\code{file.name} entirely.\n#' @param load.type .csv by default. \n#' @param layout.file Specifies the location of a layout file containing identifying information.\n#' @param single.plate Whether the GCAT is analyzing a single plate or not.\n#' @param blank.value Blank OD measurement for uninoculated wells. By default(NULL), the value of the first OD\n#'measurement in each well is used.\n#' @param start.index Which timepoint should be used as the first one after inoculation?\n#' @param time.input Either a character describing the format used to convert timestamps in the input to numbers\n#' representing number of seconds (see \\code{strptime}), or a factor to divide entries in the Time column by to get the \n#' numbers of hours.\n#' @param normalize.method Describes the method used by \\code{normalize.ODs} to normalize cell density values using blank reads.\n#' @param add.constant A value for r in the log(OD + r) transformation.\n#' @param use.log Should the analysis use log on all values.\n#' @param points.to.remove A vector of integers specifying which timepoints should be removed across all wells. \n#' By default(0) none are marked for removal.\n#' @param use.linear.param Should the linear parameter be used or not.\n#' @param use.loess Should the loess model be used or not.\n#' @param smooth.param If loess model is used, this parameter define the smoothing parameter for the loess model.\n#' @param fall.cutoff A cutoff used by \\code{check.slopes} to decide on thresholds for jumps and tanking.\n#' @param growth.cutoff A threshold used by check.growth to decide whether a well displays growth.\n#' @param remove.jumps Should jumps in OD detected by the subfunction \\code{check.slopes}?\n#' @param plate.nrow The number of rows in the input files.\n#' @param plate.ncol The number of columns in the input files.\n#' @param input.skip.lines If specified, this number of lines shall be skipped from the top when reading the input file with read.csv \n#' @param multi.column.headers The headers of the column when analyzing multiple plates.\n#' @param single.column.headers The headers of the column when analyzing a single plate.\n#' @param layour.sheet.headers The headers of the layout file.\n#' @param growth.model What growth model should be used?\n#' @param backup.growth.model If the main growth model fails, the back up model will be used.\n#' @param silent Surpress all messages.\n#' @param verbose Display all messages when analyzing each well.\n#' \n#' @return An array of well objects \ngcat.fit.main = function(file.name, input.data = NULL, load.type = \"csv\", layout.file = NULL, \n single.plate = F, blank.value = NULL, start.index = 2, time.input = NA,\n normalize.method = \"default\", add.constant = 1, use.log = T, points.to.remove = 0,\n use.linear.param=F, use.loess=F, smooth.param=0.1,\n fall.cutoff = -0.0025, growth.cutoff = 0.05, remove.jumps = F, \n plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0,\n multi.column.headers = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\"), single.column.headers = c(\"\",\"A1\"), \n layout.sheet.headers = c(\"Strain\", \"Media Definition\"),\n growth.model = NA, backup.growth.model = NA, \n silent = F, verbose = F){\n \n # Explanation of arguments:\n \n # ---File Handling---\n # file.name - full path to an excel spreadsheet, .csv or tab-delimited text file, in either the single or multiple-plate format\n # input.data - use pre-loaded data set (output from function only). will override if not NULL \n # load.type - supports \"csv.\" \n # layout.file - full path to a file containing the plate layout in the same format as . will not be used if is \"xlsx\"\n\n # ---Input file format---\n # single.plate - true denotes data in single-plate format, i.e. simple OD output. false denotes multiple-plate robotic screening output. \n # note: reading directly from excel to R results in timestamps being converted to days.\n\n # ---Normalization and Transforms---\n # blank.value - user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well. \n # start.index - which timepoint should be used as the first one after inoculation (defaults to the 2th one)\n # normalize.method - how should each growth curve be normalized? allowed values are:\n # \"first\": subtracts the first OD, assumed to be the blank, from all ODs\n # \"none\": does nothing, assumes no blank. highly recommend log(OD+1) transform in this case.\n # \"average.first\": forces all filled wells on each plate to match the average value at (after subtracting the first OD) \n # add.constant - a numeric constant that will be added to each curve before the log transform (defaults to 1)\n # use.log - should a log transform be applied to the data after normalization? \n # points.to.remove - a list of numbers referring to troublesome points that should be removed across all wells.\n \n # ---Pre-fitting processing---\n # fall.cutoff - a cutoff value for determining whether OD falls significantly between two timepoints. see in prefit.processing.R for details. \n # growth.cutoff - a cutoff value for determining whether a well contains a successfully growing culture or not.\n # remove.jumps - should the slope checking function be on the lookout for large jumps in OD?\n \n # ---Model fitting---\n # model - which parametrized growth model to use? can be richards, gompertz, or logistic. models are defined as objects of class model, see \"model.class.R\"\n # backup.model - which model should be used if fitting using fails? should ideally be simpler than the main model (less parameters)\n \n # ---Miscellanous input/output preferences---\n # silent - should messages be returned to the console?\n # verbose - should sub-functions return messages to console? (when I say verbose, I mean it!) \n # unlog - should exported graphics be transformed back to the OD scale? \n # return.fit - should the function return an array of wells? if not, it will return a list of generated files. \n \n \t\n ########################################################################\n # Read from .csv file #\n ########################################################################\n #\n # The functions used here are found in table2well.R\n \n if(!silent) cat(\"\\nReading input files...\")\n # Read from .csv or tab-delimited text file using (in load.R)\n # if is provided, it will be used here. \n \n plate.layout = NULL \n # Read layout file if it is specified. \n if(!is.null(layout.file)){\n if(load.type==\"csv\") plate.layout = read.csv(layout.file,header=T,stringsAsFactors=F)\n else plate.layout = read.table(layout.file,header=T,sep=\"\\t\",stringsAsFactors=F) \n if(!silent) cat(\"\\n\\tAdded plate layout information from\", layout.file, \"\\n\")\n }\n \n # Load the data \n\t\twell.array = try(gcat.load.data(file.name = file.name, input.data = input.data, \n plate.layout = plate.layout, plate.nrow = plate.nrow, plate.ncol = plate.ncol, \n input.skip.lines = input.skip.lines, multi.column.headers = multi.column.headers, \n single.column.headers = single.column.headers, layout.sheet.headers = layout.sheet.headers,\n blank.value = blank.value, start.index = start.index, single.plate = single.plate, \n load.type = load.type, silent=silent),silent=silent)\n\n # Return an error if there is a problem with file loading. \n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n\t\t\n\n # !---At this point, is an array of well objects, each containing raw data and media/strain information if provided--- \n \n # Attempt to apply time formatting to all wells in array \n well.array = try(aapply(well.array, gcat.start.times, start.index = start.index, time.input = time.input),silent=silent)\n \n # Return an error if there is a problem with time formatting\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \t\n ########################################################################\n # Perform normalization and transformation of raw data #\n ######################################################################## \n #\n # The functions used here are found in normalize.and.transform.R \n \n if(!silent) cat(\"\\nProcessing raw data...\")\n \n # Set all timepoints to active for now using \"points.to.remove=0\" argument with \n # adds an extra column to the \"well.array\" slot of each well specifying which points to remove when data is retrieved from the well\n well.array = aapply(well.array, remove.points, points = 0)\n \n # Normalize ODs using specified method and adding a constant if desired. \n # sets the \"norm\" slot of each well to a value to be subtracted from OD values whenever data is retrieved from the well\n well.array = try(normalize.ODs(well.array, normalize.method = normalize.method, \n start.index = start.index, blank.value = blank.value, add.constant = add.constant),silent=silent)\n \n # Return an error if there is a problem with normalization\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \n # Transform ODs on the logarithmic scale, regardless of whether is true \n # an extra column of log-transformed values is added to the \"well.array\" slot of each well \n # the \"use.log\" slot of each well is set instead to determine whether the transformed values will be returned when data is retrieved from the well.\n well.array = try(aapply(well.array, transform.ODs, start.index = start.index, blank.value = blank.value, use.log = use.log, constant.added = add.constant),silent=silent)\n \n # Return an error if there is a problem with transformation\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \n # Remove specified timepoints across wells (use \"points.to.remove=NULL\" if no points to remove) \n well.array = try(aapply(well.array, remove.points, points = points.to.remove),silent=silent)\n \n # Return an error if there is a problem with point removal\n if (class(well.array) == \"try-error\")\n stop(\"Error in : \", well.array) \n \n \n \n ########################################################################\n # Pre-fitting data processing (analysis of curve shapes) #\n ######################################################################## \n #\n # The functions used here are found in slope.analysis.R\n \n # Estimate slope at each timepoint \n # add a column to the \"well.array\" slot of each well with the local slope at each timepoint \n \n well.array = try(aapply(well.array, calculate.slopes, silent=!verbose),silent=silent)\n \n # Return an error if there is a problem with slope calculation\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \n # Check slopes for tanking and/or jumping behavior\n # fills the \"curve.par\" slot of each well with , denoting the timepoint at which tanking starts (if none, value is NA) \n # uses to remove all points after \n # It will also fill the \"jump.error\" slot with a status message, and try to use an automated process to remove the \n # erroneous points if is true (default false). \n \n well.array = try(aapply(well.array, check.slopes, fall.cutoff = fall.cutoff, remove.jumps = remove.jumps, silent=!verbose, draw = F),silent=silent)\n \n # Return an error if there is a problem with slope analysis\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \t\n \t\n # Check curves for growth above cutoff\n # fills the \"curve.par\" slot of each well with , denoting whether the well has no detectable growth. \n well.array = try(aapply(well.array, check.growth, growth.cutoff = growth.cutoff, start.index = start.index),silent=silent)\n \n # Return an error if there is a problem with growth.check\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \t\n ########################################################################\n # Fit parameterized models to data #\n ######################################################################## \n #\n # The functions used here are found in fit.model.R \n \n # Fit each well with the selected model and attempt to catch failed fittings with the backup model \n # skips wells designated as above\n # fills the \"fit.info\" slot of each well with \"success,\" \"failed,\" or \"skipped\"\n # if fit was successful:\n # fills the \"equation\" and \"model.name\" slots with the relevant info for the successful model \n # fills the \"fit.par\" slot with fitted parameters if fit is successful \n \n if(!silent) cat(\"\\nFitting models to data...\")\n well.array = aapply(well.array, fit.model, growth.model=growth.model,\n backup.growth.model = backup.growth.model, use.linear.param=use.linear.param,\n use.loess=use.loess, smooth.param=smooth.param, silent=!verbose)\n \n # Return an error if there is a problem with model fitting\n if (class(well.array) == \"try-error\")\n \tstop(\"Error in : \", well.array)\n \n if(!silent) cat(\"\\ndone!\\n\") \n return(well.array) \n } \n \n \n########################################################################\n# #\n# Output function for generating files from fitted data. #\n# #\n######################################################################## \n#' Output function for generating files from fitted data.\n#' \n#' Handles files and directories, calls \\code{table.out}, \\code{plate.overview} and \\code{view.fit} \n#' to generate output tables and graphics.\n#' \n#' @param fitted.well.array A list of fitted well objects.\n#' @param out.prefix Prefix that is in the name of output files.\n#' @param blank.value User can enter a blank OD measurement for uninoculated wells. \n#' If NULL, defaults to the value of the first OD measurement of each well.\n#' @param start.index Which timepoint should be used as the first one after inoculation (defaults to the 2th one)\n#' @param growth.cutoff Minimum threshold for curve growth.\n#' @param points.to.remove A list of numbers referring to troublesome points that should be removed across all wells.\n#' @param remove.jumps Should the slope checking function be on the lookout for large jumps in OD?\n#' @param out.dir name a directory to output the table of curve parameters to (defaults to working directory) \n#' @param graphic.dir name a directory to output the images of the fitted curves to \n#' (defaults to subdirectory \"pics\" of above)\n#' @param overview.jpgs should jpgs be generated for each plate with the overview graphic? \n#' This is for backwards compatibility with the old web server. \n#' @param silent should messages be returned to the console?\n#' @param unlog should exported graphics be transformed back to the OD scale?\n#' @param constant.added (should be the same value as add.constant above) - \n#' used to readjust for the constant added during the log transform when plotting ODs. \n#' @return A list of output files if success.\n\ngcat.output.main = function(fitted.well.array, out.prefix = \"\", source.file.list, upload.timestamp = NULL, \n add.constant, blank.value, start.index, growth.cutoff, points.to.remove, remove.jumps, \n out.dir = getwd(), graphic.dir = paste(out.dir,\"/pics\",sep = \"\"), overview.jpgs = T,\n use.linear.param=F, use.loess=F, plate.nrow = 8, plate.ncol = 12,\n unlog = F, silent = T){ \n \n # Prepare timestamp for addition to output file names. \n filename.timestamp = strftime(upload.timestamp, format=\"_%Y-%m-%d_%H.%M.%S\")\n \t\n ########################################################################\n # Prepare to write to output files #\n ######################################################################## \t\n\t if(is.null(blank.value)) blank.value = \"First timepoint in well\"\n\t \n if(!silent) cat(\"\\nFinding/creating new output directories...\")\n \n old.wd = getwd() \n # Create output directory if it doesn't exist\n if(class(try(setwd(out.dir), silent = T)) == \"try-error\"){\n \tif(!silent) cat(\"\\ncreating new output directory\")\n \tif (class(try(dir.create(out.dir))) == \"try-error\")\n \t\tstop(\"Error creating new output directory!\")\n \t}\n \n # Create graphics directory if it doesn't exist\n if(class(try(setwd(graphic.dir), silent = T)) == \"try-error\"){\n \tif(!silent) cat(\"\\ncreating new graphics directory\")\n \tif (class(try(dir.create(graphic.dir))) == \"try-error\")\n \t\tstop(\"Error creating new graphics directory!\")\n \t}\n \t\n ########################################################################\n # Populate a data table with fit results and write to file #\n ######################################################################## \n # \n # The functions used here are found in table.output.R \n \n # Creates a table with a row for each well and a column for each of various identifiers and fitted and calculated parameters. \n \n if(!silent) cat(\"\\nPopulating data table...\")\n table.fit = try(table.out(fitted.well.array, filename.timestamp=filename.timestamp,use.linear.param=use.linear.param, use.loess=use.loess, constant.added=add.constant))\n \n # Return an error if there is a problem with returning the table\n if (class(fitted.well.array) == \"try-error\")\n \tstop(\"Error in : \", fitted.well.array)\n \t\n \t\n # Set working directory to \n\tif (class(try(setwd(out.dir))) == \"try-error\")\n\t\tstop(\"Error setting directory for table output\")\n\t\t\n\t# Write output table to file in \n table.filename = paste(out.dir, \"/\", out.prefix, \"_gcat.fit\", filename.timestamp, \".txt\", sep = \"\")\n if (class(try(write.table(table.fit, table.filename, sep = \"\\t\", row.names = F))) == \"try-error\")\n\t\tstop(\"Error writing tabular output\")\n\t\t\n \t# ---If successfully written, add postscript and start a list of generated files.\t\n generated.files = table.filename\n\n ########################################################################\n # Write individual fit and overview graphics to file #\n ######################################################################## \n # \n # The functions used here are found in graphic.output.R\n \n if(!silent) cat(\"\\nDrawing graphics...\")\n \n # Set working directory to \n\tif (class(try(setwd(graphic.dir))) == \"try-error\")\n\t\tstop(\"Error setting directory for graphic output\")\n\t# Use function to write fit graphics to file. \n \n\tgraphic.files = try(pdf.by.plate(fitted.well.array, out.prefix=out.prefix, upload.timestamp = upload.timestamp, \n unlog=unlog,constant.added=add.constant,overview.jpgs=overview.jpgs, plate.ncol = plate.ncol, plate.nrow = plate.nrow),silent=silent)\n \n if (class(graphic.files) == \"try-error\")\n\t\tstop(\"Error in : \", graphic.files)\n \n # If successfully written, add to the list of generated files.\t\n generated.files = c(generated.files, graphic.files)\n \n ########################################################################\n # Add a postscript to the output table with legend and file info. #\n ######################################################################## \n # \t\n sink(table.filename, append = T)\n analysis.timestamp = strftime(Sys.time(), format=\"%Y-%m-%d %H:%M:%S\")\n cat(\"\\n# Raw OD values are adjusted and log-transformed before fitting a growth curve as follows: log.OD = log(OD - blank + const) where blank is OD of blank medium and const is specified by the user (1 by default)\",\n \"\\n# Values are reported on the above 'log.OD' scale unless otherwise specified.\",\n \"\\n# .SE columns report standard errors of those values that are estimated directly as parameters of global sigmoid models.\",\n \"\\n# .OD columns report values back-transformed to the linear 'OD - blank' scale.\",\n \"\\n\")\n \n cat(\"\\n# -- Explanation of columns --\",\n \"\\n# - model: Name of the model the well was successfully fit with (if any)\",\n \"\\n# - lag.time: Lag time estimate inferred from the fitted model\",\n \"\\n# - inflection.time: inflection time point of the growth curve when drawn on the log scale\",\n \"\\n# - max.spec.growth.rate: maximum specific growth rate estimate inferred from the fitted model. Estimated as the first derivative of the growth curve at inflection time point\", \n \"\\n# - baseline: growth curve baseline. Global sigmoid model: baseline is parameter 'b' of the model. LOESS: baseline is the same as the lowest predicted log.OD value\",\n \"\\n# - amplitude: difference between upper plateau and baseline values. Global sigmoid model: amplitude is parameter 'A' of the model. LOESS: amplitude = max.log.OD - min.log.OD\",\n \"\\n# - plateau: upper asymptote value of the fitted model. Global sigmoid model: plateau = b + A. LOESS: plateau = max.log.OD\",\n \"\\n# - inoc.log.OD: log.OD value at inoculation. Estimated value from the fitted model is used, rather than the actual measurement\", \n \"\\n# - max.log.OD: maximal log.OD value reached during the experiment. Estimated value from the fitted model is used rather than the actual measurement\",\n \"\\n# - projected.growth: maximal projected growth over inoculation value. Global sigmoid model: projected.growth = plateau - inoc.log.OD. LOESS: not reported\",\n \"\\n# - achieved.growth: maximal growth over inoculation value actually achieved during the experiment. achieved.growth = max.log.OD - inoc.log.OD\",\n \"\\n# - shape.par: shape parameter of the Richard equation\",\n \"\\n# - R.squared: goodness of fit metric. Also known as coefficient of determination. R.squared is usually between 0 and 1. A value close to 1 indicates good fit.\",\n \"\\n# - RSS: residual sum of squares. Another goodness of fit metric. Smaller values indicate better fits.\",\n \"\\n# - empty: (Well indicator)\", \n \"\\n# - an 'E' indicates that the well was empty and no growth was detected. \",\n \"\\n# - an 'I' indicates that the well was inoculated and growth was detected above the threshold. \",\n \"\\n# - an 'E*' indicates that the well was empty and growth was detected (possible contamination). \",\n \"\\n# - an '!' indicates that the well was inoculated and no growth was detected. \",\n \"\\n# - asymp.not.reached: shows “L” if the bottom asymptote (baseline) was not reached and “U” if the upper asymptote (plateau) was not reached.\",\n \"\\n# - tank: (Tanking indicator) If a number is present then the growth trend was determined to tank at that timepoint index.\", \n \"\\n# - other: Additional flag column. Displays information about whether jumps in OD were detected and what was done about them.\",\n \"\\n# - pdf.file and page.no: location of the figure for this well in the output .pdf files.\"\n )\n \n\t # Analysis information \n \n cat(\"\\n#\\n# -- Source file information--\",\n \"\\n# \", paste(source.file.list, collapse = \"\\n# \"), \n \"\\n# analyzed using GCAT v\", global.version.number, \n \"\\n# request sent: \", upload.timestamp, \n \"\\n# completed: \", analysis.timestamp, \n \"\\n#\\n# -- Parameters used in current analysis --\",\n \"\\n# - Constant added to log(OD + n) transformation:\", add.constant,\n \"\\n# - Blank OD value: \", blank.value,\n \"\\n# - Index of inoculation timepoint\", start.index,\n \"\\n# - Minimum growth threshold:\", growth.cutoff, \n \"\\n# - Removed points:\", paste(points.to.remove, collapse = \" \"),\n \"\\n# - Jump detection:\", remove.jumps) \n sink()\n \n ########################################################################\n # Return values to R #\n ######################################################################## \n # \n\n if(!silent) cat(\"\\ndone!\") \n setwd(old.wd)\n # Return list of generated files\n\treturn(generated.files)\n }\n \n",
- "created" : 1425413277808.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "1165231602",
- "id" : "3DC2D3DA",
- "lastKnownWriteTime" : 1428438949,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/GCAT.main.R",
- "project_path" : "R/GCAT.main.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4DC8D1 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4DC8D1
deleted file mode 100644
index 9eb45ad..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4DC8D1
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# class definition and functions. Objects contain raw #\n# data from screening runs on single wells from 96-well plates, and #\n# other slots for processing and model-fitting details. #\n# #\n########################################################################\n\n# Windows OS compatibility\nSys.setlocale(locale=\"C\")\n#require(RExcelXML)\n\n# Treat nls and loess as S4 classes to avoid warnings\nsetOldClass(\"nls\")\nsetOldClass(\"loess\")\nsetClass(\"well\", representation(position = \"character\",\n\t\t\t\t\t well.info = \"list\",\n\t\t\t\t\t screen.data = \"data.frame\", \n start.index = \"numeric\",\n\t\t\t\t\t use.log = \"logical\",\n\t\t\t\t\t norm = \"numeric\",\n\t\t\t\t\t curve.par = \"list\", \n\t\t\t\t\t fit.par = \"list\",\n fit.std.err = \"list\",\n\t\t\t\t\t equation = \"expression\",\n\t\t\t\t\t model.name = \"character\",\n\t\t\t\t\t fit.info = \"character\",\n\t\t\t\t\t add.info = \"character\",\n inflection.time = \"numeric\",\n rss = \"numeric\",\n loess = \"loess\",\n nls = \"nls\"))\n\n# Slots:\n# position - 3-member vector containing identifying information for the well: row (letters), column (numbers) and plate ID. \n# well.info - a list containing strain and media names if provided\n# screen.data - a data frame with Time and raw OD values. This is the only slot that is filled upon creation of a well object. \n# as different functions are run on the well the data frame gets filled with additional columns. \n# use.log - a single logical value denoting whether to return log-transformed values when data is requested from the well\n# norm - a value to subtract from all OD values before returning data. filled by (see normalize.and.transform.R)\n# curve.par - a list of parameters that denote whether the well is empty, whether it contains ODs indicating a viable culture, whether it tanks at a certain timepoint. \n\n# if model fitting using is successful:\n# fit.par - will be a list containing the fitted model parameters\n# fit.std.err - will be a list containing the standard errors for the fitted model parameters\n# equation - will contain an expression for evaluating the successfully fitted model \n# model.name - will contain the name of the successfully fit model\n\n# fit.info - a message with info about whether the fit was successful, failed, or skipped. \n# add.info - a message with info about whether jumps in OD were detected or removed, or if ODs were detected below the blank OD.\n# inflection.time - the Time value at the point where the specific growth is located. no longer a formula param NWD\n# rss - residual sum of squares\n# loess - object returned by running loess on the normalized well data\n# nls - object returned by running nls on the normalized well data\n\nsetGeneric(\"getPosition\", function(object){standeardGeneric(\"getPosition\")})\nsetMethod(\"getPosition\", \"well\", \n function(object){\n return(object@position)\n })\n\nsetGeneric(\"getWellInfo\", function(object){standeardGeneric(\"getWellInfo\")})\nsetMethod(\"getWellInfo\", \"well\", \n function(object){\n return(object@well.info)\n })\n\nsetGeneric(\"getScreenData\", function(object){standeardGeneric(\"getScreenData\")})\nsetMethod(\"getScreenData\", \"well\", \n function(object){\n return(object@screen.data)\n })\n\nsetGeneric(\"getStartIndex\", function(object){standeardGeneric(\"getStartIndex\")})\nsetMethod(\"getStartIndex\", \"well\", \n function(object){\n return(object@start.index)\n })\n\nsetGeneric(\"getUseLog\", function(object){standeardGeneric(\"getUseLog\")})\nsetMethod(\"getUseLog\", \"well\", \n function(object){\n return(object@use.log)\n })\n\nsetGeneric(\"getNorm\", function(object){standeardGeneric(\"getNorm\")})\nsetMethod(\"getNorm\", \"well\", \n function(object){\n return(object@norm)\n })\n\nsetGeneric(\"getCurPar\", function(object){standeardGeneric(\"getCurPar\")})\nsetMethod(\"getCurPar\", \"well\", \n function(object){\n return(object@curve.par)\n })\n\nsetGeneric(\"getFitErr\", function(object){standeardGeneric(\"getFitErr\")})\nsetMethod(\"getFitErr\", \"well\", \n function(object){\n return(object@fit.std.err)\n })\n\nsetGeneric(\"getEquation\", function(object){standeardGeneric(\"getEquation\")})\nsetMethod(\"getEquation\", \"well\", \n function(object){\n return(object@equation)\n })\n\nsetGeneric(\"getModelName\", function(object){standeardGeneric(\"getModelName\")})\nsetMethod(\"getModelName\", \"well\", \n function(object){\n return(object@model.name)\n })\n\nsetGeneric(\"getFitInfo\", function(object){standeardGeneric(\"getFitInfo\")})\nsetMethod(\"getFitInfo\", \"well\", \n function(object){\n return(object@fit.info)\n })\n\nsetGeneric(\"getAddInfo\", function(object){standeardGeneric(\"getAddInfo\")})\nsetMethod(\"getAddInfo\", \"well\", \n function(object){\n return(object@add.info)\n })\n\nsetGeneric(\"getInflectionTime\", function(object){standeardGeneric(\"getInflectionTime\")})\nsetMethod(\"getInflectionTime\", \"well\", \n function(object){\n return(object@inflection.time)\n })\n\nsetGeneric(\"getRSS\", function(object){standeardGeneric(\"getRSS\")})\nsetMethod(\"getRSS\", \"well\", \n function(object){\n return(object@rss)\n })\n\nsetGeneric(\"getLoess\", function(object){standeardGeneric(\"getLoess\")})\nsetMethod(\"getLoess\", \"well\", \n function(object){\n return(object@loess)\n })\n\nsetGeneric(\"getnls\", function(object){standeardGeneric(\"getnls\")})\nsetMethod(\"getnls\", \"well\", \n function(object){\n return(object@nls)\n })\n\nsetGeneric(\"getFitPar\", function(object){standeardGeneric(\"getFitPar\")})\nsetMethod(\"getFitPar\", \"well\", \n function(object){\n return(object@fit.par)\n })\n\n# --------------------------------------------------------------------\n# Function to create a new well (requires only Time and OD vectors, which will fill slot \"screen.data\")\n# slots \"nls\" and \"loess\" are initialized to empty lists\nwell = function(Time = NULL, OD = NULL){\n x = list()\n class(x) = \"loess\"\n y = list()\n class(y) = \"nls\"\n\tnew(\"well\", screen.data = data.frame(Time, OD, stringsAsFactors=F), loess=x, nls=y)\n}\n\n# -----------------------------------------------------------------------\n#### A show method for well ####\nsetMethod(\"show\", \"well\",\n function(object) {\n print(\"Object of class well\")\n print(\"@position:\")\n print(object@position)\n print(\"@well.info:\")\n print(object@well.info)\n print(\"@screen.data:\")\n print(head(object@screen.data))\n print(\"...\")\n print(paste(nrow(object@screen.data),\"rows of data\"))\n print(paste(\"@start.index:\",object@start.index))\n print(paste(\"@use.log:\",object@use.log))\n print(paste(\"@norm:\",object@norm))\n print(\"@curve.par:\")\n print(object@curve.par)\n print(\"@fit.par:\")\n print(object@fit.par)\n print(\"@fit.std.err:\")\n print(object@fit.std.err)\n print(paste(\"@equation:\",object@equation))\n print(paste(\"@model.name:\",object@model.name))\n print(paste(\"@fit.info:\",object@fit.info))\n print(paste(\"@add.info:\",object@add.info))\n print(paste(\"@inflection.time:\",object@inflection.time))\n print(paste(\"@rss:\",object@rss))\n if (length(object@nls) > 0) {\n print(\"@nls:\")\n print(object@nls)\n } else {\n print(\"no nls model\")\n }\n if (length(object@loess) > 0) {\n print(\"@loess:\")\n print(object@loess)\n } else {\n print(\"no loess model\")\n }\n }\n )\n\n#### A plot method for well ####\n# x - object of class well\n# y - not used\n# constant.added - used to readjust for the constant added during the log transform: log.OD = log(OD - blank + constant.added)\n# xlim - x axis limits, vector of length 2\n# ylim - y axis limits, vector of length 2\n# scale - determines the font scale for the entire graph. all cex values are calculated from this\n# number.points - should points be labeled with numeric indices?\n# draw.symbols - should be called on the well and markings drawn on the graph?\n# show.text - show R^2 and growth curve parameters as text on the plot\n# show.calc - draw lines that illustrate growth curve parameters\n# draw.guess - initial guess model. Drawn if specified\n# well.number - the number of the well in an array of wells\n# ... - additional arguments passed to the generic plot function\n\nsetMethod(\"plot\",\n signature(x = \"well\", y=\"missing\"),\n function (x, y, constant.added = 1.0, xlim = NULL, ylim = NULL,\n well.number = NULL, scale = 1, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, ...) \n {\n # Determine the boundaries for the axes (if user did not specify them)\n if(is.null(ylim)){\n min.y = min(data.from(x, remove = F, remove.tanking = F)[,2], na.rm = T)\n min.y = min(min.y, x@fit.par$b)\n max.y = max(data.from(x, remove = F, remove.tanking = F)[,2], na.rm = T)\n max.y = max(max.y, x@fit.par$b + x@fit.par$A)\n ylim = c(min.y, min.y + (max.y-min.y)*1.15)\n }\n if(is.null(xlim)){\n min.x = min(data.from(x, remove = F, remove.tanking = F)[,1], na.rm = T)\n max.x = max(data.from(x, remove = F, remove.tanking = F)[,1], na.rm = T)\n xlim = c(min.x - 0.05 * (max.x-min.x), max.x)\n }\n \n \n # Title of plot: [well number] plate name; well name;\n # strain name; media name\n \n main = paste(plate.name(x), \" \", well.name(x), \"\\n\",\n strain.name(x), \"; \", media.name(x), sep = \"\")\n if (!is.null(well.number)) main = paste(\"[\", well.number , \"] \", main, sep=\"\")\n \n # Draw the data and symbols if is true.\n plot.data(x, main = main, scale = scale, constant.added=constant.added, \n number.points = number.points, draw.symbols = draw.symbols, xlim = xlim, ylim = ylim, ...)\n \n # Draw the fitted model.\n plot.model(x, scale = scale, constant.added=constant.added)\n \n # Draw text info if specified. \n if(show.text)\n draw.text(x, scale = scale * 0.5, xlim = xlim, ylim = ylim,...)\n \n # Show calculated parameters if specified. \n if (show.calc)\n draw.calc.par(x, scale = scale * 0.5, constant.added = constant.added)\n \n # Draw initial guess if a model is specified. \n if (class(draw.guess) == \"model\"){\n Time = data.from(x)$Time\n guess = eval(getExpression(draw.guess), as.list(getGuess(draw.guess)(x)))\n try(lines(Time, guess, col = \"brown2\"), silent = T)\n }\n }\n)\n\n########################################################################\n# Some miscellaneous functions to extract info from well objects #\n# Most of these return a single value from the well. #\n########################################################################\n#\n# Since many of these need to be applied to all wells over an array, while conserving the dimensions of \n# that array, this file includes a wrapper function (see bottom of file).\n\nplate.name = function(well)\n\tgetPosition(well)[1]\n\n# Return the full alphanumeric well name (with leading zeros if applicable)\nwell.name = function(well){\n\trow = getPosition(well)[2]\n\tcol = as.numeric(getPosition(well)[3])\n\tif (col>9)\n\t\tcol = as.character(col)\n\telse\n\t\tcol = paste(\"0\", col, sep = \"\")\n\n\tpaste(row,col,sep = \"\")\n\t}\n\nis.empty = function(well)\n\tgetCurPar(well)$empty.well\n\nlacks.growth = function(well)\n\tgetCurPar(well)$no.growth\n\ntanking.start = function(well)\n\tgetCurPar(well)$tanking.start\n\nremoved.points = function(well)\n\t(1:length(well))[getScreenData(well)$Remove]\n\nremaining.points = function(well,...){\n\tas.numeric(rownames(data.from(well,...)))\n\t}\n\nstrain.name = function(well){\n if(is.null(getWellInfo(well)$Strain))\n return(\"\")\n else\n return(getWellInfo(well)$Strain)\n }\nmedia.name = function(well){\n if(is.null(getWellInfo(well)$Media))\n return(\"\")\n else\n return(getWellInfo(well)$Media)\n }\n\nraw.data = function(well)\n\tdata.from(well, remove.tanking = F, remove = F, na.rm = F, raw.data = T)\n\ncontains.fit = function(well)\n\tlength(getFitPar(well)) > 0\n\nsetMethod(\"length\", signature(x = \"well\"), function(x) length(x@screen.data[,1]))\n\n# The function has some options: by default it returns a two-column data frame with time and OD \n# (or log OD if the slot is true in the object), after normalization to the value specified in slot. \n# - With set to true the rows specified in the column of the slot are not returned. \n# - With set to true all the rows after the index are removed. \n# - Setting to true overrides all these settings and just returns 2 columns with Time and Raw OD.\n\ndata.from = function(well, remove = T, remove.tanking = T, raw.data = F, na.rm = F){\n\t\n\tif (length(getUseLog(well)) == 0)\n\t\tOD.column = \"OD\"\n\telse if (getUseLog(well))\n\t\tOD.column = \"log.OD\"\n\telse\n\t\tOD.column = \"OD\"\n\t\n\tif (raw.data){\n\t\tOD.column = \"OD\"\n\t\tnorm = 0\n\t\t}\n\telse if (!getUseLog(well))\n\t\tnorm = getNorm(well)\n\telse\n\t\tnorm = 0\n\n\tif(remove.tanking & is.numeric(tanking.start(well)))\n\t\twell = remove.points(well, (tanking.start(well)):length(well))\n\tif (!remove | is.null(getScreenData(well)$Remove))\n\t\toutput = getScreenData(well)[c(\"Time\", OD.column)]\n\telse\n\t\toutput = getScreenData(well)[!getScreenData(well)$Remove ,c(\"Time\", OD.column)]\n\n\toutput[,2] = output[,2] - norm\n\n\tif (!raw.data){\n\t\tif (!length(getUseLog(well)))\n\t\t\tnames(output)[2] = \"Corrected.OD\"\n\t\tif (!getUseLog(well))\n\t\t\tnames(output)[2] = \"Corrected.OD\"\n\t\t}\n\n\tif (na.rm)\n\t\toutput[!is.na(output[,2]),]\t \n\telse\n\t\toutput\n\t}\n\n\n# Functions much like but gives a single vector containing the \n# slope at each point. Has a parameter allowing removal of NA values. \n\nslopes = function(well, remove = T, remove.tanking = T, na.rm = F){\n\n\tif(remove.tanking & is.numeric(tanking.start(well)))\n\t\twell = remove.points(well, (tanking.start(well)):length(well))\n\tif (!remove | is.null(getScreenData(well)$Remove))\n\t\toutput = getScreenData(well)$Slope\n\telse\n\t\toutput = getScreenData(well)$Slope[!getScreenData(well)$Remove]\n\n\tif (na.rm)\n\t\toutput[!is.na(output)]\t \n\telse\n\t\toutput\n\t}\n\n# -----------------------------------------------------------------------\n# Well array functions: these must be used on entire arrays of well objects\n# instead of single ones. \n\nplate.names = function(well.array)\n\tdimnames(well.array)[[3]]\n\ntanking.start.values = function(well.array, array = F){\n\tif (array)\n\t\taapply(well.array, function(well) tanking.start(well))\n\telse\n\t\tsapply(well.array, function(well) tanking.start(well))\n\t}\n\n",
- "created" : 1425413290697.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2437520955",
- "id" : "4DC8D1",
- "lastKnownWriteTime" : 1425509121,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/class.well.R",
- "project_path" : "R/class.well.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4E20D973 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4E20D973
deleted file mode 100644
index ee62373..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/4E20D973
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n# Wrapper for sapply to use lapply over an array, conserving the dimensions.\naapply = function(x, FUN,...){\n dim.values = dim(x)\n\tdim.names = dimnames(x)\n\tx = lapply(x, function(x){FUN(x,...)})\n\tdim(x) = dim.values\n\tdimnames(x) = dim.names\n\treturn(x)\n\t}\n\n# A function to manually create an unchecked exception.\nexception = function(class, msg)\n{\n cond <- simpleError(msg)\n class(cond) <- c(class, \"MyException\", class(cond))\n stop(cond)\n}\n",
- "created" : 1425413273698.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "3044086993",
- "id" : "4E20D973",
- "lastKnownWriteTime" : 1424208963,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/misc.R",
- "project_path" : "R/misc.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/6FF7DC6F b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/6FF7DC6F
deleted file mode 100644
index a243cea..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/6FF7DC6F
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# Populate an output table with parameters and other useful info for #\n# each well in a fitted dataset. #\n# #\n########################################################################\n#\n# unlog - Should OD values be returned on the linear scale instead of log-transformed scale? \n# constant.added - For returning values on linear scale, what constant was added to ODs before the log transform? \n# reach.cutoff - what proportion of the plateau OD must tbe reached by the last valid timepoint for the curve to be marked as reaching its plateau OD?\n#\n\ntable.out = function(fitted.data.set, unlog = F, constant.added, reach.cutoff = 0.90, filename.timestamp = NULL,use.linear.param=F, use.loess=F){\n \n # The idea is basically to use and on the fitted data array in order \n # to get one vector for each column of the output table. \n \n # Get identifying information (plate, well, media and strain names)\n plate.ID = unlist(aapply(fitted.data.set,plate.name))\n well.ID = unlist(aapply(fitted.data.set,well.name))\n media.ID = unlist(aapply(fitted.data.set,media.name))\n strain.ID = unlist(aapply(fitted.data.set,strain.name))\n # Get fit information for each well\n # - was it marked as empty in the plate layout?\n # - did the program find it to contain no growth (\"dead\")? \n # - was the fitting procedure successful? \n # - did the curve tank? if so, at what timepoint? if not, set value to \"-\"\n \n empty = unlist(aapply(fitted.data.set, is.empty))\n dead = unlist(aapply(fitted.data.set, lacks.growth))\n fit = unlist(aapply(fitted.data.set, contains.fit))\n tanking = unlist(aapply(fitted.data.set, tanking.start))\n tanking[is.na(tanking) | tanking == 1 | dead] = \"-\"\n \n # Get calculated values for each well: specific growth, final and initial OD, fitted plateau and baseline OD, lag time, etc.\n inflection.time = unlist(aapply(fitted.data.set, inflection.time))\n max.spec.growth.rate = unlist(aapply(fitted.data.set, max.spec.growth.rate))\n max.log.OD = unlist(aapply(fitted.data.set, max.log.OD))\n inoc.log.OD = unlist(aapply(fitted.data.set, inoc.log.OD))\n projected.growth = unlist(aapply(fitted.data.set, projected.growth))\n projected.growth.OD = unlist(aapply(fitted.data.set, projected.growth.OD, constant.added))\n achieved.growth = unlist(aapply(fitted.data.set, achieved.growth))\n achieved.growth.OD = unlist(aapply(fitted.data.set, achieved.growth.OD, constant.added))\n lag.time = unlist(aapply(fitted.data.set, lag.time))\n shape.par = unlist(aapply(fitted.data.set, shape.par))\n RSS = unlist(aapply(fitted.data.set, rss))\n baseline = unlist(aapply(fitted.data.set, baseline))\n amplitude = unlist(aapply(fitted.data.set, amplitude))\n plateau = unlist(aapply(fitted.data.set, plateau))\n ########################3h#############################################\n max.spec.growth.rate.SE = unlist(aapply(fitted.data.set, max.spec.growth.rate.SE))\n shape.par.SE = unlist(aapply(fitted.data.set, shape.par.SE))\n lag.time.SE = unlist(aapply(fitted.data.set, lag.time.SE))\n amplitude.SE = unlist(aapply(fitted.data.set, amplitude.SE)) # a.k.a amplitude error\n baseline.SE = unlist(aapply(fitted.data.set, baseline.SE)) # a.k.a baseline error\n #######################################################################\n \n # If the curve falls short of 90% of plateau OD by the final timepoint.\n no.reach.plateau = !unlist(aapply(fitted.data.set, reach.plateau, cutoff = 0.9))\n # If the fitted baseline is below zero on linear scale\n no.reach.baseline = unlog(baseline,constant.added) < 0\n \n # If any of these are NA as a result of failed fits, change them to false: they don't need to be reported. \n no.reach.plateau[is.na(no.reach.plateau)] = F\n no.reach.baseline[is.na(no.reach.baseline)] = F\n # What percent of the total growth does the curve actually reach? \n # (in case of total growth being 0, change this to 100%)\n percent.reach = 100*((max.log.OD - inoc.log.OD) / (projected.growth))\n percent.reach[is.infinite(percent.reach)] = 100\n \n # Return the name of the model (if any) that was successfully fit to the well. \n model.used = unlist(aapply(fitted.data.set, function(well)well@model.name))\n \n # \"Goodness of fit\" metric\n good.fit = unlist(aapply(fitted.data.set, model.good.fit))\n \n # Code the two flags: \n flag1 = flag2 = rep(\"-\", length(tanking))\n \n for(i in 1:length(tanking)){\t\n # Flag 1 (empty/inoculated flag) possible values:\n # well was empty and no growth was found (E)\n # well was empty, but growth was found (E*)\n # well was inoculated but no growth was found (!)\n # well was inoculated and growth was found (I)\n \n if(empty[i] & !fit[i])\n flag1[i] = \"E \"\n if(empty[i] & fit[i])\n flag1[i] = \"E*\"\n if(!empty[i] & dead[i])\n flag1[i] = \"! \"\n if(!empty[i] & !dead[i])\n flag1[i] = \"I \"\n \n # Flag 2 (lower/upper asymptotes) possible values:\n # well did not reach lower asymptote (baseline OD) (L)\n # well did not reach upper asymptote (plateau OD) (U)\n # well did not reach either asymptote (L/U)\n # well reached both asymptotes (-)\n \n if(no.reach.baseline[i]){\n if (no.reach.plateau[i])\n flag2[i] = \"L/U\"\n else\n flag2[i] = \"L\"\n }\n else{\n if (no.reach.plateau[i])\n flag2[i] = \"U\"\n else\n flag2[i] = \"-\"\n }\n # Also use the and and to provie more info about why model fitting failed in some cases. \n if(dead[i])\n model.used[i] = \": skipped\"\n else if(!empty[i] & !fit[i])\n model.used[i] = \": failed\"\t\n }\n \n # Flag 3: return the additional info slot. \n flag3 = unlist(aapply(fitted.data.set, function(well){\n if (length(well@add.info) > 0) \n return(well@add.info)\n else\n return(\"\")\n }))\n \n # If something is amiss with the data table use this to check on the arguments...\n #cat(\"plate \", length(plate.ID),\" well \", length(well.ID),\" media \", length(media.ID),\" strain \", length(strain.ID),\n #\" model \", length(model.used),\" max.spec.growth.rate\", length(max.spec.growth.rate), \"projected.growth\", length(projected.growth),\n #\"lag.time\", length(lag.time), \"inoc.log.OD\", length(inoc.log.OD), \"good.fit\",\n #length(good.fit),\"empty\", length(flag1),\"asymp\", length(flag2),\" tank \", length(tanking),\" reach \", length(percent.reach),\" other \", length(flag3), sep = \"\\n\")\n # 06.28.11: Add a row number identifier for output perusal\n row.number = 1:length(plate.ID)\n \n pdf.file = page.no = c()\n # 06.29.11: Add pdf file name and page number references. Prepare timestamp for addition to output file names (for file references in last column)\n for(i in 1:length(plate.ID)){\n pdf.file[i] = paste(plate.ID[i], \"_plots\", filename.timestamp, \".pdf\", sep=\"\")\n page.no[i] = (i-1) %% 96 + 2\n }\n # Slap it all together into a data frame.\n if(use.loess){\n output.core = data.frame(row = row.number, plate = plate.ID, well = well.ID, media = media.ID, strain = strain.ID, \n model = model.used, lag.time, inflection.time, max.spec.growth.rate, \n baseline, amplitude, plateau, inoc.log.OD, max.log.OD, achieved.growth,\n baseline.OD = unlog(baseline,constant.added), amplitude.OD = unlog(amplitude,constant.added), \n plateau.OD = unlog(plateau,constant.added), inoc.OD = unlog(inoc.log.OD,constant.added), \n max.OD = unlog(max.log.OD,constant.added), achieved.growth.OD = achieved.growth.OD,\n R.squared = good.fit, RSS = RSS, empty = flag1, asymp.not.reached = flag2, tank = tanking, other = flag3, pdf.file = pdf.file, page.no = page.no)\n } else {\n output.core = data.frame(row = row.number, plate = plate.ID, well = well.ID, media = media.ID, strain = strain.ID, \n model = model.used, lag.time = lag.time, lag.time.SE, inflection.time, max.spec.growth.rate, max.spec.growth.rate.SE, \n baseline, baseline.SE, amplitude, amplitude.SE, plateau, inoc.log.OD, max.log.OD, projected.growth, achieved.growth,\n baseline.OD = unlog(baseline,constant.added), amplitude.OD = unlog(amplitude,constant.added), \n plateau.OD = unlog(plateau,constant.added), inoc.OD = unlog(inoc.log.OD,constant.added), \n max.OD = unlog(max.log.OD,constant.added), projected.growth.OD = projected.growth.OD, achieved.growth.OD = achieved.growth.OD,\n shape.par = shape.par, shape.par.SE,\n R.squared = good.fit, RSS = RSS, empty = flag1, asymp.not.reached = flag2, tank = tanking, other = flag3, pdf.file = pdf.file, page.no = page.no)\n }\n \n # Add units to column names\n names2 = names(output.core)\n names2[grep(\"time\",names2)] = sub(\"$\",\", hrs\", names2[grep(\"time\",names2)])\n names2[grep(\"rate\",names2)] = sub(\"$\",\", log.OD/hr\", names2[grep(\"rate\",names2)])\n log.OD.fields = c(\"baseline\", \"baseline.SE\", \"amplitude\", \"amplitude.SE\", \"plateau\", \"projected.growth\", \"achieved.growth\")\n names2[names2 %in% log.OD.fields] = sub(\"$\", \", log.OD\", names2[names2 %in% log.OD.fields])\n names(output.core) = names2\n \n # Add on any additional fields found in the plate layout. \n all.layout.fields = sapply(fitted.data.set, function(well) unlist(well@well.info)) \n all.layout.fields = as.data.frame(t(all.layout.fields))\n \n \n addl.info = all.layout.fields[,!(names(all.layout.fields) %in% c(\"Strain\", \"Media\"))]\n if(!is.data.frame(addl.info)){\n addl.info = data.frame(addl.info)\n names(addl.info) = names(all.layout.fields)[!(names(all.layout.fields) %in% c(\"Strain\", \"Media\"))] \n }\n \n output = cbind(output.core,addl.info)\n \n return(output)\n}\n\n\n\n\n\n\n",
- "created" : 1425413247086.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2376564637",
- "id" : "6FF7DC6F",
- "lastKnownWriteTime" : 1424208623,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/table.output.R",
- "project_path" : "R/table.output.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/74C1CB00 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/74C1CB00
deleted file mode 100644
index bcbbd2a..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/74C1CB00
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# Normalize OD readings for an entire array of well objects #\n# #\n########################################################################\n#\n# Note: This function does not write any new OD values to the well objects in the array - it only \n# fills the \"norm\" slot of each well object in the array with a value that will be subtracted \n# from all OD measurements when returning data from the wells using the function (see well.class.R) \n#\n# These functions make use of which simply returns the raw time and OD of a well (also see well.class.R)\n#\n# well.array: an array of well objects. note this is the only normalization function that acts on an entire array instead of an individual well.\n# normalize.method: \n# - (default): subtracts the blank OD (either specified by or taken from the first timepoint as default) of each well from all timepoints in that well\n# - average.blank: subtracts the mean of all first OD timepoints on a plate from all timepoints in all wells on that plate\n# - average.first: takes the mean of the difference between the OD of the specified timepoint and the first timepoint of all wells on a plate\n# and subtracts this value from all timepoints in all wells on that plate\n# - anything else: do nothing\n# blank.value - user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well. \n# start.index - which timepoint should be used as the first one after inoculation (defaults to the 2th one)\n# add.constant: add a numeric constant to all timepoints in all wells. \nnormalize.ODs = function(well.array, normalize.method = \"default\", blank.value = NULL, start.index = 2, add.constant = 1){\n \n if (normalize.method == \"default\"){\n well.array = aapply(well.array, function(well, blank.value){\n # Use the blank OD value if specified; otherwise, get it from the first OD timepoint. \n if(is.null(blank.value)) blank.value = raw.data(well)[1,2]\n # Set the blank OD (minus the constant to be added) to the \"norm\" slot of each well.\n \twell@norm = blank.value - add.constant\n \treturn(well)}, blank.value)\n }\n else if (normalize.method == \"average.blank\"){ \n # Use the blank OD value if specified; otherwise, get it from the first OD timepoint.\n\t\tblank.ODs = unlist(aapply(well.array, function(well, blank.value){\n if(is.null(blank.value)) blank.value = raw.data(well)[1,2]\n return(blank.value)}, blank.value))\n\t\tplate.IDs = unlist(aapply(well.array, plate.name))\n\t\tblank.averages = tapply(blank.ODs, plate.IDs, mean)\n\t # Set this value (minus the constant to be added) to the \"norm\" slot of each well. \n\t\twell.array = aapply(well.array, function(well){\n\t\t\twell@norm = blank.averages[plate.name(well)] - add.constant\n\t\t\treturn(well)})\n\t\t}\n\telse if (normalize.method == \"average.first\"){\n\t # Find the mean difference between starting OD (timepoint specified by ) and blank OD (first timepoint) for each plate\n # Use the blank OD value if specified; otherwise, get it from the first OD timepoint.\n blank.ODs = unlist(aapply(well.array, function(well, blank.value){\n if(is.null(blank.value)) blank.value = raw.data(well)[1,2]\n return(blank.value)}, blank.value))\n\t\tfirst.ODs = unlist(aapply(well.array, function(well) raw.data(well)[start.index,2]))\n\t\tplate.IDs = unlist(aapply(well.array, plate.name))\n\t blank.averages = tapply(first.ODs-blank.ODs,plate.IDs,mean)\n\t # Set this value (minus the constant to be added) to the \"norm\" slot of each well.\n\t\twell.array = aapply(well.array, function(well){\n\t\t\twell@norm = raw.data(well)[start,2] - blank.averages[plate.name(well)] - add.constant \n\t\t\treturn(well)})\n\t\t}\t\n\telse{\n # Simply set the negative constant to be added to the \"norm\" slot of each well. \n\t\twell.array = aapply(well.array, function(well){\n\t\t\twell@norm = - add.constant\n \t\t\treturn(well)})\n\t\t}\n if(is.null(blank.value))\n well.array = aapply(well.array, remove.points, 1)\n return(well.array)\n\t}\n\n########################################################################\n# #\n# Log-transform OD readings for a single well object #\n# #\n########################################################################\n\n# Must include this so that the checking process will not complain about\n# inconsistency S3 generic/method. Though I don't know why.\ntransform <- function(input.well, ...) {\n UseMethod(\"transform\")\n}\n\n#' Transform.Ods \n#' \n#' This function adds a \"log.OD\" column to the \"screen.data\" slot of a well object with log-transformed data. \n#' The raw data is kept intact. \n#' It also checks to see if any of the raw OD values (before a certain timepoint) is below the blank OD. \n#' This can be disastrous for the log(OD) transform. \n#' @param input.well an object of class well \n#' @param use.log gets added to the \"use.log\" slot of the well object. this will determine whether the log-transformed data \n#' or raw normalized data is returned using the function \\code{data.from}. \n#' @param blank.value user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well. \n#' @param start.index which timepoint should be used as the first one after inoculation (defaults to the 2th one) \n#' @param negative.OD.cutoff if any ODs below the specified blank value are detected before this index timepoint, the entire well is discarded. \ntransform.ODs = function(input.well, use.log = T, blank.value = NULL, start.index = 2, negative.OD.cutoff = 10, constant.added = 1.0, ...){\n \n # The default value for the log-transformed ODs will be NA. Valid values will be filled in. \n\tlog.OD = rep(NA, length(input.well))\n OD = raw.data(input.well)[,2]\n \n\t # Use the blank OD value if specified; otherwise, get it from the first OD timepoint.\n if(is.null(blank.value))\n blank.value = OD[1]\n\n # Remove any points from the analysis that weren't already removed and fall below the blank value (using below)\n OD[input.well@screen.data$Remove] = NA\n \n negative.points = which(OD + 0.2 * constant.added < blank.value)\n if(length(negative.points) > 0)\n\t input.well = remove.points(input.well, negative.points)\n \n # If any points fall below the blank value by more than 0.2 * and before the cutoff index , remove the well from analysis. \n # First adjust the cutoff to compensate for curves that don't start at timepoint 1\n negative.OD.cutoff = negative.OD.cutoff + start.index - 1\n \n if(any(negative.points <= negative.OD.cutoff)){\n input.well = remove.points(input.well, rep(T,length(input.well)))\n input.well@add.info = paste(\"ODs at timepoint(s)\", paste(negative.points[negative.points <= negative.OD.cutoff],collapse=\" \"), \"were below blank OD; well discarded\")\n }\n\n # Take the natural log of the rest of the OD values (after subtracting the normalization value)\n log.OD[which(OD > input.well@norm)] = log(OD[which(OD > input.well@norm)] - input.well@norm)\n\t\n\t# Add a column to the \"screen.data\" slot of the well\n\tinput.well@screen.data$log.OD = log.OD\t\n\t# Update the \"use.log\" slot of the well \n\tinput.well@use.log = use.log\t\n\n\treturn(input.well)\n\t}\n\n########################################################################\n# #\n# Remove timepoints from the analysis but not from the raw data #\n# #\n########################################################################\n# \n# Removes timepoints from further analysis. Does not remove them from the raw data;\n# instead, this function creates or updates the Remove column in slot \"screen.data\" of the well which dictates whether \n# individual timepoints are returned using the function. \n#\n# can be a vector containing:\n# - any combination of positive and negative integers \n# the timepoints at indices corresponding to positive integers will be set to be removed.\n# the timepoints at indices corresponding to negative integers will be be re-added if they were previously set to be removed.\n# - a single zero, which resets all timepoints (nothing will be removed)\n# - a logical vector to replace the Remove column and which will be cycled along the length of the timepoints. \n\nremove.points = function(input.well, points){\n # Copy the Remove column or create a new one if it doesn't yet exist\n\tif (is.null(input.well@screen.data$Remove))\n\t\tRemove = rep(F, length(input.well))\n\telse\n\t\tRemove = input.well@screen.data$Remove\n\n # If is a logical vector, recycle it along the length of Remove \n\tif (length(points[!is.na(points)]) != 0){\n\t# Separate positive and negative integers\n \tif (is.logical(points) & !any(is.na(points)))\n \t\tRemove = rep(points,length.out=nrow(input.well@screen.data))\t\n \t\telse{\n \tpos = points[points > 0]\n \tneg = -points[points < 0]\n \tRemove[pos] = T\n \tRemove[neg] = F \n \t\t \t\n # If contains only zeros, reset the Remove vector to all F \t\n \tif (all(points == 0))\n \t\tRemove[1:length(Remove)] = F\n \t\t}\n }\n # replace the Remove column\n\tinput.well@screen.data$Remove = Remove\n\tinput.well\n\t}\n\n\n\n\n\n",
- "created" : 1425413265888.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "3349576883",
- "id" : "74C1CB00",
- "lastKnownWriteTime" : 1428436335,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/normalize.and.transform.R",
- "project_path" : "R/normalize.and.transform.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/907B51DB b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/907B51DB
deleted file mode 100644
index 6ee699e..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/907B51DB
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "# Default NAMESPACE created by R\n# Remove the previous line if you edit this file\n\n# Export all names\nexportPattern(\".\")\n\n# Import all packages listed as Imports or Depends\nimport(\"pheatmap\", \"gplots\")\n",
- "created" : 1428439772005.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "749067415",
- "id" : "907B51DB",
- "lastKnownWriteTime" : 1423260997,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/NAMESPACE",
- "project_path" : "NAMESPACE",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_namespace"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/992E3490 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/992E3490
deleted file mode 100644
index f8ed693..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/992E3490
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead,\n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n########################################################################\n# #\n# class definition and functions. Objects contain equations #\n# and other information for parameterized growth curve models. #\n# #\n########################################################################\nsetClass(\"model\", representation(name = \"character\",\n expression = \"expression\",\n formula = \"formula\",\n guess = \"function\"))\n# Slots:\n# name - a simple description of the model.\n# expression - an object of class \"expression\" that evaluates the response (transformed OD) with respect to the variable Time.\n# formula - same as expression, but with y as the response.\n# guess - a function that computes initial guesses for the parameters given a well object with a valid \"screen.data\" slot\n# containing useable OD values and slope estimates\n# --------------------------------------------------------------------\n###################### BEGIN PROTOTYPING ACCESSOR METHODS##############\n\n# Minh: Let this code fragment be F1.\nif (!isGeneric(\"getName\")){\n if (is.function(\"getName\"))\n fun <- getName\n else\n fun <- function(object) standardGeneric(\"getName\")\n setGeneric(\"getName\", fun)\n}\n# End of F1\nsetMethod(\"getName\", \"model\", function(object) object@name)\n\n# Minh: Let this line be F2.\nsetGeneric(\"getExpression\", function(object){standardGeneric(\"getExpression\")})\n# Question: How is F1 different from F2?\n\nsetMethod(\"getExpression\", \"model\",\n function(object){\n return(object@expression)\n })\n\nsetGeneric(\"getFormula\", function(object){standeardGeneric(\"getFormula\")})\nsetMethod(\"getFormula\", \"model\", \n function(object){\n return(object@formula)\n })\n\nsetGeneric(\"getGuess\", function(object){standeardGeneric(\"getGuess\")})\nsetMethod(\"getGuess\", \"model\", \n function(object){\n return(object@guess)\n })\n######################## ENG PROTOTYPING ########################\n\n# Function to create a new model\n#' Model \n#' \n#' Function to create a new model \n#' @param name The name of the model \n#' @param expression Expression of the model \n#' @param formula The formula of this model \n#' @param guess The guess of this model \n#' @return The new model \nmodel = function(name, expression, formula, guess){\n new(\"model\", name = name, expression = expression, formula = formula, guess = guess)\n}\n\nloess.g = function(well,smooth.param=0.75){\n #data = data.from(well)\n #growth = data[,2]\n #Time = data[,1]\n Time = data.from(well)[,1]\n \n # predicted growth values to be used in estimating growth curve parameters\n loess.fit = loess(data.from(well)[,2]~Time,span=smooth.param)\n t = seq(from = min(Time), to = max(Time), by = (max(Time)-min(Time))/1000)\n y = predict(loess.fit, data.frame(Time=t))\n attr(y,\"names\") = NULL # need to remove the names to prevent them from showing up in the returned vector\n \n # Remove any data points where y has not been estimated\n filt = is.finite(y)\n t = t[filt]\n y = y[filt] # remove any NA etc\n \n # specific growth using loess to find max derivative\n delta.t = diff(t)\n dydt = diff(y)/delta.t\n u = max(dydt)\n \n # lower and upper asymptotes\n b = min(y)\n A = max(y) - min(y)\n \n # inflection point\n inflection.pt.index = which.max(dydt)\n inflection.time = t[inflection.pt.index]\n inflection.y = y[inflection.pt.index]\n \n # lag time\n lam = inflection.time - (inflection.y-b)/u\n \n # Return named array of estimates\n c(A = A, b = b, lam = lam, u = u)\n}\n\n\n",
- "created" : 1425419995370.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2225975310",
- "id" : "992E3490",
- "lastKnownWriteTime" : 1425511077,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/class.model.R",
- "project_path" : "R/class.model.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/B0B7416E b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/B0B7416E
deleted file mode 100644
index 30cb003..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/B0B7416E
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead,\n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n########################################################################\n# #\n# Fit a parameterized model to the growth data in a well object. #\n# #\n# There are now three modelling choices: #\n# 1) Sigmoid model (no linear param c) #\n# 2) Linear Sigmoid model #\n# 3) Loess (with optional smoothing parameter) #\n########################################################################\n#' fit.model \n#' \n#' This function will use the function stored in the \"guess\" slot of \\code{growth.model} to calculate initial guesses \n#' for growth.model parameters, then it will use the \"formula\" slot with \\code{nls} to fit a non-linear least squares \n#' \\code{growth.model} or Local Polynomial Regression Fitting to the data. Richards model is first fitted. \n#' If the shape parameter is statisticaly significant then Richards is used. If it is within 2 SE of 1 or Zero than \n#' a simpler model is preferred. If the Richards fit fails, then Logistic is tried. If it fails, Gompertz is tried. \n#' Model fit failure is reported if none of the models can sucessfully fit the data \n#' \n#' @param input.well The well needed to be fitted with the given model. \n#' @param growth.model What growth model should be used? \n#' @param backup.growth.model If \\code{gowth.mode} fails, this model will be used. \n#' @param fit.if.no.growth should the function attempt to fit a well even if there was no growth detected? default is F \n#' @param silent output back to R console? \n#' @param use.linear.param: Should an additional linear parameter (c) be used when fitting the data to the model? \n#' @param use.loess: Should Local Polynomial Regression Fitting (loess function) be used instead of nls? \n#' @param smooth.param: If loess is used, an optional smoothing parameter. Default is .6 \nfit.model = function(input.well, growth.model, backup.growth.model = NULL, fit.if.no.growth = F, \n use.linear.param=F, use.loess=F, smooth.param, silent = T){\n \n # Conditional breakpoint: go into debugging mode when fitting a specific well\n #if (input.well@position[\"row\"] == \"A\" && input.well@position[\"col\"] == \"12\") browser()\n \n # Change all relevant slots to or blank values\n input.well@model.name = \"\"\n input.well@fit.par = list()\n input.well@equation = expression()\n \n # Get OD vs. time data from well\n input.data = data.from(input.well, na.rm = T)\n \n # Skip well if in slot \"curve.par\" is set to true, and is false.\n if(!fit.if.no.growth & lacks.growth(input.well)){\n input.well@fit.info = \"skipped - no growth in well.\"\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n }\n # Skip well if there are fewer than 5 data points left in the analysis.\n if (length(input.data$Time) < 5){\n input.well@fit.info = \"skipped - not enough points.\"\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n }\n \n # Change column headers of input.data to the more general \"Time\" vs. \"y\"\n names(input.data) = c(\"Time\", \"y\")\n \n # Set a lower bound for nls model parameters A and b to slightly lower than min(y)\n low.y = min(input.data$y,na.rm=T)\n low.y = low.y - 0.1*abs(low.y)\n \n # Extract the model formula from (slot \"formula\")\n # Use the function from slot \"guess\" to calculate initial guesses for model parameters based on slope estimates in \n # Attempt to fit a nonlinear least squares odel using \n \n # Creating loess, logistics, richards, and gompertz model. \n loess.e = expression(\"loess\")\n loess.f = formula(Time ~ y)\n loess.model = model(\"local polynomial regression fit.\", loess.e, loess.f, loess.g)\n \n ### Testing accessor method for class model.\n #print(getName(loess.model))\n #print(getFormula(loess.model))\n #print(getGuess(loess.model))\n ### End testing ###\n \n remove(loess.e, loess.f)\n \n ########################################################################\n # Create the logistic 4-parameter model (when v ~ 1) #\n ########################################################################\n logistic.g = function(well,smooth.param=0.75) {\n loess.model@guess(well,smooth.param)\n }\n \n ########################################################################\n # Create the Richards 5-parameter model #\n ########################################################################\n richards.g = function(well,smooth.param=0.75){\n c(loess.model@guess(well,smooth.param),v=0.5)\n }\n \n ########################################################################\n # Create the Gompertz model (might be useful as a #\n # limiting case of Richards model when v ~ 0) #\n ########################################################################\n gompertz.g = function(well,smooth.param=0.75){\n loess.model@guess(well,smooth.param)\n }\n \n logistic.e = expression((A/(1+exp((4*u/A)*(lam-Time)+2)))+b)\n logistic.f = formula(y~(A/(1+exp((4*u/A)*(lam-Time)+2)))+b)\n logistic = model(\"logistic sigmoid.\", logistic.e, logistic.f, logistic.g)\n remove(logistic.e, logistic.f, logistic.g)\n \n richards.e = expression(A*(1+v*exp(1+v)*exp((u/A)*(1+v)**(1+1/v)*(lam-Time)))**(-1/v)+b)\n richards.f = formula(y~A*(1+v*exp(1+v)*exp((u/A)*(1+v)**(1+1/v)*(lam-Time)))**(-1/v)+b)\n richards = model(\"richards sigmoid\", richards.e, richards.f, richards.g)\n remove(richards.e, richards.f, richards.g)\n \n gompertz.e = expression(A*exp(-exp((u*exp(1)/A)*(lam-Time)+1))+b)\n gompertz.f = formula(y~A*exp(-exp((u*exp(1)/A)*(lam-Time)+1))+b)\n gompertz = model(\"gompertz sigmoid\", gompertz.e, gompertz.f, gompertz.g)\n remove(gompertz.e, gompertz.f, gompertz.g)\n \n \n \n # 3) Loess (with optional smoothing parameter)\n if(use.loess){\n number.of.points = nrow(input.well@screen.data)\n if (smooth.param <= 1/number.of.points)\n exception(\"Invalid input\", \"Smoothing parameter is out of range.\")\n \n fit = try(loess(y~Time, data=input.data, span=smooth.param), silent=TRUE)\n input.well@loess = fit\n if (class(fit) != \"loess\") stop(\"loess fit failed on well\", paste(input.well@position,collapse=\" \"))\n input.well@fit.info = \"Loess model fit successfully.\"\n input.well@model.name = loess.model@name\n input.well@equation = loess.model@expression\n # There are no estimated params, so just return the initial guesses\n input.well@fit.par = append(as.list(loess.model@guess(input.well,smooth.param)),list(\"smoothing parameter\"=smooth.param))\n # Note: since there are no params there are no Std. Errors either\n input.well@inflection.time = inflection.time(input.well)\n # calculate Rss for loess\n input.well@rss = sum((input.data$y-predict(fit))**2)\n } else {\n fit = fit.nls.model(input.well,richards)\n # should we use richards? Yes, unless the v param is close to 1 or Zero\n if(class(fit) == \"nls\"){\n rich.fit = fit # if v is significant or other fits consequently fail\n fit.par = as.list(coef(fit))\n # is fit similar to the Logistic?\n if(fit.par$v >= .5 && abs(fit.par$v-1) < 2*summary(fit)$parameters[\"v\",\"Std. Error\"] ){\n fit = fit.nls.model(input.well,logistic)\n input.well@fit.info = paste(\"Logistic model fit successfully.\")\n input.well@model.name = logistic@name\n input.well@equation = logistic@expression\n # is fit similar to Gompertz?\n }else if(fit.par$v < .5 && abs(fit.par$v) < 2*summary(fit)$parameters[\"v\",\"Std. Error\"]){\n fit = fit.nls.model(input.well,gompertz)\n input.well@fit.info = \"Gompertz model fit successfully.\"\n input.well@model.name = gompertz@name\n input.well@equation = gompertz@expression\n # v param is significant. stick with Richards\n }else{\n input.well@fit.info = paste(\"Richards model fit successfully.\")\n input.well@model.name = richards@name\n input.well@equation = richards@expression\n }\n # just in case logistic or gompertz failed to fit...\n if(class(fit) != \"nls\"){\n fit = rich.fit\n input.well@fit.info = paste(\"Richards model fit successfully.\")\n input.well@model.name = richards@name\n input.well@equation = richards@expression\n }\n } else{\n # Richards failed. try backup models\n fit = fit.nls.model(input.well,logistic)\n if(class(fit) != \"nls\"){\n # last resort try gompertz\n fit = fit.nls.model(input.well,gompertz)\n if(class(fit) != \"nls\"){\n input.well@fit.info = \"Model fitting failed.\" \n } else{\n input.well@fit.info = \"Gompertz model fit successfully.\"\n input.well@model.name = gompertz@name\n input.well@equation = gompertz@expression\n }\n }else{\n input.well@fit.info = paste(\"Logistic model fit successfully.\")\n input.well@model.name = logistic@name\n input.well@equation = logistic@expression\n }\n }\n }\n \n # If no error was reported by the model fitting, add coefficients to slot \"fit.par\",\n # Also add the Standard Errors for each parameter\n if (class(fit) == \"nls\"){\n input.well@nls = fit\n input.well@inflection.time = inflection.time(input.well)\n input.well@fit.par = as.list(coef(fit))\n rSs = sum(residuals(fit)**2)\n \n if (length(rSs) != 0)\n input.well@rss = rSs\n else \n input.well@rss = NA\n input.well@fit.std.err = as.list(summary(fit)$parameters[,\"Std. Error\"])\n }\n # Output to console\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n}\n\n# Fit nls model to a well using a specified model\n# Arguments:\n# input.well: object of class well\n# model: object of class model, e.g. richards, gompertz or logistic\nfit.nls.model <- function (input.well, model) {\n # Get OD vs. time data from well\n input.data = data.from(input.well, na.rm = T)\n # Change column headers of input.data to the more general \"Time\" vs. \"y\"\n names(input.data) = c(\"Time\", \"y\")\n \n # Set a lower bound for nls model parameters A and b to slightly lower than min(y)\n low.y = min(input.data$y,na.rm=T)\n low.y = low.y - 0.1*abs(low.y)\n \n # Set the initial guess\n start = model@guess(input.well)\n \n # Set lower bounds\n if (length(start) == 4) {\n lower = c(low.y,low.y,0,0)\n } else if (length(start) == 5) {\n lower = c(low.y,low.y,0,0,0.07)\n } else {\n stop(\"Unsupported model: \", model@name)\n }\n # Make sure initial guess values do not violate lower bounds\n start[start < lower] = lower[start < lower]\n \n # Fit the model\n try(nls(formula = model@formula, data = input.data, start = start, algorithm=\"port\", lower=lower), silent = TRUE)\n}\n",
- "created" : 1425413287007.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2355348307",
- "id" : "B0B7416E",
- "lastKnownWriteTime" : 1425509229,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/fit.model.R",
- "project_path" : "R/fit.model.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/E3886516 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/E3886516
deleted file mode 100644
index 6b8ffd1..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/per/t/E3886516
+++ /dev/null
@@ -1,17 +0,0 @@
-{
- "contents" : "setwd(\"~/Downloads/\")\nfile.list = file.name = \"YPDAFEXglucoseTests_2-25-10.csv\"\nlayout.file = \"YPDAFEXglucoseTests_2-25-10_Layout.csv\"\nsingle.plate = T\nout.dir = getwd()\ngraphic.dir = paste(out.dir, \"/pics\", sep = \"\")\nadd.constant = 1\nblank.value = NULL\nstart.index = 2\ngrowth.cutoff = 0.05\nuse.linear.param = F\nuse.loess = F\nsmooth.param = 0.6\npoints.to.remove = 0\nremove.jumps = F\nsilent = F\nverbose = T\nreturn.fit = F\noverview.jpgs = T\nplate.nrow = 8\nplate.ncol = 12\ninput.skip.lines = 0\nmulti.column.headers = c(\"Plate ID\", \"Well\", \"OD\", \"Time\")\nsingle.column.headers = c(\"\",\"A1\")\nlayout.sheet.headers = c(\"Strain\", \"Media Definition\")\n\nt = gcat.analysis.main(file.list, single.plate, layout.file = NULL, \n out.dir = getwd(), graphic.dir = paste(out.dir, \"/pics\", sep = \"\"), \n add.constant = 1, blank.value = NULL, start.index = 2, growth.cutoff = 0.05,\n use.linear.param=use.linear.param, use.loess=use.loess, smooth.param=0.1,\n points.to.remove = 0, remove.jumps = F, time.input = NA,\n plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0,\n multi.column.headers = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\"), single.column.headers = c(\"\",\"A1\"), \n layout.sheet.headers = c(\"Strain\", \"Media Definition\"),\n silent = F, verbose = F, return.fit = F, overview.jpgs = T)\n",
- "created" : 1425413468083.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "3533220190",
- "id" : "E3886516",
- "lastKnownWriteTime" : 1425422563,
- "path" : "~/Documents/GCAT4/trunk/R/GCAT/R/addingParams.R",
- "project_path" : "R/addingParams.R",
- "properties" : {
- "tempName" : "Untitled1"
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/11120BC5 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/11120BC5
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/11120BC5
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/246004C1 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/246004C1
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/246004C1
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2A2D32C6 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2A2D32C6
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2A2D32C6
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2E01469A b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2E01469A
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/2E01469A
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3349C922 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3349C922
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3349C922
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3E229DBB b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3E229DBB
deleted file mode 100644
index 32390ac..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3E229DBB
+++ /dev/null
@@ -1,3 +0,0 @@
-{
- "tempName" : "Untitled1"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3F8DCF42 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3F8DCF42
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/3F8DCF42
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/41E3FF1E b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/41E3FF1E
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/41E3FF1E
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/4E0B0FCA b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/4E0B0FCA
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/4E0B0FCA
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/54A894D0 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/54A894D0
deleted file mode 100644
index 32390ac..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/54A894D0
+++ /dev/null
@@ -1,3 +0,0 @@
-{
- "tempName" : "Untitled1"
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/61AC4784 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/61AC4784
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/61AC4784
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8250000E b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8250000E
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8250000E
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8341E8EB b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8341E8EB
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8341E8EB
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8CEE19EC b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8CEE19EC
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/8CEE19EC
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9127986F b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9127986F
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9127986F
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/972C75E0 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/972C75E0
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/972C75E0
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9D1E246E b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9D1E246E
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/9D1E246E
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/A3297EED b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/A3297EED
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/A3297EED
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/ADFC007A b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/ADFC007A
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/ADFC007A
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/B005EA5F b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/B005EA5F
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/B005EA5F
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/BA83F0DB b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/BA83F0DB
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/BA83F0DB
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DB5E821B b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DB5E821B
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DB5E821B
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DEDB5126 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DEDB5126
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/DEDB5126
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E533697B b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E533697B
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E533697B
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E98DE888 b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E98DE888
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/E98DE888
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/EDBD9DE b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/EDBD9DE
deleted file mode 100644
index 7a73a41..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/EDBD9DE
+++ /dev/null
@@ -1,2 +0,0 @@
-{
-}
\ No newline at end of file
diff --git a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/INDEX b/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/INDEX
deleted file mode 100644
index 0543b58..0000000
--- a/R/GCAT/.Rproj.user/7CB73FCA/sdb/prop/INDEX
+++ /dev/null
@@ -1,18 +0,0 @@
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FNAMESPACE="972C75E0"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2FGCAT.main.R="11120BC5"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2FaddingParams.R="54A894D0"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fclass.model.R="9127986F"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fclass.well.R="8250000E"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Ffit.model.R="3F8DCF42"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Ffitted.calculations.R="E533697B"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fmisc.R="8341E8EB"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fnormalize.and.transform.R="EDBD9DE"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fplot.fit.R="2E01469A"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Fslope.analysis.R="4E0B0FCA"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Ftable.output.R="B005EA5F"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2FR%2Ftable2well.R="A3297EED"
-~%2FDocuments%2FGCAT4%2Ftrunk%2FR%2FGCAT%2Fman%2Fgcat.load.data.Rd="2A2D32C6"
-~%2FDocuments%2FGCAT4_old%2Ftrunk%2FR%2FGCAT%2FR%2FGCAT.main.R="61AC4784"
-~%2FDocuments%2FGCAT4_old%2Ftrunk%2FR%2FGCAT%2FR%2FaddingParams.R="DEDB5126"
-~%2FDocuments%2FGCAT4_old%2Ftrunk%2FR%2FGCAT%2FR%2Fclass.model.R="8CEE19EC"
-~%2FDocuments%2FGCAT4_old%2Ftrunk%2FR%2FGCAT%2FR%2Ffit.model.R="9D1E246E"
diff --git a/R/GCAT/DESCRIPTION b/R/GCAT/DESCRIPTION
index bff0405..9d09fb1 100755
--- a/R/GCAT/DESCRIPTION
+++ b/R/GCAT/DESCRIPTION
@@ -9,9 +9,10 @@ Description: Imports high-throughput growth curve data from microtiter
(specific growth rate, maximum growth capacity, and lag time)
for each well in a read.
The code was written by Jason Shao (no longer at GLBRC) and Nate DiPiazza.
-Version: 5.0
-Depends: pheatmap, gplots
+Version: 5.0.2
+Depends: pheatmap, gplots, methods
Maintainer: Yury Bukhman
License: LGPL-3
-Date: 2014-02-10
+Date: 2015-04-21
Author: Jason Shao, Nate DiPiazza , Minh Duc Bui, Yury V Bukhman
+Suggests: testthat
diff --git a/R/GCAT/NAMESPACE b/R/GCAT/NAMESPACE
index 55fdde5..433990b 100755
--- a/R/GCAT/NAMESPACE
+++ b/R/GCAT/NAMESPACE
@@ -1,8 +1,24 @@
-# Default NAMESPACE created by R
-# Remove the previous line if you edit this file
+# Generated by roxygen2 (4.1.1): do not edit by hand
-# Export all names
-exportPattern(".")
-
-# Import all packages listed as Imports or Depends
-import("pheatmap", "gplots")
+export(gcat.analysis.main)
+exportClasses(well)
+exportMethods(getAddInfo)
+exportMethods(getCurPar)
+exportMethods(getEquation)
+exportMethods(getFitErr)
+exportMethods(getFitInfo)
+exportMethods(getFitPar)
+exportMethods(getInflectionTime)
+exportMethods(getLoess)
+exportMethods(getModelName)
+exportMethods(getNorm)
+exportMethods(getPosition)
+exportMethods(getRSS)
+exportMethods(getScreenData)
+exportMethods(getStartIndex)
+exportMethods(getUseLog)
+exportMethods(getWellInfo)
+exportMethods(getnls)
+import(gplots)
+import(methods)
+import(pheatmap)
diff --git a/R/GCAT/R.Rproj b/R/GCAT/R.Rproj
deleted file mode 100644
index f7442a1..0000000
--- a/R/GCAT/R.Rproj
+++ /dev/null
@@ -1,17 +0,0 @@
-Version: 1.0
-
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
-
-EnableCodeIndexing: Yes
-UseSpacesForTab: Yes
-NumSpacesForTab: 2
-Encoding: UTF-8
-
-RnwWeave: Sweave
-LaTeX: pdfLaTeX
-
-BuildType: Package
-PackageInstallArgs: --no-multiarch --with-keep.source
-PackageRoxygenize: rd
diff --git a/R/GCAT/R/.Rhistory b/R/GCAT/R/.Rhistory
deleted file mode 100644
index e69de29..0000000
diff --git a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/69756148 b/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/69756148
deleted file mode 100644
index 1305dca..0000000
--- a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/69756148
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n########################################################################\n# #\n# Functions to calculate various things about wells based on fit model #\n# #\n########################################################################\n#\n# Common arguments:\n# fitted.well - should be a well containing the results of , most functions will return NA if well has not been fit yet.\n# unlog - should the value be returned on the linear scale as opposed to the log-transformed scale?\n# constant.added - for returning values on the linear scale, what was the constant added before the log transform?\n# digits - passed to the function, default is no rounding (infinity digits)\n\nunlog = function(x, constant.added) {\n ########################################################################\n # Transform values back to OD scale #\n ########################################################################\nexp(x) - constant.added\n}\n\nwell.eval = function(fitted.well, Time = NULL){\n ########################################################################\n # Evaluate estimated OD at any timepoints using the fitted model #\n ########################################################################\n\n # If no timepoints are provided, use the ones collected in the experiment itself.\n\tif(!is.numeric(Time))\n\t\tTime = data.from(fitted.well)$Time\n\n # Use of equation is deprecated. Use nls and loess models stored in the well object instead\n # Attempt to use with the fitted equation and parameters to get estimates for OD at the given timepoints.\n\t#output = try(eval(fitted.well@equation, fitted.well@fit.par), silent = T)\n \n # Predict log.OD value(s) using nls model if present. If no nls model, try using loess.\n if (length(fitted.well@nls)>0) {\n output = try(predict(fitted.well@nls,list(Time=Time)),silent=T)\n } else if (length(fitted.well@loess)>0) {\n output = try(predict(fitted.well@loess,Time),silent=T)\n } else {\n output = NA\n }\n\n # Return values. If OD evaluation failed for any reason, return NULL.\n if (is.numeric(output)){\n return(output)\n } else {\n return(NULL)\n\t}\n}\n\nmodel.residuals = function(fitted.well, unlog = F){\n ########################################################################\n # Evaluate model residuals using the measured vs. fitted log.OD values #\n ########################################################################\n\tmeasured.OD = data.from(fitted.well)[,2]\n\n\t# Use with no Time argument to get fitted OD values at measured timepoints.\n\tpredicted.OD = well.eval(fitted.well)\n\n\t# If all values are valid, return the differences\n\tif (!is.numeric(predicted.OD))\n\t\treturn(NA)\n\telse\n return(measured.OD - predicted.OD)\n\t}\n\ndev.from.mean = function(fitted.well){\n ########################################################################\n # Evaluate deviations of log.OD values from the mean #\n ########################################################################\n measured.ODs = data.from(fitted.well,remove=T,na.rm=T)[,2]\n \n # Get the mean values of these measured ODs.\n mean.ODs = mean(measured.ODs)\n \n if (!is.numeric(mean.ODs))\n return (NA)\n else\n return (measured.ODs - mean.ODs)\n}\n\nrss = function(fitted.well){\n #######################################################################\n # Get the residual sum of square. #\n #######################################################################\n if (length(fitted.well@rss) == 0)\n return (NA)\n else\n return (fitted.well@rss)\n}\n\nmodel.good.fit = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate a metric for fit accuracy using squared residuals #\n ########################################################################\n\n # Sum of squared residuals\n\tRSS = rss(fitted.well)\n \n # Total sum of squared\n tot = sum(dev.from.mean(fitted.well)^2)\n \n # Coefficient of determination\n return (1 - RSS/tot)\n\t}\n\nparameter.text = function(fitted.well){\n ########################################################################\n # Output a string with values of fitted parameters #\n ########################################################################\n \n # Get a list of fitted parameters\n fit.par = fitted.well@fit.par\n \n # Giving the parameter text descriptive names.\n if (length(fitted.well@fit.par) != 0){\n names(fit.par)[1] = \"A\" \n names(fit.par)[2] = \"b\" \n names(fit.par)[3] = \"lambda\" \n names(fit.par)[4] = \"max.spec.growth.rate\" \n \n if (fitted.well@model.name == \"richards sigmoid\"){ \n names(fit.par)[5] = \"shape.par\" \n } \n \n if (fitted.well@model.name == \"richards sigmoid with linear par.\"){ \n names(fit.par)[5] = \"shape.param\" \n names(fit.par)[6] = \"linear term\"\n } \n \n if (fitted.well@model.name == \"logistic sigmoid with linear par.\")\n names(fit.par)[5] = \"linear.term\"\n \n # if loess, just show smoothing param\n if(fitted.well@model.name == \"local polynomial regression fit.\")\n fit.par = fitted.well@fit.par[\"smoothing parameter\"]\n }\n \n # Return nothing if the list is empty. Otherwise, concatenate the terms in the list with the parameter names.\n\tif(!is.list(fit.par))\n\t\treturn()\n else{\n \toutput = \"\"\n \ti = 1\n \twhile(i <= length(fit.par)){\n \t\toutput = paste(output, names(fit.par)[i], \"=\", round(as.numeric(fit.par[i]),3), \"; \", sep = \"\")\n \t\ti = i + 1\n if (i %% 6 == 0)\n output = paste(output, \"\\n\")\n \t\t}\n \toutput\n \t}\n\t}\n\nmax.spec.growth.rate = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate maximum specific growth rate #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n round(fitted.well@fit.par$u,digits)\n}\n\n\nplateau = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate plateau log.OD from fitted parameters #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n plat = fitted.well@fit.par$A + fitted.well@fit.par$b\n \n\tif (!is.numeric(plat)) {\n\t plat = NA\n\t} else {\n plat = round(plat, digits)\n\t}\n\treturn(plat)\n}\n\nbaseline = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate baseline log.OD from fitted parameters #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n\n base = fitted.well@fit.par$b\n\n # If A (plateau OD) is invalid, return NA.\n\tif (!is.numeric(fitted.well@fit.par$A))\n\t\tbase = NA\n # If b (baseline OD) is invalid but plateau OD was valid, return zero.\n else if (!is.numeric(base))\n\t\tbase = 0\n\telse{\n\t\t base = round(base, digits)\n\t\t}\n\treturn(base)\n\t}\n\ninoc.log.OD = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate log.OD at inoculation from fitted parameters #\n ########################################################################\n\n # Evaluated the fitted model at the inoculation timepoint (should be zero from using from table2wells.R)\n\tif (is.null(well.eval(fitted.well)))\n\t\treturn(NA)\n else{\n inoc.time = fitted.well@screen.data$Time[fitted.well@start.index]\n inoc.log.OD = well.eval(fitted.well, inoc.time)\n if (is.na(inoc.log.OD)) inoc.log.OD = fitted.well@fit.par$b # need this in a special case: loess fits with start.index = 1 \n return(round(inoc.log.OD, digits))\n }\n\t}\n\n\nmax.log.OD = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate max log.OD from model fit #\n ########################################################################\n\n # Evaluated the fitted model at the final timepoint (just the last valid timepoint in the experiment)\n\tif (is.null(well.eval(fitted.well)))\n\t\treturn(NA)\n else{\n \treturn(round(max(well.eval(fitted.well),na.rm=T), digits))\n }\n}\n\n\nprojected.growth = function(fitted.well,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n\tplateau(fitted.well,digits) - inoc.log.OD(fitted.well,digits)\n}\n\nprojected.growth.OD = function(fitted.well,constant.added,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n value = unlog(plateau(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added)\n round(value,digits)\n}\n\n\nachieved.growth = function(fitted.well,digits=Inf) {\n ########################################################################\n # Calculate achieved growth: max.log.OD minus the inoculated log.OD #\n ########################################################################\n max.log.OD(fitted.well,digits) - inoc.log.OD(fitted.well,digits)\n}\n\nachieved.growth.OD = function(fitted.well,constant.added,digits=Inf) {\n ########################################################################\n # Calculate projected growth: plateau minus the inoculated log.OD #\n ########################################################################\n value = unlog(max.log.OD(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added)\n round(value,digits)\n}\n\nreach.plateau = function(fitted.well, cutoff = 0.75){\n ########################################################################\n # Did the curve come close to the plateau OD during the experiment? #\n ########################################################################\n\n plat = plateau(fitted.well)\n inoc = inoc.log.OD(fitted.well)\n final = max.log.OD(fitted.well)\n\n\tif (!is.na(final)){\n # If the plateau is the same as the OD at inoculation, return TRUE\n if ((plat - inoc) == 0)\n return(T)\n # If the difference between the final OD and inoculation OD is at least a certain proportion\n # of the difference between the plateau and inoculated ODs, return TRUE.\n else\n return((final - inoc) / (plat - inoc) > cutoff)\n\t\t}\n\telse\n\t\treturn(T)\n\t\t# If no final OD was calculated (if curve was not fit properly) just return T.\n\t}\n\n\nlag.time = function(fitted.well, digits = Inf){\n ########################################################################\n # Calculate the lag time from the fitted OD #\n ########################################################################\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n fitted.well@fit.par$lam\n}\n\n# new params for GCAT 4.0\n\namplitude = function(fitted.well){\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n \n return(fitted.well@fit.par$A)\n}\n\nshape.par = function(fitted.well){\n if(length(fitted.well@fit.par) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.par$v), NA, fitted.well@fit.par$v)\n}\n\nmax.spec.growth.rate.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$u), NA, fitted.well@fit.std.err$u)\n}\n# The SE values are the errors.\nlag.time.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$lam), NA, fitted.well@fit.std.err$lam)\n}\n\nshape.par.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$v), NA, fitted.well@fit.std.err$v)\n}\n\namplitude.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$A), NA, fitted.well@fit.std.err$A)\n}\n\nbaseline.SE = function(fitted.well){\n if(length(fitted.well@fit.std.err) == 0)\n return(NA)\n ifelse(is.null(fitted.well@fit.std.err$b), NA, fitted.well@fit.std.err$b)\n}\n\n# used to calulate the inflection.time value\ninflection.time = function(well){\n if (length(well@loess) == 0 && length(well@nls) == 0) return(NA) # can' compute inflection time in the absence of a fit\n data = data.from(well)\n Time = data[,1]\n t = seq(from = min(Time), to = max(Time), by = (max(Time)-min(Time))/1000)\n y = well.eval(well,t)\n if (is.null(y)) return(NA)\n delta.t = diff(t)\n dydt = diff(y)/delta.t\n infl.index = which.max(dydt)\n t[infl.index]\n}\n",
- "created" : 1423868976528.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "4167103586",
- "id" : "69756148",
- "lastKnownWriteTime" : 1423868994,
- "path" : "~/Documents/GCAT4_old/trunk/R/GCAT/R/fitted.calculations.R",
- "project_path" : "fitted.calculations.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/6C03FCBD b/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/6C03FCBD
deleted file mode 100644
index 0e0a97c..0000000
--- a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/6C03FCBD
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead,\n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n########################################################################\n# #\n# Fit a parameterized model to the growth data in a well object. #\n# #\n# There are now three modelling choices: #\n# 1) Sigmoid model (no linear param c) #\n# 2) Linear Sigmoid model #\n# 3) Loess (with optional smoothing parameter) #\n########################################################################\n#' fit.model \n#' \n#' This function will use the function stored in the \"guess\" slot of \\code{growth.model} to calculate initial guesses \n#' for growth.model parameters, then it will use the \"formula\" slot with \\code{nls} to fit a non-linear least squares \n#' \\code{growth.model} or Local Polynomial Regression Fitting to the data. Richards model is first fitted. \n#' If the shape parameter is statisticaly significant then Richards is used. If it is within 2 SE of 1 or Zero than \n#' a simpler model is preferred. If the Richards fit fails, then Logistic is tried. If it fails, Gompertz is tried. \n#' Model fit failure is reported if none of the models can sucessfully fit the data \n#' \n#' @param input.well The well needed to be fitted with the given model. \n#' @param growth.model What growth model should be used? \n#' @param backup.growth.model If \\code{gowth.mode} fails, this model will be used. \n#' @param fit.if.no.growth should the function attempt to fit a well even if there was no growth detected? default is F \n#' @param silent output back to R console? \n#' @param use.linear.param: Should an additional linear parameter (c) be used when fitting the data to the model? \n#' @param use.loess: Should Local Polynomial Regression Fitting (loess function) be used instead of nls? \n#' @param smooth.param: If loess is used, an optional smoothing parameter. Default is .6 \nfit.model = function(input.well, growth.model, backup.growth.model = NULL, fit.if.no.growth = F, \n use.linear.param=F, use.loess=F, smooth.param, silent = T){\n \n # Conditional breakpoint: go into debugging mode when fitting a specific well\n #if (input.well@position[\"row\"] == \"A\" && input.well@position[\"col\"] == \"12\") browser()\n \n # Change all relevant slots to or blank values\n input.well@model.name = \"\"\n input.well@fit.par = list()\n input.well@equation = expression()\n \n # Get OD vs. time data from well\n input.data = data.from(input.well, na.rm = T)\n \n # Skip well if in slot \"curve.par\" is set to true, and is false.\n if(!fit.if.no.growth & lacks.growth(input.well)){\n input.well@fit.info = \"skipped - no growth in well.\"\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n }\n # Skip well if there are fewer than 5 data points left in the analysis.\n if (length(input.data$Time) < 5){\n input.well@fit.info = \"skipped - not enough points.\"\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n }\n \n # Change column headers of input.data to the more general \"Time\" vs. \"y\"\n names(input.data) = c(\"Time\", \"y\")\n \n # Set a lower bound for nls model parameters A and b to slightly lower than min(y)\n low.y = min(input.data$y,na.rm=T)\n low.y = low.y - 0.1*abs(low.y)\n \n # Extract the model formula from (slot \"formula\")\n # Use the function from slot \"guess\" to calculate initial guesses for model parameters based on slope estimates in \n # Attempt to fit a nonlinear least squares odel using \n \n # Creating loess, logistics, richards, and gompertz model. \n loess.e = expression(\"loess\")\n loess.f = formula(Time ~ y)\n loess.model = model(\"local polynomial regression fit.\", loess.e, loess.f, loess.g)\n remove(loess.e, loess.f, loess.g)\n \n ########################################################################\n # Create the logistic 4-parameter model (when v ~ 1) #\n ########################################################################\n logistic.g = function(well,smooth.param=0.75) {\n loess.model@guess(well,smooth.param)\n }\n \n ########################################################################\n # Create the Richards 5-parameter model #\n ########################################################################\n richards.g = function(well,smooth.param=0.75){\n c(loess.model@guess(well,smooth.param),v=0.5)\n }\n \n ########################################################################\n # Create the Gompertz model (might be useful as a #\n # limiting case of Richards model when v ~ 0) #\n ########################################################################\n gompertz.g = function(well,smooth.param=0.75){\n loess.model@guess(well,smooth.param)\n }\n \n logistic.e = expression((A/(1+exp((4*u/A)*(lam-Time)+2)))+b)\n logistic.f = formula(y~(A/(1+exp((4*u/A)*(lam-Time)+2)))+b)\n logistic = model(\"logistic sigmoid.\", logistic.e, logistic.f, logistic.g)\n remove(logistic.e, logistic.f, logistic.g)\n \n richards.e = expression(A*(1+v*exp(1+v)*exp((u/A)*(1+v)**(1+1/v)*(lam-Time)))**(-1/v)+b)\n richards.f = formula(y~A*(1+v*exp(1+v)*exp((u/A)*(1+v)**(1+1/v)*(lam-Time)))**(-1/v)+b)\n richards = model(\"richards sigmoid\", richards.e, richards.f, richards.g)\n remove(richards.e, richards.f, richards.g)\n \n gompertz.e = expression(A*exp(-exp((u*exp(1)/A)*(lam-Time)+1))+b)\n gompertz.f = formula(y~A*exp(-exp((u*exp(1)/A)*(lam-Time)+1))+b)\n gompertz = model(\"gompertz sigmoid\", gompertz.e, gompertz.f, gompertz.g)\n remove(gompertz.e, gompertz.f, gompertz.g)\n \n \n \n # 3) Loess (with optional smoothing parameter)\n if(use.loess){\n number.of.points = nrow(input.well@screen.data)\n if (smooth.param <= 1/number.of.points)\n exception(\"Invalid input\", \"Smoothing parameter is out of range.\")\n \n fit = try(loess(y~Time, data=input.data, span=smooth.param), silent=TRUE)\n input.well@loess = fit\n if (class(fit) != \"loess\") stop(\"loess fit failed on well\", paste(input.well@position,collapse=\" \"))\n input.well@fit.info = \"Loess model fit successfully.\"\n input.well@model.name = loess.model@name\n input.well@equation = loess.model@expression\n # There are no estimated params, so just return the initial guesses\n input.well@fit.par = append(as.list(loess.model@guess(input.well,smooth.param)),list(\"smoothing parameter\"=smooth.param))\n # Note: since there are no params there are no Std. Errors either\n input.well@inflection.time = inflection.time(input.well)\n # calculate Rss for loess\n input.well@rss = sum((input.data$y-predict(fit))**2)\n } else {\n fit = fit.nls.model(input.well,richards)\n # should we use richards? Yes, unless the v param is close to 1 or Zero\n if(class(fit) == \"nls\"){\n rich.fit = fit # if v is significant or other fits consequently fail\n fit.par = as.list(coef(fit))\n # is fit similar to the Logistic?\n if(fit.par$v >= .5 && abs(fit.par$v-1) < 2*summary(fit)$parameters[\"v\",\"Std. Error\"] ){\n fit = fit.nls.model(input.well,logistic)\n input.well@fit.info = paste(\"Logistic model fit successfully.\")\n input.well@model.name = logistic@name\n input.well@equation = logistic@expression\n # is fit similar to Gompertz?\n }else if(fit.par$v < .5 && abs(fit.par$v) < 2*summary(fit)$parameters[\"v\",\"Std. Error\"]){\n fit = fit.nls.model(input.well,gompertz)\n input.well@fit.info = \"Gompertz model fit successfully.\"\n input.well@model.name = gompertz@name\n input.well@equation = gompertz@expression\n # v param is significant. stick with Richards\n }else{\n input.well@fit.info = paste(\"Richards model fit successfully.\")\n input.well@model.name = richards@name\n input.well@equation = richards@expression\n }\n # just in case logistic or gompertz failed to fit...\n if(class(fit) != \"nls\"){\n fit = rich.fit\n input.well@fit.info = paste(\"Richards model fit successfully.\")\n input.well@model.name = richards@name\n input.well@equation = richards@expression\n }\n } else{\n # Richards failed. try backup models\n fit = fit.nls.model(input.well,logistic)\n if(class(fit) != \"nls\"){\n # last resort try gompertz\n fit = fit.nls.model(input.well,gompertz)\n if(class(fit) != \"nls\"){\n input.well@fit.info = \"Model fitting failed.\" \n } else{\n input.well@fit.info = \"Gompertz model fit successfully.\"\n input.well@model.name = gompertz@name\n input.well@equation = gompertz@expression\n }\n }else{\n input.well@fit.info = paste(\"Logistic model fit successfully.\")\n input.well@model.name = logistic@name\n input.well@equation = logistic@expression\n }\n }\n }\n \n # If no error was reported by the model fitting, add coefficients to slot \"fit.par\",\n # Also add the Standard Errors for each parameter\n if (class(fit) == \"nls\"){\n input.well@nls = fit\n input.well@inflection.time = inflection.time(input.well)\n input.well@fit.par = as.list(coef(fit))\n rSs = sum(residuals(fit)**2)\n \n if (length(rSs) != 0)\n input.well@rss = rSs\n else \n input.well@rss = NA\n input.well@fit.std.err = as.list(summary(fit)$parameters[,\"Std. Error\"])\n }\n # Output to console\n if (!silent)\n cat(plate.name(input.well), well.name(input.well), \":\", input.well@fit.info, \"\\n\")\n return(input.well)\n}\n\n# Fit nls model to a well using a specified model\n# Arguments:\n# input.well: object of class well\n# model: object of class model, e.g. richards, gompertz or logistic\nfit.nls.model <- function (input.well, model) {\n # Get OD vs. time data from well\n input.data = data.from(input.well, na.rm = T)\n # Change column headers of input.data to the more general \"Time\" vs. \"y\"\n names(input.data) = c(\"Time\", \"y\")\n \n # Set a lower bound for nls model parameters A and b to slightly lower than min(y)\n low.y = min(input.data$y,na.rm=T)\n low.y = low.y - 0.1*abs(low.y)\n \n # Set the initial guess\n start = model@guess(input.well)\n \n # Set lower bounds\n if (length(start) == 4) {\n lower = c(low.y,low.y,0,0)\n } else if (length(start) == 5) {\n lower = c(low.y,low.y,0,0,0.07)\n } else {\n stop(\"Unsupported model: \", model@name)\n }\n # Make sure initial guess values do not violate lower bounds\n start[start < lower] = lower[start < lower]\n \n # Fit the model\n fit = try(nls(formula = model@formula, data = input.data, start = start, algorithm=\"port\", lower=lower), silent = TRUE)\n}\n",
- "created" : 1423869161433.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "2981843703",
- "id" : "6C03FCBD",
- "lastKnownWriteTime" : 1423869188,
- "path" : "~/Documents/GCAT4_old/trunk/R/GCAT/R/fit.model.R",
- "project_path" : "fit.model.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/7189AD56 b/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/7189AD56
deleted file mode 100644
index 2accf9e..0000000
--- a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/7189AD56
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\nrequire(pheatmap)\nrequire(gplots)\n\n########################################################################\n# #\n# Graphic output functions for fitted well objects. The functions are #\n# fairly complicated and intertwined and may need revision. #\n# #\n########################################################################\n\n########################################################################\n# Basic function plots time vs. OD from a well object #\n########################################################################\n#' plot.data \n#' \n#' Basic function plots time vs. OD from a well object \n#' \n#' @param input.well The well object that need to be plottedd \n#' @param unlog should data be plotted on a linear (vs. logarithmic) scale? \n#' @param view.raw.data should the raw data be plotted? ( \n#' @param number.points should points be labeled with numeric indices? \n#' @param scale determines the font scale for the entire graph. all cex values are calculated from this. \n#' @param draw.symbols - should be called on the well and markings drawn on the graph? \n#' @param ... additional arguments passed to plot() \n\nplot.data = function(input.well, view.raw.data = F, unlog = F, scale = 1, \n main = paste(plate.name(input.well), well.name(input.well)), number.points = T, \n draw.symbols = F, constant.added,...){\n \n # Get data as well as a vector showing which points were removed. \n\tinput.data = data.from(input.well, remove = F, remove.tanking = F, raw.data=view.raw.data)\n\tremoved.points = !(rownames(input.data) %in% rownames(data.from(input.well, remove = T, remove.tanking = T)))\n point.colors = as.character(factor(removed.points,levels=c(F,T),labels=c(\"black\",\"gray80\")))\n\n # Draw the axes and all text labels first.\n\tpar(mar = c(5, 4, 4, 5)+0.1)\n\tplot(input.data, main = main, xlab = \"Time(hours)\", ylab = \"log(OD - blank + const)\",\n mex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, cex.lab = 1.2*scale, type =\"n\",...)\n \n\t# Draw a second vertical axis, showing unlogged OD scale\n\t# - Determine the range of the labels: from min.OD to max.OD\n if (class(try(ylim,silent=T)) == \"try-error\") {\n OD = unlog(input.data[,2],constant.added)\n baseline.OD = unlog(baseline(input.well),constant.added)\n min.OD = min(min(OD,na.rm=T),baseline.OD,na.rm=T)\n plateau.OD = unlog(plateau(input.well),constant.added)\n max.OD = max(max(OD,na.rm=T),plateau.OD,na.rm=T)\n } else {\n min.OD = unlog(ylim[1],constant.added)\n max.OD = unlog(ylim[2],constant.added)\n }\n\t# - Compute labels and their positions\n\tOD.labels = seq(from = min.OD, to = max.OD, length.out = 5)\n\tOD.labels = signif(OD.labels,2)\n\tOD.at = log(OD.labels+constant.added)\n # - Draw the axis\n\taxis(side=4, at=OD.at, labels=OD.labels, cex.axis = 1.2*scale, cex.lab = 1.2*scale)\n\tmtext(4, text = \"OD - blank\", line = 3, cex=1.2)\n \n\t# If is true, then label each point with the index of its timepoint and plot removed points in grey, others in black. \n if (number.points)\n\t\ttext(input.data$Time, input.data[,2], rownames(input.data), col = point.colors, cex = 0.5*scale)\n\t# Otherwise plot all points, using a different plotting character for removed points. \t\n else\n\t\tpoints(input.data$Time, input.data[,2], pch = 1 + removed.points*15)\n\n # If is set to T, then draw all the markings that makes to determine curve parameters. \n\tif (draw.symbols & !view.raw.data)\n\t\tcheck.slopes(input.well, draw = T)\n\treturn()\n\t}\n\t\n########################################################################\n# Plots the fitted model curve from a well object if it exists #\n########################################################################\n#\n# time: specify which points (in units of time) to plot fitted OD values for. if not specifies, plot all timepoints in range of well. \n\nplot.model = function(input.well, col = 1, scale = 1, lty = 1, time = NULL, unlog = F, constant.added=1){\n\t \n input.data = data.from(input.well)\n\tgrowth = input.data[,2]\n\n # If no list of timepoints is specified, get a list of 360 timepoints (should be smooth enough) from the well's range. \n\tif (is.null(time)){\n\t\ttime.fin = max(data.from(input.well, raw.data = T, remove = F, remove.tanking = F)$Time)\n\t\ttime = seq(0, time.fin, length.out = 360)\n\t}\n \n # Evaluate the predicted OD at the specified timepoints based on the fitted model. \n\tpredicted.OD = well.eval(input.well, time) \n # If any values were returned, plot the lines on the current graph. Otherwise, just return without doing anything.\n\tif (is.numeric(predicted.OD))\n\t\tlines(time, predicted.OD, col = col, lty = lty, lw = 2 * scale) \n else\n return()\t\n \n}\n\t\n########################################################################\n# Put various parameters and info in text form on the graphs #\n########################################################################\n#\t\t\t\ndraw.text = function(input.well, scale = 0.5, xlim = 0, ylim = 0,...){\n\n\tinput.data = data.from(input.well, remove = F, remove.tanking = F)\n\tfit = input.well@fit.par\n \n # - fit information (fit status, model name if available, jump detection output, fit parameters if available) from well \n # color = red if no fit, blue if fit, green if skipped\n # - empty or inoculated well.\n # color = green if empty, blue if inoculated, red if inoculated but has no growth or empty but has growth. \n \n \tcol2 = \"blue\" \n\ttext2 = paste(input.well@fit.info, input.well@model.name, input.well@add.info, \"\\n\", parameter.text(input.well))\n\n\tif (length(input.well@fit.par) == 0) # no fit\n\t\tcol2 = \"red\"\n\n\tif (is.empty(input.well)){\n\t\t text1 = \"empty well\"\n\t\tif(!lacks.growth(input.well) | length(input.well@fit.par) == 0) # growth curve fit for an empty well\n\t\t\tcol1 = \"red\"\n\t\telse\n\t\t\tcol1 = \"forestgreen\"\n\t\tif (length(input.well@model.name) == 0) # well was skipped \n\t\t\tcol1 = col2 = \"forestgreen\"\n\t\t}\n\telse{\n\t\ttext1 = \"inoculated well\"\n\t\tif(lacks.growth(input.well) | length(input.well@fit.par) == 0) # failure to fit an inoculated well\n\t\t\tcol1 = \"red\"\n\t\telse\n\t\t\tcol1 = \"forestgreen\"\n\t\t}\n\t\t\n\t# - goodness of fit metric. \n\t# color = red if below 2, yellow if between 2 and 2.72, and green if above 2.72. \n\t\n\tif(!is.na(model.good.fit(input.well))){\n\t\t#if (model.good.fit(input.well, unlog = F) > 2.72)\n\t\t#\tcol1.5 = \"forestgreen\"\n\t\t#else if (model.good.fit(input.well, unlog = F)> 2.0)\n\t\t#\tcol1.5 = \"gold2\"\n\t\t#else\n\t\t#\tcol1.5 = \"red\"\n col1.5 = \"forestgreen\"\n\t\ttext1.5 = paste(\"R squared:\", round(model.good.fit(input.well),3))\n\t\t}\n\telse\n\t col1.5 = text1.5 = NULL\n\t \n # Print all text at the top of the graph with approprate positions and scaling \t\n \ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0.025 * diff(ylim), \n \t\ttext1.5, cex = 1.5*scale, col = col1.5) \n \t\t\n\ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0 * diff(ylim), \n\t\t\ttext1, cex = 1.5*scale, col = col1)\n\n\ttext(x = xlim[1] + 0.50 * diff(xlim), y = ylim[2] - 0.03 * diff(ylim), \n\t\t\ttext2, pos = 1, cex = 1.5*scale, col = col2)\t\n\t}\t\n\n########################################################################\n# Draw lines on graph denoting calculated parameters #\n########################################################################\n#\t\n# - should curve parameters be labeled? \n\t\t\ndraw.calc.par = function(input.well, scale = 0.5, unlog = F, constant.added, show.num = T){\n\n # Don't do anything if well was not fit. \n\tif (is.null(well.eval(input.well)))\n\t\treturn()\n \n # Collect values for various curve parameters. \n\tbaseline = baseline(input.well)\n\tinoc.log.OD = inoc.log.OD(input.well)\n\tmax.log.OD = max.log.OD(input.well) \n\tplateau = plateau(input.well) \n\tinflection.time = input.well@inflection.time # was a param in model\n\tfin.time = (inflection.time+max(data.from(input.well)[,1]))/2\n\n\t# = timepoint at greatest growth \n\t# = OD measurement at , minus the constant added before taking the log (if reversing the transformation)\n\t# = slope (on log scale) at (specific growth)\n # had to add the unlog code. was calculated differently before NWD 7/21/14\n\tmax.slope = max.spec.growth.rate(input.well)\n\tmax.y = well.eval(input.well, inflection.time)\n\tlag.x = lag.time(input.well) \n lag.y = baseline\n\t\n\t# ---- Specific growth rate ---- #\n\tlines(c(lag.x, inflection.time), c(lag.y, max.y), lty = 2, col = \"red\")\n\n\n # Blue dotted line at time of maximum growth, with text label for specific growth rate. \n\tabline(v = inflection.time, lty = 2, lw = (scale^2)*2, col = \"blue\")\n if(show.num) text(inflection.time, max.y, round(max.slope,3), col = \"blue\", cex = 1.5*scale, pos = 2)\n\n # inoculation OD and baseline of the fitted model\n abline(h = inoc.log.OD, lw = scale*2, lty = 3)\n abline(h = baseline, col = \"red\", lw = (scale^2)*2, lty = 2)\n if(show.num) {\n text(fin.time, inoc.log.OD, paste(round(inoc.log.OD,3),\"\\n\",sep=\"\") , col = \"black\", cex = 1.5*scale, pos = 2)\n text(fin.time, baseline, paste(\"\\n\\n\", round(baseline,3), sep=\"\") , col = \"red\", cex = 1.5*scale, pos = 2)\n }\n\n # ---- Lag time ---- #\n # Do not draw a horizontal line to lag time if it is 0 or negative. \n # Otherwise draw a red line from the starting point to the lag time, and label with the lag time \n \tif (lag.time(input.well) == 0){\n\t\tif(show.num) text(0, inoc.log.OD, \"\\n\\n0.000\", col = \"red\", cex = 1.5*scale, pos = 4)\n\t\t}\n\telse{\n\t\tlines(c(0, lag.x), c(baseline, baseline), col = \"red\", lw = (scale^2)*2, lty = 2)\n\t\tif(show.num) text(lag.x, lag.y, paste(\"\\n\\n\", round(lag.time(input.well),3)), col = \"red\", cex = 1.5*scale, pos = 2)\n\t}\n\n # ---- Total growth ---- #\n\n # Draw horizontal lines for the max.log.OD in black, the plateau in green and the initial OD in black.\n abline(h = max.log.OD, lty = 3, lw = scale*2)\n abline(h = plateau, lty = 2, lw = (scale^2)*2, col = \"forestgreen\")\n\t\n\t# Draw a vertical line from the initial OD to the final OD in black, and then to the plateau in gray. \n\tlines(c(fin.time, fin.time), c(inoc.log.OD, max.log.OD), lw = (scale^2)*2, lty = 3)\n\tlines(c(fin.time, fin.time), c(max.log.OD, plateau), lw = (scale^2)*2, lty = 3, col = \"grey\")\n\n # Text: plateau and initial ODs (on left), difference between initial and final OD on right\n\tif(show.num){\n text(fin.time, plateau, paste(round(plateau,3),\"\\n\",sep=\"\") , col = \"forestgreen\", cex = 1.5*scale, pos = 2)\n text(fin.time, max.log.OD, paste(\"\\n\\n\\n\",round(max.log.OD,3),sep=\"\") , col = \"black\", cex = 1.5*scale, pos = 2)\n text(fin.time, .5*(max.log.OD-inoc.log.OD)+inoc.log.OD, round(max.log.OD - inoc.log.OD,3), cex = 1.5*scale, pos = 4)\n # difference between final and plateau OD (if large enough) \n if (!reach.plateau(input.well))\n\t\t text(fin.time, .5*(plateau-max.log.OD)+max.log.OD, paste(\"(\", round(plateau - max.log.OD,3), \")\", sep = \"\"), col = \"grey\", cex = 1.5*scale, pos = 2) \n }\n }\n\t\n########################################################################\n# Draw residuals from the nonlinear fit with option for lowess line #\n########################################################################\n#\t\t\nplot.residuals = function(well, xlim = NULL, lowess = T){\n\n\tdata = data.from(well, remove = F, remove.tanking = F)\n\n\tif (is.null(xlim))\n\t\txlim = c(min(data$Time, 0)-1, max(data$Time))\n\t\t\n\tplot(data.from(well)[,1], model.residuals(well), main = paste(plate.name(well), well.name(well), \"\\n[Residuals]\"),\n\t\txlab = \"Time(hours)\", ylab = paste(\"Residual\", names(data)[2]), xlim = xlim)\n\t\n\tabline(0,0, lty = 2)\n\t\n\tif (lowess)\n\t\tlines(lowess(data.from(well)[,1], model.residuals(well)), lw = 2, col = \"red\")\n\t}\n\n##############################################################################\n# This function is used to create a heatmap using: \n# specific growth, total growth, and lag time\n# for each well on a plate.\n#\n# @params\n# fitted.well.array: matrix containing well array object data\n# attribute: the data type we should use to create a heatmap\n# @returns\n# path of heatmap pdf file\n##############################################################################\ncreate.heatmap = function(fitted.well.array, attribute, unlog=NULL){\n attr.name <- deparse(substitute(attribute))\n pdf.name <- \"\"\n if(class(fitted.well.array) == \"matrix\"){\n #We may want to sub() out periods from plate.ID if it causes problems\n plate.ID = unique(unlist(aapply(fitted.well.array,plate.name)))[1]\n if(is.null(unlog)) {\n spec.growth = unlist(aapply(fitted.well.array, attribute))\n }\n # currently only total growth needs to be unlogged if unlog == T\n else {\n spec.growth = unlist(aapply(fitted.well.array, attribute))\n }\n num.dig = 3 #how many digits should be put on pdf?\n max = round(max(spec.growth, na.rm=T), digits=num.dig)\n min = round(min(spec.growth, na.rm=T), digits=num.dig)\n avg = round(mean(spec.growth, na.rm=T), digits=num.dig)\n heat.text = paste(toupper(sub(\"\\\\.\", \" \", attr.name)), \":\\n\", plate.ID, \"\\n\",\n paste(\"Max:\",max ,\"Min:\" ,min ,\"Avg:\", avg, sep=\"\"))\n \n attr.name <- sub(\"\\\\.\", \"_\", attr.name) #do not want periods in file path\n letters <- attr(fitted.well.array, \"dimnames\")[[1]]\n for(i in 1:length(letters)) letters[i] = paste(\" \", letters[i], \" \")\n nums <- attr(fitted.well.array, \"dimnames\")[[2]]\n for(i in 1:length(nums)) nums[i] = paste(\" \", nums[i], \" \")\n heat <- matrix(spec.growth, nrow=dim(fitted.well.array)[1], ncol=dim(fitted.well.array)[2], dimnames=list(letters,nums))\n pdf.name <- paste(getwd(), \"/\", plate.ID, \"_\", attr.name, \".pdf\", sep=\"\")\n pdf(pdf.name)\n #heatmap(heat, Rowv=NA, Colv=NA, revC=T, scale=\"none\", na.rm=T, main=plate.ID, col=rainbow(100), margins=c(6,6))\n #mtext(paste(\"Max:\", round(max(spec.growth, na.rm=T), digits=4),\"Min:\", round(min(spec.growth, na.rm=T), digits=4), \"Avg:\", round(mean(spec.growth, na.rm=T), digits=4)), side=1, line=3)\n pheatmap(heat, color=colorpanel(100, \"red\", \"orange\", \"yellow\"),\n border_color=\"black\", cell_width=2, cell_height=3,\n cluster_rows=F, cluster_cols=F, scale='none', main=heat.text, fontsize=16)\n dev.off()\n }\n else {\n return(\"Error\") \n }\n return(pdf.name)\n}\n\n########################################################################\n# Draw grids of 96 points as a visual representation of fit status, #\n# and other info for an array of fitted well objects, plate by plate #\n########################################################################\n#\n\t\nplate.overview = function(fitted.well.array, scale = 1){\n \n \n # Start with a list of the unique plate names in the fitted well array \n # and an appropriately-sized grid of coordinates to plot wells on.\n\tplates = unique(unlist(aapply(fitted.well.array,plate.name)))\n \n\tgrid = data.frame(x = rep(rep(1:plate.ncol, each = plate.nrow), length(plates)),\n\t\t\t y = rep( rep(-(1:plate.nrow), times = plate.ncol), length(plates)))\n\n\n # Gather information on each well to display on each of the coordinates in :\n # - was it marked as empty in the plate layout?\n # - did the program find it to contain no growth (\"dead\")? \n # - was the fitting procedure successful? \n # - did the curve tank? if so, at what timepoint? if not, or if the curve was marked as dead anyway, do not display the value. \n # - does the \"additional info\" slot indicate that any points were removed or the whole well discarded?\n \n\tempty = unlist(aapply(fitted.well.array, is.empty))\n\tdead = unlist(aapply(fitted.well.array, lacks.growth))\n\tfit = unlist(aapply(fitted.well.array, contains.fit))\n \n\ttanking = unlist(aapply(fitted.well.array, tanking.start))\n\ttanking[is.na(tanking) | tanking == 1 | dead] = \"\"\n\n\terrors = unlist(aapply(fitted.well.array, function(well){\n\t\tif (length(well@add.info) == 0)\n\t\t\t\"\"\n\t\telse if (grepl(\"removed\", well@add.info))\n\t\t\t\"-\"\n else if (grepl(\"detected\", well@add.info))\n \"+\"\n\t\telse if (grepl(\"discarded\", well@add.info))\n\t\t\t\"!\"\n\t\telse\n\t\t\t\"\"\n\t\t}))\n\t\n # Color and plotting character vectors (length = the number of wells in the array)\n # Default = 1 (open point, black)\n\tcolors = char = rep(1, length(tanking))\n\n # Desired colors\n colors[empty & dead] = \"green3\" # Empty well with no growth.\n colors[!empty & fit] = \"blue\" # Inoculated well with successfully fitted growth curve.\n \n # Undesired colors \n colors[empty & !dead] = \"darkolivegreen4\" # Inoculated well with some growth. \n colors[!empty & !fit] = \"red\" # Inoculated well with no successfully fit (either no growth or unsuccessful fit).\n \n char[!dead & fit] = 19 # Filled points for non-empty wells with successful fits \n char[!dead & !fit] = 4 # an X for non-empty wells with failed fits. \n \n char[errors == \"!\"] = 8 # Asterisk for discarded wells. \n char[errors == \"-\" & dead ] = 5 # Open diamond for empty wells (after removing points).\n char[errors == \"-\" & !dead & fit] = 23 # Filled diamond for non-empty wells with removed points and successful fits. \n char[errors == \"-\" & !dead & !fit] = 8 # Asterisk for wells with removed points and failed fits.\n \n \n\tfor (plate in 1:length(plates)){\n \n\t\tindices = (plate - 1) * plate.nrow*plate.ncol + 1:(plate.nrow*plate.ncol)\n\n # Plot the grid using colors and plotting characters determined above. \n\t\tplot(grid[indices,], col = colors[indices], bg = colors[indices], pch = char[indices], \n\t\t\t main = plates[plate], mex = scale, cex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, \n\t\t\t xaxt = \"n\", yaxt = \"n\", xlim = c(-0.5,plate.ncol + 0.5), ylim = c(-(plate.nrow + 1.5), 0.5), xlab = \"\", ylab = \"\")\n \n # Symbol legends\n \n legend.xpos = (c(-1,2.75,6.5,6.86,10.25)+0.5)*(plate.ncol+1)/13 - 0.5\n legend.ypos = -(plate.nrow + 0.5)\n \n legend(x=legend.xpos[1], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Empty, no growth\",\"Empty with growth\"),\n pch = c(1,19),\n pt.bg = c(\"green3\",\"darkolivegreen4\"),\n col = c(\"green3\",\"darkolivegreen4\")\n )\n legend(x=legend.xpos[2], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Inoculated with growth\", \"Inoculated, no growth\"),\n pch = c(19,1),\n pt.bg = c(\"blue\",\"red\"),\n col = c(\"blue\",\"red\")\n )\n legend(x=legend.xpos[3], y= legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Well tanks at specified index\", \"Some points removed\"),\n pch = c(21,23),\n pt.bg = c(\"grey\",\"grey\"),\n col = c(\"black\",\"black\")\n ) \n \n text(x=legend.xpos[4], y=legend.ypos - 0.29,\"#\",cex=0.5*scale) \n \n legend(x=legend.xpos[5], y=legend.ypos, cex = 0.7 * scale, y.intersp = 1.5, bty=\"n\", \n legend=c(\"Model fitting failed\", \"Well discarded\"),\n pch = c(4,8),\n pt.bg = c(\"black\",\"black\"),\n col = c(\"black\",\"black\")\n )\n \n # Add tanking indices if any were found. \n \ttext(grid[indices,] + 0.30, cex = 0.75*scale, \n labels = tanking[indices], col = colors[indices])\n\n # Label rows and columns\n\t\ttext(-1, -1:-plate.nrow, pos = 4, LETTERS[1:plate.nrow], cex = scale) \n\t\ttext( 1:plate.ncol, 0 , 1:plate.ncol, cex = scale) \t\n\t\t}\n \n\t}\n\n########################################################################\n# Draw each well in an array of fitted well objects in succession. #\n# Include options for adding notations, text info and fit parameters. #\n########################################################################\n#\n\nview.fit = function(fitted.data, indices = 1:length(fitted.data), \n unlog = F, constant.added, xlim = NULL, ylim = NULL, display.legend = T, \n\t\t show.text = T, show.calc = T, draw.guess = NULL, draw.symbols = F, number.points = T, \n\t\t\tuser.advance = T, show.residuals = F, scale = 1,...){\n\n if(!is.array(fitted.data))\n fitted.data = list(fitted.data)\n\n # Determine the boundaries for the axes (if user did not specify them)\n if(is.null(ylim)){\n min.y = min(unlist(aapply(fitted.data, function(well){\n \tif (unlog) well@use.log = F\n \tmin.y = min(data.from(well, remove = F, remove.tanking = F)[,2], na.rm = T)\n min(min.y, well@fit.par$b)\n \t})))\n max.y = max(unlist(aapply(fitted.data, function(well){\n if (unlog) well@use.log = F\n \tmax.y = max(data.from(well, remove = F, remove.tanking = F)[,2], na.rm = T)\n max(max.y, well@fit.par$b + well@fit.par$A)\n \t})))\n ylim = c(min.y, min.y + (max.y-min.y)*1.15) - unlog*constant.added\n }\n if(is.null(xlim)){\n min.x = min(unlist(aapply(fitted.data, function(well){\n min(data.from(well, remove = F, remove.tanking = F)[,1], na.rm = T)\n \t})))\n max.x = max(unlist(aapply(fitted.data, function(well){\n \tmax(data.from(well, remove = F, remove.tanking = F)[,1], na.rm = T)\n \t})))\n xlim = c(min.x - 0.05 * (max.x-min.x), max.x)\n }\n \n # Display a figure legend\n if(display.legend){\n well.fit.legend(xlim=xlim,ylim=ylim,scale=scale,constant.added=constant.added)\n if(user.advance){\n prompt = readline(\" to continue or Q to quit >>\")\n if (toupper(prompt) == \"Q\") break\n }\n }\n # Start to cycle through the wells \n\twell.number = 1\n\twhile (well.number <= length(fitted.data)) {\t\t\n\t\t# Only show wells specified by (default all wells)\n if (well.number %in% indices){ \n # plot the well\n fitted.well = fitted.data[[well.number]]\n plot(x=fitted.well, constant.added = constant.added, xlim = xlim, ylim = ylim,\n unlog = unlog, well.number = well.number, scale = scale, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, ...)\n \n if(user.advance)\n cat(\"\\n[\", well.number, \"] \", plate.name(fitted.well), \" \", well.name(fitted.well), \".\", sep = \"\")\n \n if (show.residuals & is.numeric(model.residuals(fitted.well))){\n if(user.advance)\n if (toupper(readline(\" for residuals >>\")) == \"Q\") break\n plot.residuals(fitted.well)\n }\n \n # Allow user to advance the currently shown well if specified. \n\t\t\tif (user.advance){\n \n\t\t\t\tprompt = readline(\" to continue, or type # of next well or Q to quit >>\")\n\t\t\t\tif (toupper(prompt) == \"Q\") break\n\n user.input = suppressWarnings(try(as.numeric(prompt),silent=T))\n \n # Go onto the next well unless input is a number. \n\t\t\t\tif (is.numeric(user.input) & !is.na(user.input) & length(user.input) > 0)\n\t\t\t\t\twell.number = user.input - 1\n\t\t\t }\n\t\t\t}\n # Advance the loop\n well.number = well.number + 1\n\t\t}\t\t\n\t}\t\n\n\nwell.fit.legend = function(xlim, ylim, scale = 1, constant.added){\n par(mar = c(5, 4, 4, 5)+0.1)\n plot(0,0, main = \"[Index] \\n; \",\n xlim = xlim, ylim = ylim, xlab = \"Time\", ylab = \"log(OD - blank + const)\", \n mex = scale, cex.main = 1.5*scale, cex.axis = 1.2*scale, cex.lab = 1.2*scale, type = \"n\")\n \n # Draw a second vertical axis, showing unlogged OD scale\n min.OD = unlog(ylim[1],constant.added)\n max.OD = unlog(ylim[2],constant.added)\n OD.labels = seq(from = min.OD, to = max.OD, length.out = 5)\n OD.labels = round(OD.labels,1)\n OD.at = log(OD.labels+constant.added)\n axis(side=4, at=OD.at, labels=OD.labels, cex.axis = 1.2*scale, cex.lab = 1.2*scale)\n mtext(4, text = \"OD - blank\", line = 3, cex=1.2)\n \n # Sample max. slope line\n abline(v=min(xlim)+0.5*max(xlim), col=\"blue\", lty=2)\n text(mean(xlim),min(ylim)+0.4*diff(ylim),labels=\"Maximum specific\\ngrowth rate\",col=\"blue\",pos=2,cex=0.75*scale)\n \n # Sample plateau line\n abline(h=min(ylim)+0.8*diff(ylim),col=\"forestgreen\",lty=2)\n text(min(xlim)+0.9*diff(xlim),ylim+0.8*diff(ylim),labels=\"Growth plateau\",col=\"forestgreen\",pos=3,cex=0.75*scale)\n\n # Sample max.log.OD line\n abline(h=min(ylim)+0.7*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.9*diff(xlim),ylim+0.7*diff(ylim),labels=\"max.log.OD\",col=\"black\",pos=1,cex=0.75*scale)\n \n # Sample inoc.log.OD\n abline(h=min(ylim)+0.1*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.1*diff(xlim),min(ylim)+0.1*diff(ylim),labels=\"Fitted growth\\nat inoculation\",col=\"black\",pos=3,cex=0.75*scale)\n \n # Sample baseline\n abline(h=min(ylim)+0.05*diff(ylim),col=\"red\",lty=2)\n text(min(xlim)+0.1*diff(xlim),min(ylim)+0.05*diff(ylim),labels=\"Baseline\",col=\"red\",pos=1,cex=0.75*scale)\n\n # Sample lag time\n lines(min(xlim)+c(0.1,0.25,0.50)*max(xlim),min(ylim)+c(0.05,0.05,0.4)*diff(ylim),col=\"red\",lty=2)\n text(min(xlim)+0.25*max(xlim),min(ylim)+0.05*diff(ylim),labels=\"Lag time\",col=\"red\",pos=1,cex=0.75*scale)\n \n # Sample achieved growth\n lines(min(xlim)+c(0.75,0.75)*max(xlim),min(ylim)+c(0.1,0.7)*diff(ylim),col=\"black\",lty=3)\n text(min(xlim)+0.75*max(xlim),min(ylim)+0.3*diff(ylim),labels=\"Achieved growth\",col=\"black\",cex=0.75*scale)\n \n # Sample plateau - achieved growth\n lines(min(xlim)+c(0.75,0.75)*max(xlim),min(ylim)+c(0.7,0.8)*diff(ylim),col=\"grey\",lty=3)\n text(min(xlim)+0.75*max(xlim),min(ylim)+0.75*diff(ylim),labels=\"Projected minus achieved growth\",col=\"grey\",cex=0.75*scale)\n \n # Symbol legend\n legend(x=\"right\", title = \"Timepoint Symbols\", legend = c(\"Normal point\", \"Ignored point\"),\n cex = 0.75*scale, pt.cex = c(0.6,0.6)*scale, pch = c(35,35), col=c(\"black\",\"gray80\"),\n x.intersp=1, xjust = 1, y.intersp=1.5)\n}\n\npdf.by.plate = function(fitted.data, out.prefix = \"\", upload.timestamp = NULL, \n out.dir = getwd(), unlog = F, constant.added, silent = T, overview.jpgs = T, ...){\n \n # Prepare timestamp for addition to output file names. \n filename.timestamp = strftime(upload.timestamp, format=\"_%Y-%m-%d_%H.%M.%S\")\n \n # Start file list with the overview pdf\n file.list.out = paste(out.dir,\"/\",out.prefix, \"_overview\", filename.timestamp, \".pdf\",sep=\"\")\n \n # Call to draw a graphic representation of each plate in this file. \n pdf(file.list.out, title = paste(out.prefix, \"plate overview\"))\n plate.overview.out = try(plate.overview(fitted.data),silent=T)\n if(class(plate.overview.out) == \"try-error\")\n stop(\"Error in : \", plate.overview.out)\n \n # Close devices\n while(dev.cur() != 1)\n dev.off() \n \n # Cycle through each plate \n for(i in 1:dim(fitted.data)[3]){\n \n # Get plate ID and position in data array.\n plate.ID = dimnames(fitted.data)[[3]][i]\n plate.indices = (i-1) * plate.nrow*plate.ncol + 1:(plate.nrow*plate.ncol)\n if(overview.jpgs){\n # most be > 1 to partition value breaks for heatmap\n well.matrix <- aapply(fitted.data[,,i], max.spec.growth.rate) \n num.wells <- length(well.matrix[!sapply(well.matrix, is.na)])\n if(num.wells > 1){\n #Heatmap block##########################################################\n #alongside the jpgs file create 3 heatmaps for each plate. NWD\n spec.heat.file = create.heatmap(fitted.data[,,i], max.spec.growth.rate)\n if(spec.heat.file == \"Error\")\n stop(\"Error in for specific growth\")\n lag.heat.file = create.heatmap(fitted.data[,,i], lag.time)\n if(lag.heat.file == \"Error\")\n stop(\"Error in for lag time\")\n total.heat.file = create.heatmap(fitted.data[,,i], achieved.growth)\n if(total.heat.file == \"Error\")\n stop(\"Error in for total growth\")\n # Add name of file if successfully written to file list output. Including heatmap files NWD\n file.list.out = c(file.list.out, spec.heat.file, lag.heat.file, total.heat.file)\n ########################################################################\n }\n jpg.name = paste(out.dir, \"/\", plate.ID, \"_overview\", \".jpg\", sep=\"\")\n jpeg(jpg.name, quality = 90, width = 600, height = 480)\n plate.overview.out = try(plate.overview(fitted.data[,,i]),silent = T)\n if(class(plate.overview.out) == \"try-error\")\n stop(\"Error in : \", plate.overview.out)\n }\n else\n jpg.name = c()\n \n # Open a separate PDF for each plate.\n if(!silent) cat(\"\\nprinting PDF for\", plate.ID)\n pdf.name = paste(out.dir, \"/\", plate.ID, \"_plots\", filename.timestamp, \".pdf\", sep=\"\")\n pdf(pdf.name, title = paste(\"R Graphics output for plate\", plate.ID))\n \n # Call to draw each well on the plate to the pdf. \n view.fit.out = try(view.fit(fitted.data, indices = plate.indices, unlog=unlog, constant.added=constant.added, user.advance=F,...),silent=T) \n \n if(class(view.fit.out) == \"try-error\")\n stop(\"Error in : \", view.fit.out)\n\n # Close all devices\n while(dev.cur() != 1)\n dev.off() \n \n if(!silent) cat(\"...done!\\n\\twritten to\", pdf.name, \"\\n\") \n file.list.out = c(file.list.out, jpg.name , pdf.name)\n }\n return(file.list.out)\n } \n",
- "created" : 1423869024830.000,
- "dirty" : false,
- "encoding" : "UTF-8",
- "folds" : "",
- "hash" : "197404081",
- "id" : "7189AD56",
- "lastKnownWriteTime" : 1423869067,
- "path" : "~/Documents/GCAT4_old/trunk/R/GCAT/R/plot.fit.R",
- "project_path" : "plot.fit.R",
- "properties" : {
- },
- "source_on_save" : false,
- "type" : "r_source"
-}
\ No newline at end of file
diff --git a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/85C00847 b/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/85C00847
deleted file mode 100644
index 1d1225a..0000000
--- a/R/GCAT/R/.Rproj.user/7CB73FCA/sdb/per/t/85C00847
+++ /dev/null
@@ -1,16 +0,0 @@
-{
- "contents" : "#Copyright 2012 The Board of Regents of the University of Wisconsin System.\n#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, \n#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman\n#\n#This file is part of GCAT.\n#\n#GCAT is free software: you can redistribute it and/or modify\n#it under the terms of the GNU Lesser General Public License as published by\n#the Free Software Foundation, either version 3 of the License, or\n#(at your option) any later version.\n#\n#GCAT is distributed in the hope that it will be useful,\n#but WITHOUT ANY WARRANTY; without even the implied warranty of\n#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n#GNU Lesser General Public License for more details.\n#\n#You should have received a copy of the GNU Lesser General Public License \n#along with GCAT. If not, see .\n\n# Notes by Jason\n# 9/07/11\n\n\n########################################################################\n# #\n# Function for loading data from tabular format into an object array #\n# #\n########################################################################\n#' Load tabular data\n#'\n#' This function handles loading data from tabular format (.csv, tab-delimited text or R data frame object)\n#' and returns an array of well objects, each filled with raw Time vs. OD data. \n#' It takes single-plate or multiple-plate format data. For single-plate data, \n#'it calls on the function \\code{gcat.reorganize.single.plate.data} to rearrange the table before creating the output object. \n#'\n#' @param file.name Complete path and file name of a comma-separated values (.csv) file containing growth curve data \n#' in the multiple-plate (long) format. \n#' @param input.data A list of tables representing input files read with \\code{read.table}. Used to save time in cases\n#' of running multiple analyses on the same dataset. If used, the function will ignore \\code{file.name} entirely.\n#' @param load.type .csv by default. \n#' @param plate.laout Specifies the layout of the given plate.\n#' @param single.plate.ID specifies a plate name for a single-plate read. If NULL, this is derived from the file name. \n#' @param blank.value Blank OD measurement for uninoculated wells. By default(NULL), the value of the first OD\n#'measurement in each well is used.\n#' @param add.constant A value for r in the log(OD + r) transformation.\n#' @param plate.nrow The number of rows in the input files.\n#' @param plate.ncol The number of columns in the input files.\n#' @param input.skip.lines If specified, this number of lines shall be skipped from the top when reading the input file with read.csv \n#' @param multi.column.headers The headers of the column when analyzing multiple plates.\n#' @param single.column.headers The headers of the column when analyzing a single plate.\n#' @param layout.sheet.headers The headers of the layout file.\n#' @param growth.model What growth model should be used?\n#' @param silent Surpress all messages.\n#' @param verbose Display all messages when analyzing each well.\n#' \n#' @return A list of well objects.\n\ngcat.load.data = function(file.name = NULL, load.type = \"csv\", input.data = NULL, single.plate.ID = NULL, \n plate.layout = NULL,plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0,\n multi.column.headers = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\"), single.column.headers = c(\"\",\"A1\"), \n layout.sheet.headers = c(\"Strain\", \"Media Definition\"),\n blank.value = NULL, start.index = 2, single.plate = F, silent = T){\n\n ########################################################################\n # Read from .csv, tab-delimited text file or data frame object #\n ########################################################################\n \n\tif(is.null(input.data)){\n\t\t# Either read from .csv.\n\t\tinput.data = read.csv(file.name, stringsAsFactors=F, skip = input.skip.lines, fileEncoding='UTF-8')\n\n\t\t# Determine single plate name if not specified. \n if (is.null(single.plate.ID)){\n # Split the file name by \".\" and discard the last member (file extension). \n single.plate.ID = strsplit(basename(file.name),\"\\\\.\")[[1]]\n single.plate.ID = paste(single.plate.ID[-length(single.plate.ID)],collapse=\".\")\n }\n\t\t}\n\n # Call to arrange data from a single plate format file\n\tif(single.plate)\n\t\tinput.data = gcat.reorganize.single.plate.data(input.data = input.data, single.column.headers,\n blank.value = blank.value, single.plate.ID = single.plate.ID,silent=silent)\n\n ########################################################################\n # Search for and standardize column headers in input data #\n ########################################################################\n \n # Go through the specified column headers, determining what their positions are in the \n # input data frame and if any are missing.\n \n # Get names of column headers in input data\n input.colnames = colnames(input.data)\n \n # Create a list denoting the column numbers in input data that match each of the specified column names, \n # and a separate list for any missing columns. \n \n column.matches = c()\n missing.list = NULL\n \n\tfor(i in 1:length(multi.column.headers)){\n\t\tif (multi.column.headers[i] %in% input.colnames)\n\t\t\tcolumn.matches[i] = min(which(input.colnames == multi.column.headers[i]))\n\t\t# Take the first column in input file that matches a specified column header.\n\t\telse{\n\t\t\tmissing.list = c(missing.list, i)\n\t\t}\n\t}\n\n # If any columns are missing, stop and report error with missing column names\n\tif (is.vector(missing.list)){\n\t\tmessage = \"The following columns:\"\n\t\tfor (i in missing.list)\n\t\t\tmessage = paste(message, paste(\" \", multi.column.headers[i]), sep = \"\\n\") \n\t\tstop(message, \"\\n were not found in the data file.\")\n\t\t}\n\n # Reorder and rename the columns, using the list of matching column numbers from above.\n\tinput.data = input.data[,column.matches]\n\tnames(input.data)[1:4] = c(\"Plate.ID\", \"Well\", \"OD\", \"Time\") \n\n # Use 'substring' to split the alphanumeric \"Well\" field into row (letters) and column (numbers)\n\tinput.data$Well.row = substring(input.data$Well, 0,1)\n\tinput.data$Well.col = as.numeric(substring(input.data$Well, 2))\n\n\n ########################################################################\n # Create an array of well objects with the Time and OD data #\n ########################################################################\n #\n # Use the by function to split up the data frame into shorter segments by well (row, column and plate)\n \n\twell.array = by(data = input.data[,c(\"OD\", \"Time\")], \n INDICES = list(input.data$Well.row,input.data$Well.col,input.data$Plate.ID), \n FUN = function(x){data.frame(Time=x$Time, OD=x$OD,stringsAsFactors = F)}, simplify = F)\n\n \n # Then apply the function (found in well.class) to each member to create a well object\n well.array = aapply(well.array,function(x){well(x$Time,x$OD)})\n\n # Differentiate any duplicate plate names in the array's dimnames \n\tnew.dimnames = dimnames(well.array)\n for (i in length(new.dimnames[[3]]):1){\n\t\tif (any(new.dimnames[[3]][-i] == new.dimnames[[3]][i]))\n\t\t\tnew.dimnames[[3]][i] = paste(\"another_\", new.dimnames[[3]][i], sep = \"\")\t\n\t\t}\n\tdimnames(well.array) = new.dimnames\n\t\n # Copy the plate/row/column names found in the dimnames into the array objects themselves (use \"position\" slot)\n\tfor(plate in unique(dimnames(well.array)[[3]])){\n\t\tfor (col in unique(dimnames(well.array)[[2]])){\n\t\t\tfor(row in unique(dimnames(well.array)[[1]])){\n\t\t\t\twell.array[[row,col,plate]]@position = c(plate=plate,row=row,col=col)\n\t\t\t\t}\n\t\t\t}\n\t\t}\n\t\t\n ########################################################################\n # Add plate layout information to well array #\n ########################################################################\n \n # Use the object to add media and strain information to the \"well.info\" slot of each well\n # Also set the value of the parameter in slot \"curve.par\" to T for empty wells. \n ########################################################################\n # Add plate layout information to well array #\n ########################################################################\n \n # Use the object to add media and strain information to the \"well.info\" slot of each well\n # Also set the value of the parameter in slot \"curve.par\" to T for empty wells. \n \n\n \n if(is.null(plate.layout)){ # If is not provided, do not add strain information, and assume all wells are inoculated. \n plate.layout = data.frame(Row=rep(PLATE.LETTERS[1:plate.nrow],plate.ncol),Column=rep(1:plate.ncol,each=plate.nrow),rep(\"Unknown Strain\",96),rep(\"Unknown Media\",96))\n colnames(plate.layout) = c(\"Row\", \"Column\", layout.sheet.headers)\n } \n else\n if(!silent) cat(\"\\n\\t\\tusing plate layout to fill well info.\")\n \n\tfor(plate in unique(dimnames(well.array)[[3]])){\n\t\tfor (col in unique(dimnames(well.array)[[2]])){\n\t\t\tfor(row in unique(dimnames(well.array)[[1]])){\n well = well.array[[row,col,plate]]\n # For each well on each plate, find the corresponding row in . \n # If refers to specific plates, then use those to find the correct row. Otherwise, generalize across all plates. \n if (\"Plate.ID\" %in% names(plate.layout)) layout.row.number = which(plate.layout$Column==well@position[\"col\"] & plate.layout$Row==well@position[\"row\"] & plate.layout$Plate.ID==well@position[\"plate\"] ) \n else layout.row.number = which(plate.layout$Column==well@position[\"col\"] & plate.layout$Row==well@position[\"row\"]) \n # Error if either no rows or more than one row matches the well\n if (length(layout.row.number) != 1)\n stop(\"incorrectly formatted plate layout! check names of columns, rows, and plates (if applicable).\")\n \n # Add any additional columns to the well's \"well.info\" slot\n well.info = plate.layout[layout.row.number,!(names(plate.layout) %in% c(\"Row\",\"Column\",\"Plate.ID\",layout.sheet.headers))]\n # Fix the column name issue if there is only one additional entry. \n if(length(well.info) == 1){\n well.info = data.frame(well.info,stringsAsFactors=F) \n names(well.info) = names(plate.layout)[!(names(plate.layout) %in% c(\"Row\",\"Column\",\"Plate.ID\",layout.sheet.headers))] \n } \n well@well.info = as.list(well.info)\n \n well@well.info$Strain = plate.layout[layout.row.number, layout.sheet.headers[1]]\n well@well.info$Media = plate.layout[layout.row.number, layout.sheet.headers[2]]\n # Set parameter in slot \"curve.par\" accordingly \n well@curve.par$empty.well = (plate.layout$Strain[layout.row.number] == \"Empty\") \n well.array[[row,col,plate]] = well\n\t\t\t\t}\n\t\t\t}\n\t\t}\n # Set start index value in each well\n well.array = aapply(well.array, function(x,start.index) { x@start.index = start.index; x }, start.index)\n \n ########################################################################\n # Return values to R #\n ######################################################################## \n # \n\t# Console output if desired, return the completed well array.\n\tif (!silent)\n\t\tcat(\"\\n\\t\", dim(well.array)[[3]], \"plates added to array from\", file.name)\n\treturn(well.array)\n\t}\n\n\n\n########################################################################\n# #\n# Reorganize data from single-plate input format before reading #\n# #\n########################################################################\n#\n# This function reorganizes the data frame from a single-plate format file. \n# input.data - data frame read straight from a single-plate format data file. \n# single.plate.ID - specifies a plate name for a single-plate read, since none is given in the single-plate format\n# The plate will be named Plate_1 unless otherwise specified. \n\ngcat.reorganize.single.plate.data = function(input.data, blank.value = NULL, single.column.headers, single.plate.ID = \"Plate_1\", silent=T){\n \n ########################################################################\n # Standardize the formatting and search for specified column names #\n ######################################################################## \n # \n # Locate the first and last rows from the table and return errors if not defined \n # Note: this only works if the time column is the first column\n \n\theader.row = min(which(input.data[,1] == single.column.headers[1])) \n if (length(header.row) != 1 | is.infinite(header.row))\n stop(\"could not locate header row in input file!\")\n \n # The last row: where column 2 starts to be blank, or the total number of rows, whichever is smaller \n extra.rows.start = min(which(input.data[-(1:header.row),2] == \"\"), which(is.na(input.data[-(1:header.row),2])), nrow(input.data[-(1:header.row),]))\n if (length(extra.rows.start) != 1 & is.infinite(extra.rows.start)) \n stop(\"could not locate last row in input file!\")\n\n # Use header row to rename the columns, then cut off extra rows (including the ones above header)\n\tnames(input.data) = as.character(unlist(input.data[header.row,]))\n input.data = input.data[(header.row+1):(header.row+extra.rows.start-1),]\n \n # Time column: allow for multiple matches to the name (since it's usually blank) but assume it's the first one\n\tTime.column = which(names(input.data) == single.column.headers[1])\n\tif (length(Time.column) != 1){\n if(!silent) cat(\"No unique time column in input.data file! Using the first one encountered.\")\n\t\tTime.column = min(Time.column)\n\t\t}\n # First well column (default A1): only allow for one match.\t\n\tWell.column.start = which(names(input.data) == single.column.headers[2])\n\tif (length(Well.column.start) != 1)\n\t\tstop(\"No unique start point for well columns in input.data file!\")\n\n # If the time column was found, rename it \"Time\" and reformat it into a numeric value\n # Adjust the blank measurement timestamp to -1 seconds if there is one\n\n names(input.data)[Time.column] = \"Time\"\n # Note: Some single plate screens have timepoints ending with \"s\" for seconds. \n # This line removes the \"s\" while maintaining general compatibility. \n\tinput.data$Time = unlist(strsplit(input.data$Time, \"s\"))\n\n # If is NULL (default - takes the first OD as the blank reading), then the first timepoint can labeled something non-numeric. \n # In that case, rename it to match the first real timepoint minus one. \n # when user input blank value, Blank timepoint i.e. input.data$Time[1] == Blank, labeled as \"Blank\" from data input file\n # It also should rename it to match the first real timepoint minus one. \n if(is.null(blank.value) || is.na(as.numeric(input.data$Time[1])))\n input.data$Time[1] = as.numeric(input.data$Time[2]) - 1\n \n ########################################################################\n # Start to fill the reformatted data frame #\n ######################################################################## \n \n # If all columns are present, make a list of all the wells.\n\twell.list = paste(rep(PLATE.LETTERS[1:plate.nrow], each = plate.ncol), rep(formatC(1:plate.ncol, digits = log(plate.ncol, 10), flag = \"0\"), plate.nrow), sep = \"\")\n\n #\tDuplicate the well names times the number of time measurements in each well\t\n Well = rep(well.list, each = length(input.data[,Time.column]))\t\n\t\t\n # Duplicate times the length of the entire output \n Plate.ID = rep(single.plate.ID, length(Well))\n\n # Duplicate the time column times the number of wells and convert to numeric\n\tTime = as.numeric(rep(input.data[,Time.column], times = length(Well.column.start:ncol(input.data))))\n\n # Append OD measurements from each well together and convert to numeric\n\tOD = c()\n\tfor (i in Well.column.start:ncol(input.data)){\n\t\tOD = as.numeric(c(OD, input.data[,i]))\n\t\tOD = unlist(OD)\n\t\t}\n\n # Fill and return the data frame containing the above four vectors.\n\toutput = data.frame(Plate.ID, Well, OD, Time)\t\n\t\n # Include any extra columns that were not Time or OD measurements?\n\tfor(i in (1:length(names(input.data)))[-c(Time.column,Well.column.start:ncol(input.data))]){\n\t\tnew.column = data.frame(rep(input.data[,i], length(Well.column.start:ncol(input.data))))\n\t\tnames(new.column) = names(input.data)[i]\n\t\toutput = cbind(output, new.column)\n\t\t}\t\n\treturn(output)\n}\n\n\n########################################################################\n# #\n# Function to combine two well array datasets by plate #\n# #\n########################################################################\n# ----------------------------------------------------------\n# This function can append together arrays created using \n# Arguments: any number of array objects as long as they are all output straight from