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

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

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

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

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

1 資料與方法

1.1 材料

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

表1 GBC 相關(guān)轉(zhuǎn)錄組數(shù)據(jù)信息和樣本分類情況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ù)處理過程：⑴去除每段測序讀段的接頭序列并過濾低質(zhì)量序列；⑵使用bowtie2 軟件[19]將序列與參考基因組hg19 進行比對；⑶在RSEM 軟件[20]中計算基因表達的FPKM（fragments per kilobase million mapped reads）值。GSE139682 數(shù)據(jù)集的預(yù)處理過程：⑴對序列進行過濾，除去堿基質(zhì)量低于13 的20%的讀段；⑵使用HISAT2.1.0 軟件[21]將序列比對到參考基因組hg38；⑶根據(jù)Chepelev等[22]報道的方案計算基因表達的RPKM（reads per kilobase per million mapped reads）值。

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

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

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

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

2 結(jié)果

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

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

圖1 GBC 表達基因的Venn 圖 A：兩個數(shù)據(jù)集的共同差異表達基因Venn 圖；B：兩個數(shù)據(jù)集共同差異表達的上調(diào)基因Venn 圖；C：兩個數(shù)據(jù)集共同差異表達的下調(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ù)集中差異表達最大的5 個共同上調(diào)和5 個共同下調(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 差異表達基因的功能富集

為研究差異表達的140個基因的相關(guān)生物學(xué)功能，對這些基因進行GO注釋，分別富集到2249條生物學(xué)過程術(shù)語、247 條細胞組分術(shù)語和295 條分子功能術(shù)語。其中符合閾值校正P＜0.05 的生物學(xué)過程術(shù)語454 條，從最顯著的5 條，可以看出這140個表達異常的基因主要影響前腦發(fā)育、端腦發(fā)育和神經(jīng)發(fā)生的正調(diào)控等生物學(xué)過程，體現(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個基因主要參與了GBC的生物學(xué)發(fā)生過程（圖2）。符合閾值的2條細胞組分術(shù)語，其中9個基因參與突觸后膜組成，4個基因參與橫小管組成。另外在閾值的設(shè)定下沒有符合的分子功能術(shù)語（圖3）。

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

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

圖2 差異表達基因的GO 生物學(xué)過程注釋（縱坐標(biāo)表示富集到的不同GO 生物學(xué)過程，橫坐標(biāo)表示注釋到某條術(shù)語的基因數(shù)目；不同色階表示為不同校正P 值，顏色越紅表示校正P 值越小，其差異越顯著，越藍表示校正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 差異表達基因的GO 細胞組分注釋（縱坐標(biāo)表示富集到的不同GO 細胞組分，橫坐標(biāo)表示注釋到某條術(shù)語的基因數(shù)目；圖中展示的兩條術(shù)語校正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 差異表達基因的PPI 網(wǎng)絡(luò)圖Figure 4 PPI network diagram of differentially expressed genes

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

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

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

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

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

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

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

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

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

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