@我的博客：有味
导读
液相色谱 - 质谱（LC-MS）因为可以检测的代谢物的覆盖范围广，灵敏度高和样品制备简单快捷而广泛运用于非靶向代谢组学研究中。如
血清、尿液和脑脊液等样本。然而，在大样本多批次数据采集的背景下，仪器本身随着时间的迁移，系统污染以及色谱柱活化再生过程的不当等一系列原因导致仪器对数据反映不稳定。这些不稳定性会直接影响测定的精确分子量、保留时间以及信号响应强度。尤其是在大样本，多批次样本采集背景下，数据产生的信号漂变会比较严重，一般来说，批次间的系统信号偏差会大于批次内样本间的信号偏差。这些测量误差降低了数据的重复性和再现性，因此可能降低检测生物反应和模糊解释的能力。
batchCorr包介绍
那么，R软件包
batchCorr就是在这样的一个背景下孕育而生的。其主要目的就是为了解决样本在测定过程中随着时间的推移所产生的信号漂变，从而为我们带来更加清新的数据呈现。该软件包的工作模式是将批次间和批次内的变异分开来单独处理，其工作路线图如下：
1 | 该软件包主要包含三个模块：
- 批次间一致性信号提取和比对；
- 批次内离子信号强度偏移的校正；
- 批次间数据的归一化。
2 | 该教程拟解决的常见问题
- 只有一个批次数据信号的校正：这个是比较容易实现的，适用于只有一个批次的数据，但是又想通过去除原始数据中仪器响应偏离的那部分而提高数据质量的情况。
- 多个批次间和批次内数据的校正：这个相对来说是个比较难啃的骨头，因为这涉及数据一致性，批次间或者批次内数据的变异识别，信号漂变校正和归一化处理，步骤比较繁琐。
针对以上这两种情况，该教程将讨论有关特定任务的问题以及使用程序包随附的数据显示工作代码。因此，您应该能够按照本教程生成数据输出和数据，以检查算法及其使用是否正确。
使用流程
包的安装和加载
install.packages("devtools") # 先安装devtools包，后面马上要用到
install_git("https://gitlab.com/CarlBrunius/batchCorr.git")
示例数据
- batchCorr算法需要表格格式的峰数据（peakTable），行表型进样个体，列表示变量。另外还需要一个包含样品信息，批次信息，进样分类，进样顺序的metadata表格。值得注意是，metadata表和peadTable表中的进样顺序要一致。如果不确定的话，可以通过以下指令进行确证。
identical(rownames(peakTable),metadata$sampleName)
- 另外一个在使用前需要注意的地方就是：对于信号偏移的校正，需要
peakTable不能含有缺失值；对于批次间的比对和数据一致性分析，需要用到没有进行缺失值处理的peakTable表。那么在今天介绍的数据流程分析中，作者采用xcms包进行峰提取，函数是fillPeaks()，接下来就是通过内部的多变量缺失值填充的方法。
单个批次信号校正
1 | 首先加载分析包以及示例数据
library(batchCorr)
data('OneBatchData')
2 | 数据组成
- B_PT表含peak信息的表（没有确实数据，如经过fillPeaks处理和imputation填充）
- B_meta表包含批次信息（所有来自批次B的样本）、样品分组（QC或者Reference）以及进样顺序（进样序列的编号）
3 | 单个批次信号漂变的校正
3.1 | 第一条命令
batchBCorr <- correctDrift(peakTable = B_PT, injections = B_meta$inj,
sampleGroups = B_meta$grp, QCID = 'QC',
modelNames = 'VVE', G = 17:22)
运行这条命令需要大概几分钟的时间，取决于你的PC电脑的配置。因为batchCorr采用的是mclust包进行分类，分析的算法采用的是
VVE，聚类数目从17开始试到22。
- peakTable：样本信息在行，特征值信息在列，特征值id需要唯一；
- injection：样品或者QC的进样顺序；
- sampleGroup：后接一个向量，如分组信息；
- QCID：将QC在sampleGroup的标识名称写出，如这里就是
QC； - modelNames和G参数：可以参考mclust包的解释
3.1.1 | 结果解释
该命令运行结束后，会得到几个PDF的图，先看下cluster_BIC图：
- 从该图得到的信息有，BIC（BIC值越大则说明所选取的变量集合拟合效果越好）值取最高时对应的cluster数目，这里是22。
- 当然这几个参数都是可以根据现实需要进行调节的，如modelNames = c('VVV','VVE','VEV','VEE','VEI','VVI','VII') 和G = seq(5,35,by=10)，只不过运行时间也会延长。
- 最终变量集合拟合的结果依赖于原始数据的质量（因为仪器设备的质量和信号漂变也是随机的），如果是比较正常且质量较高的数据，聚类簇的数目一般在6~25。如果是质量比较差的数据，那么G的设置可能会>50。但是实际上，作者并不认为设置如此大的聚类数是一个合理的选择，因此，40就是上限了。上述的这些一般性程序可以帮助你在有限的时间内找到合适你数据的参数。
- 上述构建的信号漂变的校正用到了QC样本，通过在整个进样序列中QC信号强度的漂变来构建模型。作者认为batchCorr之所以和其它信号校正工具相比显得独树一帜，是因为运用了上述提及的那些聚类方法来将特征值聚成一个个信号漂变的结合从而避免将仪器噪声也聚类进去的可能。
- batchCorr校正信号漂变的方法有QC和外来的参考样本。当然，QC样本是目前大多数实验室采用的一种方法。主要原因有两个：1）监控仪器性能和稳定性；2）对信号校正是否提高了数据质量的五篇估计——如果质量并没有提高，那么就不会对该cluster进行信号校正，因为这可能会带来其他的数据偏差成分。
3.2 | 第二条命令
当然，如果在数据采集的批次中间隔插入了参考样本信息，那么
correctDrift()函数也可以进行适当的调整来适应。可以通过增加参数'refID'来指定那些样品是参考样品（在这个数据集中就是"Ref"）:
batchBCorr <- correctDrift(peakTable = B_PT, injections = B_meta$inj,
sampleGroups = B_meta$grp, QCID = 'QC',
RefID='Ref', modelNames = 'VVE', G = 17:22)
testFeatsFinal <- batchBCorr$TestFeatsFinal # 提取信号校正后，在QC样本中CV 值小于0.3的数据，后续的统计分析用到的就是这张校正好且变异系数在合理范围内的数据
testFeatsFinal <- as.data.frame(testFeatsFinal)
3.2.1 | 结果解释
- 这条命令运行完后得到的结果存储在
batchCorr对象中，那么相关的信息查看可以通过$字符来提取。如batchCorr$actionInfo可以查看每个集合发生了什么变化；batchCorr$testFeatsCorr可以查看提取经过信号校正的数据；batchCorr$testFeatsFinal可以查看提取通过评判标准QC CV（变异系数）值<所设置的阈值（默认是0.3）。
多批次数据信号校正和归一化
4 | 多批次间数据信号漂变的校正
4.1 | 加载含多个批次的数据集合
library(batchCorr)
data('ThreeBatchData') # 加载数据
加载的这个数据集包含三个对象：
- PTnofill是含有缺失数据的质谱峰表，比如没有进行fillPeak和imputation；
- PTfill则是不含缺失数据，也就是经过fillPeak和imputation的质谱峰强度信息表；
- meta如上所述，是包含批次信息（如样品是来自三个批次，B、F和H）、样品分组和进样顺序信息。
4.2 | 首先来执行批次间的比对
## Perform batch alignment
# Extract peakinfo (i.e. m/z and rt of features)
peakIn <- peakInfo(PT = PTnofill, sep = '@', start = 3) # These column names have 2 leading characters describing LC-MS mode -> start at 3
dim(PTnofill)
# [1] 90 11815 没有比对之前是共有11815个features值
# Perform multi-batch alignment
alignBat <- alignBatches(peakInfo = peakIn, PeakTabNoFill = PTnofill, PeakTabFilled = PTfill, batches = meta$batch, sampleGroups = meta$grp, selectGroup = 'QC')
# Extract new peak table
PT=alignBat$PTalign
dim(PT)
# [1] 90 11284 比对后的feature个数减少到11284个。
4.2.1 | 结果解释
- 这一步就是以PTfill表为参考，因为PTfill是经过人为进行缺失值处理和填充处理的表。而PTnofill是原始的未经过处理的表（含有很多NA值）。这里的比对主要是将每个批次的缺失特征值集中起来并在不同批次之间进行比较：也就说如果两个特征值的mz值和rt都很接近，并且在不同批次间的缺失/存在特性一致，那么就将这两个特征值合并。
4.3 | 分批进行批次内信号校正
# Batch B
batchB <- getBatch(peakTable = PT, meta = meta, batch = meta$batch, select = 'B')
BCorr <- correctDrift(peakTable = batchB$peakTable, injections = batchB$meta$inj,
sampleGroups = batchB$meta$grp, QCID = 'QC',
G = seq(5,35,by=3), modelNames = c('VVE', 'VEE'))
# Batch F
batchF <- getBatch(peakTable = PT, meta = meta, batch = meta$batch, select = 'F')
FCorr <- correctDrift(peakTable = batchF$peakTable, injections = batchF$meta$inj,
sampleGroups = batchF$meta$grp, QCID = 'QC',
G = seq(5,35,by=3), modelNames = c('VVE', 'VEE'))
# Batch H
batchH <- getBatch(peakTable = PT, meta = meta, batch = meta$batch, select = 'H')
HCorr <- correctDrift(peakTable = batchH$peakTable, injections = batchH$meta$inj,
sampleGroups = batchH$meta$grp, QCID = 'QC',
G = seq(5,35,by=3),modelNames = c('VVE', 'VEE'))
4.3.1 | 结果解释
- 关于这几步的结果其实和前面讲的单个批次信号校正的参数和结果是一致的，可参考前面的说明
4.4 | 最后一步：不同批次校正结果合并及归一化
mergedData <- mergeBatches(list(BCorr,FCorr,HCorr))
normData <- normalizeBatches(peakTable = mergedData$peakTable,
batches = meta$batch, sampleGroup = meta$grp,
refGroup = 'Ref', population = 'sample') # refGroup也可以设置为"QC"
PTnorm <- normData$peakTable # 该数据结果可以用于后续的统计分析
4.1 | 结果解释
- mergeBatches函数可以将QC 中CV值<30%且在50%样本中存在的features合并，最终得到包含9,42个特征值的代谢物特征值丰度表。当然这个质控的比例可以通过参数qualRatio进行调整；
- normalizeBatches函数将会检查参考样本是否如何一定的质控标准，如果符合，那么就对其进行归一化；如果不符合标准，那么就一整个群体的中位数进行归一化。
结语
基本上的流程上面已经介绍完了，由于本人知识储备不是很足，所以在个人理解上会有瑕疵，请大家酌情而论。如果我有任何错误，请您及时发送邮件jnzd_hemaozhang@hotmail.com告知我，不胜感激！最后祝大家科研顺利！
参考
[1] Large-scale untargeted LC-MS metabolomics data correction using between-batch feature alignment and cluster-based within-batch signal intensity drift correction
[2] A Brief Tutorial on batchCorr: An R package for between- and within-batch drift correction of high-resolution mass spectrometry-based data.链接: https://pan.baidu.com/s/1rW2LTuBKolMKFfox34L3PQ
提取码：lrxi
