基于TCGA数据库探究可溶性因子对乳腺癌患者的预后作用

刘孟晨 孙贻安 李静蔚

乳腺癌是以乳房组织中乳腺上皮细胞的增殖失控为主要表现的恶行异质性疾病。有统计表明，仅2018年在全球范围内便有60多万人死于乳腺癌[1]，且存活的人也因焦虑[2]、放疗[3]、乳房切除[4]等因素极易诱发抑郁症等精神类疾病。因此，乳腺癌对女性群体的生命及身心健康产生了极大的威胁。值得注意的是，乳腺癌患者的死亡原因大多可归咎于转移及相关并发症[5]。在转移的过程，肿瘤细胞从原发肿瘤脱离并进入血流[6]。这些循环肿瘤细胞最终滞留在远处器官的毛细血管床中，导致继发部位产生转移性集落[5]。因此，可溶性因子在转移等过程中发挥重要作用，它们的变化往往能加剧乳腺癌患者病情恶化甚至死亡的风险[7]。

可溶性因子是机体内可溶于血清的各类细胞因子，可通过增加炎性损伤、干预信号传导、调节免疫等方式参与各类疾病的发展[8]。研究表明，骨桥蛋白、表皮生长因子和IL-6等可溶性因子会参与乳腺癌细胞和巨噬细胞之间的细胞相互作用[9]。成纤维细胞分泌的可溶性因子可募集单核细胞进而促进三阴性乳腺癌患者反应性基质的激活使其更具备侵袭性[10]。乳腺癌微环境下所分泌的可溶性因子甚至可触发异常信号传导，刺激骨微环境中的破骨细胞分化，进而加大骨转移的风险[11]。这些研究表明，可溶性因子极大程度参与了乳腺癌的发病进展及预后，但大多研究并未基于整体性及系统性对可溶性因子的作用进行筛选和分析。基于此，本研究拟基于癌症基因图谱（the cancer genome Atlas,TCGA）数据库中乳腺癌患者的基因表达等数据，探究可溶性因子在乳腺癌中的变化与影响。

1 资料和方法

1.1 数据的采集 通过TCGA数据库（https://portal.gdc.cancer.gov/projects）下载乳腺癌患者的转录组数据及临床数据。删除生存数据缺失的样本后，通过R语言中的limma包对肿瘤组织与正常组织进行差异分析，以FDR值＜0.05为条件筛选差异基因。通过基因集富集分析（gene set enrichment analysis,GSEA）数据库（http://www.gsea-msigdb.org/gsea/msigdb/cards/NABA_SECRETED_FACTORS.html）下载可溶性因子数据集。并通过survival包对样本的生存时间和生存状态与可溶性因子基因相关表达数据进行生存分析得到可溶性因子预后相关基因。将预后相关可溶性因子数据集与TCGA中差异基因取交集后得到可溶性因子预后相关基因表达数据。以热图和森林图的方式对可溶性因子预后相关基因的差异表达及对乳腺癌患者预后的影响进行展示。

1.2 蛋白质-蛋白质相互作用（protein-protein inter‐action,PPI）、基因本体论（gene ontology,GO）与京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes,KEGG）分析 通过基因搜索工具（search tool for recurring instances of neighbouring genes,String）数据库（https://string-db.org/cgi/input.pl）对交集基因进行PPI分析，并根据TCGA数据中乳腺癌患者转录组数据的表达计算可溶性因子预后相关基因的相关性，最后通过R语言对可溶性因子预后相关基因进行GO、KEGG分析。

1.3 风险模型的构建及主成分分析、随机临近嵌入分析 Glmnet包是构建Lasso回归模型的重要工具，通过样本中所有可溶性因子预后相关基因的表达及危险系数的乘积作为风险值为基础建模。以所有样本风险值的中位数为界限，划分高低风险组。通过survival、survminer、timeROC包比较高低风险组患者的生存差异并获取受试者工作特征曲线。使用Rtsne包对高低风险组进行主成分分析和T分布随机临近嵌入分析展现高低风险组基因的分布情况并进行降维可视化。

1.4 高低风险组免疫浸润模式分析 基于单样本基因富集分析（single sample GSEA,ssGSEA）的方法对每个样本进行免疫打分，并根据模型中高低风险组中可溶性因子集内核心基因的表达进行免疫细胞及免疫功能的差异分析。

2 结果

2.1 数据的分析及处理 通过对TCGA数据集中的数据进行整理，发现数据库中共有1 222个乳腺癌相关样本，其中正常样本有113个，乳腺癌样本1 109个。GSEA数据库中可溶性因子基因共有343个，其中有42个基因与乳腺癌患者的预后相关（图1a）。将TCGA数据库中乳腺癌患者体内的差异基因与预后相关的可溶性因子取交集后共得到27个预后相关的可溶性因子的差异基因（图1b），并通过森林图（图1c）及热图（图1d）的形式展现交集基因的风险因素及表达。从图1c的结果中可以看出，除FGF8及WNT11为高风险致病基因外，其他可溶性因子均对患者预后影响较小，这可能是由于众多可溶性基因综合作用下单基因产生的生物学效应相互之间存在一定的抵消。从图1d的结果来看，FGF8、WNT11、CCL23、FLT3LG、TNFSF12、IL-7、CXCL1、FGFBP1、IL-17B、S-100B、SFRP1、NRG1、IFNE、IL-4、CXCL14、IL-16、XCL1在正常组中表达明显，但在肿瘤组中表达较低，故这些基因可能有抑癌的作用。CCL25、CCL5、CXCL13、CXCL9、IL-12B、IL-18、IL-24、LTA、LTB、TNFSF12、WNT7B等基因在正常组中表达较低但在肿瘤组中表达较高，这表明这些基因可能有抑癌的作用。

图1 数据的分析与处理（a：42个可溶性因子预后相关基因；b：可溶性因子预后相关基因与差异基因的韦恩图；c：预后相关可溶性因子森林图；d：预后相关可溶性因子热图）

2.2 PPI网络图及GO、KEGG分析 PPI网络图显示，在可溶性因子预后相关基因相互间存在较强的相互作用（图2a），其中相关度≥10的靶点为IL-7、IL-4、IL-16、CCL5、CXCL13、CXCL9、CXCL1、LTA，这表明这些靶点与其他基因之间具有较强的相互作用。相关性分析发 现，FLT3LG、CXCL9、IL-12B、IL-18、IL-24、LTA、LTB、TNFSF12的表达之间存在明显的正相关（图2b）。GO分析表明,可溶性因子预后相关基因在生物过程中可以参与补体的激活，循环免疫球蛋白介导的体液免疫反应，体液免疫反应等（图2c）。在细胞组分中涉及免疫球蛋白复合物、免疫球蛋白复合物，循环的浆膜外侧血液微粒子等。在分子功能中多涉及抗原结合、免疫球蛋白受体结合、趋化因子活性等。KEGG分析（图2d）表明，可溶性因子预后相关基因与病毒蛋白与细胞因子和细胞因子受体的相互作用、细胞因子-细胞因子受体的相互作用、趋化因子信号传导途径、造血细胞谱系、细胞粘附分子、NF-κB信号传导途径、Th1和Th2细胞分化、TNF信号传导途径、IL-17信号传导途径、癌症中的PD-L1表达和PD-1检查点通路、JAK-STATd等信号传导途径密切相关。

图2 PPI、GO、KEGG分析（a：PPI分析图；b：相关性分析图；c：GO分析图；d：KEGG分析图）

2.3 风险模型的构建及主成分分析、随机临近嵌入分析 使用Lasso回归分析构建风险模型，最终筛选出14个基因参与构建风险模型，风险模型计算公式为：风险值=CCL25×（-0.267 761 762 190 859）+CXCL1×（-0.0343324083608361）+CXCL13×（-0.0427683660862547）+CXCL14×（-0.0498671701605376）+FGF8×1.2580095296421+FGFBP1×（-0.0758992553296801）+FLT3LG×（-0.159 083 190 747 296）+IL-12B×（-0.088 335 107 649 380 9）+IL-24×（-0.0514706044485947）+IL-4×（-1.22386864779545）+NRG1× （-0.252 446 906 132 768 ）+S-100B×（-0.0277742304519763）+WNT11×（-0.0736641709388857）+WNT7B×0.123 111 236 484 031。

以风险值的中位数为界，将患者划分为高、低风险两组，依据患者的分组及风险值绘制患者的生存状态（图3a）及风险值热图（图3b）。使用主成分分析、随机临近嵌入分析将患者的分组信息进行可视化（图3c、3d）。通过比较两组间患者的总生存率，发现高风险患者的生存时间显著低于低风险组（图3e，P＜0.01）。通过R语言绘制的ROC曲线发现，患者1～3年的AUC均＞0.65，这表明本模型较为可靠（图3f）。最后通过单因素及多因素分析，发现可溶性因子预后相关基因的致病性与乳腺癌患者的年龄、性别、分期均显著相关（图3g、3h）。

图3 风险模型的构建及主成分分析、随机临近嵌入分析（a：高、低风险组生存状态图；b：高、低风险组风险值热图；c：高、低风险组主成分分析图；d：高、低风险组随机临近嵌入分析；e：高、低风险组生存分析；f：高、低风险组ROC曲线下面积图；g：高、低风险组单因素Cox回归；h：高、低风险组多因素Cox回归）

2.4 高低风险组免疫浸润模式分析 通过ssGSEA对高低风险组中免疫细胞及免疫功能进行了差异分析，研究表明可溶性因子可以很大程度上改变免疫细胞在高低风险组的评分。在高风险组中所有的免疫细胞及免疫功能评分均发生了降低。这表明aDCs、B细胞、CD8+T细胞、DCs、iDCs、巨噬细胞、肥大细胞、中性粒细胞、NK细胞、pDCs、辅助T细胞等免疫细胞在高风险患者中自身的功能有明显的抑制，见图4。

图4 高低风险组免疫浸润模式分析（a：高、低风险组免疫细胞差异分析；b：高、低风险组免疫功能分析；**P＜0.01）

3 讨论

可溶性因子是人体内各类免疫细胞、内分泌细胞、神经细胞分泌的具有高度诱导性细胞因子，深度参与了乳腺癌发展、转移、恶化等过程[12]。

本研究表明，可溶性因子能通过多种途径影响乳腺癌的发生、发展。如CXCL1、CXCL9、CXCL13、CX‐CL14等可通过影响机体的免疫反应在乳腺癌中发挥作用。其中CXCL1是肿瘤相关巨噬细胞分泌的最丰富的趋化因子，可以将各种基质细胞募集到肿瘤环境中，构建促进肿瘤细胞生长、血管生成和转移的环境[13]。CXCL9可通过调节乳腺癌中免疫细胞浸润的丰度，进而干预乳腺癌的发展[14]。CXCL13是重要的趋化因子之一，B淋巴细胞聚集的标志物，在B淋巴细胞的归巢、迁移和积累中起关键作用[15]。中性粒细胞被认为是提供抵御入侵病原体的主要免疫细胞，它们的募集在很大程度上被CXCL1和CXCL2等可溶性因子所诱导[16-17]。CCL5可募集T细胞、B细胞、嗜酸性粒细胞、嗜碱性粒细胞、中性粒细胞、巨噬细胞和成纤维细胞等各类免疫细胞[18-19]，这与本研究后续免疫浸润分析的结果高度重合。CCL5还可通过PI3K/Akt通路激活αvβ3整合素以促进细胞迁移，而αvβ3整合素又能激活IKKα/β和NF-κB通路形成联机反应[20]。当然，其他白介素家族所导致的炎症损伤也不能被忽视。有研究表明，IL-4由巨噬细胞和T淋巴细胞产生，可控制T辅助细胞的成熟，IL-4诱导的谷氨酰胺代谢对乳腺癌的生长有较大的刺激作用[21]。IL-7最近参与了产生IL-17A的先天性样T细胞的选择性扩增和功能，包括自然杀伤T细胞和Th17细胞[22-23]。此外，IL-7水平与前列腺癌患者的不良预后相关[24-25]，与乳腺癌患者的肿瘤侵袭性正相关[26]。还有研究表明，FGF8作为一种可调节乳腺癌生长的一种可溶性生长因子，可通过调控类固醇激素的分泌影响乳腺癌的预后[27]。这些研究均显示可溶性因子自身具有较高的致病性，而本研究基于这些关键基因作为预后模型在后续的验证中也表现出极强的特异性。

更重要的是，可溶性因子不仅自身具有一定致病性，还可以与其他因子的相互作用下形成联级反应。KEGG结果中，可溶性因子所涉及的NF-κB、Th1/Th2、PD-L1/PD-1等通路均与乳腺癌的发生发展密切相关。据报道，乳腺癌中NF-κB的激活可上调细胞周期蛋白D1、细胞周期蛋白依赖性激酶2（CDK2）和c-Myc[28-30]的表达，从而加速细胞周期进程并导致细胞增殖失控。NF-κB还调节IL-1β、TNFα、EGFR和HER2等可溶性因子的分泌以促进肿瘤细胞的生长[31]。NF-κB 活性的下降能改变Bcl-2[32]、IAP蛋白（XIAP、cIAP-1/2）[32]和 TNF 受体相关因子（TRAF）1/2[33]等细胞死亡调节基因的表达，进而导致抗凋亡和促存活基因的上调并抑制对化疗药物的细胞凋亡反应。Th1/Th2的平衡也影响着乳腺癌的发生、发展。人体内幼稚的CD4+T细胞在不同的细胞因子刺激下分化为Th1和Th2细胞。Th1细胞产生具有抗癌活性的IFN-γ和IL-12。Th2细胞介导体液免疫并产生IL-4和IL-10。与Th1反应相反，Th2反应实际上可以促进癌症进程。Th2细胞分泌的IL-4抑制IFN-γ的分泌[34]。IL-10可下调Th1细胞的增殖和IFN-γ的产生[35]。因此Th1/Th2的平衡会极大影响着乳腺癌患者肿瘤微环境的状况。PD-1/PD-L1是一个免疫检查点，被认为是乳腺癌治疗的重要靶点。PD-1与其配体PD-L1结合，限制T细胞活性，并使乳腺癌细胞逃避免疫系统[36-37]。在正常情况下，PD-1/PD-L1通路通过负性调节免疫反应防止过度刺激并维持对自身抗原的免疫耐受，相反，PD-L1在乳腺癌中过度表达，导致肿瘤微环境中的免疫功能受损[38]。故这些可溶性因子可能是通过影响这些信号传导干预乳腺癌的进展。

综上所述，可溶性因子影响着乳腺癌患者的预后。它们自身对乳腺癌组织的损伤及相互间的联级反应是推动乳腺癌发展的重要因素。此外，可溶性因子还对乳腺癌患者的免疫功能产生了重要影响。但这此过程所涉及的因素较多，故本文并未详细论述。后续课题组将在严格执行伦理审查的情况下广泛收集临床病理样本，通过单细胞测序等手段明确可溶性因子对各类免疫细胞的影响。