基于生物信息學(xué)的膽囊癌差異表達(dá)譜中關(guān)鍵蛋白調(diào)控基因分析

韓偉光，莘瑋，蘇水霞，王青

（空軍軍醫(yī)大學(xué)第二附屬醫(yī)院 普通外科，陜西 西安 710038）

膽囊癌（gallbladder cancer，GBC）作為膽道系統(tǒng)中最常見(jiàn)的惡性腫瘤，具有高度侵襲性，其病死率高且預(yù)后不良[1-2]。GBC占膽道惡性腫瘤的80%～95%，其總發(fā)生率為0.8%～1.2%[3]。以往研究[4]表明，超過(guò)85%的GBC患者被診斷時(shí)患有膽結(jié)石，證實(shí)膽結(jié)石的存在與GBC發(fā)生密切相關(guān)，是GBC的一個(gè)重要危險(xiǎn)因素。然而，膽結(jié)石患者中少于3%會(huì)發(fā)展為GBC，這表明需要額外的因素來(lái)驅(qū)動(dòng)GBC的進(jìn)展。因此，解析GBC發(fā)生的潛在分子生物學(xué)機(jī)制對(duì)GBC的診斷及治療具有重要的臨床意義。

在組織學(xué)上，GBC主要分為腺癌和鱗癌兩種，分別占大約90%和1%～12%[5-6]。手術(shù)切除是GBC患者可能治愈的一線選擇，但由于其在早期發(fā)病階段通常無(wú)明顯癥狀，因此往往難以進(jìn)行早期診斷[7-8]。當(dāng)出現(xiàn)黃疸和疼痛等明顯癥狀時(shí)，患者大多已經(jīng)處于晚期，并伴有浸潤(rùn)和轉(zhuǎn)移，普遍錯(cuò)過(guò)了接受根治性手術(shù)的機(jī)會(huì)[9-10]。一些非手術(shù)方法，如化學(xué)療法、放射療法和分子靶向療法等可為GBC患者提供姑息性緩解，減少腫瘤生長(zhǎng)和抑制腫瘤轉(zhuǎn)移[11]。然而，在過(guò)去的幾十年中，GBC的總生存期并沒(méi)有得到明顯的改善[12-14]。因此，對(duì)GBC潛在發(fā)病機(jī)制的分子理解以及研究顯得尤為重要，同時(shí)發(fā)現(xiàn)特異的標(biāo)記物用于GBC的早期診斷和靶向治療也面臨著巨大的挑戰(zhàn)。

隨著高通量測(cè)序技術(shù)的發(fā)展，基因組與轉(zhuǎn)錄組層面的數(shù)據(jù)為研究GBC疾病的分子生物學(xué)機(jī)制提供了契機(jī)。目前已有研究[15-16]表明基因組的不穩(wěn)定性和變異性可能影響GBC疾病的發(fā)生發(fā)展。通過(guò)二代測(cè)序技術(shù)分析已鑒定出與GBC疾病發(fā)生有關(guān)的許多基因突變，包括TP53和ERBB3[17-18]。本研究通過(guò)轉(zhuǎn)錄組層面的GBC數(shù)據(jù)來(lái)篩選差異表達(dá)基因，旨在發(fā)現(xiàn)GBC中異常表達(dá)的基因與相關(guān)生物學(xué)過(guò)程。此外對(duì)出現(xiàn)的關(guān)鍵蛋白調(diào)控基因進(jìn)行評(píng)估，為GBC患者的治療提供新思路，有助于對(duì)GBC疾病進(jìn)行更深入認(rèn)識(shí)。

1 資料與方法

1.1 材料

為研究GBC患者差異基因的表達(dá)特征，本文在美國(guó)國(guó)家生物技術(shù)信息中心（National Center for Biotechnology Information，NCBI）的基因表達(dá)匯編（Gene Expression Omnibus，GEO）數(shù)據(jù)庫(kù)中下載人類(lèi)GBC相關(guān)的轉(zhuǎn)錄組公共數(shù)據(jù)，數(shù)據(jù)集編號(hào)為GSE100363和GSE139682，數(shù)據(jù)信息和樣本分類(lèi)見(jiàn)表1。

表1 GBC 相關(guān)轉(zhuǎn)錄組數(shù)據(jù)信息和樣本分類(lèi)情況Table 1 Data information of GBC-related transcriptome and sample classification

1.2 方法

1.2.1 轉(zhuǎn)錄組數(shù)據(jù)預(yù)處理GSE100363 數(shù)據(jù)集的預(yù)處理過(guò)程：⑴去除每段測(cè)序讀段的接頭序列并過(guò)濾低質(zhì)量序列；⑵使用bowtie2 軟件[19]將序列與參考基因組hg19 進(jìn)行比對(duì)；⑶在RSEM 軟件[20]中計(jì)算基因表達(dá)的FPKM（fragments per kilobase million mapped reads）值。GSE139682 數(shù)據(jù)集的預(yù)處理過(guò)程：⑴對(duì)序列進(jìn)行過(guò)濾，除去堿基質(zhì)量低于13 的20%的讀段；⑵使用HISAT2.1.0 軟件[21]將序列比對(duì)到參考基因組hg38；⑶根據(jù)Chepelev等[22]報(bào)道的方案計(jì)算基因表達(dá)的RPKM（reads per kilobase per million mapped reads）值。

1.2.2 篩選轉(zhuǎn)錄組數(shù)據(jù)差異表達(dá)基因在GBC 和正常樣本中，同時(shí)使用雙總體Student t 檢驗(yàn)和差異倍數(shù)FC（fold change）法來(lái)識(shí)別差異表達(dá)基因，本文定義符合閾值為P＜0.05 且|log2FC| ≥1的基因?yàn)椴町惐磉_(dá)基因（differentially expressed genes，DEGs）。

1.2.3 差異表達(dá)基因的GO 注釋為深入了解GBC 中差異表達(dá)基因的相關(guān)生物學(xué)過(guò)程，本文使用R軟件中clusterProfiler、org.Hs.eg.db、DOSE 等程序包對(duì)差異表達(dá)基因進(jìn)行GO（gene ontology）功能注釋[23]，包括基因參與的生物學(xué)過(guò)程（biological process，BP）、基因所處的細(xì)胞組分（cellular component，CC）和基因執(zhí)行的分子功能（molecular function，MF）。在此采用多重檢驗(yàn)校正的方法Benjamini-Hochberg 對(duì)P 值進(jìn)行校正，設(shè)置滿足校正P＜0.05 的GO 術(shù)語(yǔ)具有顯著性。

1.2.4 差異表達(dá)基因的蛋白質(zhì)互作網(wǎng)絡(luò)構(gòu)建與模塊挖掘隨后采用STRING[24]蛋白質(zhì)互作數(shù)據(jù)庫(kù)中的蛋白質(zhì)互作信息對(duì)GBC 差異表達(dá)基因構(gòu)建蛋白質(zhì)互作（protein-protein interaction，PPI）網(wǎng)絡(luò)，并用Cytoscape 軟件對(duì)蛋白質(zhì)互作網(wǎng)絡(luò)進(jìn)行可視化，設(shè)置最低互作分?jǐn)?shù)為0.4。同時(shí)使用Cytoscape 中的圖論聚類(lèi)算法MCODE（molecular complex detection）在構(gòu)建的蛋白互作網(wǎng)絡(luò)中進(jìn)行模塊挖掘，設(shè)置度的閾值為2，節(jié)點(diǎn)閾值為0.2。

1.2.5 關(guān)鍵蛋白調(diào)控基因的表達(dá)分析和預(yù)測(cè)能力評(píng)估采用Graphpad Prism 8.3.0 軟件對(duì)最終保留的關(guān)鍵蛋白調(diào)控基因表達(dá)量進(jìn)行兩獨(dú)立樣本t 檢驗(yàn)，P＜0.05 為差異具有統(tǒng)計(jì)學(xué)意義。并且對(duì)關(guān)鍵蛋白調(diào)控基因繪制受試者工作特征曲線（receiver operating characteristic curve,ROC）。ROC 曲線下面積（area under curve，AUC）越接近1.0，其診斷效果越好。

2 結(jié)果

2.1 GBC 差異表達(dá)基因分析結(jié)果

在兩個(gè)轉(zhuǎn)錄組數(shù)據(jù)中，同時(shí)采用t檢驗(yàn)和FC的差異基因分析方法，設(shè)置P＜0.05且|log2FC|≥1的閾值條件。GSE100363數(shù)據(jù)集篩選出644個(gè)（2.54%）差異表達(dá)基因，其中上調(diào)基因282 個(gè)（43.79%），下調(diào)基因362 個(gè)（56.21%）。GSE139682數(shù)據(jù)集篩選出2206個(gè)（5.28%）差異表達(dá)基因，上調(diào)基因829個(gè)（37.58%），下調(diào)基因1377個(gè)（62.42%）。兩個(gè)數(shù)據(jù)集中可重復(fù)差異基因144 個(gè)（圖1 A）。共同上調(diào)基因20 個(gè)（圖1B），共同下調(diào)基因120個(gè)（圖1C），表2列舉了在各自數(shù)據(jù)集中差異表達(dá)最大的5個(gè)共同上調(diào)基因和5個(gè)共同下調(diào)基因。由于這140個(gè)共同差異表達(dá)基因可能在GBC發(fā)生發(fā)展中發(fā)揮一定作用，故隨后的分析都基于此。

圖1 GBC 表達(dá)基因的Venn 圖 A：兩個(gè)數(shù)據(jù)集的共同差異表達(dá)基因Venn 圖；B：兩個(gè)數(shù)據(jù)集共同差異表達(dá)的上調(diào)基因Venn 圖；C：兩個(gè)數(shù)據(jù)集共同差異表達(dá)的下調(diào)基因Venn 圖Figure1 The Venn diagram of differentially expressed genes in GBC A:The Venn diagram of differentially expressed gene in GBC in both data sets;B:The Venn diagram of up-regulated genes in both data sets;C:The Venn diagram of down-regulated genes in both data sets

表2 GBC 組織樣本在各自數(shù)據(jù)集中差異表達(dá)最大的5 個(gè)共同上調(diào)和5 個(gè)共同下調(diào)基因Table 2 The 5 common up-regulated and 5 common down-regulated genes with the largest differential expression in GBC tissue samples in each dataset

2.2 差異表達(dá)基因的功能富集

為研究差異表達(dá)的140個(gè)基因的相關(guān)生物學(xué)功能，對(duì)這些基因進(jìn)行GO注釋?zhuān)謩e富集到2249條生物學(xué)過(guò)程術(shù)語(yǔ)、247 條細(xì)胞組分術(shù)語(yǔ)和295 條分子功能術(shù)語(yǔ)。其中符合閾值校正P＜0.05 的生物學(xué)過(guò)程術(shù)語(yǔ)454 條，從最顯著的5 條，可以看出這140個(gè)表達(dá)異常的基因主要影響前腦發(fā)育、端腦發(fā)育和神經(jīng)發(fā)生的正調(diào)控等生物學(xué)過(guò)程，體現(xiàn)了GBC相關(guān)的通路出現(xiàn)了不同程度干擾。在這5條通路中，CDON（cell adhesion associated,oncogene regulatged）、FEZ1（fasciculation and elongation protein zata 1）、GLI3（GLI family zinc finger 3）、SLIT2（slit homolog 2）、SOCS2（suppressor of cytokine signaling 2）基因出現(xiàn)頻率最高，表明這5個(gè)基因主要參與了GBC的生物學(xué)發(fā)生過(guò)程（圖2）。符合閾值的2條細(xì)胞組分術(shù)語(yǔ)，其中9個(gè)基因參與突觸后膜組成，4個(gè)基因參與橫小管組成。另外在閾值的設(shè)定下沒(méi)有符合的分子功能術(shù)語(yǔ)（圖3）。

2.3 差異表達(dá)基因PPI 網(wǎng)絡(luò)分析及模塊挖掘

對(duì)這140個(gè)差異表達(dá)基因在STRING數(shù)據(jù)庫(kù)中篩選互作值≥0.4的互作對(duì)，之后用Cytoscape軟件對(duì)蛋白互作網(wǎng)絡(luò)進(jìn)行可視化（圖4）；用Cytoscape軟件中的MCODE插件在PPI網(wǎng)絡(luò)中進(jìn)行模塊挖掘（表3），PPI網(wǎng)絡(luò)圖中不同的顏色對(duì)應(yīng)表3中不同的模塊。結(jié)果表明SFRP1基因與GLI3、CCND2互作形成模塊1；SEMA3D基因與SLIT2、NTN1互作形成模塊2；ANK2基因與GNAO1、PTPRD互作形成模塊3。在不同基因表達(dá)的蛋白互作下，3個(gè)種子基因?qū)τ贕BC的生物學(xué)過(guò)程可能顯得更為重要。

圖2 差異表達(dá)基因的GO 生物學(xué)過(guò)程注釋?zhuān)v坐標(biāo)表示富集到的不同GO 生物學(xué)過(guò)程，橫坐標(biāo)表示注釋到某條術(shù)語(yǔ)的基因數(shù)目；不同色階表示為不同校正P 值，顏色越紅表示校正P 值越小，其差異越顯著，越藍(lán)表示校正P 值越大越不顯著）Figure 2 Annotation of GO biological process of differentially expressed genes (the vertical axis showing enrichment of different biological processes,the horizontal axis showing the number of genes annotated to a certain term;different color gradations standing for different adjusted P values,namely the darker the red color,the smaller the adjusted P value and the more significant the difference,and the darker the blue color,the bigger the adjusted P value and the less significant the difference )

圖3 差異表達(dá)基因的GO 細(xì)胞組分注釋?zhuān)v坐標(biāo)表示富集到的不同GO 細(xì)胞組分，橫坐標(biāo)表示注釋到某條術(shù)語(yǔ)的基因數(shù)目；圖中展示的兩條術(shù)語(yǔ)校正P 值在0.036附近）Figure 3 Annotation of GO cellular component of differentially expressed genes (the vertical axis showing enrichment of different cellular components,the horizontal axis showing the number of genes annotated to a certain term;the adjusted P values of the two terms displayed in above picture are around 0.036)

圖4 差異表達(dá)基因的PPI 網(wǎng)絡(luò)圖Figure 4 PPI network diagram of differentially expressed genes

表3 PPI 網(wǎng)絡(luò)中對(duì)應(yīng)的模塊信息Table 3 The corresponding module information in PPI network

2.4 關(guān)鍵蛋白調(diào)控基因表達(dá)分析情況

將GBC組織和正常膽囊組織相比，采用兩獨(dú)立樣本t檢驗(yàn)，發(fā)現(xiàn)GBC組織的關(guān)鍵蛋白調(diào)控基因SFRP1、SEMA3D、ANK2的表達(dá)均明顯下降，其差異均具有統(tǒng)計(jì)學(xué)意義（均P＜0.05）（圖5）。這一結(jié)果表明調(diào)控這3種基因的高表達(dá)可能降低GBC的發(fā)生。SFRP1和SEMA3D基因的差異情況略?xún)?yōu)于ANK2基因。

圖5 關(guān)鍵蛋白調(diào)控基因表達(dá)情況分析Figure 5 Analysis of expressions of the key protein regulating genes

2.5 評(píng)估關(guān)鍵蛋白調(diào)控基因?qū)BC 的預(yù)測(cè)能力

通過(guò)繪制ROC曲線來(lái)評(píng)估關(guān)鍵蛋白調(diào)控基因診斷和預(yù)測(cè)GBC的效果。結(jié)果發(fā)現(xiàn)SFRP1和ANK2基因?qū)BC具有一定的預(yù)測(cè)能力，其AUC值均＞0.9。SFRP1基因的預(yù)測(cè)效果略?xún)?yōu)。SFRP1、ANK23、SEMA3D聯(lián)合應(yīng)用的AUC值最大，為0.9643（圖6）。

圖6 關(guān)鍵蛋白調(diào)控基因診斷能力分析Figure 6 Analysis of diagnostic ability of the key protein regulatory genes

3 討論

GBC是最常見(jiàn)的膽道惡性腫瘤，由于缺乏特異性癥狀和可靠敏感的疾病標(biāo)志，大多數(shù)GBC患者被診斷即為晚期，因此迫切面臨對(duì)分子機(jī)制的確切理解與研究。本研究旨在利用轉(zhuǎn)錄組層面的高通量數(shù)據(jù)，使用T檢驗(yàn)和倍數(shù)變換法篩選差異表達(dá)基因，結(jié)果發(fā)現(xiàn)在兩個(gè)數(shù)據(jù)集中符合閾值為P＜0.05且|log2FC|≥1的共同出現(xiàn)的差異基因140個(gè)（上調(diào)20個(gè)，下調(diào)120個(gè)）。通過(guò)GO中3種類(lèi)型的功能注釋?zhuān)陂撝禐樾Ｕ齈＜0.05的條件下，發(fā)現(xiàn)差異基因主要與前腦發(fā)育、端腦發(fā)育和神經(jīng)發(fā)生的正調(diào)控等生物學(xué)過(guò)程有關(guān)，與突觸后膜和橫小管這兩種細(xì)胞組分密切相關(guān)，而沒(méi)有發(fā)現(xiàn)與分子功能顯著相關(guān)的注釋信息。隨即對(duì)差異表達(dá)基因構(gòu)建蛋白互作網(wǎng)絡(luò)以發(fā)現(xiàn)存在明顯相互作用的基因。值得注意的是3個(gè)種子節(jié)點(diǎn)的基因SFRP1、SEMA3D、ANK2，它們與正常膽囊組織相比，在GBC組織中表達(dá)量均為明顯下降趨勢(shì)，可能在GBC發(fā)生中起著尤為重要的作用。其中，SFRP1基因相比于其他兩個(gè)基因的差異表達(dá)量最顯著且診斷效果也最明顯。

既往研究證實(shí)由SFRP1基因所編碼的SFRP1蛋白作為Wnt信號(hào)蛋白的拮抗劑，在Wnt/catenin信號(hào)通路的調(diào)控中發(fā)揮重要作用。此外，SFRP1基因在膽道系統(tǒng)中也被報(bào)道是高度相關(guān)的甲基化基因，基于其甲基化的分子標(biāo)記物已成功地用作常規(guī)、預(yù)后及預(yù)測(cè)指標(biāo)，以更好地管理各種癌癥患者，一些研究還描述了SFRP1作為癌癥的預(yù)后和預(yù)測(cè)生物標(biāo)志物的價(jià)值[25-26]。在人類(lèi)中，高度相關(guān)甲基化基因的生物標(biāo)記，允許從相應(yīng)的正常組織中區(qū)分出腫瘤組織，這些標(biāo)志物還與獨(dú)特的癌癥表型有關(guān)，通過(guò)進(jìn)一步優(yōu)化分析SFRP1基因可能會(huì)增加對(duì)GBC早期的發(fā)現(xiàn)與診斷，并在人體液體中表現(xiàn)相對(duì)可靠的監(jiān)測(cè)潛力[27-28]。

同時(shí)，針對(duì)目前發(fā)現(xiàn)的在差異最顯著的5 條GO生物學(xué)過(guò)程中出現(xiàn)的頻率最多的基因也有研究證實(shí)它們的發(fā)生與GBC有關(guān)。有研究表明，GBC的許多抑癌基因中已經(jīng)描述了雜合子丟失，包括染色體1p34-36（p73），5q21（APC），8p21-23（PRLTS和FEZ1）和17p13（p53基因）等[29]。最新報(bào)道[30]發(fā)現(xiàn)，參與GBC發(fā)生發(fā)展的長(zhǎng)鏈非編碼RNA（long-chain non-coding RNA，lncRNA）SNHG6的沉默和miR-26b-5p的上調(diào)可促進(jìn)細(xì)胞凋亡、抑制細(xì)胞生長(zhǎng)和上皮間質(zhì)轉(zhuǎn)化，并促進(jìn)GLI3和E-鈣黏著蛋白表達(dá)的上調(diào)，進(jìn)而影響GBC的生物學(xué)過(guò)程。另外轉(zhuǎn)錄因子GLI3是Hedgehog（Hh/HH）信號(hào)通路的成員，可以全長(zhǎng)（Gli3-FL/GLI3-FL）或阻遏物（Gli3-R/GLI3-R）的形式存在。GLI3調(diào)節(jié)各種對(duì)于癌細(xì)胞生長(zhǎng)和發(fā)展至關(guān)重要的生物過(guò)程，錨定非依賴(lài)性生長(zhǎng)、增殖和遷移的腫瘤細(xì)胞系[31]。另外SLIT2基因的甲基化在多種腫瘤中也有不同程度的研究[32-33]。

此外，本文也存在如下不足之處。首先，由于是基于生物信息學(xué)方法進(jìn)行的研究，所以部分結(jié)果還需要進(jìn)一步的功能實(shí)驗(yàn)進(jìn)行驗(yàn)證，因而可以為之后的臨床應(yīng)用提供更可靠的依據(jù)。 其次，由于GBC疾病本身發(fā)生機(jī)制的異質(zhì)性和研究樣本數(shù)量的有限性，快速發(fā)展的高通量單細(xì)胞層面測(cè)序技術(shù)可能會(huì)解決這一問(wèn)題。最后，因?yàn)镚BC疾病的發(fā)生是受多方面因素影響的，只從轉(zhuǎn)錄組層面看存在一定的片面性，本文尚未結(jié)合基因組、蛋白質(zhì)組、代謝組等多組學(xué)數(shù)據(jù)信息，多組學(xué)數(shù)據(jù)之間的相關(guān)性有待研究。同時(shí)環(huán)境的因素也可能是GBC發(fā)生與進(jìn)展的重要原因。因此，本研究結(jié)果可能適用于部分GBC患者，并非適用于所有GBC患者。

綜上，本研究通過(guò)分析兩個(gè)轉(zhuǎn)錄組層面數(shù)據(jù)集取交集的相關(guān)差異基因的GO生物學(xué)過(guò)程和細(xì)胞組分，關(guān)鍵蛋白調(diào)控基因SFRP1可能在GBC中發(fā)揮不可忽視的作用。研究結(jié)果具有一定的可靠性，可以為更深刻的認(rèn)識(shí)GBC分子機(jī)制提供新的見(jiàn)解與方向。