周謝海，張 婷，謝 剛，崔志華，許杜娟,,3，夏 泉

原發(fā)性肝癌主要是肝細胞癌 (hepatocellular carcinoma, HCC)，簡稱肝癌，是常見的惡性腫瘤之一，臨床發(fā)病率和死亡率均較高?；熓侵型砥诟伟┲匾闹委熓侄?，然而，患者往往因為多藥耐藥的產(chǎn)生而導致化療失敗。肝癌耐藥是一個多進程多種因素相互作用的結(jié)果，其中包括多藥耐藥蛋白的高表達[1]，凋亡抑制途徑異常[2]，以及細胞內(nèi)在信號通路的調(diào)節(jié)[3]等。但迄今為止仍尚未完全闡明肝癌耐藥的機制。RNA-seq，即對經(jīng)RNA反轉(zhuǎn)錄以及PCR擴增得到的cDNA進行高通量測序分析，具有數(shù)據(jù)量大、定量準確、可重復性高、分析可靠等特點[4]，近些年來已廣泛運用于生命科學領(lǐng)域。該研究基于RNA-seq技術(shù)，篩選出肝癌耐藥細胞與親本細胞的差異表達基因(differentially expressed genes, DEGs)，以期探討與肝癌耐藥相關(guān)的基因及信號傳導通路。

1 材料與方法

1.1材料人肝癌細胞株BEL-7402購自中科院上海細胞庫，人肝癌耐藥細胞株BEL-7402/FU購自南京凱基生物有限公司，由本實驗室液氮凍存保種；5-氟尿嘧啶(5-FU)購自美國Sigma公司(貨號：V900394)；TRIzol Reagent 購自美國 Invitrogen公司。

1.2方法

1.2.1細胞培養(yǎng) 使用含10% 胎牛血清、雙抗(青霉素100 IU/ml、鏈霉素100 mg/L)的RPMI-1640培養(yǎng)基進行培養(yǎng)。置于 37 ℃、5% CO2培養(yǎng)箱中貼壁培養(yǎng)，細胞在對數(shù)生長期時用胰酶消化傳代。肝癌耐藥細胞株BEL-7402/FU每48 h使用含有20 mg/L 5-FU的培養(yǎng)基刺激，以維持細胞的耐藥性。為了避免5-FU的影響，在樣本處理前將BEL-7402/FU細胞在正常培養(yǎng)基中培養(yǎng)48 h后再進行取樣。

1.2.2Total RNA提取與鑒定 使用TRIzol試劑分別提取兩組細胞的總RNA，每組細胞重復3個生物學樣品，使用Agilent 2100 Bioanalyzer 檢測總RNA的濃度、RIN值、28S/18S和片斷大小。

1.2.3轉(zhuǎn)錄組測序及DEGs篩選 將提取的 RNA 樣品送至華大基因進行轉(zhuǎn)錄組文庫構(gòu)建，并通過BGISEQ-500平臺進行測序。運用DEseq2 和 PossionDis 算法進行DEGs檢測，DEseq2方法根據(jù)差異倍數(shù)在兩倍及以上且校正后的P值≤0.05來篩選DEGs。PossionDis 方法計算得到差異基因后，對差異檢驗的P值作多重假設(shè)檢驗校正，并通過控制錯誤發(fā)現(xiàn)率(false discovery rate, FDR)決定P值的域值。本研究中，DEGs默認定義為FDR≤0.001且倍數(shù)差異在2倍以上的基因。

1.2.4DEGs分析 根據(jù)GO富集、KEGG注釋結(jié)果以及官方分類，將DEGs進行功能分類、生物通路等分類，同時使用R軟件中的Phyper函數(shù)進行富集分析。使用DIAMOND軟件將DEGs比對至STRING數(shù)據(jù)庫，利用與已知蛋白的同源性獲得DEGs編碼蛋白間的互作關(guān)系，即蛋白互作網(wǎng)絡(protein protein interaction network，PPI)分析。取得分最高的前1 000個關(guān)系畫示例圖。

1.3統(tǒng)計學處理本項目使用R軟件里的Cor函數(shù)對每兩個樣品之間進行Pearson相關(guān)性分析，P<0.05為差異有統(tǒng)計學意義。樣品基因表達量服從泊松分布，本項目P值采用Wald test方法計算并通過多重假設(shè)檢驗校正得出，P<0.05為差異具有統(tǒng)計學意義。

2 結(jié)果

2.1轉(zhuǎn)錄組測序結(jié)果分析本項目使用BGISEQ-500平臺一共測了兩組細胞的樣品，每組重復3個樣品，樣品比對基因組的平均比對率為81.25%，比對基因集的平均比對率為69.47%；共檢測到基因數(shù)為17 738個，其中已知的基因為16 897個，預測的新基因為928 個；共檢測出12 995個新轉(zhuǎn)錄本，其中11 390個屬于已知蛋白編碼基因的新的可變剪接亞型，928個屬于新的蛋白編碼基因的轉(zhuǎn)錄本，剩下的677個屬于長鏈非編碼RNA。

為了反映樣本間基因表達的相關(guān)性，本研究計算了每兩個樣品之間所有基因表達量的Pearson相關(guān)系數(shù)，并將這些系數(shù)以熱圖的形式反映出來，圖中顏色越深代表相關(guān)性越高，顏色越淺代表相關(guān)性越低，見圖1。

圖1 樣品間相關(guān)性分析熱圖X、Y軸均代表每個樣品；顏色代表相關(guān)性系數(shù)

2.2DEGs表達水平分析根據(jù)各個樣品基因表達水平不同，篩選出樣品組之間的DEGs。統(tǒng)計結(jié)果顯示，共篩選出樣品組之間的DEGs數(shù)為6 212個，其中上調(diào)數(shù)為2 899個，下調(diào)數(shù)為3 313個，使用MA plot圖展示DEGs的分布，見圖2。為了篩選表達差異最為顯著的基因，項目對表達量上、下調(diào)10倍以上的基因進行分析。結(jié)果顯示，上、下調(diào)10倍以上的DEGs共有83個，其中72個屬于已知基因，11個log2FoldChange指經(jīng)過log2轉(zhuǎn)換后的BEL7402和BEL7402/FU細胞樣品組之間的差異表達倍數(shù)；NA指檢測到的新基因為預測的新基因。本次選擇上、下調(diào)的前10名分別列表展示，見表1。

表1 BEL7402和BEL7402/FU細胞的部分DEGs表達

圖2 DEGs的MA-plot分布圖

2.3DEGs的GO富集分析根據(jù)DEGs檢測的結(jié)果，本項目利用GO富集分析對轉(zhuǎn)錄組數(shù)據(jù)進行功能注釋。GO分析結(jié)果表明，分別有4 138條注釋到生物學進程(biological process, P)，有4 743條注釋到細胞組分(cellular component, C)，有4 465條基因注釋到分子功能(molecular function, F)。之后又將其細分為62個小類。在P分類中, 70.5%的DEGs可被富集到單生物體過程類別(single-oranism process)條目中，富集的基因數(shù)最多，其次為單生物細胞過程類別(single-oranism cellular process)，約占61.3%。在C分類中，分別有約82.5%的DEGs可被富集到細胞類別(cell)和細胞部分類別(cell part)條目，約75.5%的DEGs可富集到細胞內(nèi)類別條目(intracellular)。在F分類中，結(jié)合分子功能類別(binding)條目所占比例最多(84.3%)，其次為蛋白結(jié)合類別(protein binding)條目，約占42.9%。

2.4DEGs的KEGG生物通路富集分析KEGG 分析表明，共有5 718條DEGs涉及共312個通路(pathways)，主要涉及的通路見圖3，其中代謝通路 (metabolic pathways) 富集的基因數(shù)目最多，約占DEGs總數(shù)的10.84%；其次為腫瘤中的蛋白聚糖通路(proteoglycans in cancer)以及人類嗜T細胞病毒感染通路(HTLV-I infection)。

與此同時，通過對上、下調(diào)10倍以上的基因進一步進行富集分析發(fā)現(xiàn)，篩選出的83條DEGs中，共有57條可被富集到KEGG通路中，其中有11條基因涉及代謝通路，約占19.30%；5條涉及絲裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)信號通路，約占8.77%；其余的主要涉及RNA降解(RNA degradation)，磷脂酰肌醇-3-激酶/蛋白激酶B信號通路(PI3K-Akt signaling pathway)等。接下來，對代謝通路中上、下調(diào)10倍以上的基因進行匯總，上調(diào)的基因有UGT1A1、GALNT6、CDS1、MOCS1、HSD3B1、UGT8、DGKG、B3GALNT1、TFF1，下調(diào)的基因為GALNT3、LINC00689。

圖3 差異表達基因的Pathway顯著性富集分析

2.5DEGs編碼蛋白質(zhì)之間的PPI分析蛋白之間通常通過相互作用結(jié)合成復合物之后行使相應的功能，具有相互作用的DEGs通常具有相似的功能。根據(jù)STRING蛋白互作數(shù)據(jù)庫，對每組DEGs進行PPI分析。采用可信度最高的前1000個互作關(guān)系做圖，結(jié)果見圖4。如圖所示，與各蛋白聯(lián)系最密切的為TP53和表皮生長因子受體(epidermal growth factor receptor, EGFR)，蛋白質(zhì)相互作用主要集中在CDKN1A、BCL2、ABALON、BCL2L1、MAPK8、FYN、ESR1、STAT5B及STAT5A等基因之間。以上基因平均涉及到23條信號通路，基因的表達情況見表2。

3 討論

本項目運用RNA-seq技術(shù)，從整體水平研究肝癌親本細胞BEL-7402與耐藥細胞BEL-7402/FU的基因轉(zhuǎn)錄水平，篩選出DEGs并通過生物信息學分析尋找與肝癌耐藥關(guān)系密切的基因和通路。研究顯示，BEL-7402和BEL-7402/FU細胞之間存在著多個DEGs的上調(diào)和下調(diào)，涉及到多種細胞信號通路和生物學過程。

GO富集分析發(fā)現(xiàn)，兩組細胞在各分類中均富集了大量的DEGs，說明肝癌耐藥的過程不是單一因素引發(fā)的，而是多條途徑相互作用的結(jié)果。例如在P分類中，Single-organism process和Single-organism cellular process顯著富集，它們是最基礎(chǔ)的生命活動類別；在C分類中，Cell占比最多, 說明細胞適應化療藥物的過程是整體的和全局性的；在F分類中，Binding和Protein binding富集最為顯著，說明肝癌耐藥涉及到細胞中各種蛋白和分子的結(jié)合作用。

圖4 差異表達基因的蛋白互作關(guān)系網(wǎng)絡圖圈的大小表示關(guān)系密集程度，顏色深淺表示基因在網(wǎng)絡中的重要程度

表2 PPI分析中主要DEGs的表達

KEGG分析結(jié)果顯示，代謝通路富集的基因數(shù)目最多，并且在上、下調(diào)10倍的DEGs中，富集到代謝通路的基因占比最大，多數(shù)基因涉及到腫瘤能量代謝和糖蛋白的合成。腫瘤在適應化療和放療帶來的損傷時，需要通過適應環(huán)境并開啟防御機制，其中包括藥物的外排、抵抗凋亡、DNA的損傷修復以及信號通路的激活等，這些生物過程都需要足量的能量供給。早在20世紀初，科學家就發(fā)現(xiàn)腫瘤細胞具有高水平的糖酵解作用，并通過代謝重組維持細胞內(nèi)的ATP和NADH水平的正常，以滿足細胞生存和大分子合成的需要[5]。本課題組前期研究[6]發(fā)現(xiàn)，在腫瘤酸性微環(huán)境中，ASIC1a通過調(diào)控Ca2+/PI3K/AKT信號通路誘導人肝癌耐藥。腫瘤酸性微環(huán)境的產(chǎn)生原因之一正是腫瘤缺氧、腫瘤細胞高水平的糖酵解的發(fā)生[7]，此外，Zhou et al[8]研究也發(fā)現(xiàn)，結(jié)直腸癌細胞ATP損耗可以提高化療藥物的敏感性。因此，有理由認為腫瘤代謝通路與肝癌耐藥關(guān)系密切。

在代謝通路中，尿苷二磷酸葡萄糖醛酸轉(zhuǎn)移酶1A1(Uridine diphosphoglucu-ronosyl transferase 1A1, UGTlA1)表達差異顯著。UGT1A1是參與膽紅素代謝的一種酶，也是伊立替康代謝的關(guān)鍵酶，與結(jié)直腸癌耐藥有關(guān)[9]，并且特異性沉默UGT1A1基因可顯著提高伊立替康對結(jié)直腸癌細胞的化療敏感性[10]，但UGT1A1與肝癌耐藥的研究卻鮮有報道。腫瘤耐藥的機制之一為腫瘤細胞對化療藥物的代謝增強，因此，在BEL-7402/FU細胞中，可能存在由于UGT1A1的過度表達而提高了對5-FU的代謝，從而導致肝癌耐藥的產(chǎn)生。本次研究中還發(fā)現(xiàn)GALNT6、B3GALNT1、GALNT3參與到O-糖基化過程，與細胞膜上的糖蛋白合成和癌癥進展息息相關(guān)[11]。以上DEGs可能在肝癌耐藥細胞的代謝調(diào)節(jié)中起到不可或缺的作用，它們均有可能是調(diào)控肝癌耐藥的潛在基因。

PPI分析顯示，蛋白質(zhì)相互作用主要集中在CDKN1A、BCL2、ABALON等基因之間，與各蛋白聯(lián)系最密切的為TP53和EGFR基因。以上基因平均可富集到23條信號通路，它們多數(shù)參與了細胞周期、凋亡以及增殖等信號通路的調(diào)節(jié)，在腫瘤耐藥的調(diào)節(jié)中起到了關(guān)鍵作用。TP53是重要的抑癌基因之一，主要通過轉(zhuǎn)錄合成細胞周期相關(guān)蛋白進而調(diào)控細胞增殖和凋亡，TP53的突變與化療耐藥有關(guān)[12]。本研究發(fā)現(xiàn)TP53表達輕度下調(diào)(-1.56倍)，但其位點是否存在突變?nèi)孕枰M一步的檢測。EGFR是表皮生長因子受體(HER)家族成員之一，在大部分肝癌患者標本中呈高表達[13]，目前認為，EGFR主要通過調(diào)控激活RAS/MAPK信號通路誘導原癌基因的表達促進肝癌的發(fā)展[14]。本次檢測發(fā)現(xiàn)EGFR在肝癌親本細胞中呈現(xiàn)高表達，而在耐藥細胞中表達輕度下調(diào)(-1.57倍)，考慮到EGFR的調(diào)節(jié)需要磷酸化后發(fā)揮功能，因此僅僅只觀察受體表達水平并不能代表功能的變化，仍須進一步檢驗其磷酸化水平。

本項研究共檢測了17 738個基因，測序范圍廣，并且在本次測序中發(fā)現(xiàn)了928個屬于新的蛋白編碼基因的轉(zhuǎn)錄本。這些新基因的發(fā)現(xiàn)，為以往測序的結(jié)果做了補充，在這些基因中有可能存在尚未被發(fā)現(xiàn)的引起肝癌耐藥的潛在位點。

綜上所述，本研究篩選出代謝通路為肝癌耐藥最為密切的通路，其中UGT1A1、GALNT6、B3GALNT1和GALNT3在代謝通路中高表達，可能是與肝癌耐藥最為密切的基因，而TP53和EGFR等基因可能是蛋白調(diào)控網(wǎng)絡中最為關(guān)鍵的基因。這些DEGs的篩選為解決肝癌耐藥問題提供了新思路，課題組后期將針對DEGs的篩選結(jié)果進行驗證，以期進一步探討肝癌耐藥的機制。