@我的博客：有味
導(dǎo)讀
液相色譜 - 質(zhì)譜（LC-MS）因?yàn)榭梢詸z測的代謝物的覆蓋范圍廣倡缠，靈敏度高和樣品制備簡單快捷而廣泛運(yùn)用于非靶向代謝組學(xué)研究中昙沦。如
血清载荔、尿液和腦脊液等樣本懒熙。然而工扎，在大樣本多批次數(shù)據(jù)采集的背景下肢娘，儀器本身隨著時(shí)間的遷移橱健，系統(tǒng)污染以及色譜柱活化再生過程的不當(dāng)?shù)纫幌盗性驅(qū)е聝x器對數(shù)據(jù)反映不穩(wěn)定畴博。這些不穩(wěn)定性會(huì)直接影響測定的精確分子量蓝仲、保留時(shí)間以及信號(hào)響應(yīng)強(qiáng)度袱结。尤其是在大樣本垢夹，多批次樣本采集背景下，數(shù)據(jù)產(chǎn)生的信號(hào)漂變會(huì)比較嚴(yán)重犀盟，一般來說阅畴，批次間的系統(tǒng)信號(hào)偏差會(huì)大于批次內(nèi)樣本間的信號(hào)偏差贱枣。這些測量誤差降低了數(shù)據(jù)的重復(fù)性和再現(xiàn)性纽哥，因此可能降低檢測生物反應(yīng)和模糊解釋的能力栖秕。
batchCorr包介紹
那么簇捍，R軟件包
batchCorr就是在這樣的一個(gè)背景下孕育而生的垦写。其主要目的就是為了解決樣本在測定過程中隨著時(shí)間的推移所產(chǎn)生的信號(hào)漂變梯投，從而為我們帶來更加清新的數(shù)據(jù)呈現(xiàn)分蓖。該軟件包的工作模式是將批次間和批次內(nèi)的變異分開來單獨(dú)處理么鹤，其工作路線圖如下：
1 | 該軟件包主要包含三個(gè)模塊：
- 批次間一致性信號(hào)提取和比對蒸甜；
- 批次內(nèi)離子信號(hào)強(qiáng)度偏移的校正柠新；
- 批次間數(shù)據(jù)的歸一化恨憎。
2 | 該教程擬解決的常見問題
- 只有一個(gè)批次數(shù)據(jù)信號(hào)的校正：這個(gè)是比較容易實(shí)現(xiàn)的净蚤，適用于只有一個(gè)批次的數(shù)據(jù)输硝，但是又想通過去除原始數(shù)據(jù)中儀器響應(yīng)偏離的那部分而提高數(shù)據(jù)質(zhì)量的情況塞栅。
- 多個(gè)批次間和批次內(nèi)數(shù)據(jù)的校正：這個(gè)相對來說是個(gè)比較難啃的骨頭穆端，因?yàn)檫@涉及數(shù)據(jù)一致性赖草，批次間或者批次內(nèi)數(shù)據(jù)的變異識(shí)別盗温，信號(hào)漂變校正和歸一化處理袄琳，步驟比較繁瑣亡脸。
針對以上這兩種情況怔软，該教程將討論有關(guān)特定任務(wù)的問題以及使用程序包隨附的數(shù)據(jù)顯示工作代碼衣厘。因此，您應(yīng)該能夠按照本教程生成數(shù)據(jù)輸出和數(shù)據(jù)影暴，以檢查算法及其使用是否正確错邦。
使用流程
包的安裝和加載
install.packages("devtools") # 先安裝devtools包，后面馬上要用到
install_git("https://gitlab.com/CarlBrunius/batchCorr.git")
示例數(shù)據(jù)
- batchCorr算法需要表格格式的峰數(shù)據(jù)（peakTable）型宙，行表型進(jìn)樣個(gè)體撬呢，列表示變量。另外還需要一個(gè)包含樣品信息妆兑，批次信息魂拦，進(jìn)樣分類，進(jìn)樣順序的metadata表格搁嗓。值得注意是芯勘，metadata表和peadTable表中的進(jìn)樣順序要一致。如果不確定的話腺逛，可以通過以下指令進(jìn)行確證荷愕。
identical(rownames(peakTable),metadata$sampleName)
- 另外一個(gè)在使用前需要注意的地方就是：對于信號(hào)偏移的校正，需要
peakTable不能含有缺失值棍矛；對于批次間的比對和數(shù)據(jù)一致性分析安疗，需要用到?jīng)]有進(jìn)行缺失值處理的peakTable表。那么在今天介紹的數(shù)據(jù)流程分析中茄靠，作者采用xcms包進(jìn)行峰提取茂契，函數(shù)是fillPeaks()，接下來就是通過內(nèi)部的多變量缺失值填充的方法慨绳。
單個(gè)批次信號(hào)校正
1 | 首先加載分析包以及示例數(shù)據(jù)
library(batchCorr)
data('OneBatchData')
2 | 數(shù)據(jù)組成
- B_PT表含peak信息的表（沒有確實(shí)數(shù)據(jù)，如經(jīng)過fillPeaks處理和imputation填充）
- B_meta表包含批次信息（所有來自批次B的樣本）、樣品分組（QC或者Reference）以及進(jìn)樣順序（進(jìn)樣序列的編號(hào)）
3 | 單個(gè)批次信號(hào)漂變的校正
3.1 | 第一條命令
batchBCorr <- correctDrift(peakTable = B_PT, injections = B_meta$inj,
sampleGroups = B_meta$grp, QCID = 'QC',
modelNames = 'VVE', G = 17:22)
運(yùn)行這條命令需要大概幾分鐘的時(shí)間脐雪，取決于你的PC電腦的配置厌小。因?yàn)閎atchCorr采用的是mclust包進(jìn)行分類，分析的算法采用的是
VVE战秋，聚類數(shù)目從17開始試到22璧亚。
- peakTable：樣本信息在行，特征值信息在列脂信，特征值id需要唯一癣蟋；
- injection：樣品或者QC的進(jìn)樣順序；
- sampleGroup：后接一個(gè)向量狰闪，如分組信息疯搅；
- QCID：將QC在sampleGroup的標(biāo)識(shí)名稱寫出，如這里就是
QC埋泵； - modelNames和G參數(shù)：可以參考mclust包的解釋
3.1.1 | 結(jié)果解釋
該命令運(yùn)行結(jié)束后幔欧，會(huì)得到幾個(gè)PDF的圖，先看下cluster_BIC圖：
- 從該圖得到的信息有丽声，BIC（BIC值越大則說明所選取的變量集合擬合效果越好）值取最高時(shí)對應(yīng)的cluster數(shù)目礁蔗，這里是22。
- 當(dāng)然這幾個(gè)參數(shù)都是可以根據(jù)現(xiàn)實(shí)需要進(jìn)行調(diào)節(jié)的雁社，如modelNames = c('VVV','VVE','VEV','VEE','VEI','VVI','VII') 和G = seq(5,35,by=10)浴井，只不過運(yùn)行時(shí)間也會(huì)延長。
- 最終變量集合擬合的結(jié)果依賴于原始數(shù)據(jù)的質(zhì)量（因?yàn)閮x器設(shè)備的質(zhì)量和信號(hào)漂變也是隨機(jī)的）霉撵，如果是比較正常且質(zhì)量較高的數(shù)據(jù)滋饲，聚類簇的數(shù)目一般在6~25。如果是質(zhì)量比較差的數(shù)據(jù)喊巍，那么G的設(shè)置可能會(huì)>50屠缭。但是實(shí)際上，作者并不認(rèn)為設(shè)置如此大的聚類數(shù)是一個(gè)合理的選擇崭参，因此呵曹，40就是上限了。上述的這些一般性程序可以幫助你在有限的時(shí)間內(nèi)找到合適你數(shù)據(jù)的參數(shù)何暮。
- 上述構(gòu)建的信號(hào)漂變的校正用到了QC樣本奄喂，通過在整個(gè)進(jìn)樣序列中QC信號(hào)強(qiáng)度的漂變來構(gòu)建模型。作者認(rèn)為batchCorr之所以和其它信號(hào)校正工具相比顯得獨(dú)樹一幟海洼，是因?yàn)檫\(yùn)用了上述提及的那些聚類方法來將特征值聚成一個(gè)個(gè)信號(hào)漂變的結(jié)合從而避免將儀器噪聲也聚類進(jìn)去的可能跨新。
- batchCorr校正信號(hào)漂變的方法有QC和外來的參考樣本。當(dāng)然坏逢，QC樣本是目前大多數(shù)實(shí)驗(yàn)室采用的一種方法域帐。主要原因有兩個(gè)：1）監(jiān)控儀器性能和穩(wěn)定性赘被；2）對信號(hào)校正是否提高了數(shù)據(jù)質(zhì)量的五篇估計(jì)——如果質(zhì)量并沒有提高，那么就不會(huì)對該cluster進(jìn)行信號(hào)校正肖揣，因?yàn)檫@可能會(huì)帶來其他的數(shù)據(jù)偏差成分民假。
3.2 | 第二條命令
當(dāng)然，如果在數(shù)據(jù)采集的批次中間隔插入了參考樣本信息龙优，那么
correctDrift()函數(shù)也可以進(jìn)行適當(dāng)?shù)恼{(diào)整來適應(yīng)羊异。可以通過增加參數(shù)'refID'來指定那些樣品是參考樣品（在這個(gè)數(shù)據(jù)集中就是"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 # 提取信號(hào)校正后彤断，在QC樣本中CV 值小于0.3的數(shù)據(jù)野舶，后續(xù)的統(tǒng)計(jì)分析用到的就是這張校正好且變異系數(shù)在合理范圍內(nèi)的數(shù)據(jù)
testFeatsFinal <- as.data.frame(testFeatsFinal)
3.2.1 | 結(jié)果解釋
- 這條命令運(yùn)行完后得到的結(jié)果存儲(chǔ)在
batchCorr對象中，那么相關(guān)的信息查看可以通過$字符來提取宰衙。如batchCorr$actionInfo可以查看每個(gè)集合發(fā)生了什么變化平道；batchCorr$testFeatsCorr可以查看提取經(jīng)過信號(hào)校正的數(shù)據(jù)；batchCorr$testFeatsFinal可以查看提取通過評(píng)判標(biāo)準(zhǔn)QC CV（變異系數(shù)）值<所設(shè)置的閾值（默認(rèn)是0.3）菩浙。
多批次數(shù)據(jù)信號(hào)校正和歸一化
4 | 多批次間數(shù)據(jù)信號(hào)漂變的校正
4.1 | 加載含多個(gè)批次的數(shù)據(jù)集合
library(batchCorr)
data('ThreeBatchData') # 加載數(shù)據(jù)
加載的這個(gè)數(shù)據(jù)集包含三個(gè)對象：
- PTnofill是含有缺失數(shù)據(jù)的質(zhì)譜峰表巢掺，比如沒有進(jìn)行fillPeak和imputation；
- PTfill則是不含缺失數(shù)據(jù)劲蜻，也就是經(jīng)過fillPeak和imputation的質(zhì)譜峰強(qiáng)度信息表陆淀；
- meta如上所述，是包含批次信息（如樣品是來自三個(gè)批次先嬉，B轧苫、F和H）、樣品分組和進(jìn)樣順序信息疫蔓。
4.2 | 首先來執(zhí)行批次間的比對
## 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個(gè)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個(gè)數(shù)減少到11284個(gè)含懊。
4.2.1 | 結(jié)果解釋
- 這一步就是以PTfill表為參考，因?yàn)镻Tfill是經(jīng)過人為進(jìn)行缺失值處理和填充處理的表衅胀。而PTnofill是原始的未經(jīng)過處理的表（含有很多NA值）岔乔。這里的比對主要是將每個(gè)批次的缺失特征值集中起來并在不同批次之間進(jìn)行比較：也就說如果兩個(gè)特征值的mz值和rt都很接近，并且在不同批次間的缺失/存在特性一致滚躯，那么就將這兩個(gè)特征值合并雏门。
4.3 | 分批進(jìn)行批次內(nèi)信號(hào)校正
# 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 | 結(jié)果解釋
- 關(guān)于這幾步的結(jié)果其實(shí)和前面講的單個(gè)批次信號(hào)校正的參數(shù)和結(jié)果是一致的，可參考前面的說明
4.4 | 最后一步：不同批次校正結(jié)果合并及歸一化
mergedData <- mergeBatches(list(BCorr,FCorr,HCorr))
normData <- normalizeBatches(peakTable = mergedData$peakTable,
batches = meta$batch, sampleGroup = meta$grp,
refGroup = 'Ref', population = 'sample') # refGroup也可以設(shè)置為"QC"
PTnorm <- normData$peakTable # 該數(shù)據(jù)結(jié)果可以用于后續(xù)的統(tǒng)計(jì)分析
4.1 | 結(jié)果解釋
- mergeBatches函數(shù)可以將QC 中CV值<30%且在50%樣本中存在的features合并掸掏，最終得到包含9,42個(gè)特征值的代謝物特征值豐度表茁影。當(dāng)然這個(gè)質(zhì)控的比例可以通過參數(shù)qualRatio進(jìn)行調(diào)整；
- normalizeBatches函數(shù)將會(huì)檢查參考樣本是否如何一定的質(zhì)控標(biāo)準(zhǔn)丧凤，如果符合募闲，那么就對其進(jìn)行歸一化；如果不符合標(biāo)準(zhǔn)愿待，那么就一整個(gè)群體的中位數(shù)進(jìn)行歸一化浩螺。
結(jié)語
基本上的流程上面已經(jīng)介紹完了靴患，由于本人知識(shí)儲(chǔ)備不是很足，所以在個(gè)人理解上會(huì)有瑕疵年扩，請大家酌情而論蚁廓。如果我有任何錯(cuò)誤访圃，請您及時(shí)發(fā)送郵件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
