像我这种做Enzyme kinetic的砖工，每天都有大把时间耗在处理数据+作图上。过去的几年里，懒惰的我都是先从运行软件（UVProbe）上导出数据。由于该软件及其不好用，只能导出由机器二次拟合的Activity Data和在软件自带的处理界面上得到的斜率数据即Main Table。而对于原始数据，我捣鼓了很久都没办法找到直接从.kin文件中读取的方法，只能打开每个文件，分几组（一次run几个cuvette就有几组）打包得到Data Print数据，然后手动复制粘贴。得到斜率数据后，还得在Excel上打开，Separate into columns，计算真实速率。然后再打开SigmaPlot，使用其Enzyme Kinetic模块，通过手动输入许多组fixed Substrate / varied Inhibitor 的速率值后，得到拟合不同inhibition model的图片。如果要在图片上进行标记，那还要花费更多的时间。 每一组数据，大概都需要花费将近半小时的时间进行如上处理，感觉自己宝贵的生命就一点点地耗费在这毫无意义的重复劳动上了。

而前段时间，由于组里的HPLC无法导出图片，我只能转去研究出一套让软件每次运行后自动导出原始数据的方法，然后自己作图。期间对比了Matlab，R，Python等工具，觉得R还算是一个轻量化比较好用的工具，便大致看了下它的初学者文档，写了一个将当天run的每次HPLC分别做图的小script。每周把跑的HPLC数据放在一个文件夹里，然后打开script，在RStudio里轻轻一点击，所有图便出现在文件夹里。然后选出有意义的图，在imageJ中montage一下，就能打印出来跟老板talk了，感觉自己萌萌哒。

从那时以后，一个念头就在我脑中生出：为什么不用R做一个一键做kinetic 图的软件？构思如下：

先读取Main Table，因为我有给每个数据做标记的习惯，所以可以直接从Main Table中将速率与底物浓度对应起来，排除掉错误/ 失败数据并进行数据处理。同时根据其从打包的Data Print中提取每一组原始数据并导出。

然后将处理过的数据fit进不同的模型中，利用AIC值将各模型进行排序，选取最优值。

最后在不同Inhibitor浓度下作v0-S图， 并将得到的各种Parameter标记在图上。

由于只是既HPLC图之后第二次使用R，自己又很多年没有编程了，script写得磕磕绊绊，过程几乎就是不断试错。途中查阅了R的基本官方手册以及cookbook，作图时研究了好几个关于ggplot2包的博客，还看了许多stat.ethz.ch与stackoverflow上的案例，进展依旧十分缓慢，断断续续大概花了半个月才搞出一个运行结果比较好看，code也相对来说比较简洁的方案。

这其中，我在两个地方卡得时间最久：均超过一天。

一处是模型的拟合。有三种常用的简化模型，分别是

[code language=“r”] formulaComp <- as.formula(“v0 ~ (Vmax * S) / (Km * (1 + I / Kis) + S)“) formulaNoncomp <- as.formula(“v0 ~ (Vmax * S) / (Km * (1 + I / Kis ) + S * ( 1 + I / Kii)) “) formulaUncomp <- as.formula(“v0 ~ (Vmax * S) / (Km + S * (1 + I / Kii))“) [/code]

其中，仅有v0，S和I已知，其余均需要进行拟合。对于这种非线性模型，R中自带的有nls()函数：

Nonlinear Least Squares Description Determine the nonlinear (weighted) least-squares estimates of the parameters of a nonlinear model.

Usage nls(formula, data, start, control, algorithm, trace, subset, weights, na.action, model, lower, upper, …)

这个函数可以通过三种算法进行拟合：默认的Gauss-Newton算法， 对类线性模型适用的Golub-Pereyra算法以及Port包中采用的nl2sol有边界算法。实际使用中我发现默认算法可以轻易地拟合competitive model，但这个函数对另两种模型无效。于是只能另寻它法。

我先加载了nls2包，包内的nls2()函数适用与nls相同的算法，但是可以使用 “brute-force”, “random-search”, “plinear-brute”和”plinear-random”等附加选项，不明白这些分别有什么卵用的我在粗暴地试验下，发现这个函数并不管用。

Nonlinear Least Squares with Brute Force Description Determine the nonlinear least-squares estimates of the parameters of a nonlinear model.

Usage nls2(formula, data = parent.frame(), start, control = nls.control(), algorithm = c(“default”, “plinear”, “port”, “brute-force”, “grid-search”, “random-search”, “plinear-brute”, “plinear-random”), trace = FALSE, weights, …, all = FALSE)

后来又发现了minpack.lm包，包内的nlsLM()函数采用Levenberg-Marquardts算法，居然一次成功。

Standard ‘nls’ framework that uses ‘nls.lm’ for fitting Description nlsLM is a modified version of nls that uses nls.lm for fitting. Since an object of class ‘nls’ is returned, all generic functions such as anova, coef, confint, deviance, df.residual, fitted, formula, logLik, predict, print, profile, residuals, summary, update, vcov and weights are applicable.

Usage nlsLM(formula, data = parent.frame(), start, jac = NULL, algorithm = “LM”, control = nls.lm.control(), lower = NULL, upper = NULL, trace = FALSE, subset, weights, na.action, model = FALSE, …)

另一处则是各种Parameter的标记。为了得到类似Km = 1.11 ± 0.11 µM一开始的试验环节，我简单粗暴地使用手动重复的annotation，效果拔群

[code language=“r”] results.comp <- summary(comp)$parameters[,1:2] annotation.comp <- annotation_custom(grobTree( textGrob(bquote(paste(K[m] == .(format(results.comp[1,1], digits = 3)) %+-% .(format(results.comp[1,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.2, hjust=0), textGrob(bquote(paste(K[is] == .(format(results.comp[2,1], digits = 3)) %+-% .(format(results.comp[2,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.15, hjust=0), textGrob(bquote(paste(V[max] == .(format(results.comp[3,1], digits = 3)) %+-% .(format(results.comp[3,2], digits = 2)) ,’ ‘,mu,M/s)), x=0.5, y=0.1, hjust=0), gp=gpar(col=“black”, fontsize=13, fontface=“italic”))) # Results

results.noncomp <- summary(noncomp)$parameters[,1:2] annotation.noncomp <- annotation_custom(grobTree( textGrob(bquote(paste(K[m] == .(format(results.noncomp[1,1], digits = 3)) %+-% .(format(results.noncomp[1,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.25, hjust=0), textGrob(bquote(paste(K[is] == .(format(results.noncomp[2,1], digits = 3)) %+-% .(format(results.noncomp[2,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.2, hjust=0), textGrob(bquote(paste(K[ii] == .(format(results.noncomp[3,1], digits = 3)) %+-% .(format(results.noncomp[3,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.15, hjust=0), textGrob(bquote(paste(V[max] == .(format(results.noncomp[4,1], digits = 3)) %+-% .(format(results.noncomp[4,2], digits = 2)) ,’ ‘,mu,M/s)), x=0.5, y=0.1, hjust=0), gp=gpar(col=“black”, fontsize=13, fontface=“italic”))) # Results

results.uncomp <- summary(uncomp)$parameters[,1:2] textGrob(bquote(paste(K[m] == .(format(results.uncomp[1,1], digits = 3)) %+-% .(format(results.uncomp[1,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.2, hjust=0), textGrob(bquote(paste(K[ii] == .(format(results.uncomp[2,1], digits = 3)) %+-% .(format(results.uncomp[2,2], digits = 2)) ,’ ‘,mu,M)), x=0.5, y=0.15, hjust=0), textGrob(bquote(paste(V[max] == .(format(results.uncomp[3,1], digits = 3)) %+-% .(format(results.uncomp[3,2], digits = 2)) ,’ ‘,mu,M/s)), x=0.5, y=0.1, hjust=0), gp=gpar(col=“black”, fontsize=13, fontface=“italic”))) # Results [/code]

但是，当我想将Km等作为变量调用，这个方法就一直报错。后来研究了quote, bquote, expression, substitue等得用法后，发现可以将简单地表达式以字符串形式存储在变量内，再通过parse调用变量，呼出表达式，再使用bquote可以直接进行变换，简化后的代码如下：

[code language=“r”] paraMatrix <- data.frame(para = c(“K[m]“, “K[i]“, “K[ii]“, “K[is]“, “V[max]“), unit = c(“mu*M”, “mu*M”, “mu*M”, “mu*M”, “mu*M/s”), row.names = c(“Km”, “Ki”, “Kii”, “Kis”, “Vmax”)) results <- summary(fittingModel)$parameters[,1:2] numPara <- length(results[,1])

for (j in 1:numPara) { t <- row.names(results)[j] para <- parse(text = as.character(paraMatrix[t,1]))[[1]] unit <- parse(text = as.character(paraMatrix[t,2]))[[1]] p <- p + annotation_custom(grobTree(textGrob(bquote(.(para) == .(format(results[j,1], digits = 3)) %+-% .(format(results[j,2], digits = 2)) ~ .(unit)), x=0.5, y=0.05 * (numPara+2-j), hjust=0), gp=gpar(col=“black”, fontsize=13, fontface=“italic”))) } [/code]

在fail & trial 的过程中，我发现R在数据拟合方面，模型还是很多的。但是在一些处理细节方面，则有一点混乱。比如bquote和substitute两个函数，虽然调用方式有一些区别，但效果基本上是一样的。究其愿意，大抵是bquote是LISP的遗留物云云。如果能统一一下就更好了。不过我现在对bquote等的适用场景还有一点茫然，改天完全搞明白了再予总结。