縱向組學數(shù)據(jù)統(tǒng)計分析方法和研究策略*

呂嘉麗 范冰冰 魏夢珂 張 濤△

近年來，隨著高通量檢測技術與機器學習方法的發(fā)展，以前瞻性隊列設計為基礎的縱向組學研究已經(jīng)成為系統(tǒng)生物學研究新趨勢[1]。統(tǒng)計學上，縱向數(shù)據(jù)統(tǒng)計分析方法已經(jīng)形成了系統(tǒng)的統(tǒng)計分析框架，但這些方法主要集中于變量數(shù)目小于觀察單位數(shù)的低維醫(yī)學縱向數(shù)據(jù)[2]。而針對縱向組學數(shù)據(jù)，目前仍缺乏成熟的統(tǒng)計分析策略。本文擬對近年來國內(nèi)外研究者提出的縱向組學數(shù)據(jù)統(tǒng)計分析方法進行介紹，并系統(tǒng)地總結(jié)各個方法的核心思想及優(yōu)缺點，給出縱向組學數(shù)據(jù)統(tǒng)計分析策略。

縱向組學研究設計與數(shù)據(jù)特點

在不同的生命及疾病狀態(tài)下，機體組學標記物濃度處于連續(xù)動態(tài)變化過程?？v向組學研究設計是指在疾病發(fā)生發(fā)展過程中或采取干預措施后采集多個時間點的生物標本，進行高通量組學檢測。該研究設計可以分析組學標記物的動態(tài)變化規(guī)律，探討生物體對危險因素累積的反應過程及疾病發(fā)生發(fā)展機制。

在數(shù)據(jù)特點上，縱向組學數(shù)據(jù)包含個體、時間、組學標記物、相關暴露因素及結(jié)局測量信息(如表1所示)?？v向組學數(shù)據(jù)具有以下特征：(1)非平衡性：不同觀察單位隨訪時間點不一樣，同一觀察單位間的隨訪間隔也不一樣；(2)自相關性：同一觀察單位的不同測量變量之間具有復雜的相關與因果調(diào)控關系；(3)測量誤差：組學研究常常需要借助高通量檢測儀器完成生物樣本測量，數(shù)據(jù)集難免存在由于測量誤差引起的變異；(4)時依混雜：在縱向組學研究中，組學測量數(shù)據(jù)及環(huán)境暴露因素均隨時間而動態(tài)變化，會對因果效應估計產(chǎn)生影響；(5)高維災難：縱向組學數(shù)據(jù)除具有一般縱向數(shù)據(jù)特點外，還存在一般高維小樣本組學數(shù)據(jù)的高維災難問題。

表1 縱向組學數(shù)據(jù)結(jié)構

縱向組學數(shù)據(jù)常用統(tǒng)計分析方法

縱向組學研究設計的研究目的包括：(1)研究不同組別間組學標記物的動態(tài)輪廓差異，發(fā)現(xiàn)不同組別間的動態(tài)差異組學標記物；(2)研究組學標記物隨時間變化的動態(tài)趨勢；(3)基于動態(tài)差異組學標記物，建立預測模型。根據(jù)上述研究目標及數(shù)據(jù)特點，常見縱向組學設計分析方法主要有以下三類[3-5]：(1)單變量統(tǒng)計分析：用于識別隨時間變化而改變的差異動態(tài)標記物；(2)聚類方法：對標記物的動態(tài)變化進行趨勢分析，對變化趨勢一致的標記物進行聚類；(3)降維方法：考慮到變量間的復雜相關性，利用多變量統(tǒng)計分析方法對高維小樣本數(shù)據(jù)進行降維，發(fā)現(xiàn)組學標記物在不同組別間的組學輪廓差異。

1.單變量分析

重復測量方差分析(repeated measures ANOVA) 是早期用于縱向數(shù)據(jù)分析的方法[6]。目前，混合效應模型(mixed effects models)及廣義估計方程(generalized estimating equations)是縱向組學數(shù)據(jù)分析的常用單變量統(tǒng)計方法[2]。曹紅艷等[7]在廣義估計方程基礎上提出了一種懲罰廣義估計方程(penalized generalized estimating equations)，并運用該方法對小鼠進行糖尿病發(fā)病關聯(lián)基因位點篩選。該方法的核心思想是基于LASSO或SCAD等懲罰方法進行廣義估計方程建模，不僅保持了廣義估計方程的重要特性，同時將該方法推廣到高維數(shù)據(jù)分析，適用于協(xié)變量個數(shù)p隨樣本例數(shù)n同階變化的情況。

2.聚類分析

聚類分析(clustering analysis)能同時考察所有變量，識別變化趨勢一致、功能相似的組學標記物，對于生物機制的研究具有重要意義[8]。時間序列聚類(time series clustering)根據(jù)時間序列相似度對研究對象進行聚類，從而使不同聚類的類間距離最大，類內(nèi)距離最小。模糊C均值聚類(fuzzy c-means clustering) 是動態(tài)組學研究設計中應用最為廣泛的一種算法[9-10]。其核心思想是對j個觀察單位X={X1，X2…Xj}尋找c個模糊組，并求每組的聚類中心，使得目標函數(shù)達到最小。模糊C聚類的優(yōu)點在于能適應分離性不好的數(shù)據(jù)集，允許數(shù)據(jù)性質(zhì)的模糊性，為數(shù)據(jù)結(jié)構描述提供了詳細信息[11]。一項模擬研究表明，相較于K均值聚類算法，模糊聚類算法具有相對較高的聚類效能[12]。但該算法的性能尚依賴于聚類個數(shù)和初始隸屬度矩陣。

3.降維分析

常用于組學數(shù)據(jù)的降維方法包括主成分分析(principal component analysis，PCA)、偏最小二乘判別分析(partial least squares-discriminant analysis，PLSDA)與平行因子分析(parallel factor analysis)等[3]。然而，這些降維方法均未將縱向組學數(shù)據(jù)集的時間順序信息納入模型，即打亂時間順序后仍然能得到相同的結(jié)果，且并不適用于非線性組學數(shù)據(jù)[13]。目前，已開發(fā)的降維方法主要包括基于多水平思想的線性降維方法及基于核函數(shù)的非線性降維方法。

(1)線性降維方法

縱向數(shù)據(jù)具有多水平結(jié)構資料的特征，觀察單位為1水平單位，該觀察單位的重復測量資料為2水平單位。多水平模型的核心思想是通過估計兩個水平上的方差，并考慮解釋變量對方差的影響，充分利用各水平內(nèi)的聚集信息，從而獲得回歸系數(shù)的有效估計，提供正確的標準誤與置信區(qū)間[14-15]。

多水平同步成分分析(multilevel simultaneous component analysis，MSCA) 是結(jié)合多水平思想與主成分分析思想的降維方法[16-17]。多水平同步成分分析模型將數(shù)據(jù)集總變異分為個體間和個體內(nèi)兩個水平上的變異。在相同的成分數(shù)下，PCA與MSCA相比，解釋的變異相同或更多，對MSCA施加限制越多則解釋變異越小。但MSCA相較于PCA的可解釋度更好，其中不同的亞模型能夠解釋數(shù)據(jù)中不同的變異，其個體內(nèi)模型比PCA能更好地展示數(shù)據(jù)中的動態(tài)變異，而個體間模型又比PCA更好地展示個體間的非動態(tài)變異。

多水平偏最小二乘判別分析(multilevel PLS-DA，ML-PLS-DA) 的核心思想是將個體間和個體內(nèi)兩水平的變異分開解釋[18-19]。數(shù)據(jù)分析時，首先將個體間變異和個體內(nèi)變異的部分分開，其中個體間變異是對個體兩次測量的均值分析得到，而個體內(nèi)變異是對個體前后兩次測量的差值進行分析。當使用ML-PLS-DA來描述個體內(nèi)變異時，主要是關注個體間相同的處理效應。因此，ML-PLS-DA的第一主成分描述主要效應，其后的成分描述個體間不同的處理效應。

(2)非線性降維方法

線性降維方法計算簡便，原理簡單，易于解釋。但該類方法在處理非線性組學數(shù)據(jù)時，仍然存在一定局限性。為更精確挖掘非線性組學數(shù)據(jù)信息，研究者在傳統(tǒng)降維分析方法中引入非線性。其中，應用最廣泛的是基于核函數(shù)的主成分方法(kernel principal component analysis)[20]，核主成分分析方法將原始數(shù)據(jù)空間RM中的樣本x映射至特征空間F(Φ：RM→F；x→X)，在特征空間內(nèi)對樣本X實現(xiàn)主成分分析[20-21]。該方法較傳統(tǒng)主成分分析方法有以下優(yōu)點[22]：① 引入了非線性映射函數(shù)Φ，將原始數(shù)據(jù)映射至特征空間，能夠更好地解釋原始數(shù)據(jù)中非線性變異部分；②可使用不同核函數(shù)，對不同種類非線性組學數(shù)據(jù)實現(xiàn)降維；③能提供比主成分分析更多的特征數(shù)目，可以最大限度地提取特征信息。與核主成分分析類似，核偏最小二乘判別分析(kernel partial least squares-discriminant analysis)也是利用核函數(shù)方法將原始數(shù)據(jù)映射至特征空間，以描述特征間的非線性關系[23]。除基于核函數(shù)的分析方法外，非線性降維方法還包括局部線性嵌入、等距映射、多尺度變換等流形學習方法。

縱向組學數(shù)據(jù)統(tǒng)計分析策略

縱向組學研究設計在成為研究生命體功能變化有力手段的同時，也給統(tǒng)計分析帶來了新的機遇與挑戰(zhàn)。研究者們在進行縱向組學數(shù)據(jù)分析過程中，常常忽略了數(shù)據(jù)集的時序性及相應統(tǒng)計分析方法原理與前提假設，降低了研究結(jié)論的可靠性。本文針對縱向組學數(shù)據(jù)特點，探索了組學統(tǒng)計分析策略，具體總結(jié)如圖1所示。

高通量測量技術中的實驗環(huán)境、儀器性能及人工操作等均會對數(shù)據(jù)質(zhì)量產(chǎn)生影響。因此，組學測量數(shù)據(jù)變異來源廣泛，除生物變異外，還包括環(huán)境影響及測量誤差等。目前常用數(shù)據(jù)預處理方法包括噪聲濾除、基線校正及標準化處理。

圖1 縱向組學數(shù)據(jù)統(tǒng)計分析策略流程圖

單變量統(tǒng)計分析思想簡單，易于理解，常用于快速篩選組學研究中隨時間動態(tài)變化的組學標記物[24-25]。在單變量統(tǒng)計分析過程中，重復測量方差分析對數(shù)據(jù)資料要求極為嚴格，須同時滿足方差分析條件及協(xié)方差陣球?qū)ΨQ性，混合效應模型及廣義估計方程能考慮到縱向數(shù)據(jù)的相關性，處理缺失值問題，但后者無法處理高維縱向數(shù)據(jù)的非平衡性問題，存在一定局限性。此外，高維情境常涉及多重比較問題，需著重考慮對假設檢驗的檢驗水準α進行校正，目前常用校正方法包括Bonferroni校正法及FDR(false discovery rate) 校正法[26-27]。

聚類方法常用于組學標記物時序變化趨勢分析。在經(jīng)單變量分析后，可通過聚類算法篩選出規(guī)律變化的不同類組學標記物，對不同類的組學標記物選擇不同的模型進行研究。K均值聚類算法計算簡便快捷，適應性廣，但在聚類過程中未考慮到縱向數(shù)據(jù)的時間序列信息；有序樣品聚類大多用于樣品聚類，是一種特殊條件系統(tǒng)聚類算法，但難以直接得出相關序列特征的結(jié)論；模糊聚類算法用于時間序列數(shù)據(jù)時，以時間為維度，計算隸屬度，同時允許了數(shù)據(jù)性質(zhì)的模糊性，為數(shù)據(jù)結(jié)構的描述提供了詳細的信息。

主成分分析與偏最小二乘法判別分析是組學研究中常用的降維方法。線性降維方法計算簡便，原理簡單，易于解釋。在降維的同時，考慮到數(shù)據(jù)的時序信息，能更好地展示數(shù)據(jù)集的動態(tài)變異。非線性降維方法在降維分析中進一步引入非線性，更精確構建判別模型，挖掘非線性數(shù)據(jù)信息。最終，使用外部測試集評價潛在組學標記物預測效果，探索潛在組學標記物的生物學功能，為分析結(jié)果提供合理的生物學解釋。

由于縱向組學數(shù)據(jù)的復雜特性，上述分析手段能在一定程度上解決組學數(shù)據(jù)統(tǒng)計分析問題，但仍存在局限性。目前組學標記物變量篩選的方法主要依靠單變量統(tǒng)計分析及后續(xù)改進的偏最小二乘法判別分析等方法。數(shù)據(jù)發(fā)生輕微變化時，變量篩選效果也會受到影響，因此建立穩(wěn)定有效的縱向高維數(shù)據(jù)變量篩選方法仍然值得研究者們進行探討[28]。此外，利用縱向組學數(shù)據(jù)進行因果推斷分析時，如何校正時依混雜因素對因果效應估計的影響也需要進一步研究[29-30]。以上兩個關鍵科學問題的解決將會對縱向組學數(shù)據(jù)提供新的思路與契機。