融合生物學通路的變分自編碼器在肺癌蛋白與代謝組學數據中的應用研究*

劉芝霖 榮志煒 俞軼培 邱滿堂 侯艷,4△

【提 要】 目的 本研究提出了融合生物學通路的變分自編碼器(variational auto encoder,VAE),對蛋白質與代謝組學數據進行整合分析,并應用于肺腺癌患者探索可能的病理機制。方法 為VAE的隱變量層節點賦予通路的實際意義,解碼器按通路包含的生物學分子信息構建稀疏神經網絡,使通路節點只與自身包含的分子連接,將隱變量作為提取的高級特征。對隱變量進行Kmeans聚類分析并使用調整蘭德系數評估效果,引入基因差異表達分析方法limma探索差異表達通路,在北京大學人民醫院胸外科肺腺癌患者的蛋白質與代謝組學數據中進行實例分析。結果 融合生物學通路的VAE提取的高級特征不僅將聚類準確度提高了38%,還通過差異表達分析鑒別出了實性與亞實性結節肺腺癌間的差異通路。結論 融合生物學通路的VAE可用于組學數據整合分析,其提取的高級特征具有通路表達活性的實際生物學意義。

隨著各種組學測序技術的發展,產生了基因組、轉錄組、蛋白組、代謝組等各種組學數據,這些數據從不同組學層面描述了癌癥的異質性,有利于癌癥亞型分類研究及生物標志物探索,對于癌癥的早期診斷、臨床治療和預后都有重要意義[1]。變分自編碼器(variational auto encoder,VAE)[2]作為一種深度學習模型,具有組合低級特征為抽象的高級特征的功能,其強大的建模非線性數據關系的能力在組學數據分析中顯示出極大的優越性,但缺乏解釋性的問題也成為其在醫學領域應用的一大阻礙[3],需要新的設計以提高模型的解釋性。因此,本文提出融合生物學通路的VAE[4],將通路信息內置于神經網絡的架構中,使模型提取的高級特征同時具有表征通路整體狀態的實際意義。在早期肺腺癌患者的蛋白質與代謝組數據的實際應用中,揭示了亞實性結節與實性結節肺腺癌及癌旁樣本間的生物學關系,為癌癥機制研究提供了有生物學意義的分析結果。

原 理

ELBO=Eq(z|x)(logp(x|z))-KL(q(z|x)‖p(z))

(1)

其中q(z|x)用于近似后驗分布的變分分布,p(x|z)是隱變量生成樣本的條件分布,KL表示KL散度(Kullback-Leibler divergence)[6],是衡量兩個分布差異的一種度量,這里用于量化分布q(z|x)與先驗分布p(z)間的差異。ELBO的第一項可以通過小批量數據的蒙特卡羅抽樣得到,若訓練神經網絡的批次大小為M,即抽樣M次,則該項計算如下:

(2)

當樣本生成的p(x|z)服從正態分布時,最大化該項等價于最小化重構數據的均方誤差。為了能夠在訓練中應用標準的反向傳播,需要使用重參數化技巧,假設q(z|x)服從正態分布:

q(z|x)=N(μ(x),Σ(x))

(3)

其中Σ(x)為單位對角矩陣,μ(x)和Σ(x)由編碼器學習得到,然后在q(z|x)中采樣,計算:

(4)

隱變量作為高級特征通常沒有實際意義,這也是VAE等深度學習模型被稱為“黑箱”的原因之一,缺乏生物學上的解釋性限制了其應用場景。本研究為隱變量賦予了通路整體表達狀態的實際意義,設定每個隱變量節點代表一個通路,使隱變量節點與生物學實體相對應,根據通路中是否包含某一生物學分子決定隱變量是否連接某一特征。這種稀疏連接的方式,能夠讓每個隱變量節點只匯聚指定特征的信息,強制其表征一組分子共同作用的情況,提高了模型的解釋性(圖1)。

圖1 變分自編碼器結合生物學通路

AE也是由編碼器解碼器組成,但AE得到的隱變量分散在多塊不連續的低維流形上,而VAE假設隱變量服從先驗分布p(z),其隱變量空間平滑連續[2],更符合生物學通路表達的實際情況。解碼器作為數據生成過程的條件分布p(x|z)能夠解析數據的結構,而編碼器用于近似q(z|x),需要足夠復雜的神經網絡的強大擬合能力,因此本模型只在解碼器上使用了稀疏連接的結構設計。此外,由于生物學通路信息涉及多個層級的組學分子相互作用,在有合適通路信息的情況下,本模型能適用于多個組學的數據的整合分析。

對于稀疏連接神經網絡的具體構建,本研究構造了一個01掩碼矩陣M來實現,如果通路i中包含分子j,則Mi,j為1,否則為0。在網絡的前向傳播與反向傳播中,網絡的權重都會乘以這個掩碼矩陣,使得掩碼為0的權重無論是多少都強制為0,而且梯度也只在通路中的指定特征上計算,從而實現代表通路的隱變量只與通路內分子特征連接的效果。在生物學通路信息中,有的分子可能具有多種功能,從而包含在不同的通路中,為了使其在各通路中都被學習,采用神經網絡的暫退法(dropout)與隨機梯度下降,每次訓練只使用少量擋本計算梯度并隨機丟棄一些節點,在提高神經網絡泛化能力[7-8]的同時也使每個包含該分子的通路都充分學習了其信息。在實際數據中訓練時,由于在生物學上的認知有限,存在一些不屬于任何已知通路的特征分子,因此設定了少量稱為“缺失通路”的隱變量節點來連接這些特征,保證了數據能夠較好的重建,同時也作為額外的通路信息以解釋更多的數據變異。此外,為了提高模型的泛化能力,減少對高維組學數據的過擬合,除了使用dropout與隨機梯度下降,也需要根據實際數據控制編碼器解碼器網絡層數與節點數以調整模型復雜度。

當探索不同標簽的樣本間差異時,特征分子或通路表達上的變化值得關注,它們是癌癥機制研究的重要線索,需要使用差異表達分析的方法。由于隱變量一般被假設為正態分布,其學習到的參數中μ作為均值代表了該通路的平均表達情況,本研究從分析差異基因表達的線性回歸方法R包limma[9]得到啟發,將其運用到通路表達的差異分析上,判斷通路表達狀態在兩組間是否存在差異。此外,提取的隱變量亦可以用于其他的下游分析,如降維可視化、分析與生存時間的關聯等。本研究使用的模型由python 3.9及pytorch深度學習模塊搭建,隱變量的統計分析及作圖使用R 4.1.3完成。

實例分析

1.實驗數據集來源與整理

本研究使用來自北京大學人民醫院胸外科早期肺腺癌患者的蛋白組與代謝組數據,其標簽是醫生通過CT影像表現確定的,分為實性結節肺腺癌、亞實性結節肺腺癌與癌旁正常組織共三類樣本。組學數據由非標記定量蛋白質組學和非靶向代謝組學方法測定,所有測定均經過患者的知情同意。獲得數據后,首先將蛋白與代謝物的名稱轉換為通路數據庫中的ID,然后對數據進行標準化。排除只有一個組學信息的樣本,將兩個組學矩陣拼接,得到85個患者與3240分子特征的矩陣,其中實性結節肺腺癌樣本37個,亞實性結節肺腺癌樣本18個,癌旁正常組織樣本30個;蛋白質分子特征2946個,代謝物分子特征294個。通路信息來自京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes,KEGG)數據庫[10],通過R包KEGGREST獲取人類的通路,其名稱由hsa加五位的數字編碼構成,排除無關的通路后,得到77條通路信息用于設定通路節點與特征分子的連接情況。

2.實例分析結果

為了可視化隱變量對數據的分類效果,使用統一流形逼近和投影(uniform manifold approximation and projection,UMAP)方法對隱變量進行降維,該方法能夠將高維數據降維并保留點之間的相互關系,本研究使用R包umap采用默認設置完成這一過程。在降至2維后,由圖2可以發現正常的癌旁組織與癌癥組織的差異較大,亞實性與實性的肺腺癌差異雖小,但還是明顯分為了兩部分,表明該方法提取的特征即使是對于早期癌癥也能顯示出差異。

圖2 顯示樣本分型的隱變量UMAP圖

使用Kmeans聚類方法進行無監督聚類,比較特征提取前后的變化,以評估特征提取效果。評估指標為調整蘭德系數(adjusted rand index,ARI),其用于比較聚類結果與真實結果的差異,取值范圍為0～1,越接近1則表明聚類效果越好。當分類為三類時,ARI在原始數據上為0.45,在隱變量上為0.62,提高了38%。當分類僅為兩類,即癌與非癌時,ARI由原始數據上的0.77,提升為隱變量上的0.86,增幅為12%?？梢婋[變量對特征進行了優化,能夠更好地對不同組織類型進行表征。

使用基于線性回歸的差異表達分析方法limma對較難區分的實性與亞實性樣本進行差異表達分析。對于原始數據,雖然也能分析出差異分子特征如SCEL、UBA6等,但這些分子未能在KEGG中顯著富集,在生物學上也沒有文獻表明與癌癥有明顯關聯,該結果解釋性較差。隱變量差異分析結果的火山圖(圖3)顯示hsa00010、hsa00630等代謝通路的表達狀態上呈現顯著差異,P值最小的hsa00010、hsa00630,與其他通路相比更可能與癌癥代謝變化相關,而且該結果由合并通路中所有分子狀態得到,具有更高的可信度。

圖3 隱變量差異表達分析火山圖

hsa00010是與糖酵解與糖異生相關的通路,其調控關系如圖4所示,實性結節肺腺癌與亞實性結節肺腺癌相比,紅色表示上調,綠色表示下調,白色表示因目前技術局限性未測定的分子。由于通路中大部分分子的表達差異在統計學上并不顯著,但將它們當作整體分析時若呈現一致的變化可能會使得通路整體活性表現出統計學差異,為了更好地呈現這種趨勢,上調和下調僅根據表達均值計算。從圖中可以看出大部分分子呈現上調的趨勢,而且在通路調控關系中距離較近,能夠相互影響作用。該結果表明,在分析單個分子的差異時難以發現的一系列分子微小同步變化,累積起來會有不小的效應,導致通路整體表達水平的變化,直接將通路作為整體分析能夠較好地發現這種趨勢的差異。

圖4 實性結節肺腺癌相較于亞實性結節肺腺癌在糖酵解/糖異生通路的變化

hsa00630是乙醛酸和二羧酸代謝通路(圖5),通路里的大量分子也顯示出同步上調,帶來通路整體水平的變化,然而其在火山圖里呈現出與hsa00010不同方向的變化,這是因為作為通路狀態的高級特征值的上升不一定代表通路活性的上調,高級特征的提取經過了神經網絡的非流形轉換,值的高低不等于通路活性的高低,但值的差異卻可以反映通路活性的差異。

討 論

本研究提出了一種新的VAE架構,采用已知的生物學通路信息構建稀疏解碼器,以此分析不同樣本在通路水平上的活動情況。融合生物學通路信息的VAE能夠將高維的組學數據按照通路信息有序整合為各個通路整體表達水平,提高了表征樣本的能力,同時便于在通路的層級上進行差異表達分析,探索疾病機制。對于這一思路衍生的其他模型,如使用AE代替VAE或在編碼器中使用稀疏連接,實際測試表明它們的效果不如當前的模型,與理論分析的結果一致,最終本研究采用了稀疏連接編碼器的VAE。犧牲小部分數據重建能力,以獲得更有意義的高級特征,這對于醫學方面的應用更加重要,而且未來更全面的通路信息也許能夠實現比全連接更好的數據重建。此外,本研究使用了通路信息作為生物學先驗知識構建神經網絡,其可以推廣到調控網絡、生物學模塊等多種生物學概念上,根據包含的生物學分子構建多種多樣的連接方式,具有廣泛的應用前景。

在癌癥的發生發展中,為了滿足腫瘤細胞的生物能量、生物合成和氧化還原需求,會對代謝通路的調控機制進行重編程,導致異常的通路活性[11],所以代謝通路的改變在癌癥機制分析中十分重要,能夠發現重要的生物標志物,這一點在肺腺癌中也得到了證實[12]。以前的研究也發現亞實性結節的肺腺癌相對于實性結節的具有更好的生存率[13],本研究也得到了相似的結論:基于現有數據的通路活性的差異分析表明,實性結節的樣本中能量相關代謝通路比亞實性結節樣本有更高的活性,可能是導致其進展更快的原因。目前在蛋白質與代謝組學上對實性與亞實性肺腺癌的研究還較少,本文能夠為相關的機制研究提供可能的線索。

本研究使用的數據來自早期癌癥的患者,所以數據本身的區分度還不夠大,特別是在實性與亞實性間,但本方法還是能夠發現一些與致癌機制相關的線索,提升了樣本間的區分度,證明了方法的有效性。傳統的差異分析方法先分析出差異特征,再根據這些特征做通路富集,然而這種做法存在一定缺陷。本研究中,分析原始數據得到的差異分子不是每一個都得到富集;同時,從隱變量中得到的差異通路包含的分子本身并沒有顯示出統計學上的差異,這些問題都降低了結果的解釋性。本研究發現。一系列同步的“不顯著”變化能夠導致通路總體表達水平的顯著變化[14],因此以通路為整體分析在生物學上更合理,錯誤率也更低[15]。在癌癥等疾病早期患者中,相關的分子表達還未能顯現出差異,但通路內的分子若出現同步微小變化,就能在通路水平上被捕獲到,這對于早期診斷具有重要意義。

組學數據分析方法不僅需要適應數據的復雜性與異質性[16],還要有良好的生物學解釋性。深度學習方法擬合能力強但解釋性較差,基于通路的分析錯誤率低且生物學解釋性好,融合生物學通路的VAE結合了兩者的優勢,對多組學數據的整合分析和精準醫學的發展都能起到重要作用[17]。