利用TCGA数据库筛选食管鳞状细胞癌免疫相关mRNA

孟玉娟， 亓 民， 马 军， 袁 翎， 任凯凯， 马佳康， 侯明宇， 张中冕

1.郑州大学附属洛阳中心医院医学检验中心，河南 洛阳 471009； 2.郑州大学第二附属医院 郑州大学消化疾病研究所； 3.郑州大学附属肿瘤医院放疗科； 4.郑州大学第二附属医院肿瘤科

食管鳞状细胞癌(esophageal squamous cell carcinoma，ESCC)是我国常见恶性肿瘤之一，尽管内镜技术的发展提高了诊疗效果，但根据2019年中国肿瘤登记年报统计，其发病率排第6位，死亡率排第4位[1]。

免疫学的发展揭示了机体免疫功能平衡和肿瘤发生的相关性，免疫相关基因对恶性肿瘤的预后有预测作用，一些调节抗肿瘤免疫功能的药物如PD-1、CTLA4单抗及免疫细胞治疗技术，已在临床上取得了较好的效果[2]。利用免疫细胞浸润分析方法发现，ESCC中9个免疫相关基因被认为是预测预后的因素，与ESCC患者的总生存率密切相关[3]，低级别胶质瘤患者中发现8个预后关键基因[4]；免疫细胞浸润的程度和分布也与乳腺癌、前列腺癌和结直肠癌的总生存率相关[5-7]。然而，关于ESCC免疫浸润之间关系的研究还很有限。本研究利用TCGA数据库中的ESCC测序数据，分析免疫细胞浸润和免疫相关基因，筛选用于预测ESCC患者预后的免疫相关基因。

1 方法

1.1 数据下载从TCGA数据库(http://tcga-data.nci.nih.gov/tcga)中下载ESCC转录组FPKM(fragments per kilobase of transcript per million mapped reads)、Counts数据以及临床数据。

1.2 单样本基因集富集分析(single sample GSEA，ssGSEA)利用R软件中“GSVA”包对ESCC的FPKM数据进行处理得到每个样品中29种免疫相关基因集(免疫基因集包括免疫细胞类型、功能和通路)评分，根据ssGSEA的结果，使用“hclust”(R包)将ESCC肿瘤样本进行免疫浸润分组(分为低、高和中等免疫组)。

1.3 分组效果的验证基于ESCC的FPKM转录组表达数据，通过ESTIMATE算法得到每个样本的基质、免疫、ESTIMATE评分和肿瘤纯度，验证ssGSEA分组的效果。同时分析各免疫浸润组中人白细胞相关抗原(HLA)和程序性死亡配体1(PD-L1)的基因表达水平用于分组效果再验证。

1.4 加权共表达网络分析(weighted correlation network analysis，WGCNA)根据|log2FC|>1和错误发现率(FDR)<0.05的标准，利用R软件“edger”软件包对ESCC的mRNA Counts数据进行mRNA差异分析。结合免疫浸润分组，对差异表达mRNA，利用R软件“WGCNA”软件包进行加权共表达网络分析筛选免疫最相关模块，根据性状与模块基因的相关性(gene significance，GS)>0.3以及基因与模块的相关性(module membership，MM)>0.8筛选模块中的免疫相关的hub基因以进行进一步分析。

1.5 生存分析利用ESCC患者的生存时间以及状态数据结合hub基因的表达情况，通过R软件“survminer”包中的surv_cutpoint函数选择最佳的cut-off值并进行Kaplan-Meier生存分析以比较hub基因表达的总生存期(OS)差异，筛选出与生存相关的hub基因。

1.6 GO(Gene Ontology)富集和KEGG(Kyoto Encyclopedia of Genes and Genomes)分析对最相关模块中的基因利用R软件“Cluster Profiler”包进行GO及KEGG富集分析，筛选出模块相关的功能及通路。

2 结果

2.1 ESCC样本的免疫分组及验证使用ssGSEA算法分析了TCGA ESCC样本在29个免疫浸润数据集中的评分，用层次聚类的方法将81例肿瘤样本分为高免疫组(14例)、中等免疫组(46例)及低免疫组(21例)(见图1)。为了验证上述分组策略的可行性，基于ESCC的FPKM表达谱数据使用ESTIMATE算法计算基质、免疫、ESTIMATE评分和肿瘤纯度。随着样本免疫性的增加，肿瘤纯度降低，基质、免疫和ESTIMATE评分逐渐升高(见图2A)。4项评分在3组中差异均有统计学意义(见图2B)。此外，ESCC中HLA相关基因和PD-L1基因表达水平在高免疫组与中等免疫组中差异有统计学意义(见图2C～D)。

2.2 免疫相关hub基因的筛选利用ESCC转录组mRNA Counts数据筛选差异mRNA，下载的数据包括ESCC样本81例及癌旁样本11例(非完全为ESCC配对癌旁样本，10例为食管腺癌配对样本)。共筛选到4 867个差异mRNAs，根据81例样本的免疫浸润分组，结合差异基因表达量行WGCNA分析，将基因分为17个模块，其中与免疫分组最相关的为brown模块(cor=-0.73，P=5e-14)(见图3A)。根据MM>0.8及GS>0.3的标准对模块中的基因进行过滤，共筛选到6个hub基因，分别为HAVCR2、CCL4、FPR3、SLAMF8、SAMSN1、SLC15A3(见图3B)。

2.3 hub基因的生存分析基于ESCC临床数据以及6个hub基因的表达量，通过R软件“survminer”包中的surv_cutpoint函数选择最佳的cut-off值，发现HAVCR2与患者OS相关，HAVCR2的表达量越高，患者的OS越短(见图4)。

注：聚类1(绿色)、聚类2(红色)和聚类3(蓝色)分别代表高免疫组、中等免疫组和低免疫组。

注：Immunity-L:低免疫组；Immunity-M: 中等免疫组；Immunity-H:高免疫组。* *P<0.01，* * *P<0.001，ns：差异无统计学意义。

图2 分组效果的验证A：热图显示使用ESTIMATE算法得到的基质、免疫、ESTIMATE评分和肿瘤纯度在免疫浸润分组之间的比较；B：小提琴图表示其在免疫浸润分组之间的统计学差异；C～D：箱形图显示了免疫浸润分组中HLA和PD-L1的差异表达

Fig 2 Verification of grouping effectA: the thermogram showed the comparison of matrix, immunity, ESTIMATE score and tumor purity between immune infiltration groups using the ESTIMATE algorithm; B: Violin chart showed the statistical difference between groups with immune infiltration; C-D: the box diagram showed the differential expressions of HLA and PD-L1 in the immune infiltration groups

图3 hub基因的筛选 A：不同模块与免疫表型之间相关性的热图；B：brown模块中GS和MM相关性分析散点图Fig 3 Screening of hub gene A: thermogram of correlation between different modules and immunophenotypes; B: scatter diagram of GS and MM correlation analysis in brown module

图4 HAVCR2的生存分析图Fig 4 Survival analysis chart of HAVCR2

2.4 GO富集和KEGG分析对最相关模块中的基因进行富集分析，GO富集结果显示，该模块与趋化因子受体结合、趋化因子活性、G蛋白偶联受体结合、免疫球蛋白结合、细胞因子受体结合、CCR趋化因子受体结合、IgG结合、受体配体活性、CXCR趋化因子受体结合、细胞因子活性、糖胺聚糖结合、碳水化合物结合等功能相关(见图5A)。KEGG富集结果提示，该模块中基因可能参与病毒蛋白与细胞因子和细胞因子受体的相互作用、细胞因子与细胞因子受体的相互作用、Toll样受体信号通路等(见图5B)。

图5 Brown模块的GO富集结果及KEGG富集结果 A：GO富集结果；B：KEGG富集结果Fig 5 GO enrichment results and KEGG results of brown module A: GO enrichment results; B: KEGG enrichment results

3 讨论

早期的研究已经发现，肿瘤灶浸润淋巴细胞的多少与预后有关，随着肿瘤免疫相关机制的深入研究，一些靶向免疫机制的药物已经在临床应用，并取得了一定的疗效，但是临床上期待发现更多特异性诊疗靶点。虽然ESCC的诊断和治疗技术已经在全球范围内得到了改进，但还未发现可靠的免疫治疗方案，越来越多的研究报道免疫浸润与ESCC的预后密切相关[7-12]。

文献报道，利用基因芯片技术，发现9个免疫相关基因与ESCC患者的总生存率相关[3]。但是，利用生物信息学技术可以更有效地利用现有的数据库信息，是探索肿瘤发病机制的常用方法之一，在肿瘤免疫浸润和免疫分型方面的应用逐渐增多。在乳腺癌中，基于ssGSEA和WGCNA分析，发现4个免疫相关基因被鉴定为预后的生物标志物[13]。本研究中，利用ssGSEA分析，构建ESCC免疫分型，将TCGA数据库中的ESCC标本分为低、中、高免疫浸润组。并发现免疫相关基因HAVCR2、CCL4、FPR3、SLAMF8、SAMSN1、SLC15A3可能与ESCC免疫调控有关。进一步的预后分析，从上述基因中找到和预后相关的HAVCR2，该基因可能通过细胞因子相互作用，参与ESCC的发生和发展。

HAVCR2基因，也叫TIM3，是TIM家族中具有共同基序的Ⅰ型跨膜糖蛋白，参与机体天然免疫和适应性免疫反应，在多种实体肿瘤细胞表面表达。TIM3可抑制免疫细胞功能，诱导T细胞、NK细胞的凋亡，参与肿瘤免疫抑制。还可以可溶性方式分泌到细胞外[14]。国内文献也已经报道，TIM3基因的过表达与ESCC的预后呈负相关[15-16]。

随着计算机技术的飞速发展，多种生物信息学方法被用来筛选预测预后的指标。筛选出的6个免疫相关基因可能在ESCC发生、发展中起重要作用，尤其是TIM3基因。针对这些基因的进一步研究，探明其在ESCC中的重要作用，可能会为ESCC增添新的诊疗靶点。