Identification of epigenetic dysregulation gene markers and immune landscape in kidney renal clear cell carcinoma by comprehensive genomic analysis

通过综合基因组分析识别肾脏透明细胞癌的表观遗传失调基因标志物和免疫景观

发表期刊：Front Immunol

发表日期：2022 Aug 18

影响因子：8.786

DOI:  10.3389/fimmu.2022.901662

一、研究背景

        肾癌是泌尿系统中最致命的癌症，而且近年来发病率呈上升趋势。由于肾癌缺乏特殊的临床表现，大约20-25%的患者在诊断时已经发生了远处转移。

        肿瘤发生和发展的本质是抑癌基因的失活和促癌基因的激活。研究表明，表观遗传学变化可以调节基因表达。常见的表观遗传学修饰包括DNA甲基化和乙酰化、组蛋白甲基化和乙酰化。特别是组蛋白修饰，通常指的是甲基化和乙酰化，在基因的异常表达中起着重要作用。目前的研究表明，组蛋白甲基化异常是肾癌的独立预后标志，也是肾癌的潜在临床目标。

二、材料与方法

1、数据来源

1） 从TCGA数据库下载了肾脏透明细胞癌（KIRC）的基因表达谱和表达谱数据（n=526）

2） 从国际癌症基因组联盟（ICGC）数据库中下载了带有生存时间的RECA-EU数据集（n=91）

3） 从GEO数据库下载了GSE86091数据集的hg19版本，包括癌旁样本和肿瘤样本，该数据集包含三种组蛋白信息，即H3K4me1、H3K4me3和H3K27ac

4） HK2细胞系（正常人肾小管上皮细胞系）和所有四种人类RCC细胞系（786-O、A498、Caki-1和ACHN）

2、分析流程

流程图

三、实验结果

01 - 识别具有表观遗传失调的PCG

        作者使用limma来识别显著差异表达的基因（共2755个PCG），结合组蛋白修饰数据和甲基化数据，最终发现了872个表观PCG（Epi-PCG）和18629个非表观PCG（non-epi-PCG）。Epi-PCG仅占所有PCG的4.47%。

        对表观PCG和非表观PCG的基因外显子和转录物的数量和长度进行了比较，以显示表观遗传失调的PCG的基因组特征。epi-PCG转录本的数量多于非epi-PCG，而epi-PCG转录本的长度则短于非epi-PCG（图2A，B）。同时，与非epi-PCG相比，epi-PCG有更多的外显子和更长的外显子长度（图2C，D）。此外，还系统地分析了KIRC中的epi-PCG，并揭示了以不同组蛋白修饰和不同甲基化区域为特征的epi-PCG景观。数据显示，大多数表观遗传失调的PCG都伴随着各种组蛋白修饰的异常，这些异常的组蛋白修饰主要集中在启动子区域（图2F）。

图2    基因组特征比较

        为了确定组蛋白修饰异常引起的PCG失调的潜在功能，作者系统地分析了KIRC中epi-PCG的表达与通路之间的关系。提取了不同组蛋白修饰引起的PCG的表达谱，并使用ssGSEA计算每个样本在这些PCG中的富集分数。结果发现，6种组蛋白失调的GSEA得分在肿瘤样本中明显高于正常样本，说明这些组蛋白失调有促癌作用（图3A）。此外，还评估了每个样本的KEGG Pathway得分，并分析了每种类型的表观遗传因子的富集得分与KEGG Pathway之间的关系，得到每种类型表观遗传因子的相关KEGG Pathway。图3B显示了与6种类型的epi-PCG最相关的KEGG Pathway，共24条途径。结果表明，不同类型的epi-PCG相关途径具有一定的一致性。这些数据表明，epi-PCG与肿瘤的发生、发展和免疫力密切相关。

图3 外显性PCG的功能富集分析

        作者还分析了6种表观基因组与m6A和m5C RNA修饰之间的关系。从TCGA中提取了KIRC中m6A、m5C和m1A的表达谱，并分析了6种表观遗传因子与m6A、m5C和m1A的富集分数之间的相关性（图4A）。发现这些富集分数与m6A、m5C和m1A显著相关。进一步使用R软件包clusterProfiler对epi-PCG进行KEGG通路分析和GO功能富集分析。在基因的GO功能注释方面，注释了519个差异显著的BPs；注释了79个差异显著的CCs；注释了48个差异显著的MFs；KEGG通路富集分析被注释为38条重要的通路（图4B -E）。

图4 Epi-PCG和RNA修饰的相关性和功能富集分析

02 - 基于epi-PCG的3种预后差异的分子亚型的鉴定

        在TCGA和ICGC的数据集中，对epi-PCG进行了单因素生存分析，两个数据集中的生存相关基因被选为分子亚型的聚类基因，最后51个交叉基因被纳入其中（图5A）。对这51个基因在正常和肿瘤样本中的表达差异分析表明，这些基因的表达有明显的差异。基于51个Epi-PCG，用ConsensusClusterPlus对两个数据集进行聚类，得到三种Epi-PCG相关的亚型（图5B，C）。KM分析表明，在TCGA数据集中，C2的预后较差，而C1的预后较好（图5D）。在ICGC的数据集中也观察到类似的结果（图5E）。

图5 识别与epi-PCG相关的分子亚型

        作者分析了三种分子亚型之间的趋化因子是否存在表达差异。在TCGA数据集中，41个趋化因子中有26个在不同亚型中表现出明显的表达差异（图6A），这表明不同亚型的免疫细胞浸润程度可能不同，这些差异可能导致肿瘤进展和免疫治疗效果的差异。此外，18个趋化因子受体基因中的17个在三种分子亚型的表达中存在明显差异（图6B）。使用ssGSEA计算每位患者的IFNγ评分、CYT评分和血管生成评分。可以看出，各亚组的IFNγ评分存在明显差异（图6C）。其中，C1和C3的免疫T细胞细胞溶解活性最高，而C2的细胞溶解活性最低（图6D）。C2的血管生成得分最低（图6E）。在对47个免疫检查点相关基因的相关分析中，43个基因在三个亚组中有明显的差异（图6F）。这些结果表明，不同亚组对免疫治疗的反应可能不同。

图6 不同亚组中的免疫差异

        在TCGA数据集中，使用CIBERSORT方法评估每个样本中22个免疫细胞的得分，并观察这些免疫细胞得分在三个亚组中的分布，如图7A所示。16个免疫细胞在不同亚型中表现出明显的差异（图7B）。用ssGSEA的方法计算了28个免疫细胞的分数，然后比较它们在亚型中的差异，发现28个免疫浸润分数在亚型中有显著差异（图7C）。

图7 评估不同分子亚型的免疫和路径特征

        作者还分析了不同分子亚型对免疫治疗和化疗反应的差异。TIDE软件被用来评估免疫疗法对我们定义的分子亚型的潜在临床效果。TIDE预测得分越高，表明免疫逃逸的可能性越大，这表明患者从免疫治疗中获益的可能性越小。如补充图2所示 ，在TCGA的数据集中，C2的TIDE得分最低（补充图2A）。同时，还比较了不同分子亚型的预测T细胞功能障碍评分(补充图2B )和T细胞排斥评分(补充图2C )，不同组别之间也有明显差异。

补充图2 分子亚型的TIDE评分、T细胞功能紊乱评分和T细胞排斥评分的比较

03 - 建立基于epi-PCG相关基因的预后风险模型

        最终训练集数据共316个样本，测试集数据共210个样本。利用训练集数据，对每一个表观PCG进行单变量Cox分析，最后纳入46个预后基因。使用R软件包glmnet进行lasso cox回归分析。首先分析了各自变量的变化轨迹，如图8A所示。可以看出，随着lambda的逐渐增大，自变量系数接近0的数量也逐渐增多。采用10倍交叉验证法建立模型，并确定每个lambda下的置信区间，如图8B所示。可以看出，当lambda=0.0316时，该模型是最佳的。因此，当lambda=0.0316的10个基因被认为是进一步分析的目标基因。为了减少基因的数量，使用了MASS软件包中的stepAIC方法，最终将10个基因减少到8个基因。最终的RiskScore公式如下：RiskScore=0.28*ETV4+0.631*SH2B3-0.338*FATE1+0.363*GRK5-0.42*MALL-0.196*HRH2-0.354*SEMA3G-0.431*SLC10A6。

        根据样品的表达水平，计算出每个样品的RiskScore，样品的RiskScore如图8C所示。此外，使用R软件包timeROC分析RiskScore的预后分类ROC，并分别确定了1年、3年和5年预后的分类效率。如图8D所示，该模型具有较高的AUC区域。最后，对Riskscore进行了z-score，将Riskscore大于0的样本分为高风险组，而Riskscore小于0的样本为低风险组，并绘制KM曲线，如图8E所示。

图8 在TCGA数据集中构建一个预后模型

        为了评估模型的稳健性，根据样本的表达水平，使用与训练集相同的模型和相同的系数计算TCGA验证集、TCGA整个数据集和ICGC数据集中每个样本的RiskScore。应用R软件包timeROC来分析TCGA验证集的RiskScore的预后分类。1年、3年和5年的ROC效率分别为0.73、0.69和0.63（补充图3A）。最后对Riskscore进行了zscore。将Riskscore大于0的样本分为高危组，而低于0的样本为低危组，并绘制KM曲线。结果显示，高危组患者的预后明显差于低危组（补充图3B）。具体来说，99个样本被归类为高风险，111个样本被归类为低风险。在所有TCGA数据集中，1年、3年和5年的ROC效率分别为0.77、0.73和0.70（补充图3C）。高危组患者的预后明显比低危组差（补充图3D）。这里，242个样本被归类为高风险组，284个样本被归类为低风险组。使用独立的验证集ICGC来验证该模型的适用性。采用TimeROC来评估RiskScore对ICGC的预后分类。1年、3年和5年的ROC效率分别为0.77、0.73和0.70（补充图3E）。然后对Riskscore进行Z-score，将Riskscore大于0的样本分为高风险组，而低于0的样本为低风险组，并绘制KM曲线。结果表明，高危组患者的预后明显差于低危组（补充图3F ）。其中，49个样本被列为高危组，42个样本被列为低危组。

补充图3

        通过Riskscore将临床亚组特征分为高风险组和低风险组。结果表明，Riskscore可以显著区分年龄、性别、TMN分期和分级两组的预后差异（图9A-M）。此外，RiskScore与临床亚组特征的相关性比较也显示Riskscore在T期、N期、M期、阶段、等级和性别方面有显著差异（图9N -S）。

图9 基于Riskscore的临床亚组生存分析和相关性分析

        计算了每个样本在不同功能上的ssGSEA得分，并进一步分析了这些功能与RiskScore之间的相关性。挑选出相关性大于0.4的功能，以发现一个功能与RiskScore呈正相关，而其余22个功能与RiskScore呈负相关。挑选出最相关的23条KEGG途径，并根据其富集分数进行聚类分析。在这23条通路中，例如P53信号通路，随着RiskScore的增加而增加，而代谢通路如脂肪酸代谢、糖酵解葡萄糖生成、半乳糖代谢随着RiskScore的增加而减少。此外，利用STRING在线工具表征了这8个预后基因之间的蛋白-蛋白相互作用（PPI）。结果显示，SEMA3G、ETV4和SH2B3有密切的相互作用，GRK5和HRH2有密切的相互作用，表明它们可能对影响KIRC的预后有协同作用。

        此外，评估了预后基因的表达与免疫细胞浸润之间是否存在相关性。Pearson相关分析显示，M0、M1和M2巨噬细胞以及调节性T细胞的富集与预后基因明显相关；特别是ETV4和活化的CD4记忆T细胞之间观察到相对较强的相关性；SEMA3G、SLC10A6和SH2B3的表达与调节性T细胞明显相关。此外，发现八个预后基因的表达在三个分子亚型中的分布。总生存期最差的C2亚型在TCGA和ICGC数据集中所有8个基因的表达量都是最低的，这与之前的结果一致，即高危组的这些基因的表达量相对较低(图8C) 。

04 - 8个基因签名是KIRC的独立预后危险因素

        为了验证8基因签名模型在临床应用中的独立性，作者对TCGA数据集进行了单因素和多因素的COX分析。单因素COX回归分析表明，RiskType与患者的生存期有明显关系。相应的多变量COX回归分析显示，RiskType仍与生存密切相关。上述结果表明，8个基因特征是KIRC的一个独立的预后危险因素（补充图8A，B）。

        将年龄、M期和RiskScore等重要的临床特征结合在多因素考克斯分析中，构建了一个nomogram模型（补充图8C）。结果表明，RiskScore特征对生存率预测的影响最大，说明基于5个基因建立的风险模型可以更好地预测预后。此外，对1年、3年和5年生存率的nomogram数据进行了校正（补充图8D），证明该方法具有很强的预测性能。绘制了年龄、M分期、RiskScore和nomogram的DCA图，结果显示nomogram有很高的净收益（补充图8E）。

补充图8 对8个基因特征的单因素和多因素cox分析

        作者最终选择了4个与预后相关的风险模型，即9基因签名（Zhong）、7基因签名（Jiang）、7基因签名（Chen）和6基因签名（Ren），与8基因模型进行预测性能比较。为了使模型具有可比性，根据4个模型中的相应基因，用同样的方法计算TCGA中每个KIRC样本的风险分数。对RiskScore进行Z-score，RiskScore大于0的样本被归入高风险组，而RiskScore小于0的样本则归入低风险组。计算了两组之间KIRC样本的预后差异。四个模型的ROC和KIRC-KM曲线见图10。可以看出，9基因签名（Zhong）模型的1、3、5年AUC值均低于本研究模型（图10A）；7基因签名（Jiang）（图10C）和6基因签名（Ren）（图10G）模型的1、3年AUC值均低于本研究模型，但5年AUC值却高于本研究模型。7-基因签名（陈）模型的1年AUC值高于本研究模型，但3年和5年的AUC值低于本研究模型（图10E）。这五个模型预测的高组和低组样本的KIRC预后也不同（图10B，D，F，H）。

图10 风险模型与其他模型的比较

05 - 在体外对8个外显基因的表达水平进行验证

        检测了4个人类肾癌细胞系（786-O、A498、Caki-1和ACHN）和正常人类肾小管上皮细胞系HK2中8个表观基因（ETV4、SH2B3、FATE1、GRK5、MALL、HRH2、SEMA3G和SLC10A6）的mRNA和蛋白表达水平。如图11A所示，观察到与HK2细胞系相比，ETV4的mRNA表达水平明显增加，SH2B3、FATE1、GRK5、MALL、HRH2、SEMA3G和SLC10A6的表达水平在肾癌细胞中明显下降。8个epi-PCG的蛋白表达水平与mRNA表达水平相似（图11B）。这些发现与生物信息学的结果一致，表明在多组学数据分析中发现的差异表达的表观PCG在癌细胞中表现出明显的变化。

图11 8种表观基因在体外的表达水平的验证

四、结论

        在这项研究中，作者系统地分析了KIRC中PCG的异常表达，并结合组蛋白修饰，筛选出872个epi-PCG和18629个非epi-PCG。根据差异表达的epi-PCG相关基因，KIRC样本被分为三个亚型，这些亚型在预后方面表现出明显的差异。基于epi-PCG，构建了一个8个基因的预后风险模型，该模型在训练集和验证集都有很强的稳定性和预测性能，不同的RiskScores可以充分反映患者的临床特征。