陳治宏，羅鳳媛，馬麗媛，李紅東

（1. 贛南醫(yī)學(xué)院2019級碩士研究生；2. 贛南醫(yī)學(xué)院2020級碩士研究生；3. 贛南醫(yī)學(xué)院2021級碩士研究生；4. 贛南醫(yī)學(xué)院信息工程學(xué)院，江西 贛州 341000）

膽道閉鎖（Biliary atresia，BA）是新生兒常見的一種進行性阻塞性膽道疾病，臨床表現(xiàn)為梗阻性黃疸，如果BA 患兒得不到及時治療，易發(fā)展為膽汁淤積性肝硬化、門靜脈高壓，最終導(dǎo)致肝衰竭而死亡［1］。膽道閉鎖具有地區(qū)和種族差異，一般認為亞洲人發(fā)病率較高，尤以日本和中國的發(fā)病率最高［2］。目前對BA 患兒應(yīng)用最廣泛且療效較好的治療手段是Kasai 手術(shù)，但由于進行性肝纖維化、膽汁淤積性肝硬化和門脈高壓的出現(xiàn)，多數(shù)患兒最終須選擇肝臟移植救治［3］。BA 因其早期診斷困難、病理改變特殊以及治療效果不理想，已成為威脅新生兒生命健康的重要危險因素。

目前，研究者提出了許多理論來解釋BA的病因，包括遺傳變異、病毒感染和免疫介導(dǎo)等。CHENG G等［4］對89 例BA 患者開展全基因組拷貝數(shù)突變分析發(fā)現(xiàn)了29 個BA 相關(guān)的拷貝數(shù)變異，其中包括富集在炎癥相關(guān)通路的基因，提示了炎癥在BA 中的重要作用。SHEN C 等［5］研究表明，60%的BA 患兒在診斷時檢測到巨細胞病毒（CMV）的DNA，同時感染CMV 的BA 患兒患有膽管炎和較高的膽管纖維化程度。SQUIRES J E 等［6］通過對輪狀病毒（RV）誘導(dǎo)的BA小鼠模型研究發(fā)現(xiàn)，BA小鼠的NK細胞數(shù)量增加會促進慢性膽管炎癥。SAXENA V 等[7]發(fā)現(xiàn)RV 誘導(dǎo)的BA 小鼠模型中的漿細胞樣樹突狀細胞（pDCs）顯著增加，同時pDCs產(chǎn)生白細胞介素15(IL-15)激活CD80，加速膽管纖維化進程。BA 的復(fù)雜致病因素，提示BA中可能存在不同的分子亞型。

目前，基于高通量數(shù)據(jù)對疾病進行分子亞型的研究主要是根據(jù)樣本基因表達水平的相似性。通過無監(jiān)督聚類的方法，例如非負矩陣分解（NMF）、層次聚類和潛在因子分析（LF）等挖掘潛在的分子亞型［8-9］。但這類方法往往基于較為復(fù)雜的數(shù)量統(tǒng)計模型，難以理解和進行生物學(xué)解釋，且不易進行臨床實踐和驗證。有研究表明，基因間相對表達秩次（大?。╆P(guān)系在正常組織樣本內(nèi)存在廣泛的穩(wěn)定性，在疾病組織樣本中受到了廣泛的擾動。基于這一現(xiàn)象，本文提出了一種基于樣本內(nèi)基因間相對表達秩次關(guān)系的方法識別潛在疾病分子的亞型方法。這種方法不僅易于理解和解釋，還可整合來自不同實驗室的樣本進行生物信息學(xué)分析，對基因表達水平的個體間生物變異具有魯棒性[10-11]。本研究將該方法應(yīng)用于挖掘BA 潛在的分子亞型，進而識別亞型相關(guān)的核心基因。

1 材料與方法

1.1 數(shù)據(jù)來源及其預(yù)處理本文所分析的基因表達譜數(shù)據(jù)集，GSE46960 和GSE15235 均來自基因表達綜合數(shù)據(jù)庫（GEO，https：//www.ncbi.nlm.nih.gov/geo/）。GSE46960數(shù)據(jù)集檢測平臺為GPL6244，共檢測了64 例BA 患兒的肝組織、7例正常嬰兒的肝組織。GSE15235 數(shù)據(jù)集檢測平臺為GPL570，共檢測了43例BA患兒的肝組織。

各數(shù)據(jù)集均從GEO 下載原始表達譜數(shù)據(jù)（. CEL 文件），利用R 軟件Affy 軟件包提供的RMA（Robust Multichip Analysis）算法對其進行背景校正［12］，并利用相應(yīng)的平臺注釋文件將探針I(yè)D 對應(yīng)到基因ID 上。如果多個探針映射到同一基因，則取探針的平均表達值作為該基因的表達水平。最終，數(shù)據(jù)集GSE46960 共檢測18 859 個基因，數(shù)據(jù)集GSE15235 共檢測20 684 個基因。兩套數(shù)據(jù)集共同檢測的17 179個基因作為背景基因用作后續(xù)分析。

1.2 基于基因間相對表達秩次關(guān)系識別BA 分子亞型假設(shè)樣本內(nèi)基因間相對表達秩次關(guān)系與疾病的分子亞型相關(guān)，那么應(yīng)存在以下幾種現(xiàn)象：⑴在不同分子亞型樣本內(nèi)，存在基因間相對表達秩次關(guān)系不同的現(xiàn)象。⑵具有潛在鑒別疾病分子亞型能力的基因?qū)?，按其相對表達秩次關(guān)系劃分的兩組樣本間存在差異表達的基因，同時涉及生物學(xué)功能的改變。⑶識別具有鑒別疾病分子亞型能力的潛在基因?qū)Υ嬖诰奂F(xiàn)象。

基于上述假設(shè)，本文識別疾病分子亞型的具體步驟如下。

1.2.1 篩選正常組織樣本中基因間相對表達秩次關(guān)系穩(wěn)定的基因?qū)⒈尘盎騼蓛山M合，如果一對基因（基因i 和基因j）在正常樣本內(nèi)滿足P（Gi＞Gj）≥0.9［公式⑴]，那么該基因?qū)Χx為穩(wěn)定基因?qū)Α?/p>

其中n表示所有的正常樣本數(shù)，t代表第t個樣本，Gi＞Gj表示基因i和基因j的表達秩次關(guān)系。

1.2.2 基于穩(wěn)定基因?qū)Y選疾病組織樣本中的顯著逆轉(zhuǎn)基因?qū)谡颖局械玫降姆€(wěn)定基因?qū)Γu估其在疾病中相對表達秩次關(guān)系改變情況。對于任一穩(wěn)定對（Gi，Gj），假設(shè)在正常樣本中觀察到Gi＞Gj的樣本有n1個、Gi≤Gj的樣本有n2個，在疾病樣本中滿足Gi＞Gj的樣本有m1個、Gi≤Gj的樣本有m2個。利用Fisher 精確檢驗，可檢驗Gi、Gj間相對表達秩次關(guān)系在正常和疾病樣本間分布是否存在差異。經(jīng)多重假設(shè)檢驗校正后，滿足假陽性發(fā)現(xiàn)率（FDR）＜5%的基因?qū)Χx為逆轉(zhuǎn)基因?qū)Α?/p>

1.2.3 篩選具有潛在疾病分子亞型鑒別能力的基因?qū)τ谠诩膊≈泻Y選得到的逆轉(zhuǎn)基因?qū)?，按如下步驟進一步篩選具有潛在疾病分子亞型鑒別能力的基因?qū)Γ孩拍孓D(zhuǎn)對的秩次關(guān)系的改變只在20%～80%的疾病患者之中發(fā)生。這是為了確?；?qū)哂需b別疾病分子亞型的能力。⑵對于任意逆轉(zhuǎn)對（Gi，Gj），根據(jù)其相對表達秩次關(guān)系，可將疾病樣本分為兩組，Gi＞Gj為1 組，Gi≤Gj為1 組。然后利用t檢驗，控制FDR 為5%，識別兩組間的差異基因，并基于KEGG（京都基因與基因組百科全書，https：//www. kegg. jp/）通路數(shù)據(jù)庫對兩組間識別的差異基因進行功能富集分析。若一個生物學(xué)通路顯著富集了差異基因，則認為該通路發(fā)生了生物學(xué)功能的擾動。對于一個逆轉(zhuǎn)基因?qū)Γ砂雌湓诩膊颖局械南鄬Ρ磉_秩次關(guān)系將疾病樣本分成兩組。如果這兩組樣本間存在差異表達基因，同時存在生物學(xué)功能的擾動，則認為該基因?qū)哂需b別潛在疾病分子亞型的能力，稱其為分子亞型候選基因?qū)Α?/p>

1.2.4 分子亞型識別首先，按照分子亞型候選基因?qū)_動通路的數(shù)目，將基因?qū)M行降序排列。然后，構(gòu)建一個分子亞型候選基因?qū)Α良膊颖镜木仃?，其中元素rtij為矩陣中第t個樣本中基因?qū)Γ╥，j）的相對秩次關(guān)系，其值為1 或0（如果Gi＞Gj，則為1，如果Gi≤Gj，則為0）。以1 000 為梯度，分別取前k個基因?qū)Γ╧=1 000*n，n=1，2，3…）對應(yīng)的矩陣，利用歐式距離進行聚類分析，識別疾病分子亞型。為了保證結(jié)果的可靠性，對k個基因?qū)垲惤Y(jié)果，利用隨機實驗進行驗證：隨機從所有分子亞型候選基因?qū)χ刑暨xk 個基因?qū)M行聚類分析，重復(fù)100 次。如果隨機分類結(jié)果與前k個基因?qū)Φ姆诸惤Y(jié)果一致率＜95%，則重新選擇標記。反之，將其作為分子亞型預(yù)測標記，用作后續(xù)分析。

1.3 基于加權(quán)基因共表達網(wǎng)絡(luò)識別疾病分子亞型相關(guān)基因模塊基于加權(quán)基因共表達網(wǎng)絡(luò)（WGCNA）分析提取與疾病亞型顯著相關(guān)的模塊（用R 包WGCNA 實現(xiàn)）［13］。首先，對背景基因表達矩陣相關(guān)系數(shù)進行加權(quán)，使基因間的相互作用關(guān)系符合無標度分布。然后對基因進行分類，并將具有相似表達模式的基因分成一個模塊，最小模塊大?。╩in-ModuleSize）設(shè)置為100，其他參數(shù)設(shè)置為默認值。同一模塊的基因往往表現(xiàn)出相似的表達模式和功能［14］。

1.4 模塊基因功能富集分析基于京都基因和基因組百科全書（KEGG）數(shù)據(jù)庫對模塊基因進行通路富集分析，用R包clusterProfile實現(xiàn)［15］。

1.5 PPI 網(wǎng)絡(luò)構(gòu)建和核心基因篩選基于蛋白質(zhì)互作在線數(shù)據(jù)庫STRING（www. string-db. org）對候選模塊內(nèi)的基因進行蛋白與蛋白之間的互作分析（PPI），利用網(wǎng)絡(luò)最大集團度（MCC）評分算法［公式⑵］［16］，識別PPI網(wǎng)絡(luò)中的關(guān)鍵基因。

其中，集合S表示集合中元素的數(shù)量，S（ν）是包含ν 的最大集團的集合，（|C|-1）！是所有小于|C|的正整數(shù)的乘積。

2 結(jié) 果

2.1 具有潛在疾病分子亞型鑒別能力的基因?qū)Π凑辗椒▽W(xué)篩選具有潛在疾病分子亞型鑒別能力的基因?qū)α鞒??；跀?shù)據(jù)集GSE46960，在每個樣本中分別對基因的表達水平進行兩兩比較，形成基因?qū)Α翗颖镜?-1 矩陣X。Xij表示第i 個基因?qū)Γ℅a，Gb）在第j 個樣本中的相對表達秩次關(guān)系，其值為1或0（1代表Ga＞Gb，0代表Ga≤Gb）。在7例正常嬰兒樣本中，有144 822 185個基因?qū)υ?0%以上的正常樣本組織中具有Ga＞Gb的表達模式，其中有5 105 085個基因?qū)Φ南鄬Ρ磉_秩次關(guān)系在BA 樣本中發(fā)生顯著逆轉(zhuǎn)Ga≤Gb（Fisher 精確檢驗，F(xiàn)DR＜0.05），暗示了這些基因?qū)Φ南鄬Ρ磉_秩次關(guān)系可能與BA 的發(fā)生發(fā)展相關(guān)。對每一個逆轉(zhuǎn)基因?qū)ΓM一步根據(jù)其在疾病樣本中的相對表達秩次關(guān)系，將疾病樣本分為兩組（Ga＞Gb為一組，Ga≤Gb為一組），然后識別兩組間的差異基因并基于KEGG數(shù)據(jù)庫進行功能富集分析，發(fā)現(xiàn)有4 836 202 個基因?qū)M足條件，用作后續(xù)分析。

2.2 識別BA 的分子亞型基于亞型候選基因?qū)Γ紫葮?gòu)建亞型候選基因?qū)Α良膊颖镜?-1 矩陣X。矩陣中xij表示第i 個基因?qū)υ诘趈 個樣本中的相對表達秩次關(guān)系，其值為1 或0（1 代表＞，0 代表≤）。按照亞型候選基因?qū)Ω患降纳飳W(xué)功能通路個數(shù)由大到小進行排序，然后分別選取前1 000 對、2 000 對、3 000 對、4 000 對和5 000 對基因?qū)Γ瑢?4例BA 樣本分別進行聚類分析。結(jié)果顯示前3 000 2A），采用動態(tài)剪切法劃分模塊，合并相似度＞75%的模塊，最終構(gòu)建了17 個共表達模塊（圖2B）。進一步分析了BA 不同亞型與各種基因模塊的相關(guān)性，發(fā)現(xiàn)洋紅色（Magenta）模塊與BA 亞型最相關(guān)（|r|=0.58，P=6×10-7）（圖2C）。

圖2 BA亞型基因共表達網(wǎng)絡(luò)模塊的構(gòu)建

基于KEGG數(shù)據(jù)庫，對Magenta模塊中涉及的基因進行富集分析，發(fā)現(xiàn)這些基因主要涉及PI3K-Akt信號通路、黏附斑激酶通路（FAK）和ECM 受體相互作用通路等（圖3A）。對、4 000 對、5 000 對的聚類結(jié)果與隨機情況一致。因此，本文將前3 000 對基因?qū)Φ木垲惤Y(jié)果作為BA分子亞型識別依據(jù)。聚類結(jié)果顯示64例BA 樣本明顯聚集為兩類：一類樣本個數(shù)為54 個，另一類樣本個數(shù)為10個（圖1A）。

為了進一步驗證BA 中存在的不同分子亞型，我們把得到的前3 000 個亞型候選基因?qū)?yīng)用于GSE15235 數(shù) 據(jù)集中，GSE15235 數(shù)據(jù)集中43 例BA樣本中依然明顯聚集為兩大類：一類樣本個數(shù)15個，另一類樣本個數(shù)為28 個（圖1B）。通過數(shù)據(jù)集GSE46960、GSE15235 間兩類樣本間的差異基因比較發(fā)現(xiàn)，在GSE46960 兩類間識別了436 個差異基因，GSE15235 兩類間識別了584 個差異基因，共交疊了75個差異基因（超幾何分布P=4.8×10-32），并且這75 個差異基因有69 個基因在兩組間失調(diào)方向一致（二項分布P=5.77×10-15）。

圖1 BA樣本分子亞型分類

2.3 BA 疾病亞型共表達模塊構(gòu)建及核心基因識別為了進一步分析BA 中不同亞型的基因表達模式，利用WGCNA 方法對數(shù)據(jù)集GSE46960構(gòu)建加權(quán)共表達網(wǎng)絡(luò)。使用無標度拓撲標準選擇β=8（圖

將所有Magenta 模塊內(nèi)的基因?qū)隨TRING 數(shù)據(jù)庫以構(gòu)建PPI 網(wǎng)絡(luò)，通過MCC 算法發(fā)現(xiàn)了10 個BA 亞型相關(guān)的核心基因（LUM、COL6A3、FBN1、SPARC、DCN、LAMA4、FAP、ANTXR1、LAMA2 和COL1A2）。LAMA2、LAMA4 和LUM 基因（圖中紅色基因）與其他蛋白質(zhì)的相互作用更頻繁（圖3B）。

圖3 Magenta模塊中基因的KEGG通路富集分析和PPI網(wǎng)絡(luò)構(gòu)建

對10個核心基因在GSE46960數(shù)據(jù)集的正常組織、BA亞型Ⅰ和亞型Ⅱ三組間的表達情況進行分析發(fā)現(xiàn)，相對于正常組織，10 個基因在亞型Ⅰ中顯著高表達，而在亞型Ⅱ中顯著低表達（圖4A）。在GSE15235 數(shù)據(jù)集中，我們還發(fā)現(xiàn)這10 個核心基因的表達水平在纖維化和炎癥兩類樣本間有顯著差異，且纖維化組的表達水平顯著高于炎癥組（圖4B），與BA亞型Ⅰ和亞型Ⅱ一致。

圖4 10個核心基因在不同亞型間的表達分布情況

3 討 論

本研究基于樣本內(nèi)基因間相對表達秩次關(guān)系，開發(fā)了疾病分子亞型識別算法，確定了在BA 疾病中存在兩組不同的疾病亞型，并利用WGCNA 分析篩選出了與亞型相關(guān)共表達基因模塊。通路富集分析顯示，Magenta 模塊內(nèi)基因主要涉及PI3K-Akt信號通路、黏附斑激酶通路（FAK）和ECM 受體相互作用通路等。有研究表明,PI3K-Akt 信號通路可能抑制NF-κB 和NLRP3 炎癥通路的蛋白表達，從而減少炎癥因子的分泌[17]。FAK 和ECM 受體相互作用通路抑制TGF-β1信號傳導(dǎo)，降低成纖維細胞分化[18]。最后，本文還識別了與BA 分子亞型相關(guān)的10 個核心基因。其中，LAMA2 和LAMA4 可促進轉(zhuǎn)化生長因子（TGF-β）信號傳導(dǎo)，而TGF-β 信號會阻止慢性損傷期間促纖維化FAP 的細胞凋亡，從而促進纖維化環(huán)境［19-20］。KRISHNAN A 等［21］研究表明，LUM 可在體外誘導(dǎo)生長因子β1（TGF-β1）轉(zhuǎn)化為促纖維化細胞因子。COL1A2和COL6A3編碼膠原蛋白，而膠原蛋白的過度表達是肝纖維化的基本特征。當肝星形細胞（HSCs）被激活時，基質(zhì)金屬蛋白（MMPs）的表達增加，同時MMPs 的抑制劑基質(zhì)金屬蛋白酶（TIMPs）的表達增加。如果TIMPs 增加過快，MMPs/TIMPs 的比值會發(fā)生變化，使細胞外基質(zhì)（ECM）的合成和降解失衡，促進纖維化的發(fā)展［22-23］。其余5個關(guān)鍵基因（FBN1、SPARC、DCN、FAP 和ANTXR1）雖然還未有其與BA 相關(guān)報道，但其參與了蛋白質(zhì)消化吸收、ECM-受體相互作用和PI3K-Akt 信號通路。其中ECM 的變化由多種介質(zhì)和生長因子誘導(dǎo)，從而調(diào)節(jié)各種效應(yīng)，如刺激血管生成、炎癥反應(yīng)和促進基質(zhì)侵襲，可導(dǎo)致成纖維細胞的遷移或分化，最終使細胞纖維化［24-25］。提示這10 個核心基因與BA 進行性肝炎、肝纖維化有關(guān)。

在本研究中，正常樣本相對較少，在穩(wěn)定基因?qū)Y選過程中，可能引入一定的假陽性基因?qū)?，后續(xù)需要更多的正常樣本來驗證。同時，本研究篩選的亞型候選基因?qū)?shù)量較多，不利于后續(xù)聚類鑒別疾病亞型。因此，如何篩選出更少、更具有生物學(xué)意義的亞型候選基因?qū)κ俏磥淼难芯恐攸c。后續(xù)研究我們將考慮如基因在候選對中出現(xiàn)的頻率、基因間的表達相關(guān)性等因素進一步改進算法。