肖 雯，牛芊芊，孫智勇，楊 琴*，吳本清

（1．長(zhǎng)江大學(xué) 物理與光電工程學(xué)院，湖北 荊州 434023；2．深圳愛(ài)灣醫(yī)學(xué)檢驗(yàn)實(shí)驗(yàn)室 深圳罕見(jiàn)病代謝組學(xué)精準(zhǔn)醫(yī)學(xué)工程研究中心，廣東 深圳 518000；3．中國(guó)科學(xué)院大學(xué) 深圳醫(yī)院，廣東 深圳 518000）

戊二酸血癥Ⅰ型（GA-Ⅰ）作為一種常染色體隱性遺傳病，因編碼基因發(fā)生突變導(dǎo)致對(duì)應(yīng)的戊二酰輔酶A脫氫酶（GCDH）活性降低或缺失，有毒代謝產(chǎn)物（如戊二酰肉堿、戊二酸、3-羥基戊二酸、戊烯二酸等）在體液與組織中異常蓄積，形成酸中毒而造成嚴(yán)重的不可逆神經(jīng)系統(tǒng)受損［1-2］?；颊吲R床表現(xiàn)復(fù)雜，易與其他神經(jīng)系統(tǒng)疾病混淆，缺乏特異性［3］。因此，建立高效的GA-Ⅰ早期篩查模型，對(duì)緩解新生兒遺傳代謝病公共健康難題具有積極的促進(jìn)作用。隨著新生兒篩查項(xiàng)目的推廣和普及，基于氣相色譜-質(zhì)譜（GC-MS）的尿液代謝組學(xué)技術(shù)在GA-Ⅰ早期檢測(cè)中展示出獨(dú)特優(yōu)勢(shì)，能全面揭示代謝產(chǎn)物的濃度差異［4-6］?；瘜W(xué)計(jì)量學(xué)中各種建模算法的挖掘能力可將數(shù)據(jù)內(nèi)在信息有效轉(zhuǎn)化為可理解的知識(shí)，已廣泛應(yīng)用于疾病的早期檢測(cè)和輔助診斷［7-9］。然而，質(zhì)譜技術(shù)的不斷發(fā)展產(chǎn)生高維小樣本的疾病代謝組學(xué)圖譜，使得應(yīng)用單個(gè)數(shù)據(jù)分析學(xué)習(xí)模型進(jìn)行特征篩選不夠穩(wěn)健，對(duì)樣本變異（增加或刪除一些樣本）敏感，不僅影響疾病早期檢測(cè)模型的建模性能，也會(huì)降低各領(lǐng)域?qū)＜覍?duì)分析結(jié)果的可信度［10-11］。

自助抽樣法（Bootstrap）是集成特征篩選中一種常用的重抽樣技術(shù)，先通過(guò)有放回的抽樣方式產(chǎn)生大量樣本變異較小且與原始數(shù)據(jù)結(jié)構(gòu)相同的擾動(dòng)數(shù)據(jù)集，再使用相同的數(shù)據(jù)分析學(xué)習(xí)算法建立多個(gè)特征選擇器［12-13］。該法能夠在多個(gè)擾動(dòng)數(shù)據(jù)集中挑選出持續(xù)被篩選的變量，其穩(wěn)健性和解釋樣本組別間差異的優(yōu)勢(shì)已得到實(shí)驗(yàn)和理論上的證實(shí)［14］。偏最小二乘判別分析（PLS-DA）是處理高維數(shù)據(jù)中共線性和信息冗余問(wèn)題的有力統(tǒng)計(jì)分析工具，其數(shù)據(jù)解釋能力強(qiáng)，通過(guò)計(jì)算系列信息向量如載荷（LW）、變量投影重要性（VIP）、顯著性多元相關(guān)（sMC）等直接刻畫各變量對(duì)模型的貢獻(xiàn)大?。?5-16］。本文應(yīng)用GC-MS技術(shù)采集新生兒尿液的代謝組學(xué)圖譜信息，結(jié)合自助抽樣法和PLS-DA實(shí)現(xiàn)真正解釋組別間差異且與疾病機(jī)理密切相關(guān)的穩(wěn)健特征變量篩選，建立了高性能的GA-Ⅰ早期檢測(cè)模型，并進(jìn)行分析和驗(yàn)證。

1 實(shí)驗(yàn)部分

1.1 樣本材料

尿液樣本采集于2015~2020年，由中國(guó)科學(xué)院大學(xué)深圳醫(yī)院提供。疾病組由30例GA-Ⅰ患者（＜28天）組成，對(duì)照組由32例年齡匹配的健康新生兒組成。本研究經(jīng)中國(guó)科學(xué)院大學(xué)深圳醫(yī)院倫理委員會(huì)同意，尿液代謝組學(xué)圖譜檢測(cè)均獲得參與者監(jiān)護(hù)人的同意。尿液樣本保存于-20℃，不添加任何防腐劑。

1.2 儀器與試劑

GCMS-QP2010型氣相色譜-質(zhì)譜聯(lián)用儀（日本島津公司）。乙酸乙酯（色譜純）、甲基硅烷化衍生試劑（99∶1，上海安譜實(shí)驗(yàn)科技股份有限公司），十七烷酸、尿素酶、鹽酸羥胺（美國(guó)Sigma-Aldrich公司），鹽酸（優(yōu)級(jí)純，湖南凱信公司），苦味酸（分析純，成都西亞化學(xué)工業(yè)有限公司），氫氧化鈉（分析純，美國(guó)Sigma-Aldrich公司）。以十七烷酸作為內(nèi)標(biāo)，用乙酸乙酯配制質(zhì)量濃度為0．5 mg/mL的十七烷酸溶液。

1.3 GC-MS代謝組學(xué)圖譜數(shù)據(jù)采集

尿液樣本用蒸餾水稀釋定容至2 mL，其中肌酐質(zhì)量為0．2 mg。樣本中加入20 μL尿素酶分解尿素，并與0．02 mg十七烷酸混合，然后經(jīng)鹽酸羥胺、氫氧化鈉和鹽酸處理，再加入3 mL乙酸乙酯對(duì)有機(jī)酸萃取2次。將離心（4 000 r/min，5 min）后的有機(jī)酸層移至干凈離心管中，用氮?dú)庠?0℃下吹干。加入100 μL甲基硅烷化衍生試劑（99∶1），在70℃衍生反應(yīng)30 min后，取1 μL樣本采用GC-MS測(cè)定尿液中的有機(jī)酸。測(cè)試條件為：進(jìn)樣口溫度280℃，色譜柱溫度100~280℃，以4℃/min逐步升溫，質(zhì)譜接口溫度280℃，離子源溫度200℃。電子束能量70 eV，掃描范圍m/z50~500，掃描速率1 000 Da/s。采集的原始質(zhì)譜數(shù)據(jù)利用配套的GCMSsolution軟件和代謝物質(zhì)譜數(shù)據(jù)庫(kù)進(jìn)行處理，鑒定出132種代謝物。對(duì)于同一個(gè)色譜峰，鑒別出的代謝物濃度通過(guò)其峰面積與內(nèi)標(biāo)物十七烷酸峰面積的比值確定。

2 GA-Ⅰ早期檢測(cè)模型

GA-Ⅰ早期檢測(cè)模型的技術(shù)原理如圖1所示。首先，將62例樣本隨機(jī)劃分為訓(xùn)練集和測(cè)試集。然后，對(duì)劃分的訓(xùn)練集采用自助抽樣法進(jìn)行有放回的抽樣，生成多個(gè)樣本大小與原始訓(xùn)練集相同的擾動(dòng)數(shù)據(jù)集?；赑LS-DA分別建模，計(jì)算LW、VIP和sMC的信息向量，并根據(jù)指標(biāo)的絕對(duì)值大小對(duì)變量進(jìn)行排序?？紤]到各變量對(duì)模型的貢獻(xiàn)，只取排序后的變量序列中前10%的變量，作為真正表征組別間代謝差異的特征變量。最后，集成所有擾動(dòng)數(shù)據(jù)集對(duì)應(yīng)的排序變量序列，統(tǒng)計(jì)各變量跨越模型間的篩選頻率，并根據(jù)篩選頻率再次對(duì)變量進(jìn)行排序，設(shè)置頻率篩選閾值（如0．1，0．15，0．2，…，1）確定最終的特征變量序列［13］。頻率篩選閾值設(shè)置得越高，最終特征變量序列的確定越缺乏偶然性，篩選越嚴(yán)格。對(duì)于測(cè)試集，根據(jù)確定的最終特征變量序列挑選出相應(yīng)的變量組成數(shù)據(jù)集，利用PLS-DA建模，計(jì)算y?對(duì)樣本類別進(jìn)行預(yù)測(cè)。

圖1 GA-Ⅰ早期檢測(cè)模型流程圖Fig．1 The flowchart of early detection model for GA-Ⅰ

62例樣本分別按照7∶3和6∶4比例劃分訓(xùn)練集和測(cè)試集。自助抽樣法的抽樣次數(shù)設(shè)為100次，每個(gè)擾動(dòng)數(shù)據(jù)集對(duì)應(yīng)的PLS-DA模型的隱變量個(gè)數(shù)利用10折交叉驗(yàn)證法確定。另外，為驗(yàn)證GA-Ⅰ早期檢測(cè)模型的特征篩選穩(wěn)健性，對(duì)62例樣本隨機(jī)劃分50次，生成的50個(gè)訓(xùn)練集分別采用圖1的流程篩選最終的特征變量序列，計(jì)算Kuncheva指數(shù)（KI）［17］：

其中fi表示第i個(gè)訓(xùn)練集篩選的最終特征變量序列；h=|fi|=|fj|，表示所有篩選的最終特征變量序列中包含元素個(gè)數(shù)的最小值；r=|fi∩fj|，表示fi和fj中共同元素的個(gè)數(shù)；N表示尿液中鑒定出的代謝物總數(shù)；修正項(xiàng)h2/N表示在fi和fj中選擇一個(gè)共同特征其偶然性的概率。KI取值范圍為-1~1，負(fù)值表明特征重疊在很大程度上是由于偶然性。KI指數(shù)越大，fi和fj的特征重疊程度越高。最后，總穩(wěn)定性KItot取所有兩兩比較相似度值KI的平均：

為方便表述，各種方法的總穩(wěn)定性用KI簡(jiǎn)化表示。對(duì)于劃分50次生成的50組訓(xùn)練集和測(cè)試集，分別采用受試者工作特征曲線下面積（AUC）、正確率（ACC）、靈敏度（SEN）、特異性（SPE）以及馬修斯相關(guān)系數(shù)（MCC）表征模型性能。

3 結(jié)果與討論

3.1 GC-MS代謝組學(xué)圖譜

圖2展示了GA-Ⅰ和對(duì)照組樣本的尿液代謝總離子流色譜圖（TIC），在指定的保留時(shí)間范圍內(nèi)觀察到豐富的代謝產(chǎn)物圖譜，顯示出GC-MS能夠同時(shí)檢測(cè)尿液中多種有機(jī)酸的能力。對(duì)鑒定的132種代謝物進(jìn)行單變量統(tǒng)計(jì)分析，發(fā)現(xiàn)多種代謝物在GA-Ⅰ和對(duì)照組之間均出現(xiàn)顯著的濃度變化（p＜0．01），揭示了GA-Ⅰ的發(fā)病機(jī)制導(dǎo)致復(fù)雜的組別間差異。此外，代謝物之間還存在共線性問(wèn)題，如3-羥基-3-甲基戊二酸和3-羥基異戊酸的濃度相關(guān)系數(shù)為0．977 3，異枸櫞酸和烏頭酸的濃度相關(guān)系數(shù)為0．922 0。因此，采用多元統(tǒng)計(jì)分析方法PLS-DA結(jié)合自助抽樣法處理變量間的高相關(guān)性，以提取GCMS尿液代謝組學(xué)圖譜中更細(xì)微的代謝信息變化，選擇能夠很好解釋疾病組與對(duì)照組間差異的穩(wěn)健特征。

圖2 GA-Ⅰ和對(duì)照組樣本的尿液代謝總離子流色譜圖Fig．2 Total ion chromatograms（TIC）of urine metabolic profiling for GA-Ⅰand controls

3.2 特征變量穩(wěn)健篩選

基于LW、VIP和sMC信息向量，對(duì)每個(gè)訓(xùn)練集利用單個(gè)PLS-DA建立的模型分別為L(zhǎng)W-PLSDA、VIP-PLSDA和sMC-PLSDA，結(jié)合自助抽樣法生成多個(gè)擾動(dòng)數(shù)據(jù)集建立的模型分別為BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA。當(dāng)樣本按照7∶3比例（訓(xùn)練集43例，測(cè)試集19例）劃分時(shí)，對(duì)50個(gè)訓(xùn)練集分別建模后，公式（1）中的h=18，計(jì)算得到LW-PLSDA、VIP-PLSDA和sMCPLSDA的KI值均低于0．4（圖3）?？梢?jiàn)，數(shù)據(jù)集的高維小樣本特點(diǎn)大大影響了基于單個(gè)PLS-DA建模的性能，導(dǎo)致LW-PLSDA、VIP-PLSDA和sMCPLSDA篩選的特征變量序列在各訓(xùn)練集之間差異過(guò)大，即使預(yù)測(cè)性能優(yōu)異也大大降低了各領(lǐng)域?qū)＜覍?duì)分析結(jié)果的可信度。引入自助抽樣法，BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA的KI值均顯著增加，特別是BS-VIP-PLSDA在篩選特征變量個(gè)數(shù)為12時(shí)的KI值高達(dá)0．807 5，表明基于樣本擾動(dòng)的集成特征選擇策略能夠?qū)Ｗ⑻暨x被多個(gè)擾動(dòng)數(shù)據(jù)集持續(xù)篩選的變量，有效提高特征選擇的穩(wěn)健性。進(jìn)一步按照6∶4比例劃分（訓(xùn)練集37例，測(cè)試集25例）增大50個(gè)訓(xùn)練集間的樣本差異，以避免不同訓(xùn)練集始終篩選相同的變量。此時(shí)，公式（1）中的h=16，BS-PLSDA各模型仍然比其單獨(dú)建模的PLS-DA展示出更優(yōu)異的特征變量篩選穩(wěn)健性（圖4）。

圖3 各模型篩選的特征變量序列在50個(gè)訓(xùn)練集之間的篩選穩(wěn)健性比較Fig．3 Comparison of selection stability across 50 training sets among various techniques training set：43；test set：19

圖4 各模型篩選的特征變量序列在50個(gè)訓(xùn)練集之間的篩選穩(wěn)健性比較Fig．4 Comparison of selection stability across 50 training sets among various techniques training set：37；test set：25

3.3 穩(wěn)健篩選特征變量的建模性能

為嚴(yán)格篩選特征變量，結(jié)合自助抽樣法的3種模型BS-LW-PLSDA、BS-VIP-PLSDA和BSsMC-PLSDA對(duì)應(yīng)的頻率篩選閾值最高為0．3，即在100次重抽樣中至少被30個(gè)擾動(dòng)數(shù)據(jù)集建立的模型篩選。在頻率篩選閾值設(shè)置為0．3時(shí)，3個(gè)模型由50個(gè)訓(xùn)練集確定的最終特征變量序列分別輸入到對(duì)應(yīng)的測(cè)試集中。同時(shí)，由132種代謝物組成的全變量序列也分別輸入50個(gè)測(cè)試集中采用PLS-DA建模。當(dāng)樣本按照7∶3比例（訓(xùn)練集43例，測(cè)試集19例）劃分時(shí)，建模結(jié)果如表1所示。全變量PLS-DA模型對(duì)50個(gè)測(cè)試集的AUC平均值為0．837 9，ACC平均值為0．821 5，MCC平均值為0．724 2。BSPLSDA 3個(gè)模型中，BS-VIP-PLSDA和BS-sMC-PLSDA對(duì)50個(gè)測(cè)試集的AUC平均值分別為0．854 8和0．847 1，ACC平均值分別為0．835 3和0．850 5，MCC平均值分別為0．783 8和0．801 3。結(jié)果表明相比全變量PLS-DA模型，BS-VIP-PLSDA和BS-sMC-PLSDA采用更少的信息變量可獲得更好的建模性能，顯示了特征篩選的重要性。進(jìn)一步，分別統(tǒng)計(jì)3個(gè)模型確定的最終特征變量序列在50個(gè)訓(xùn)練集之間各變量的篩選頻率，再次排序后的特征變量序列如表1所示。已有文獻(xiàn)報(bào)道，在GA-Ⅰ確診患兒的尿液中戊二酸（Glutaric acid）、3-羥基戊二酸（3-Hydroxyglutaric acid）、戊烯二酸（Glutaconic acid）的濃度增高，具有一定的診斷意義［3，18］。由表1可見(jiàn)，sMC信息向量的數(shù)據(jù)解釋能力優(yōu)于LW和VIP信息向量，上述3種具有一定診斷意義的代謝物均被篩選出，且在對(duì)應(yīng)的熱圖分析（圖5）中，GA-Ⅰ和對(duì)照組之間具有顯著的濃度差異（p＜0．05）。另外，LW和VIP信息向量雖然均篩選出戊二酸和戊烯二酸，但對(duì)于LW戊烯二酸僅被50個(gè)訓(xùn)練集中的2個(gè)篩選，而對(duì)于VIP戊烯二酸被50個(gè)訓(xùn)練集中的13個(gè)篩選，體現(xiàn)了兩者不同的特征變量搜索能力。需要注意的是，BS-LW-PLSDA的模型性能指標(biāo)雖然低于全變量PLS-DA模型，但分類效果令人滿意，最重要的是能夠提供表征疾病組與對(duì)照組之間代謝信息差異的特征變量信息，顯著提升模型的解釋能力和臨床應(yīng)用價(jià)值。其他被篩選出的特征代謝物應(yīng)進(jìn)行實(shí)驗(yàn)驗(yàn)證，有助于挖掘更多與GA-Ⅰ代謝機(jī)理相關(guān)的知識(shí)。

圖5 BS-sMC-PLSDA篩選的特征變量組別間差異熱圖分析Fig．5 Heatmap of the top-ranked significant metabolites selected by BS-sMC-PLSDA

表1 結(jié)合自助抽樣法各種模型的篩選特征變量及建模結(jié)果Table 1 The top-ranked significant metabolites and classification parameters selected by various techniques combined bootstrap

當(dāng)樣本按照6∶4比例（訓(xùn)練集37例，測(cè)試集25例）劃分時(shí)，BS-PLSDA 3個(gè)模型對(duì)應(yīng)的頻率篩選閾值最高也為0．3，建模結(jié)果如表2所示。雖然訓(xùn)練集樣本減少，但經(jīng)特征篩選后，BS-VIP-PLSDA和BS-sMC-PLSDA對(duì)50個(gè)測(cè)試集的AUC平均值分別為0．857 9和0．807 6，ACC平均值分別為0．838 4和0．813 6，MCC平均值分別為0．790 9和0．740 6，模型性能指標(biāo)仍高于全變量PLS-DA模型。另外，相比于BS-VIP-PLSDA，BS-sMC-PLSDA的建模性能略低，但仍展示出最優(yōu)的數(shù)據(jù)解釋能力，能同時(shí)篩選出3種具有一定診斷意義的代謝物。對(duì)比兩表中BS-VIP-PLSDA和BS-sMC-PLSDA分別確定的特征變量序列，發(fā)現(xiàn)對(duì)于不同樣本劃分比例，除了排序稍有不同，其組成大致相同。對(duì)于BSLW-PLSDA，相比7∶3比例樣本劃分，盡管在6∶4比例劃分時(shí)篩選的特征變量個(gè)數(shù)降為2個(gè)，但代謝物分別為與GA-Ⅰ代謝機(jī)理密切相關(guān)的戊二酸和戊烯二酸。表明了基于自助抽樣法的集成特征選擇策略對(duì)提高變量篩選穩(wěn)健性和建模性能的有效性。

表2 結(jié)合自助抽樣法各種模型的篩選特征變量及建模結(jié)果Table 2 The top-ranked significant metabolites and classification parameters selected by various techniques combined bootstrap

為進(jìn)一步展示本文GA-Ⅰ早期檢測(cè)模型的可行性，在基于樣本擾動(dòng)的集成特征選擇策略基礎(chǔ)上，對(duì)每個(gè)擾動(dòng)數(shù)據(jù)集采用支持向量機(jī)遞歸特征消除法（SVM-RFE）［19］進(jìn)行建模。在特征變量篩選時(shí)，考慮保留其對(duì)模型貢獻(xiàn)的信息，SVM首選線性核函數(shù)，對(duì)應(yīng)模型分別為L(zhǎng)IN-SVMRFE和BS-LINSVMRFE。另外，考慮SVM對(duì)非線性數(shù)據(jù)的優(yōu)異分析能力，利用徑向基核函數(shù)分別建模RBFSVMRFE和BS-RBF-SVMRFE。2種核函數(shù)對(duì)應(yīng)的參數(shù)（C和gamma）利用10折交叉驗(yàn)證法確定。在遞歸循環(huán)中，SVM-RFE每次刪除20%的變量個(gè)數(shù)。與BS-PLSDA 3個(gè)模型采用相同的自助抽樣法配置，隨著篩選變量個(gè)數(shù)的增加，BS-SVMRFE兩個(gè)模型的KI值約為單獨(dú)建模SVM-RFE模型KI值的2倍（圖3和圖4）。同樣，頻率篩選閾值設(shè)置為0．3，確定的最終特征變量序列分別輸入到50組訓(xùn)練集和測(cè)試集中，BS-SVMRFE 2個(gè)模型的建模性能均優(yōu)于BS-PLSDA 3個(gè)模型，特別是樣本按照7∶3比例劃分時(shí)BS-RBF-SVMRFE模型的性能指標(biāo)平均值均超過(guò)0．900 0（表1和表2）。但在數(shù)據(jù)解釋能力方面，線性核函數(shù)保留了變量對(duì)模型貢獻(xiàn)的信息，BS-LIN-SVMRFE篩選出戊二酸和3-羥基戊二酸兩個(gè)具有診斷意義的代謝物；而徑向基核函數(shù)丟失了這部分信息，導(dǎo)致BS-RBF-SVMRFE只篩選出戊二酸代謝物，且其他特征變量與BS-LIN-SVMRFE和BS-PLSDA 3個(gè)模型也明顯不同。表明GA-Ⅰ早期檢測(cè)模型采用基于樣本擾動(dòng)的集成特征選擇策略是合理的，對(duì)每個(gè)擾動(dòng)數(shù)據(jù)集采用PLS-DA建?？赏瑫r(shí)兼顧建模性能和模型解釋能力，符合實(shí)際臨床需求；采用SVM-RFE建模雖然可獲得較好的建模性能，但解釋能力略低。

4 結(jié)論

為提高對(duì)高維小樣本疾病數(shù)據(jù)的建模能力，基于GC-MS聯(lián)用技術(shù)，結(jié)合自助抽樣法和PLS-DA建立了GA-Ⅰ早期檢測(cè)模型。通過(guò)對(duì)尿液的代謝組學(xué)圖譜進(jìn)行分析驗(yàn)證，相比單個(gè)模型建模，自助抽樣法通過(guò)重復(fù)抽樣的方式生成多個(gè)擾動(dòng)數(shù)據(jù)集，結(jié)合PLS-DA的特征選擇能力，建立的GA-Ⅰ早期檢測(cè)模型專注于持續(xù)被篩選的變量，可有效提升特征篩選的穩(wěn)健性。進(jìn)一步，基于LW、VIP和sMC在特征空間的搜索能力，篩選出的穩(wěn)健特征變量不僅分類效果令人滿意，而且真正解釋了組別間的差異與GA-Ⅰ的代謝機(jī)理密切相關(guān)，展示了豐富的生物學(xué)意義以及優(yōu)異的數(shù)據(jù)解釋能力。由此可見(jiàn)，本研究提出的模型在GA-Ⅰ的早期檢測(cè)、輔助診斷以及疾病機(jī)理研究中具有一定的潛力。