胰腺导管腺癌差异基因表达及临床预后分析

李自昂，王帆，王晓月，王倩，程洁，林军

(武汉大学中南医院消化内科，武汉 430072)

胰腺导管腺癌(pancreatic ductal adenocarcinoma，PDAC)是一种极度恶性肿瘤，5年存活率<7%，中位生存期约为6个月，预后差[1-2]。PDAC能够显示出很大的细胞异质性，具有明显的上皮和间充质癌细胞群，在PDAC小鼠模型中，白血病抑制因子可以通过调节癌细胞分化和上皮-间充质转化状态促进癌症进展，而在PDAC患者中，间皮素高表达患者的中位生存期明显短于低表达者，表明间皮素能够促进胰腺癌的进展[3-5]。由于早期广泛转移和对化疗药物的耐药作用，PDAC的死亡率较高[6]。有研究表明，免疫调节细胞(如髓样来源的抑制细胞)在肿瘤发展和化学耐药中起重要作用[7]。预计到2030年，PDAC将成为所有恶性肿瘤中的第二大常见癌症类型[8]，但目前对该病的研究较少。因此，迫切需要寻找新的生物标志物，能够预测PDAC的预后并作为PDAC潜在的治疗靶点及检测手段。本研究主要筛选与PDAC发展高度相关并影响其预后的基因，以期为未来研究胰腺癌进展途径及靶向治疗提供新思路。

1 资料与方法

1.1微阵列数据集来源 从GEO(Gene expression omnibus)在线数据库(https://www.ncbi.nlm.nih.gov/geo/)获取4个基因数据集[GSE(gene series expression)28735、GSE62165、GSE62452、GSE71989]和2个微RNA(microRNA,miRNA)数据集(GSE25820、GSE41372)。使用R 3.6.1中的程序包“limma”读取，标准化和筛选差异表达的基因[9]。筛选|logFC|≥1且调整后P<0.05的基因作为差异基因。利用FunRich软件筛选均在各数据集差异表达的基因和miRNA[10]。

1.2基因本体论(gene ontology，GO)富集通路和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)通路分析 使用注释、可视化和集成发现数据库(the Database for Annotation，Visualization and Integrated Discovery，DAVID)(https://david.ncifcrf.gov/)进行差异基因GO富集和KEGG通路分析[11]。高度富集于信号通路[伪发现率(false discovery rate，FDR)<0.05]的基因被认为参与PDAC的生物代谢途径。R3.6.1的“GOplot”包用于可视化富集通路[12]。

1.3蛋白质相互作用(protein-protein interaction,PPI)网络和模块的构建 采用String数据库(https://string-db.org/)分析基因间的相互作用[13-14]。Cytoscape及其包含的功能组件MCODE(Molecular complex detection)用于构建可视化PPI网络并分析不同的生物模块[15-16]。

1.4miRNA靶基因的预测和预后分析 采用miRwalk网站(http://mirwalk.umm.uni-heidelberg.de/)和miRNet网站(https://www.mirnet.ca/)构建miRNA-基因互作网络[17-18]。综合差异基因、富集分析、miRNA靶基因筛选基因，并使用UCSC Xena在线工具(http://xena.ucsc.edu/)构建基因的表达聚类图[19]。最后，使用Kaplan-Meier plot网站(http://kmplot.com/)评估这些基因对PDAC的临床预后[20]。该网站根据基因中位表达水平分为高、低表达组，进而分析两组5年生存率的差异。

2 结 果

2.1差异表达基因和差异表达miRNA GSE28735数据集包括45个肿瘤组织和45个相邻正常组织。GSE62165数据集包含118个PDAC肿瘤组织和13个非肿瘤组织。GSE62452数据集包括69个胰腺肿瘤和61个相邻的正常组织。GSE71989数据集包含8个正常胰腺组织和14个PDAC组织。GSE25820数据集的miRNA表达谱包括4个正常胰腺导管细胞，4个腺泡细胞，5个慢性胰腺炎和原发性PDAC的转移组织。GSE41372数据集包含15个非肿瘤样品和15个肿瘤样品。从4个数据集中共获得217个差异基因，见图1；从2个miRNA数据集中获得13个miRNA(hsa-miR-199b-5p、hsa-miR-216a、hsa-miR-217、hsa-miR-145、hsa-miR-302c、hsa-miR-100、hsa-miR-337-3p、hsa-miR-770-5p、hsa-miR-193a-5p、hsa-miR-1236、hsa-miR-455-3p、hsa-miR-342-3p、hsa-miR-196a)，见图2。

图1 4个基因数据集差异表达基因

图2 2个miRNA数据集差异表达miRNA

2.2GO富集和KEGG通路分析 DAVID数据库中发现217个基因高度富集于8个生物学过程：细胞外基质组织、胶原分解代谢过程、细胞黏附、细胞外基质分解、蛋白水解、骨骼系统发育、胶原纤维组织、内胚层细胞分化，见图3；8个细胞成分通路：细胞外空间、细胞外区域、细胞外基质、细胞外泌体、蛋白质细胞外基质、细胞表面、内质网腔和胶原三聚体，见图4；5个分子功能通路：细胞外基质的结构成分、丝氨酸型内肽酶活性、整联蛋白结合、钙离子结合、胶原蛋白结合，见图5；4个KEGG通路：蛋白质的消化吸收、细胞外基质-受体相互作用、胰腺分泌、黏着斑，见图6。

图3 8个重要的生物学过程

图4 8个重要的细胞成分

图5 5种重要的分子功能

KEGG：京都基因与基因组百科全书

2.3PPI基因网络的构建和预后分析 通过GEO数据库共获得150个上调基因和67个下调基因。连接系数高于0.4的167个差异基因被用来构建PPI基因网络互作图。从这些互作基因中分析出10个重要的模块，见图7。在miRwalk数据库中分析了6 711个miRNA靶基因。通过miRNA-基因网络图，发现在胰腺癌中镁依赖性蛋白磷酸酶1δ(protein phosphatase magnesium-dependent 1 delta,PPM1D)是3个miRNA(hsa-miR-342-3p、hsa-miR-217、hsa-miR-770-5p)的共同靶点，见图8。筛选出10个与PDAC显著相关的基因，在UCSC xena数据库中，从癌症基因组数据集提取出4例正常样本和177例肿瘤样本，筛选出的10个基因均差异表达于肿瘤组织与非肿瘤组织(P<0.05,|logFC|>1)，见图9。Kaplan-Meier plot网站收集了172例PDAC患者的临床信息用于进行生存分析，总生存期设定为7 642 d，生存曲线显示，纤连蛋白Ⅰ(fibronectin 1，FN1)、血小板反应蛋白2型(thrombospondin 2,THBS2)、Ⅻ 型胶原蛋白α-1链(collagen type Ⅻ alpha 1 chain，COL12A1)、Ⅵ 型胶原蛋白α-3链(collagen type Ⅵ alpha 3 chain，COL6A3)和载脂蛋白L1型(apolipoprotein L1，APOL1)高表达组患者的生存率较低，而PPM1D高表达组患者的生存率较高，见图10～15。

PPI:蛋白质相关作用；红色节点代表上调基因，蓝色节点代表下调基因，线条代表每个基因之间相互作用

黄色节点代表miRNA，蓝色节点代表基因，红色节点代表PPM1D基因

FN1：纤连蛋白Ⅰ；IGFBP5:胰岛素样生长因子结合蛋白5；CDH11:钙黏蛋白11；THBS2：血小板反应蛋白2型；COL12A1：Ⅻ型胶原蛋白α-1链；COL5A2：Ⅴ型胶原蛋白α-2链；CP:外壳蛋白；APOL1：载脂蛋白L1型；COL6A3：Ⅵ型胶原蛋白α-3链；PPM1D：镁依赖性蛋白磷酸酶1δ；红色表示该基因在样本中高表达，蓝色表示该基因在样本中低表达；在样本组中，灰色表示正常组，黄色表示肿瘤组

PPM1D：镁依赖性蛋白磷酸酶1δ

FN1：纤连蛋白Ⅰ

THBS2：血小板反应蛋白2型

COL12A1：Ⅻ型胶原蛋白α-1链

COL6A3：Ⅵ型胶原蛋白α-3链

APOL1：载脂蛋白L1型

3 讨 论

miRNA是一个小的非编码调控RNA的大家族[21]。它们控制人类60%以上的蛋白质编码和基因表达[22]。miRNA参与PDAC发展过程中的大多数遗传突变，如几乎所有PDAC病例中，miRNA介导的KRAS突变激活是胰腺导管上皮细胞转化为腺癌的最早的遗传变化[23]。这说明miRNA能够通过调控基因活动参与癌症的发展。我国胰腺癌导致的死亡在癌症相关死亡中占比从2009—2019年增加了9%，因为人口老龄化及生活饮食等因素影响，这一比例还会继续上升[24]。因此，对胰腺癌的机制进行研究将成为今后研究的重点。

本研究基于公众数据库，发现miR-342-3p、miR-217、miR-199b-5p,miR-770-5p调控的基因THBS2、FN1、COL12A1在胰腺癌组织中高表达，PPM1D、COL6A3、APOL1在胰腺癌组织中低表达。在功能富集方面，这些基因显著富集于胶原分解代谢途径、细胞黏附、细胞外泌体、丝氨酸型内肽酶活性和细胞外基质受体相互作用等通路。Kaplan-Meier曲线显示，PPM1D、FN1、THBS2、COL12A1、COL6A3、APOL1表达与PDAC患者的预后显著相关。PPM1D编码的蛋白质是玉米2C型丝氨酸/苏氨酸蛋白磷酸酶家族的成员。该基因的表达以p53依赖的方式响应各种环境压力而被诱导[25-26]。PPM1D通过下调凋亡刺激p53蛋白2增强Wnt/β联蛋白途径，促进胰腺癌细胞迁移和侵袭[27]。FN1编码纤连蛋白，是糖蛋白的一种，在血浆中以可溶性二聚体的形式存在，在细胞表面和细胞外基质中以二聚体或多聚体形式存在。纤连蛋白参与细胞黏附和迁移过程，包括胚胎发生、伤口愈合、血液凝固、宿主防御和转移[28-29]。在切除的PDAC组织中，FN1高表达于肿瘤瘤体更大且远端淋巴结转移的肿瘤组织[30]。THBS2编码的蛋白质属于血小板反应蛋白家族，是二硫键连接的同型三聚体糖蛋白，介导细胞与细胞之间以及细胞与基质之间的相互作用[31]。临床研究表明，PDAC的血清THBS2水平显著高于高危人群，并且该基因在预后较差的患者中表达水平更高[32]。COL12A1编码Ⅻ型胶原蛋白的α链，是具有三重螺旋间断的原纤维相关胶原蛋白。Ⅻ型胶原蛋白是与Ⅰ型胶原蛋白结合的同型三聚体[33]。目前关于COL12A1与PDAC进展相关的研究较少。但有研究表明，COL12A1通过促分裂原活化的蛋白激酶通路维持的正反馈促进胃癌细胞的转移[34]。COL6A3是Ⅵ型胶原蛋白编码的α-3链，这是Ⅵ型胶原蛋白(在大多数结缔组织中发现的珠状细丝胶原蛋白)的三个α链之一。Ⅵ型胶原基因突变与Bethlem肌病有关，Bethlem肌病是一种罕见的常染色体显性遗传性近端肌病，伴有儿童早期发作[35]。Kang等[36]研究表明，COL6A3在胰腺癌组织及胰腺癌患者血清中显著高表达，表明该基因参与PDAC的发生发展。APOL1编码与载脂蛋白A-1结合的高密度脂蛋白。载脂蛋白A-1是一种相对丰富的血浆蛋白，是高密度脂蛋白的主要载脂蛋白。它参与血浆中大多数胆固醇酯的形成，并促进胆固醇从细胞中流出。APOL1能够导致肾损伤并且在肾脏足细胞中高表达，但损伤的具体机制仍不明确[37]。胰腺癌血清中APOL1基因表达水平显著高于普通人群，受试者工作特征曲线结果显示其可能成为PDAC的诊断标志物[38]。

综上所述，10个miRNA靶基因在PDAC与正常组织中差异表达并参与其发生发展，其中PPM1D、FN1、THBS2、COL12A1、COL6A3和APOL1与PDAC预后相关。因此，进一步研究这些基因在胰腺癌中的作用，可能是治疗胰腺癌，提高患者生存率的重要方向。