基于差异基因网络特性探究乳腺癌他莫昔芬耐药的关键标志物

马 军,王雪庭,陈梦杰,孟 磊,王 乐

(1.陕西科技大学 食品与生物工程学院，陕西 西安 710021；2.广东省科学院 广东省微生物研究所 省部共建华南应用微生物国家重点实验室 广东省菌种保藏与应用重点实验室 广东省微生物应用新技术公共实验室,广东 广州 510070；3.西安交通大学第一附属医院 肿瘤外科，陕西 西安 710061；4.西安交通大学第一附属医院 肿瘤内科，陕西 西安 710061)

0 引言

近些年来，中国的乳腺癌发病率呈持续上升的趋势[1].前期研究表明，大约有70%～80%的乳腺癌患者表达雌激素受体(estrogen receptor，ER)[2].他莫昔芬作为治疗乳腺癌雌激素受体阳性的常用药，虽然能够对雌激素受体阳性乳腺癌患者的死亡率和复发率起到控制和降低的作用，但临床资料显示约有40%左右的雌激素受体阳性的乳腺癌患者对他莫昔芬是原始耐药，还有些患者是获得性耐药[3].伴随乳腺癌临床研究的开展，尽管内分泌治疗药物层出不穷，但他莫昔芬的基础地位仍无法动摇，由于其作用机制，无论患者的月经状态均可应用，价格低廉，在基层医疗单位有广泛的患者群体.因此有关于TAM耐药的研究仍然具有十分重要的临床意义.鉴定可以用于预测TAM 疗效的可靠标记物或作为耐药逆转的新靶点，将有利于优化乳腺癌患者的个体化用药方案，进而提高肿瘤的治疗效果，具有重要的临床应用价值.

本课题将利用网络分析方法，在系统水平探究与乳腺癌他莫昔芬的耐相关的表达差异基因.最后结合文献证据与预后数据对鉴定出的他莫昔芬耐药相关的关键基因进行验证.本研究结果将即为探究乳腺癌他莫昔芬耐药的机制研究提供了新思路，也为他莫昔芬的临床用药提供参考.

1 材料与方法

1.1 公共数据

人类蛋白相互作用(protei-protein interaction,PPI)网络来自STRING数据库9606.protein.links(http://string-db.org/cgi/download.pl?sessionId=fkKHAL4ksNUb).STRING数据库收集了2 000多个物种的蛋白互作信息.综合考虑了蛋白对接、文献、基因共表达等方面的数据对蛋白与蛋白之间的互作关系进行赋值；人类乳腺癌细胞系表达谱来自于CCLE(Cancer Cell Line Encyclopedia).CCLE 数据库存储了不同肿瘤的1000多种细胞系的转录组、基因突变、拷贝数以及24种药物处理细胞后药理数据；从GEO数据库中下载获得他莫昔芬耐药相关得基因芯片数据(GSE67916).这个数据集采用GPL570芯片平台，包含了8例他莫昔芬敏感的乳腺癌细胞(MCF-7/TAM)同时还包含了10例他莫昔芬耐药的乳腺癌细胞(MCF-7/TAMR).

1.2 构建乳腺癌蛋白互作网络

从STRING中下载得到人类蛋白质相互作用网络数据导入R中，为了去除假阳性连接的存在，保留蛋白网络中连接分数大于800的相互连接.利用R的 PharmacoGx包提取CCLE 数据集存储的56种人乳腺癌细胞系的RNAseq转录组数据[4].乳腺癌细胞系中的基因reads数经标准化后，得到基因的FPKM(The number of fragments per kilobase per million)数值，再对FPKM数值进行log2(FPKM+1)转换后，作为基因的表达值.基因表达值满足一下条件将被用于从PPI网络中筛选乳腺癌蛋白互作(B_PPI)网络：即基因的表达值大于1且至少在30种乳腺癌细胞中都表达[5].

1.3 网络拓扑结构分析

通过用igraph R包中集合多种进行生物网络功能群落分析的算法.之前的实验表明Walktrap算法能更好探究蛋白互作网络中的拓扑结构[6].随机游走算法(Walktrap，CW)中提到的短随机漫游有停留在同一个模块的倾向.因此，本实验使用Walktrap算法对B_PPI 网络中的模块聚类分析后得到大小不一的群落，被划归为同一个群落中的基因，被认为参与了相似或相关的生物学功能.节点数为10以上的称为大群落，将小于10同时大于等于2的称为小群落.

1.4 表达差异基因的筛选

通过使用R语言GEOquery软件包中的getGEO函数对GSE67916数据集里芯片的原始数据进行质量控制，对样本的探针值进行log2转换.调用R中的Limma包对TAMR组和TAM组中的基因差异表达基因进行筛选，筛选所采用的标准是log2foldchange大于1.5倍且错误发现率((false discovery rate，FDR)小于等于 0.05.删除未被注释的探针或探针对应多个基因名.最终得到了459个表达差异基因，用于后续实验.

1.5 TAMR核心基因表达与预后分析

将TAMR核心基因输入到PrognoScan(http://dna00.bio.kyutech.ac.jp/PrognoScan/index.html)网站用于收集乳腺癌患者中，核心基因表达与预后的关系.然后从GEO GSE9893数据集中得到乳腺癌患者的临床信息及CDH1和STAT1基因表达谱.GSE9893数据集使用GPL5049平台对147例ER阳性患者，2例ER，PR双阴性患者，6例ER阴性PR阳性患者的芯片数据采集.依据耐药关键基因在患者中的表达量进行递增排序，选取前30%和后30%患者的表达值为基因高、低表达的阈值.使用Cox风险比例回归模型(Cox proportional hazards models)判定基因表达对预后的影响是否具有显著性(即P<0.05)，从而确定预后标记物.最后，用Kaplan-Meier分析以生成相应基因的生存曲线.

1.6 耐药相关群落的检验

Reciprocal best-hit方法被用于评估不同细胞系和大鼠肝脏的群落相似性[7].乳腺癌的蛋白互作网络经聚类后产生了多个大群落.大量的表达差异基因同样也可以被认为是一个大群落.因此，相同的方法计算表达差异基因与每个群落相似性的统计学意义，从而确定与他莫昔芬耐药相关的大群落.利用超几何分布算法分别评估重合基因在表达差异基因和蛋白群落的统计学意义.然后，采用同样的方法进一步计算蛋白群落与预后标记物之间相似性，从而验证耐药相似群落与他莫昔芬的预后密切相关.

1.7 生物学功能分析

从基因富集分析(Gene Set Enrichment Analysis，GESA)数据库中下载功能分析所需的 Gene ontology (GO)和KEGG 数据集.GO注释包含了三种功能信息：基因参与的生物学过程、细胞位置和分子功能、KEGG数据库中信号通路数据集.本文使用R功能包piano计算超几何分布，从而检测群落所含基因和表达差异基因在不同功能数据集中的富集情况，FDR≤0.05的功能集表明所测基因在此生物学功能显著聚集.

2 结果与讨论

2.1 乳腺癌蛋白互作网络

STRING数据库中的蛋白互作网络包含了人的所有表达基因及其互作关系.但是肿瘤组织是具有异质性的，这意味着其肿瘤组织的表达谱与正常组织不同.为了能准确探究乳腺癌组织内的蛋白间的互作关系，有必要依据其表达情况对源蛋白网络中的基因以及连接关系进行筛选.细胞是被广泛用于药效，药理研究的体外模型，其相比于体内肿瘤组织而言，可以有效避免在采样过程混入正常组织造成测序数据分析的干扰[8].CCLE数据库是研究药物基因组的重要数据库，包含了56个不同乳腺癌细胞系.基因在不同乳腺癌细胞系的表达情况也不同，因此需要保留在多数细胞内都表达的基因用于后续实验，从而避免特异性基因的存在.

从CCLE中下载提取得到人乳腺癌细胞的细胞系，取log2(FPKM+1)处理后结果大于1的数值为基因的有效表达，得到13 096个在30种乳腺癌细胞中表达的基因.图1(a)反映的是乳腺癌细胞表达的基因数量，红线表示所有细胞内表达基因的数目.用这些基因对STRING数据库中得到的人类蛋白相互作用网络进行过滤.当连接度大于800时，该连接为可靠连接，因此最终得到了基于乳腺癌细胞表达谱产生的蛋白互作网络包含了7 904个基因，和213 422个连接数(如图1(b)所示).

在网络分析中，度(degree)是刻画节点中心性的最直观的度量指标.一个节点的degree越大，那么也就说明这个节点的degree centrally(度中心性)越高，也就是说明这个节点在该网络中越重要.图1(c)表明原始蛋白互作网络与乳腺癌蛋白互作网络中的基因节点度具有统计学意义上的差别(Mann-Whitney test，P<0.05).结果表明随着网络中基因组成的改变，其网络拓扑结构也随之改变.

(a)56个乳腺癌细胞内表达基因的数目

2.2 他莫昔芬耐药相关群落及其功能分析

前期研究较多关注的是基因的差异倍数，或差异方向(如上调或下调).考虑了基因本身，却忽略了基因与基因的内在联系，缺乏对差异基因在系统水平上的分析.机体内的各个基因都是参与特性的生物学功能，之间存在着内在联系.而对蛋白互作网络进行聚类分析后，得到的群落可以将参与相关功能的基因划归为一组，提供了一种系统探究基因间互作关系的新视角.然后这些群落反映的是所有细胞功能.因此需再结合表达差异基因在群落存在是否具有统计学意义，确定与耐药相关的群落.最后通过计算耐药相关群落所含的预后标记物是否具有统计学意义，进而证明群落与他莫昔芬耐药相关.

采用CW算法对乳腺癌的PPI网络进行聚类分析，将网络中的基因进行分组，有利于探究参与特定生物学功能的基因组成.该算法依据网络的拓扑结构特性，将7 904个基因划归成大小不一的群落，其中包含了77个大群落和250个小群落(如图2(b)所示).最大的群落CW5中含有1 948个基因，占了将近25%的网络基因.接下来是含有702个基因的CW4、CW33、CW21、CW13分别含有大约400个基因.还有多个群落分别含有200～300个基因，例如由299个基因组成的CW30和包含了192 个基因的CW35.图2(a)进一步展示了群落所含基因数目的分布情况，大多数的群落是由小于100个基因组成的，整个乳腺癌蛋白互作网络中的群落大小分布呈现拖尾形.

(a)乳腺癌蛋白网络群落大小分布

利用GEO数据库中GSE67916的数据集进行分析，得到了459个他莫昔芬耐药与敏感的表达差异基因.接下来对表达差异基因及群落基因做富集分析，筛选得到了CW4、CW5、CW13、CW21和CW33群落(共3 863个基因)，结果表明这些群落与他莫昔芬耐药相关.再将他莫昔芬治疗的预后标记物与耐药相关群落基因进行Fisher′s exact 检测，结果具有显著性(P<0.05)的群落被进一步验证为与耐药相关.

五个耐药相关群落涉及不同的生物学过程(图3)，其中CW5群落与表达差异基因、预后标记物的相似性最高.在图4(a)中，CW5中基因参与的功能有：跨膜受体蛋白酪氨酸激酶超分子纤维组织、凋亡信号的传导途径、生物合成过程的负调节等.由图4(b)可知，通过对CW5基因进行KEGG信号通路分析，发现在丝裂原活化蛋白激酶 (mitogen-activated protein kinase,MAPK) 信号通路具有较高的显著性.MAPK信号通路不仅能促进细胞的增殖、分裂和分化，还在乳腺癌内分泌治疗药物耐药机制中扮演着重要角色.有研究结果显示，阻断MAPK信号通路可以抑制乳腺癌细胞的增殖和迁移，而激活该信号通路会增强ER转录效率并导致[9].

图3 TAM耐药相关群落功能富集分析

(a)GO 生物过程分析

2.3 TAM耐药核心基因

基因节点度是生物网络特征指数之一，其大小反应了基因在网络中所起的作用大小.利用差异基因在耐药相关群落中的节点度，确定他莫昔芬耐药核心基因.基因的节点度越高指其在耐药模块中与更多基因存在着相互影响，表明这些基因表达量的改变对他莫昔芬的药效产生重要影响.CW5包含了1 948个基因，其中有74个基因为表达差异基因.依据表达差异基因在CW5群落中的节点度，选取前5个具有高节点度的基因STAT3、STAT1、CDH1、EGR1和PRKCA作为TAM耐药的核心基因，并通过患者经他莫昔芬治疗后预后数据及文献证据验证核心基因表达与他莫昔芬耐药密切相关.前期研究也表明，5个核心基因会对他莫昔芬的药效产生重要影响.同时除了PRKC外，核心基因也可以作为他莫昔芬治疗ER阳性乳腺癌的预后标记物,如表1所示.

表1 耐药核心基因与他莫昔芬治疗效果的关系

信号转导和转录激活因子3(STAT3)在CW5群落中是具有最高节点度的表达差异基因，且其表达对他莫昔芬的疗效密切相关.有文献报道，降低STAT3在TAM耐药细胞中的表达量，可以增强细胞对他莫昔芬的敏感性[10].因此，若将STAT3抑制剂与他莫昔芬联合使用，具有逆转他莫昔芬耐药的潜力.转录因子Egr1是包括TGFβ1、PTEN、p53和纤维连接蛋白等多种肿瘤抑制因子的直接调节因子.这些因子的下游路径显示出相互作用的多个节点，这表明存在一个功能性的抑制因子网络，用于维持正常的生长调节和抵抗转化变体的出现.在侵袭性乳腺肿瘤中，p130Cas/BCAR1(Crk相关底物/乳腺癌抗雌激素抵抗1)水平升高，并且与乳腺癌的他莫昔芬耐药有关[13].在MCF-7细胞中存在着一个正反馈回路，p130Cas正调控EGR1和NAB2，进而诱导p130Cas的表达.重要的是，与MCF-7相比，TAM耐药株中NAB2表达的增强并增加EGR1与bcar15′区结合的可能性，导致TAMR细胞中p130Cas/BCAR1水平的持续升高.PRKCA 基因编码蛋白激酶Cα蛋白.Kim等[15]研究了通过抑制蛋白激酶Cα(PKCα)的活性来逆转人类多药耐药乳腺癌细胞的耐药性，这与P-糖蛋白(Pgp)介导的抗癌药物外排有关.Pgp是一种170kda的跨膜蛋白，可介导抗癌药物从细胞外排.Pgp过表达在多药耐药(MDR)细胞中有明显作用.耐药MCF-7/ADR细胞的PKCα基础活性高于药物敏感的MCF-7细胞.PKCα抑制剂Ro-31-7549能有效抑制PKCα的活性，并通过抑制MCF-7/ADR细胞内PKCα活性，观察到其逆转阿霉素在细胞内的积累，显著提高阿霉素对MDR细胞的治疗效果.这些结果表明，传统抗癌药物通过抑制PKCα活性来克服MDR是有潜力的.

2.4 STAT1和CDH1与他莫昔芬治疗预后的关系

选取STAT1(信号转导与转录激活因子1)和CDH1(钙黏蛋白1)作为示例，进一步探究核心基因与ER+乳腺癌经他莫昔芬治疗预后的关系.从GEO GSE9893数据集中得到乳腺癌患者他莫昔芬耐药的表达谱和生存情况数据，进行分析后，绘制成图5.STAT1和CDH1基因的表达量是具有个体化差异的.将基因的表达值按照递增顺序进行排序，选取前39位乳腺癌患者与后39位乳腺癌患者的基因表达值作为划分低表达组和高表达组的阈值(如图5(a)所示).如图5(b)所示，分别绘制基因高表达组和低表达组的生存曲线.图5表明STAT1和CDH1的高表达与他莫昔芬治疗ER+乳腺癌的预后差相关(STAT1,P=0.02；CDH,P=0.01).

(a)STAT1基因高低表达组的阈值

信号转导与转录激活因子(STAT)家族分子是正常细胞生理的关键调节因子，在许多病理条件下都有重要作用它与肿瘤发生、细胞周期进展、细胞存活、转化和血管生成有关.图6所示为从KEGG中提取与STAT1有关的信号通路.目前，STAT1被普遍认为是一种肿瘤标志物,具有干扰抗肿瘤免疫反应作用[16].STAT1的表达先前已经被证明与对DNA损伤和基因毒性药物的抵抗有关，这可能与不良结果有关.在长期暴露于雌激素的MCF7细胞中，STAT1表达水平增加，这再次被认为与放疗和化疗抵抗力增加有关.因此，在内分泌耐药的乳腺癌中，STAT1的表达和激活可能增加，STAT1抑制剂可能对耐药细胞有效[12].这些结果表明，STAT1可能是内分泌抵抗疾病的一个可行靶点，值得进一步研究STAT1靶向抑制剂.

钙黏蛋白1(cadherin 1,CDH1)基因位于16q22.1染色体上的编码一种钙依赖性跨膜糖蛋白E-cadherin，它是上皮细胞标志物，主要在上皮细胞中表达.参与的信号通路如图6所示.E-cadherin是上皮细胞中不可或缺的跨膜黏附分子，所以它对于成熟细胞间形成黏附连接时所必品，通过与联蛋白相互作用，调节细胞的分化、增殖、凋亡以及信号转导，有助于成熟细胞间黏附连接的形成[17].如果E-cadherin出现低表达，那么就可能导致Ras (一种小的GTP酶)、丝裂原活化蛋白激酶 (mitogen-activated protein kinase,MAPK)、以及与Ras相关的C3肉毒素底物1 (Rac1) 等多条致癌信号通路激活.而造成E-cadherin表达的降低或缺失的影响因素中就有CDH1基因突变.一旦E-cadherin出现表达异常后就会发生干扰细胞与细胞间黏附连接的强度，使细胞分离加快，细胞从原发肿瘤部位逃逸的现象.前期研究表明，相比于他莫昔芬敏感细胞，CDH1基因在他莫昔芬耐药细胞中被高甲基化，利用5-azacytidine处理细胞后，可实现耐药逆转.

(a)STAT1相关信号通路

3 结论

本论文主要是在基于网络分析的情况下对他莫昔芬在乳腺癌治疗中的耐药机制进行研究，从中筛选出他莫昔芬耐药的关键基因.通过对蛋白互作网络的处理，将网络拓扑结构运用于分析中，将基因进行聚类后对关键的群落进行分析之后得到信号通路，接着从信号通路中对预后标记物进行筛选，最终选出STAT1和CDH1两个基因.用基因表达和TAM治疗的预后数据对STAT1和CDH1进行生存分析，发现了STAT1和CDH1低表达的患者具有更好的预后情况.

综上所述，耐药核心基因具有作为他莫昔芬治疗的预后生物标记物，以及实现他莫昔芬耐药逆转的可行靶点，值得进一步研究针对核心基因的靶向抑制剂.