融合生物學通路的變分自編碼器在肺癌蛋白與代謝組學數(shù)據(jù)中的應(yīng)用研究*

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

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

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

原 理

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的第一項可以通過小批量數(shù)據(jù)的蒙特卡羅抽樣得到,若訓練神經(jīng)網(wǎng)絡(luò)的批次大小為M,即抽樣M次,則該項計算如下:

(2)

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

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

(3)

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

(4)

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

圖1 變分自編碼器結(jié)合生物學通路

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

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

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

實例分析

1.實驗數(shù)據(jù)集來源與整理

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

2.實例分析結(jié)果

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

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

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

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

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

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

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

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

討 論

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

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

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

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