肝细胞癌基因组学免疫分型的鉴定及临床意义

张博超，浦春，2，胡蝶，张艳珍，朱萍

(1. 皖南医学院第一附属医院检验科，安徽 芜湖 241001；2. 皖南医学院检验学院，安徽 芜湖 241002)

肝癌是世界上第六大癌症，在病情发现时已到中晚期，因此死亡率较高[1]。据估计，2012年度死亡人数已达到70万[2]。肝细胞癌(HCC)占所有原发性肝癌的90%，由肝脏的慢性炎症引起，是一种典型的免疫原性癌症[3]。免疫逃逸在肝癌发生发展中起到重要作用[4]。肝癌肿瘤微环境包括众多免疫抑制细胞，如调节性T细胞、肿瘤相关巨噬细胞，可干扰免疫检测，进而导致肿瘤免疫逃逸[5-6]。近年来，免疫检查点抑制剂(ICI)疗法可以重新激活人体T淋巴细胞的相关调节信号，并恢复肿瘤患者的免疫系统，进而杀死肿瘤细胞，显著延长了患有不同实体瘤患者的预期寿命[7]。免疫检查点基因的高表达是预测ICI疗法的常用指标，如存在高T细胞浸润、干扰素-γ(IFN-γ)信号、检查点基因(如PD-1和PD-L1)表达或高肿瘤突变负荷(TML)可能有利于治疗反应[8-10]。但是，关于免疫亚型的分类及基于免疫检查点基因找出最佳疗效的证据却很少。

癌症基因组图谱(TCGA)是一个大型的癌症研究参考数据库，收集整理了癌症相关的各种组学数据。本研究通过TCGA-HCC的表达谱获取预后相关基因，并将HCC样本分为两个亚型，随后对来自两个外部验证集数据库的肝细胞癌样本的数据集进行了独立验证。结合肿瘤突变负荷、药物敏感性和免疫治疗等多组学数据，进一步研究每个免疫表型的潜在免疫逃逸机制，以及哪些基因组改变可能导致这些不同表型的形成。研究发现这两种免疫亚型具有不同的细胞、分子和临床特征，这可能为肿瘤学家提供有价值的信息。

1 材料和方法

1.1 数据获取 TCGA数据库(https：//portal.gdc.cancer.gov/)作为目前最大的癌症基因信息数据库，保存包括基因表达数据、miRNA表达数据、拷贝数变异、DNA甲基化、SNP等。本课题组下载了已处理的HCC原始的mRNA表达数据，包括正常组(n=50)和肿瘤组(n=373)。癌症基因组学cBioPortal(cBioPortal，http：//www.cbioportal.org，v3.2.11版)是一个开放访问的在线工具，集成了来自大规模基因组项目的原始数据，包括但不限于TCGA和国际癌症基因组联盟(ICGC)。从NCBI GEO 公共数据库下载GSE14520的Series Matrix File 数据文件，注释平台为GPL3921，下载221例带有完整表达谱及生存信息的HCC患者数据。从ICGC数据库下载了202例带有完整表达谱及生存信息的LIHC患者数据，用于本次分析的免疫基因集(共包含1 811个免疫相关基因)通过immport数据库获取。

1.2 免疫亚型的分类 对免疫相关基因的表达谱通过NMF包执行无监督NMF聚类，使用R软件包“survival”进行Cox回归分析，评估所有候选基因与总生存率(OS)的相关性。使用NMF包执行无监督NMF聚类方法，并使用相同的候选基因将该方法应用于两个外部验证集。选择相关系数开始下降的k值作为最佳聚类数。基于T-SNE的方法，利用上述免疫基因的mRNA表达数据验证亚型分配。

1.3 药物敏感性分析 基于最大的药物基因组学数据库(GDSC癌症药物敏感性基因组学数据库，https：//www.cancerrxgene.org/)，使用R软件包“pRRophetic”来预测每个肿瘤样本的化疗敏感性；用回归的方法得到每种特定化疗药物治疗的IC50估计值，并用GDSC训练集进行10次交叉验证检验回归和预测精度。所有参数都选择了默认值，包括去除批处理效应的“combat”以及取重复基因表达的平均值。

1.4 免疫细胞浸润分析 使用CIBERSORT算法对不同亚组的HCC患者RNA-seq数据进行分析，用来推断22种免疫浸润细胞的相对比例，并对基因表达量以及免疫细胞含量进行Spearman 相关性分析，P<0.05被认为有统计学差异。

1.5 基因集差异分析 基因集变异分析(GSVA)是评估转录组基因集富集情况的一种非参数无监督方法。GSVA通过对感兴趣的基因集合进行综合打分，将基因水平变化转变为通路水平变化，进而判断样本的生物学功能。本研究将从Molecular signatures database数据库(v7.0版)下载基因集合，采用GSVA算法对每个基因集合进行综合打分，评估不同样本潜在的生物学功能变化。

1.6 WGCNA分析 通过构建加权基因共表达网络，寻找协同表达的基因模块，并探索基因网络与表型之间的关联关系，以及网络中的核心基因。利用WGCNA-R包分别构建数据集中所有基因的共表达网络，筛选方差前5 000的基因用该算法进行筛选，以便进一步分析，其中软阈值为3。将加权邻接矩阵转化为拓扑重叠度矩阵(TOM)来估计网络连接度，运用层次聚类的方法来构建TOM矩阵的聚类树结构。聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块，将几万个基因通过基因表达模式被分成了多个模块。

1.7 基因本体论(GO)和京都基因与基因组百科全书(KEGG)功能注析 使用clusterProfiler (R3.6)对关键基因进行功能注释，以全面探讨这些候选基因的功能相关性。GO和KEGG被用来评估相关的功能类别。P值和q值均<0.05的GO和KEGG富集通路被认为是显著性类别。

1.8 亚型GSEA分析 通过R包GSVA对KEGG通过做GSEA分析，最后挑选在亚型间显著性最高的前50个结果用于展示。其中基因集来自MSigDB数据库(http：//www.gsea-msigdb.org/gsea/downloads.jsp)。

1.9 肿瘤免疫评估(Timer) tumor immune estimation resource(Timer，https：//cistrome.shinyapps.io/Timer/)是一个全面的资源，用于系统分析不同癌症类型的免疫浸润。在这项研究中，使用Timer观察抗原提呈细胞(APC，antigen-presenting cell)浸润与已鉴定的有效抗原表达之间的关系。

1.10 统计学方法 采用Cox比例风险模型进行单变量和多变量分析。所有统计分析均采用R语言(version 3.6)进行。所有统计检验均为双侧，P<0.05具有统计学意义。

2 结果

2.1 获取预后相关基因及免疫亚型分析 收集HCC患者的临床信息，经Cox单因素回归(P<0.05)分析共筛选出267个预后相关基因，通过NMF共识聚类法根据上述267个候选基因的表达谱对包含肝癌样本的TCGA数据集进行聚类，综合考虑后最终选择k=2作为最佳聚类数(见图1A)。当k=2时，通过t-SNE来减小特征的维数，发现亚型间在很大程度上与二维t-SNE分布模式一致(见图1B)。在TCGA数据集中观察到显着的预后差异，与C1相比，C2具有更佳的生存可能(见图1C)。此外，两个验证集的亚型差异同样也观察到相似的差异，C1的OS时间明显短于C2亚型(见图1D～1E)。分析了亚型的病理特征，发现亚型与HCC患者的肿瘤分级和分期密切相关(P<0.05)，见表1。

图1 免疫亚型

表1 肝癌患者亚型的临床病理特征

表1(续) 肝癌患者亚型的临床病理特征

2.2 亚型的免疫特征 分析发现B cells memory、T cells CD4 memory activated、T cells CD4 memory resting、Macrophages M1、Neutrophils等肿瘤微环境因子在亚型间存在显著性(见图2A)，C2亚型免疫细胞浸润程度高于C1亚型。免疫检查点基因CTLA-4和PD-L1及人类白细胞抗原的表达在亚型间存在显著性(见图2B～图2D)，C1亚型的表达水平明显高于C2亚型，提示两个亚型之间的免疫调节通路的异常扰动是两组患者预后差异的潜在机制。

图2 免疫特征

2.3 免疫亚型与常见化疗药物及免疫治疗的相关性 早期HCC手术结合化疗治疗效果明确。研究基于GDSC数据库的药物敏感性数据，通过R包“pRRophetic”来预测每个肿瘤样本的化疗敏感性，进一步探讨免疫亚型与常见化疗药物的敏感性。研究结果表明，免疫亚型与患者对Bleomycin、Bosutinib、Cisplatin、Cytarabine、Docetaxel、Doxorubicin的敏感性显著相关(见图3A)。同时本研究进一步基于黑色素瘤免疫治疗的数据集，预测两个亚组对抗肿瘤免疫治疗的敏感性，结果显示， C1 亚组对免疫治疗更为敏感(见图3B)。

图3 免疫亚型与常见化疗药物及免疫治疗的相关性

2.4 免疫亚型患者的突变图谱 结果表明，高风险组患者TP53、CTNNB1、TTN等多个基因的突变比例在亚型间显著不同，其中C1亚型患者的TP53突变比例高于C2亚型，见图4。

图4 亚型突变图谱

2.5 HCC免疫基因共表达模块 为了确定HCC队列中免疫相关基因的共表达网络，本课题组进行了WGCNA分析。以免疫亚型C1和C2作为样本的临床性状，进一步用于构建WGCNA网络，探讨LIHC中的biomarkers。软阈值β由函数“sft$powerEstimate”确定，软阈值设置为3。然后基于tom矩阵检测基因模块，在本次分析中共检测到3个基因模块，分别为blue(n=427)、turquoise(n=568)、grey(n=66)模块。进一步通过模块与性状之间的分析，发现MEblue模块与样本类别(免疫亚型)的相关性最高(cor=-0.59，p=6e-36) (见图5A、图5B)，因此，将选择MEblue模块用于后续的相关验证分析。对两种免疫亚型的模块特征基因进行了分析，除grey外的两个模块中特征基因的表达在亚型间均具有显著性(见图5C、图5D)。

2.6 免疫基因共表达模块的功能富集以及蛋白互作网络构建 通过GO和KEGG富集分析发现MEblue模块的基因在大量的通路中均显著富集(见图6A、图6B)。比如在GO富集中有regulation of innate immune response、proteasome accessory complex、growth factor binding等通路均有大量基因被富集(见图7A)；比如在KEGG富集中有Cytokine-cytokine receptor interaction、T cell receptor signaling pathway、Lipid and atherosclerosis等通路均有大量基因被富集(见图6B)。接下来研究与MEblue模块的基因可能相关的具体信号通路，探讨可能影响HCC病变及进展的潜在分子机制。GSVA结果表明，在免疫亚型组间的差异通路主要富集到了FATTY_ACID_METABOLISM、MYC_TARGETS_V1、NOTCH_SIGNALING、COAGULATION、EPITHELIAL_MESENCHYMAL_TRANSITION等信号通路(见图6C)。此外，通过预后相关分析显示，其中blue和turquoise模块中的基因表达与HCC患者的预后显著相关(见图7A)。模块MEblue和MEturquoise中的基因表达与PCA1均显著相关(见图7B、图7C)；模块MEblue中的基因表达与PCA2的相关性最高(见图7D、图7E)。

图6 免疫基因共表达模块的功能富集

图7 免疫基因共表达模块的主成分分析

2.7 探讨在LIHC 队列的表达谱中免疫关键基因 通过在线工具cBioPortal获得了TCGA-HCC中共有20 887个基因发生了拷贝数变异(见图8A)，有13 640个基因发生了突变(见图8B)。有129个基因的OS生存分析具有高显著性(P<0.0001)，有169个基因的RFS生存分析具有高显著性(P<0.0001)(见图8C)。对前10位核心基因进行单因素(见图8D)和多因素(见图8E)分析，单因素联合多因素分析筛选出EIF5B作为肝癌的免疫关键基因。对EIF5B进行预后分析和TIMER分析，发现EIF5B基因高表达与患者不良预后相关，并与免疫细胞高浸润水平相关(见图8F、图8G)，P<0.05。

图8 免疫关键基因

3 讨论

肿瘤的发生发展与宿主的免疫系统联系密切，免疫系统可以杀死和清除肿瘤细胞，然而肿瘤细胞可以通过免疫逃逸躲避免疫系统的监视，进而诱发癌症。目前的免疫治疗主要克服肿瘤的免疫耐受及加强效应细胞的杀伤效应。既往研究[11-14]发现，在头颈部癌、结肠癌、三阴性乳腺癌中通过对免疫特征进行相应评估从而确定的免疫分型，有助于指导治疗及改善患者的预后。本研究基于从TCGA、GEO和ICGC数据库筛选的267个关键预后相关基因并建立了HCC分类，鉴定出HCC的两个亚型(C1和C2)，结果表明，C1表现出较差的预后特征。C1和C2都与免疫信号相关，C1具有免疫检查点基因的高表达。众所周知，固有免疫逃逸有两个主要方面，包括免疫原性和免疫检查点分子的表达[15]。之前的研究报告里，在HCC中，免疫检查点基因通常被激活，并伴随着高免疫应答和基因表达上调[16]。研究显示，与C2亚型相比，C1亚型中的检查点基因CTLA-4和PD-L1及人类白细胞抗原显著过度表达。这些结果表明，C1亚型可能与HCC的固有免疫逃逸密切相关，这可能导致患者的预后较差，并为使用免疫检查点阻断剂进行HCC免疫治疗提供新的见解。此外，本课题组分析了亚型的病理特征，发现亚型与HCC患者的肿瘤分级和分期密切相关。总的说本研究鉴定了HCC的两种免疫亚型。

通过药敏数据库分析发现两种免疫亚型与患者对博来霉素、博苏替尼、顺铂、阿糖孢苷、多烯紫杉醇、阿霉素的敏感性显著相关，并且C1亚组对免疫治疗更为敏感，这可能为HCC患者的治疗提供思路。通过突变数据库分析，验证了HCC亚型中TP53、CTNNB1和TTN基因的突变，其中C1亚型患者的TP53突变比例高于C2亚型。据报道TP53突变与免疫检查点的高表达、活跃的IFN-γ信号和效应T细胞信号有关[17]。本研究证明C1亚型患者具有较高的TP53突变率，进一步验证了该亚型在基于免疫检查点的治疗中的预测作用。

本研究还利用WGCNA分析探索了HCC免疫微环境中两种亚型的潜在靶点和途径。免疫相关基因富集在3个不同的模块中，结果显示固有免疫、T细胞信号调节通路和T细胞受体信号通路等免疫相关途径主要富集在蓝色模块中。在所有成对比较模块中，蓝色模块代表最常见的功能富集路径，这表明这个模块中的基因可能在HCC免疫微环境中发挥一定的作用。与目前的结果一致，以前的研究报告显示肿瘤免疫微环境在调节HCC发生、肿瘤侵袭和转移过程中起着关键作用[18]。固有免疫在调节HCC的发生和发展中起着关键作用[19]。也有人假设在HCC进展过程中可能通过T细胞信号通路诱导适应性免疫反应[20]。此外，在蓝色模块还显示了MYC突变、NOTCH信号通路、上皮间充质转化等癌变相关信号通路，进一步探明HCC发生发展的分子机制。因此，阐明HCC免疫微环境的分子机制可能有助于确定HCC新的治疗靶点。此外，本课题组还筛选出EIF5B作为HCC的免疫关键基因，可为HCC的mRNA疫苗研制提供思路。

目前的研究存在一些局限性。首先，为了全面反映影响HCC微环境表型的因素和影响，亚组分析应包括更多HCC患者的临床特征。其次，每个亚型的样本量都相对较小，只有IGCG和GEO队列用于外部验证，这可能导致片面结果和高假阳性率。在内部验证中进行交叉验证，并增加外部验证的样本量，将有助于HCC免疫微环境的未来研究。