汪圣毅，程 彦，李旭升，闫亚飞，张尚鑫，闫 强，李永翔
摘要 目的基因集富集分析（geneset enrichment analysis, GSEA）方法寻找胃癌抵抗顺铂的免疫标记基因（immunologicsignature genes, ISGs）。方法 用GEO数据库的GSE94 714数据集，GEO2 R分析差异基因，观察条件筛选对基因数的作用，GSEA纳入耐药、未耐药组胃癌细胞的全部差异表达基因，与分子标签数据库比较，获取ISGs，交集筛选，KaplanMeier Plotter分析交集基因对胃癌预后的影响。结果 差异表达基因共34183个，其中上调12452个、下调17381个，筛选差异倍数的增加使排除基因数增加。GSEA富集到标准化富集评分（NES）排序前6的条目（P< 0.01），其中的交集基因包括线粒体核糖体蛋白L12 、富含脯氨酸蛋白13、毛状蛋白样F-肌动蛋白结合蛋白1、聚（RC）结合蛋白1、艾杜糖2-硫酸酯酶、LIM结构域2、富含嘌呤元素结合蛋白A、小视觉叶同源物、CCR4-非转录复合物亚单位3、转化生长因子β1、二酰甘油激酶ζ、接头蛋白2，12个基因与胃癌的总生存时间有关，均具有统计学意义（P< 0.05）。结论GSEA方法可有效获取胃癌顺铂抵抗的ISGs，新发现的基因作为潜在靶点，可促进胃癌化疗抵抗机制的研究。
关键词 胃癌；化疗抵抗；基因集富集分析；基因表达数据库；免疫标记基因；预后
中图分类号 R735.2
文献标志码 A 文章编号
Theenrichment analysis of immunologic signature genes associated with cisplatinresistance in gastric cancer
WangShengyi, Cheng Yan, Li XuSheng, et al
(Deptof Gastrointestinal Surgery, Dept of General Surgery, The First AffiliatedHospital of Anhui Medical University, Hefei 230022)
Abstract Objective To find theimmunologic signature genes (ISGs) associated with cisplatin resistance ingastric cancer based on gene set enrichment analysis (GSEA). Methods GSE 94 714 dataset from GEO database was used, and the differentiallyexpressed genes (DEGs) were analyzed with GEO 2 R. The effects ofconditional gene screening on the gene numbers were observed. GSEA included allDEGs from drug resistant and non-resistant groups of gastric cancer cells. The DEGswere compared with the molecular signatures database (MSigDB), and the obtainedISGs were screened for intersection, and the effects of ISGs on the prognosisof gastric cancer were analyzed by Kaplan Meier Plotter method. Results A total of 34 183 DEGs included 12 452 up-regulated and 17 381down-regulated genes. The increased fold change (FC) added the number ofexcluded genes. Six entries with the top normalized enrichment score (NES) wereidentified by GSEA (P<0.01).The intersection ISGs included mitochondrial ribosomal protein L 12 (MRPL 12),proline-rich protein 13 (PRR 13), coactosin like F-actin binding protein 1 (COTL1) , poly (RC) binding protein 1 (PCBP 1), Iduronate 2-sulfatase (IDS), LIMdomain containing 2 (LIMD 2), purine rich element binding protein A (PURA), smalloptic lobes homolog (SOLH), CCR4-NOT transcription complex subunit 3 (CNOT 3), transforminggrowth factor beta 1 (TGFB 1), diacylglycerol kinase zeta (DGKZ), and dockingprotein 2 (DOK 2). Twelve ISGs were associated with the overall survival ofgastric cancer (P<0.05).Conclusion The GSEA method can effectively extract the ISGsof cisplatin resistance in gastric cancer. Newly discovered ISGs, as potentialtargets, can enhance the study of chemoresistant mechanisms in gastric cancer.
Keywords gastric cancer；chemoresistance； gene set enrichment analysis；GEO；immunologic signature gene； prognosis
胃癌在我国高发，手术是主要治疗方法，辅助化疗可改善预后[1]，但化疗抵抗制约疗效，与基因、通路的作用有关[2,3]，免疫标记基因（immunologicsignature genes, ISGs）与化疗抵抗有关[4]。外周髓磷脂蛋白22（Peripheralmyelin protein 22，PMP 22）增强胃癌顺铂抵抗[3]，是否通过ISGs调控尚不清楚。常规富集分析可发现PMP22相关的化疗抵抗通路和基因[2]，但预设差异基因的倍数和统计P值，不能较为全面地反映整体信息[5]。基因集富集分析（Geneset enrichment analysis, GSEA）纳入全部的差异基因，具有不预筛组间差异基因的优点[6]，可发现细微变化。为明确ISGs与PMP22相关胃癌顺铂抵抗的关系，用GSEA法分析GSE94 714数据集，寻找PMP22下游的ISGs，验证新发现ISGs对胃癌总生存的影响，为深入揭示胃癌化疗抵抗的机制提供参考。
1资料与方法
1.1 表达谱数据 用基因表达数据库（geneexpression omnibus，GEO）数据库的GSE94714数据集，包含6个样本，3个样本（GSM2 481 329~GSM 2 481 331）的胃癌细胞MGC-803进行了PMP22基因敲减，另外3个样本（GSM2 481 332~GSM 2 481 334）的PMP 22未行敲减处理。GEO2 R分析差异基因，设置分组：PMP22基因敲减3个样本为A组，未敲减3个样本为B组，保存所有差异基因至excel表，得到差异基因表。
1.2 芯片数据及探针基因名转换 下载GSE 94 714的soft文件，导入excel，复制探针和对应基因名称（GENE_SYMBOL），形成新的excel表即芯片数据表，探针名称升序排序。差异基因表中的探针名称进行相同的升序排序，复制ID、adj.P.Val、logFC共3列到芯片数据表，excel数据高级筛选，匹配出探针对应的基因名称，包含差异基因表，形成探针基因名转换表，定位GENE_SYMBOL列的空值，删除没有GENE_SYMBOL的所有行，另存为包含GENE_SYMBOL的差异基因表，用于后续分析。
1.3 差异表达分析 包含GENE_SYMBOL的差异基因表，导入在线分析工具（EasyTools v 2.2 by Chris Lou）中（www.chrislifescience.club:3838/R/AnnoE2）的Volcano Online，绘制火山图。差异倍数（Fold Change, FC）以2为底的对数（Log 2 FC或Log FC）分别设置为1、1.5，P值<0.05，计算2种LogFC条件下，有差异（上调、下调）和无差异的基因数目，并计算无差异基因数目的差值。
1.4 基因集富集分析 未筛选的所有差异基因，用含GENE_SYMBOL的差异基因表，导入在线工具（www.chrislifescience.club:3838/R/AnnoE2）的GSEA中，列名中分别填写导入表格中对应各列的名称，分子标签数据库（MSigDB，MolecularSignatures Database） ID选择c7。标准化富集评分（NormalizedEnrichment Score，NES)降序排列，取NES前6位的基因集进一步分析。
1.5 免疫标记基因集的交集基因 NES前6位的ISGs条目，导入MSigDB（http://software.broadinstitute.org/gsea/msigdb/search.jsp），查询并下载基因，在线工具分析交集基因（http://bioinformatics.psb.ugent.be/cgi-bin/liste/Venn/），绘制韦恩图。
1.6 统计学处理 在线工具的R软件进行统计，富集评分标准化的值（NES）进行排列检验（permutationtest），错误发现率（Falsediscovery rate，FDR）方法判断NES的统计学意义；组间基因表达差异进行t检验；交集基因导入KM-plotter数据库（http://kmplot.com/analysis/），根据基因表达水平自动选择临界值，分割为高、低表达组，Kaplan-Meier法分析基因表达水平与胃癌患者总生存（Overallsurvival，OS）的关系，Log-rank检验比较组间差异，给出风险比（Hazardratio，HR）、95%的可信区间（Confidenceinterval，CI），P< 0.05为差异具有统计学意义。
2 结果
2.1 差异表达矩阵GEO 2 R分析共获取差异基因34183个，用芯片数据表对表达数据进行整理，转换基因名，删除无GENE_SYMBOL的行，获差异表达基因共29833个。COUNTIF函数计算LogFC > 0的上调基因共12 452个，LogFC < 0下调基因共17 381个。
2.2 差异基因表达 差异倍数的对数（Log FC）为1时，结果表格中，COUNTIF函数计算上调（UP）、下调（DOWN）、无差异（NOT）基因数目，分别为3819个、3730个、22284个；LogFC为1.5时，UP、DOWN、NOT基因数，分别为2300个、1492个、26041个。2种条件下的火山图见图1。2种条件下，UP、DOWN、NOT基因数的差值分别为1519个、2238个、3757个。上下调基因数目与筛选条件的关系见图2。
图1 差异表达基因的火山图
A：Log FC为1，P<0.05的火山图；B：Log FC为1.5，P<0.05的火山图
图2 上调和下调基因数与筛选条件的关系
2.3 富集分析 MSigDB中的c7作为富集分析的基因集，富集得到条目4872个，NES降序排序，排序前6的结果见表1，前6位的条目作GSEA图，见图3。
表1 基因集富集分析NES前6的条目
Code
ID
NES
P值
A
GSE9006_HEALTHY_VS_TYPE_2_DIABETES_PBMC_AT_DX_UP
2.684
0.004
B
GSE9006_TYPE_1_VS_TYPE_2_DIABETES_PBMC_AT_DX_UP
2.626
0.004
C
GSE9006_TYPE_1_DIABETES_AT_DX_VS_1MONTH_POST_DX_PBMC_UP
2.532
0.004
D
GSE42724_NAIVE_VS_B1_BCELL_DN
2.426
0.004
E
GSE6269_FLU_VS_E_COLI_INF_PBMC_UP
2.399
0.003
F
GSE9988_LPS_VS_VEHICLE_TREATED_MONOCYTE_DN
2.368
0.004
图3前6位免疫标记相关条目的基因集富集分析
A、B、C、D、E、F为表1中Code对应条目
2.4 免疫标记基因集的交集基因NES前6位的ISGs条目，导入MSigDB，查询、下载基因，在线工具分析交集，A、C、D基因集得到的交集基因1个，为MRPL12，见图4A。E、F基因集的交集基因2个：PRR13、PURA；A、E的交集基因5个：IDS、LIMD2、SOLH、CNOT3、TGFB1；A、F的交集基因4个：COTL1、PCBP1、DGKZ、DOK2，见图4B。
图4交集基因韦恩图
2.5 相关基因对胃癌总生存的影响 上述12个交集基因，导入KM-plotter，分析发现：MRPL12（Mitochondrialribosomal protein L 12）、PRR 13（Prolinerich 13）、COTL1（Coactosinlike F-actin binding protein 1）、PCBP 1（Poly(rC)binding protein 1）的高表达与胃癌OS延长关联（P<0.05），IDS（Iduronate2-sulfatase）、LIMD 2（LIMdomain containing 2）、PURA（Purinerich element binding protein A）、SOLH（Smalloptic lobes homolog）、CNOT 3（CCR4-NOT transcription complex subunit 3）、TGFB 1（Transforminggrowth factor beta 1）、DGKZ（Diacylglycerolkinase zeta）、DOK 2（Docking protein2）的高表达与胃癌OS降低关联（P<0.05），见表2，部分基因的生存曲线见图5。
基因名
基因中文名
HR
95% CI
P 值
MRPL 12
线粒体核糖体蛋白L 12
0.80
0.67-0.94
8.4 e-03
PRR 13
富含脯氨酸蛋白13
0.81
0.69-0.97
1.8 e-02
COTL 1
毛状蛋白样F-肌动蛋白结合蛋白1
0.70
0.59-0.83
4.3 e-05
PCBP 1
聚（RC）结合蛋白1
0.56
0.47-0.67
6.0 e-11
IDS
艾杜糖2-硫酸酯酶
1.36
1.15-1.62
4.1e-04
LIMD 2
LIM 结构域2
1.55
1.30-1.84
4.7e-07
PURA
富含嘌呤元素结合蛋白A
1.72
1.41-2.10
5.1 e-08
SOLH
小视觉叶同源物
2.27
1.84-2.80
2.9 e-15
CNOT 3
CCR4-非转录复合物亚单位3
2.41
1.95-2.98
1.0 e-16
TGFB 1
转化生长因子β 1
1.53
1.24-1.90
7.8 e-05
DGKZ
二酰甘油激酶ζ
2.18
1.80-2.65
4.4 e-16
DOK 2
接头蛋白2
1.28
1.07-1.52
6.5 e-03
注：HR：hazardratio；CI：confidenceintervals
A：MRPL 12高、低组的OS生存曲线；B：PRR 13高、低组的OS生存曲线；C：IDS高、低组的OS生存曲线；D：LIMD 2高、低组的OS生存曲线
3 讨论
胃癌的综合治疗可提高疗效，作为综合治疗重要组成的化疗，在胃癌手术前后的辅助治疗中发挥重要作用，明确化疗耐药机制、降低化疗抵抗、提高化疗敏感性，具有重要的临床意义。顺铂通过多种机制抑制胃癌细胞，其耐药性的产生可由PMP22介导，但其下游是否与ISGs存在关联，尚无文献报告。为探索PMP22下调后减弱胃癌顺铂抵抗是否通过下游的ISGs发挥作用，本研究对GEO数据库的GSE94 714数据集进行了GSEA分析，发现PMP22可通过多种ISGs相关的下游机制，对胃癌顺铂抵抗的表型进行调节。
本研究发现，根据差异倍数（Foldchange, FC）和统计P值进行筛选，PMP22基因敲减、未敲减组差异表达基因的数目减少，因此采用未筛选的差异基因，GSEA分析发现了NES前6位的ISGs条目，其组成基因用韦恩图方法取交集，发现了12个ISGs，均与胃癌的预后生存（OS）有关。
随着差异倍数FC或logFC的增大，统计分析无差异的基因数逐渐增加，被排除的基因增多，筛选出有差异的基因数目逐渐减少，即上下调基因的数目均减少，与Xiao等人进行条件筛选后、使组间差异表达基因数减少的研究结果一致[7]。预筛选可遗漏部分基因，造成后续分析无法找到有意义、但是改变细微的重要基因。GSEA方法可避免遗漏，较为全面地分析差异基因，通过预定义基因集，可获取有分类的基因标记信息。我们通过GSEA分析获得的12个基因，既与免疫有关，又与胃癌的预后有关，并且在前期进行的常规富集分析中未能被发现[2]，表明GSEA方法是常规富集分析方法的重要补充。
新发现的ISGs通过不同机制参与肿瘤细胞的化疗抵抗和发生发展，部分基因与胃癌化疗抵抗的关系尚未见文献报告。MRPL12基因可通过改变线粒体代谢、调节细胞的氧耗量，影响肿瘤生长，MRPL12敲减后，胰腺肿瘤细胞内的线粒体活性降低，糖酵解加强，肿瘤细胞体外生长减慢，但体内生长加速[8]。本研究发现PRR13高表达胃癌的OS延长，但Pontikakis等发现：卡铂、紫杉醇治疗后，PRR13高表达卵巢上皮性癌的OS缩短，与化疗抵抗有关[9]，表明PRR13表达对不同肿瘤生存时间的作用存在差异，其机制有待进一步阐明。COTL1基因可增强乳腺癌对化疗的敏感性，与白细胞介素24、非经典转化生长因子β1通路有关[10]。PCBP1可降低豆荚蛋白表达、增强胃癌腹膜转移细胞对多西他赛的敏感性[11]。LIMD2与甲状腺乳头状癌的转移有关[12]。TGFB1可能通过细胞粘附、转移相关分子通路影响胃癌预后[13]，与本研究结果相似。DGKZ作为关键基因，与非小细胞肺癌的进展和预后有关[14]。上调因启动子甲基化而沉默的DOK2表达，可增强卵巢癌对铂类药物的敏性[15]。IDS、PURA、SOLH、CNOT3与胃癌化疗抵抗的关系，目前尚未见文献报告，是进一步研究的方向。
本研究的局限性：新发现的ISGs与胃癌顺铂抵抗的关系尚需进一步的体内外实验验证，其作用方向和机制依然存在不确定性，初步结果仅为继续研究的基础和线索。
总之，用PMP22敲减致胃癌顺铂抵抗减弱与未敲减比较的全部差异表达谱数据，进行GSEA研究，首次发现了与顺铂抵抗有关的ISGs，作为潜在靶点，有望增强化疗敏感性，改善预后。
参考文献
[1] HuSB, Liu CH, Wang X, et al. Pathological evaluation of neoadjuvant chemotherapyin advanced gastric cancer[J]. World J Surg Oncol, 2019, 17(1): 3.
[2] 汪圣毅,张永红,闫亚飞,等.胃癌化疗抵抗基因PMP22下游的生物信息学机制[J].安徽医科大学学报,2019, 54(4): 509-514.
[3] CaiW, Chen G, Luo Q, et al. PMP22 Regulates self-renewal and chemoresistance ofgastric cancer cells[J]. Mol Cancer Ther, 2017, 16(6): 1187-98.
[4] KlapdorR, Wang S, Hacker U, et al. Improved Killing of Ovarian Cancer Stem Cells byCombining a Novel Chimeric Antigen Receptor-Based Immunotherapy andChemotherapy[J]. Hum Gene Ther, 2017, 28(10): 886-896.
[5] 王潇,尹天舒,李柏逸,等.基因功能富集分析的研究进展[J].中国科学:生命科学,2016, 46(04): 363-373.
[6] PowersRK, Goodspeed A, Pielke-Lombardo H, et al. GSEA-InContext: identifying noveland common patterns in expression experiments[J]. Bioinformatics. 2018, 34(13):i555–i564.
[7] XiaoY, Feng M, Ran H, et al. Identification of key differentially expressed genesassociated with nonsmall cell lung cancer by bioinformatics analyses[J]. MolMed Rep, 2018, 17(5): 6379-6386.
[8] ChenY, Cairns R, Papandreou I, et al. Oxygen consumption can regulate the growth oftumors, a new perspective on the Warburg effect[J]. PLoS One, 2009, 4(9):e7033.
[9] PontikakisS, Papadaki C, Tzardi M, et al. Predictive value of ATP7b, BRCA1, BRCA2, PARP1,UIMC1 (RAP80), HOXA9, DAXX, TXN (TRX1), THBS1 (TSP1) and PRR13 (TXR1) genes inpatients with epithelial ovarian cancer who received platinum-taxane first-linetherapy[J]. Pharmacogenomics J, 2017, 17(6): 506-514.
[10] XiaL, Xiao X, Liu WL, et al. Coactosin-like protein CLP/Cotl1 suppresses breastcancer growth through activation of IL-24/PERP and inhibition of non-canonicalTGFβ signaling[J]. Oncogene, 2017, 37(3): 323–331.
[11] JiFJ, Wu YY, An Z, et al. Expression of both poly r(C) binding protein 1 (PCBP1)and miRNA-3978 is suppressed in peritoneal gastric cancer metastasis[J]. SciRep, 2017, 7(1): 15488.
[12] PinheiroDos Santos MJC, Bastos AU, da Costa VR, et al. LIMD2 Is Overexpressed in BRAFV600E-Positive Papillary Thyroid Carcinomas and Matched Lymph NodeMetastases[J]. Endocr Pathol, 2018, 29(3): 222-230.
[13] 刘梦园,吴瑛,孙明军.公共数据库GEO分析TGFB1I1在胃癌中表达及临床意义[J].临床军医杂志,2016, 44(8): 862-865.
[14] FengA, Tu Z, Yin B. The effect of HMGB1 on the clinicopathological and prognosticfeatures of non-small cell lung cancer[J]. Oncotarget, 2016, 7(15):20507-20519.
[15] FangF, Cardenas H, Huang H, et al. Genomic and Epigenomic Signatures in OvarianCancer Associated with Resensitization to Platinum Drugs[J]. Cancer Res, 2018,78(3): 631-644.
