介绍

SERRF是 U.C Davis Fiehn Lab的S.L Fan开发的一款QC-based数据标准化工具，用于校正非靶代谢组学数据的系统误差，具有出色的能力，校正后features RSD显著降低。目前使用方式有3种：

通过网页工具https://slfan2013.github.io/SERRF-online/# 进行上传数据分析；
通过本地运行shiny app，需要在将仓库克隆到本地，然后用Rstudio打开脚本运行run app;
通过本地运行normalize function运行。

emmm... 经过尝试，第三种运行比较顺利。但是文件准备设置输出等步骤还是太麻烦，所以打算对SERRF方法单独拿出来进行重构。并且接入Tidymass的massdataset，加入分析流程。

脚本结构分析

Tidymass本身自带了很多数据标准化的方法，所以像SERRF脚本中的其他标准化算法先不重构。主要对SERRF这种基于随机森林的标准化算法进行重构，当然方向也是趋于tidyverse的风格。首先看看这个脚本的步骤和逻辑。

数据清洗

首先对样品信息的数据进行了数据清洗，最后转换为data.table格式。p
然后对feature的信息进行了清洗，同样转换为data.table形式。f
最后对代谢物积累矩阵的标准化, 最后转换为matrix形式。e

我们观察了SERRF的数据格式，和tidymass的有很大的相似之处。所以打算将输入数据改成tidymass的mass_dataset

首先尝试都做data.frame处理

SERRF算法实现过程

整体将不同的batch分开计算。所有的迭代也是建立在batch个数上的

首先第一步把dataset按照batch拆开, 按batch计算训练集和目标集的spearman相关。
第二步根据feature个数进行迭代，对每个feature分析建模做标准化。
- 首先从第一步计算的训练集和目标集相关性矩阵中分别将对应feature的相关性结果提取出来；
- 对该结果进行排序，不考虑相关性的正负属性，只要是高度相关就排在前面corr_train_order = order(abs(corr_train[,j]),decreasing = TRUE)；
- 第三步，从整体表达矩阵中筛选测试集，简单讲就是求第二步两个相关性向量的交集，但是有条件，具体方法如下：
  - 首先设定一个交集数量的上限值num，脚本这里默认10，然后引入一个空置的参数筛选出的变异值sel_var对该值进行迭代计算，迭代的条件是sel_var中元素的数量 >= 0的时候终止, 也就是说比对训练集和目标集相关性向量的排序表，从各自的top10个开始做交集，每次迭代增加一个材料（第二次1:11，第三次1:12 ....）直到交集的数量大于了我们刚设定的10这个阈值就终止。
  - 这样我们就确定了用作训练数据的材料（最终的交集）。
- 首先训练集会对该batch下的所有qc做一个中心化scale(x,center = T,scale = F), 然后对交集材料做标准化；
- 对应测试集只做标准化，这样我们就得到了这10个材料在qc和sample各自的标准化表；后续会检查有是否存在这些feature在sample是否都缺失了，只有在训练集和测试集中都不含NA的feature才会被保留；
- 最终将去中心化的train_data_x和标准化的train_data_y合并为新的训练集；
- 用ranger构建模型；
- 吐槽一下，这baseR的代码直接给我看吐了。
- 分解下下面鬼见愁就的代码（针对训练集的qc）具体标准化过程如下：

norm[train.index_current_batch=='qc'] = e_current_batch[j, train.index_current_batch=='qc']/((predict(model, data = train_data)$prediction+mean(e_current_batch[j,train.index_current_batch=='qc'],na.rm=TRUE))/mean(all[j,sampleType.=='qc'],na.rm=TRUE))
norm[train.index_current_batch=='qc'] = norm[train.index_current_batch=='qc']/(median(norm[train.index_current_batch=='qc'],na.rm=TRUE)/median(all[j,sampleType.=='qc'],na.rm=TRUE))

  - 首先通过`predict`用构建的模型对`train_data`进行预测，取得观测值后加上该`batch`中`feature`所有QC样本峰面积的均值，然后用该数值除以所有QC样本峰面积的均值，得到一个比值，最后用该`batch`中`feature`所有QC样本峰面积对应的除去该比值实现对`QC`样本的的标准化。

$QC_{norm} = \frac {QC_{raw}}{\frac {QC_{model}+QC_{mean}}{ALL_{mean}}}$

  - 得到标准化后的$QC_{norm}$后，用该值的中位值除以所有`QC`样品的中位值，然后每个值再除以这个比值

$QC_{norm} = \frac {QC_{norm}} {\frac {QC_{norm.median}}{ALL_{median}}}$

下面分解对训练集中sample标准化的部分

norm[!train.index_current_batch=='qc'] =(e_current_batch[j,!train.index_current_batch=='qc'])/((predict(model,data = test_data)$predictions  + mean(e_current_batch[j, !train.index_current_batch=='qc'],na.rm=TRUE))/(median(all[j,!sampleType.=='qc'],na.rm = TRUE)))
norm[!train.index_current_batch=='qc'][norm[!train.index_current_batch=='qc']<0]=e_current_batch[j,!train.index_current_batch=='qc'][norm[!train.index_current_batch=='qc']<0]
norm[!train.index_current_batch=='qc'] = norm[!train.index_current_batch=='qc']/(median(norm[!train.index_current_batch=='qc'],na.rm=TRUE)/median(all[j,!sampleType.=='qc'],na.rm=TRUE))
norm[!is.finite(norm)] = rnorm(length(norm[!is.finite(norm)]),sd = sd(norm[is.finite(norm)],na.rm=TRUE)*0.01)
out = boxplot.stats(norm, coef = 3)$out
norm[!train.index_current_batch=='qc'][norm[!train.index_current_batch=='qc']%in%out] = ((e_current_batch[j,!train.index_current_batch=='qc'])-((predict(model,data = test_data)$predictions  + mean(e_current_batch[j, !train.index_current_batch=='qc'],na.rm=TRUE))-(median(all[j,!sampleType.=='qc'],na.rm = TRUE))))[norm[!train.index_current_batch=='qc']%in%out];
norm[!train.index_current_batch=='qc'][norm[!train.index_current_batch=='qc']<0]=e_current_batch[j,!train.index_current_batch=='qc'][norm[!train.index_current_batch=='qc']<0]

  - 第一步和上面一样同样的算法；第二步是探测了标准化后的`sample`值，如果有`＜0`的值，用原始的值填充；
  - 接下来和上面一样用标准化的值除以两中位值的比值；
  - 下一步检查结果中是否有无穷值，如果有的话用随机数填充，标准差选择非无穷值的标准差，最终的值`*0.01`；
  - 筛选出离群的值，然后填充离群值；

$Sample_{raw}-((QC_{model}+QC_{mean}) - ALL_{all.median})$

  - 先用`QC`模型+`QC`均值减去所有矩阵中非`QC`的中位值，然后用sample的原始峰面积减去该值。然后用该值填充离群值
  - 如果有引入的负数，再用上面的公式填充一下负数。

最终标准化完毕。
完成所有feature的标准化后会对NA和负数进行重新填充；
- 对于NA进行随机数填充，均值为其他非NA值的最小值，sd为其他非NA值的sd，最终得到的随机数*0.1；
- 对负数的填充，1*非负数值的最小值。
- 至此SERRF标准化过程完毕
交叉检验
- 使用的数据集是原始的qc峰面积数据。默认5-fold。为啥要做5折？
- 这里发现个错误，我们看下any函数的用法
- ，很明显作者在这里是想看所有batch中是否存在某个batch中的qc样本数量小于7，按照目前的代码any(table(p_qc$batch))总是会返回一个逻辑值TRUE，any本来就是个做判断的函数，所以把<7移动到括号内。这... 之前所有工作中ratio都被默认成了0.7. 0.7和0.8差异大吗？

## 原始的
if(any(table(p_qc$batch))<7){
  ratio = 0.7
}else{
  ratio = 0.8
}
## 改造后的
if(any(table(p_qc$batch)<7)){
  ratio = 0.7
}else{
  ratio = 0.8
}