基于RNA-seq分析人肝癌耐藥細胞BEL-7402/FU的差異表達基因

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

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

1 材料與方法

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

1.2方法

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

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

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

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

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

2 結果

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

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

圖1 樣品間相關性分析熱圖X、Y軸均代表每個樣品；顏色代表相關性系數

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

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

圖2 DEGs的MA-plot分布圖

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

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

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

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

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

3 討論

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

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

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

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

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

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

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

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

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