基于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的筛选结果进行验证，以期进一步探讨肝癌耐药的机制。