基于氣相色譜-質譜的尿液代謝組學技術結合化學計量學用于戊二酸血癥Ⅰ型早期檢測研究

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

（1．長江大學 物理與光電工程學院，湖北 荊州 434023；2．深圳愛灣醫學檢驗實驗室 深圳罕見病代謝組學精準醫學工程研究中心，廣東 深圳 518000；3．中國科學院大學 深圳醫院，廣東 深圳 518000）

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

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

1 實驗部分

1.1 樣本材料

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

1.2 儀器與試劑

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

1.3 GC-MS代謝組學圖譜數據采集

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

2 GA-Ⅰ早期檢測模型

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

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

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

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

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

3 結果與討論

3.1 GC-MS代謝組學圖譜

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

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

3.2 特征變量穩健篩選

基于LW、VIP和sMC信息向量，對每個訓練集利用單個PLS-DA建立的模型分別為LW-PLSDA、VIP-PLSDA和sMC-PLSDA，結合自助抽樣法生成多個擾動數據集建立的模型分別為BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA。當樣本按照7∶3比例（訓練集43例，測試集19例）劃分時，對50個訓練集分別建模后，公式（1）中的h=18，計算得到LW-PLSDA、VIP-PLSDA和sMCPLSDA的KI值均低于0．4（圖3）。可見，數據集的高維小樣本特點大大影響了基于單個PLS-DA建模的性能，導致LW-PLSDA、VIP-PLSDA和sMCPLSDA篩選的特征變量序列在各訓練集之間差異過大，即使預測性能優異也大大降低了各領域專家對分析結果的可信度。引入自助抽樣法，BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA的KI值均顯著增加，特別是BS-VIP-PLSDA在篩選特征變量個數為12時的KI值高達0．807 5，表明基于樣本擾動的集成特征選擇策略能夠專注挑選被多個擾動數據集持續篩選的變量，有效提高特征選擇的穩健性。進一步按照6∶4比例劃分（訓練集37例，測試集25例）增大50個訓練集間的樣本差異，以避免不同訓練集始終篩選相同的變量。此時，公式（1）中的h=16，BS-PLSDA各模型仍然比其單獨建模的PLS-DA展示出更優異的特征變量篩選穩健性（圖4）。

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

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

3.3 穩健篩選特征變量的建模性能

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

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

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

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

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

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

4 結論

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