综合生物信息学分析和组织学验证探究m6A 修饰模式对 PCa 进

Construction of a risk prediction model using m6A RNA methylation regulators in prostate cancer: comprehensive bioinformatic analysis and histological validation

利用m6A RNA甲基化调节因子构建前列腺癌风险预测模型：综合生物信息学分析和组织学验证

发表期刊：Cancer Cell Int

发表日期：2022 Jan 19

DOI:  10.1186/s12935-021-02438-1

期刊相关信息

一、背景

        前列腺癌（PCa）是男性中最主要的恶性肿瘤。PCa的治疗主要是通过外科前列腺切除术或雄性激素剥夺疗法（ADT）。然而，它可以成为去势抵抗性前列腺癌（CRPC），并且在传统治疗过程中可能出现生化复发或转移。

        表观遗传重新编程在PCa的进展中起着关键作用。其中，N6-甲基腺苷（m6A）广泛存在于整个转录组；事实上，m6A约占RNA化学修饰的60%，存在于真核生物中总腺苷残基的0.1%至0.4%。m6A调控基因参与了各种致癌和肿瘤进展过程。据报道，METTL3通过m6A甲基化稳定MYC、LEF1和整合素β1（ITGB1）的mRNA，从而推进PCa的进展，并与不良预后相关。

二、材料与方法

1.数据来源

1）TCGA中前列腺腺癌的转录组分析和单核苷酸变异数据、拷贝数和临床表型数据、基因表达矩阵

2）临床PCa样本：从北京同仁医院和北京朝阳医院招募了45名经病理诊断的PCa患者（15名Gleason评分GS<7，15名GS=7，15名GS>7）

2.实验流程

1）差异表达基因（DEG）分析：使用R软件包"limma "和 "ggpurb "来确定正常组和肿瘤组之间的DEGs；R软件包"clusterProfiler "和 "enrichplot "被用来分析GO和KEGG的路径富集

2）存活率和相关性分析：生存分析；利用m6A聚类的DEGs建立PCa预后模型，进行单变量Cox回归分析

3）共识聚类分析和主成分分析（PCA）：R软件包"ConsensusClusterPlus "进行聚类分析

4）基因组变异分析（GSVA）：GSVA用于评估KEGG基因集富集度

5）肿瘤突变负担（TMB）的估计：分析了m6A评分（以下显示为m6Ascore）与TMB之间的关系，然后进行生存分析，比较高TMB组和低TMB组的预后

6）肿瘤微环境（TME）细胞浸润：单样本基因集富集分析（ssGSEA）计算TME浸润水平；对m6Ascore和免疫相关基因进行了相关性分析

7）实验：定量实时PCR（qPCR）分析、免疫分数（IPS）分析

三、实验结果

01 - PCa中m6A调节因子的特点

        在TCGA-PRAD数据集中，通过与正常样本的比较，分析了PCa的m6A调节因子的拷贝数变异（CNV）分析改变、DEGs和突变频率。对于CNV事件，大约76%的m6A调节因子失去了DNA拷贝数，其中ZC3H14的拷贝数损失程度最高（图1A）。六个m6A调节因子获得了拷贝数，其中VIRMA的拷贝数最高（图1A）。m6A调节因子CNV的改变和在染色体上的位置显示在图1A的下部。在TCGA数据集中的499个肿瘤和52个正常样本中，使用TPM数据统计了m6A调节因子的DEGs。METTL3、HNRNPA2B1、RBM15B和IGFBP2的水平较高，而ZC3H13、FTO和IGFBP3的水平在PCa组织中低于正常组织（图1B）。在突变频率分析中，484个样本中的19个与m6A相关的基因发生了突变；在一个样本中检测到VIRMA（KIAA1429）、YTHDC2、RBM15B、YTHDF2和IGFBP1的突变（图1C）。ZC3H14表现出最高的突变率，所有25个m6A调节因子在TCGA样本中表现出低突变率。这表明在PCa进展过程中，m6A调节因子的表达相对保守和稳定。

图1    前列腺癌中m6A调节因子的特点

02 - m6A调节因子在PCa预后和不同临床病理特征中的表达情况

        利用TCGA-PARD TPM数据估计m6A调节因子在不同GS和病理T（pT）期的表达。由于GS≥7的PCa与预后较差有关，因此将PCa样本分为三组：GS <7，GS = 7和GS >7。与GS<7组相比，IGFBP3、HNRNPA2B、RBMX、RBM15B、YTHDF1、HNRNPC、VIRMA和FMR1在GS>7组明显高表达（P<0.001）（图2A）。此外，还比较了T2（188）、T3（295）和T4（10）pT阶段。在25个m6A调节因子中，IGFBP3、HNRNPA2B1、VIRMA和RBMX在T3期比T2期表达得更高（图2B）。根据生存分析中的K-M曲线，HNRNPA2B1、IGFBP1和ELAVL1的高表达与不良的RFS有关（图2C）。此外，PCa中的TP53突变与较短的放射学无进展生存期（rPFS）和CRPC时间有关。与TP53野生型PCa样本组相比，TP53突变组的VIRMA和IGFBP3的水平较高，IGFBP2的水平较低（图2D）。图2E所示的预后网络中总结了m6A调节因子的风险因素和关系。

        作者分析了m6A调节因子在GS（GS>7与GS<7）、pT（T3与T2）、RFS（单变量Cox回归分析的P值<0.05）、TN（肿瘤与正常PCa组织）和TP53（突变与野生型）中的DEGs比较。通过维恩图确定交叉基因，HNRNPA2B1和IGFBP3在所有的比较中都有差异表达（图2F）。HNRNPA2B1的水平在PCa组高于正常组，并且与晚期参数有关：即GS>7，pT3，HR>1，和TP53突变。然而，IGFBP3的水平在PCa组织中较低，但在存在上述晚期指标时却较高。这表明，PCa中可能存在不同的分子机制和表达模式，并贯穿于整个进展过程。然后，评估了45个PCa组织和15个邻近的正常前列腺组织中HNRNPA2B1和IGFBP3的表达。HNRNPA2B1和IGFBP3的水平在66.7%和40.0%的PCa组织中分别高于邻近的正常组织，与GS＜7组相比，GS＞7组表现出HNRNPA2B1和IGFBP3的表达升高（图2G）。

图2    m6A调节因子的表达与不同的PCa临床病理特征和预后的关系

03 - 根据m6A调节因子的共识聚类分析

        为了进一步探索m6A调节因子在TCGA-PRAD队列中不同预后和临床病理特征中的表型，进行了共识聚类分析，根据25个m6A调节因子（m6Acluster）确定495个PARD病例的亚组。关于群集数的CDF曲线下面积的相对变化，k=3被确定为最佳的群集类别数（图3A）。通过PCA评估子类，所有三个簇都有明显的区分，特别是在簇A和簇C（图3B），表明对m6Acluster的预测是正确的（这些簇显示为 "簇1-3或A-C"）。接下来，比较了集群间PRAD的不同临床病理特征，结果显示集群A的病理N期、GS 8-10和生化复发率相对较高（图3C）。通过GSVA分析了集群A和C之间不同的KEGG途径（图3D）。还分析了不同m6Aclusters中免疫细胞浸润的变化，发现23个免疫细胞亚群中有12个在三个集群中表达不同（图3E） 基于m6Aclusters的RFS的Kaplan-Meier生存曲线显示各组之间没有明显差异（图3F）。

图3    基于m6A调节因子的共识聚类分析

04 - 基于m6Aclusters的DEGs的共识聚类分析

        根据m6A调节因子的三个集群，在每一个配对比较中分析DEGs，并通过Venn图确定74个相交的基因（图4A）。然后分析了这74个基因的生物学功能，只有1个特定的CC（质子运输的V型ATP酶，V0结构域）被富集（图S1A-B）。此外，对74个基因的KEGG信号通路分析表明在氧化磷酸化KEGG通路中明显富集（图S1C-D）。

图S1    对m6Aclusters的74个相交的DEGs的富集分析

        对74个交叉基因进行了单变量Cox回归分析，选出了6个（BAIAP2、TEX264、MMAB、JAGN1、TIMM8AP1和IMP）与PCa复发有关的基因（图4B）。与正常样本相比，分析了6个或5个DEGs（TIMM8AP1被排除在外，因为它是一个经过处理的假基因，并且缺乏拷贝数数据）关于CNV改变、DEGs和突变频率的特点。发现6个基因中的4个（BAIAP2, IMP3, JAGN1, MMAB）在PCa组织中的水平明显上调。有趣的是，JAGN1表现出最高程度的拷贝数丢失，但其在PCa组织中的表达高于正常前列腺组织。当评估这6个基因在临床病理特征GS、pT和TP53方面的表达时，发现大多数基因在晚期（GS>7、pT3和TP53突变）下调。生存分析也表明，这6个基因的低表达与不良RFS明显相关（图S3A-D）。图S3E所示的预后网络中总结了这6个基因的风险因素和关系。

图S3    m6Aclusters的6个DEGs在PCa预后和不同临床病理特征中的表达情况

        对这6个基因进行了共识聚类分析（geneCluster），并确定k=3的结果是关于delta区域结果的最佳分类（聚类1-3或A-C）（图4C）。PCA验证了这三个亚群的意义（图S4A）。三个基因簇的临床病理特征热图显示，基因簇A的pT期、N期、GS和生物复发率都明显高于其他簇（图4D）。K-M生存分析显示三个m6A群之间存在明显差异，其中基因群A的RFS最差（图4E）。接下来，分析了不同基因群中m6A调节因子的表达，发现25个m6A调节因子中有21个在它们之间有差异表达（图4F）。此外，23个免疫细胞亚群中的16个在三个群中有差异表达（图S4B）。在不同KEGG途径的GSVA中，一个与小细胞肺癌相关的途径在基因群A和C之间出现失调（图S4C）。

图S4    基于m6Aclusters的6个DEG的共识聚类分析

        作者还分析了关于GS、pT、RFS、TN和TP53的6个基因的DEGs，如图2G所示。通过维恩图确定了交叉基因MMAB和BAIAP2（图4G），两者在PCa组织中与正常组织相比表达量高，但与四个参数相比表达量低，表明是晚期。这表明MMAB和PAIAP2在PCa发生和发展中的表达模式相反。如上所述，利用45个PCa组织和15个邻近的正常前列腺组织，通过qPCR分析了MMAB的mRNA表达水平。与邻近的正常组织相比，它在60.0%的PCa组织中上调，与GS＜7组相比，GS＞7组的表达量减少（图4H）。

图4    根据m6Aclusters的DEGs进行的共识聚类分析

05 - 基于6个基因的低m6Ascore与RFS的不良预后有关

        上面的聚类分析是基于患者的群体，不能用来定量评估m6A调节因子的模式。考虑到m6A调节因子的复杂性和个体异质性，通过对TCGA-PARD中6个基因水平的PCA计算m6Ascore，并进行统计分析，考察m6Ascore与geneCluster/m6Acluster之间的关系。结果发现两个簇都有明显的差异：基因簇A的m6Ascore最低，基因簇C最高（图5A）。根据RFS分析中m6Ascore的最佳阈值分割，将PRAD患者分为两组，根据K-M生存分析，不良的RFS与低m6Ascore组有关（图5B）。冲积图显示了m6Acluster、geneCluster、m6Ascore和复发状态（fustat）的属性变化（图5C）。所有低m6Ascore组的患者都被归入基因簇A，这与较差的RFS结果有关。作者评估了m6Ascore和免疫相关基因的相关性，结果发现23个免疫相关细胞浸润中有16个与m6Ascore呈负相关（图5D）。这表明在m6Ascore升高过程中对免疫反应的负向调节。以前的研究表明，低的免疫浸润与延长生存期有潜在的关系，与本研究结果一致，高m6Ascore组在RFS方面有良好的预后倾向（图5B）。

图5    基于6个基因的低m6Ascore组与RFS中的不良预后有关

06 - TCGA-PRAD肿瘤突变和亚型的m6Ascore状态的特点

        如图5B所示，将TCGA-RAD患者分为两组（低m6Ascore组58人，高m6Ascore组437人），评估m6Ascore和TMB的关联。结果显示，m6Ascore组的TMB得分没有明显差异（图6A，左），Spearman相关分析也显示m6Ascore和TMB的表达水平之间没有明显的相关性（图6A，右）。根据RFS分析中TMB值的最佳阈值分割，将PRAD患者分为两组（低TMB组242人，高TMB组231人）。在K-M分析中，发现高TMB组与低TMB组相比有不良RFS的趋势，但差异不明显（图6B，左）。这表明PCa患者较高的恶性程度与高TMB评分组有潜在的相关性。当把m6A和TMB一起分析时，发现低TMB与高m6Ascore的组合有很大的RFS优势（图6B，右）。低TMB和高m6Ascore组在RFS方面有更好的预后（图5B和图6B，左）。

        作者分析了m6Ascore组之间体细胞突变的分布差异，结果显示低m6Ascore组的TP53突变频率更广泛，但SPOP的突变频率降低（图6C）。对于m6Ascore和复发状态（fustat）之间的关系，低m6Ascore组的复发率较高（图6D左），这与图5B的RFS分析一致。生化复发也与低m6Ascore有关（图6D，右）。然后通过RFS分析来分析m6Ascore组在不同临床病理特征方面的预后（图6E和图S5A-B）。在GS 8-10和pT3阶段，低m6Ascore组与不良RFS有关（图6E）。

图S5    TCGA-PRAD肿瘤亚型中m6Ascore状态的特征

        检测PD-L1的表达以阐明对免疫疗法的潜在反应，低m6Ascore组显示出相对较高的表达水平图6F）。进一步预测了免疫检查点阻断（ICB）的风险分数的价值，PD1和CTLA4被纳入IPS分析，并进一步分为四个部分：ips_ctla4_neg_pd1_neg（CTLA4负反应和PD1负反应），ips_ctla4_neg_pd1_pos（CTLA4负反应和PD1正反应），ips_ctla4_pos_pd1_neg和ips_ctla4_pos_pd1_pos。在低和高m6Ascore组的比较中，平均IPS显示在PD1和CTLA4的阴性或阳性反应的四个部分没有明显的差异（图6G、图S5C）。这些结果表明，m6Ascore在预测PD1和CTLA4治疗反应的风险评分模型中可能缺乏功效。

图6    TCGA-PRAD肿瘤突变和亚型中m6Ascore状态的特征

四、结论

        在这项研究中，作者通过生物信息学分析和组织验证，全面分析了TCGA-PRAD中的m6A调节因子。首先，筛选了m6Aclusters的DEGs，发现了6个基因（BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1和 IMP3），将PCa患者分为三个亚组，并计算m6Ascore，构建了一个对复发有高预测价值的风险模型。这项研究可能有助于确定m6A信号对PCa的进展和预后的影响。