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

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

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

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

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

1 材料與方法

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

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

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

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

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

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

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

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

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

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

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

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

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

2 結(jié) 果

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

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

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

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

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

圖1 BA樣本分子亞型分類

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

將所有Magenta 模塊內(nèi)的基因?qū)隨TRING 數(shù)據(jù)庫(kù)以構(gòu)建PPI 網(wǎng)絡(luò)，通過(guò)MCC 算法發(fā)現(xiàn)了10 個(gè)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)建

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

圖4 10個(gè)核心基因在不同亞型間的表達(dá)分布情況

3 討 論

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

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