今天小编和大家分析一篇23年12月发表在iScience（IF:4.6）杂志的文章《Experimental prognostic model integrating N6-methyladenosine-related programmed cell death genes in colorectal cancer》作者在本研究中，利用 TCGA-COADREAD/CRC 队列，鉴定了 854 个 m6A 相关的 PCD 基因，构成了通过 LASSO Cox 回归建立的稳健 10 基因风险模型 （CDRS） 的基础。使用CRC细胞系和新鲜组织进行qPCR实验进行验证。揭示了独特的基因组特征、通路激活以及与肿瘤微环境的关联，为CRC风险分层和个性化治疗途径提供了新的范式。
背景：
结直肠癌 （CRC） 是全球第三大最普遍的癌症，是癌症相关死亡的第二大原因。越来越多的证据表明，基因表达谱可用于估计癌症患者的生存率。因此，确定额外的预后生物标志物，可以根据风险对CRC患者进行分层，从而能够更精确地预测患者的预后，甚至促进个性化治疗计划的制定，最终提高CRC患者的生存率。最近的研究已经确定失调的 RNA 修饰，是癌症进展和发展的重要因素。m6A在正常生理学以及包括癌症在内的疾病中起着关键作用，其失调与肿瘤微环境（TME）重塑有关。在结直肠癌中，异常的 m6A 修饰与肿瘤生长、转移和预后不良有关。程序性细胞死亡 （PCD） 是一个关键过程，通过消除不必要或受损的细胞来维持组织稳态。m6A修饰可以调节PCD基因的表达，从而在癌症的发生和发展中起着至关重要的作用，m6A RNA修饰在CRC发育过程中调节PCD的作用在很大程度上仍未得到探索。
方法：
1. 数据下载。
从癌症基因组图谱（TCGA）数据集和基因表达综合（GEO）中获得了CRC的基因表达谱和相关临床信息。总生存期信息不完整或生存时间少于一个月的患者被排除在外。对于 TCGA-COAD（READ/CRC），使用 R 包“TCGAbiolinks”来访问 RNA 测序数据（TPM 值）、突变谱、拷贝数变异 （CNV） 数据和临床数据。使用“maftools”R 包来可视化体细胞突变数据，并使用 GISTIC2.0（基因模式）来检测体细胞拷贝数变异。CMS亚型是使用“CMScaller”包确定的。整合了meta-GEO数据集，该数据集由749个样本组成，包括来自GSE17537（54个样本）、GSE38832（122个样本）和GSE39582（573个样本）的数据。在将不同数据集合并到meta-GEO数据集中之前，批量效应已得到校正。此外，互联网数据库中获得了 IMvigor210 队列的数据。
2. M6A相关的PCD基因获取
从文献中获取了23 个 m6A RNA 甲基化调节因子和 1078 个 PCD 相关基因的整理和分析，采用Spearman相关性分析评估PCD相关基因和m6A调节因子的表达。
3. 构建风险评分模型
使用 "limma "R软件包中的鉴定正常样本和CRC样本中m6A相关PCD基因的差异表达，显著性标准为调整后P值<0.05且|log2FC|≥0.5。随后，对确定的基因进行了单变量 Cox 回归分析和 Kaplan-Meier 生存分析，以确定与 m6A 相关的 PCD 预后基因。为减少遗漏，将 Cox 分析中的 P 值截断值调整为 0.1，将 Log-rank 检验中的 P 值截断值调整为 0.01。进一步应用 LASSO Cox 回归算法来减少 m6A 相关 PCD 候选基因的数量，并构建最合适的风险评分模型：
βi 表示风险系数，Ei 表示每个基因的表达。根据中位数将结直肠癌患者分为高 CDRS 组和低 CDRS 组。此外，采用 ROC曲线曲线、曲线下面积(AUC)和一致性(c)指数来评估 CDRS 模型的预测准确性。此外，根据临床特征和 CDRS 建立了 CRC 患者的预后列线图。用校正图和决策曲线分析(DCA)评价列线图的疗效。
1. 功能富集分析
为了探索高CDRS组和低CDRS组之间生物学过程的潜在差异，使用 "GSVA "R包进行了GSVA富集分析。还使用 "clusterProfiler "R软件包进行了基因组富集分析（GSEA）。基因组 "h.all.v7.5.1.symbols "和 "c2.cp.kegg.v7.5.1.symbols "是从MSigDB数据库（https://www.gsea-msigdb.org/gsea/msigdb/index.jsp）下载的。
2. 肿瘤微环境分析
为了确定DElncRNAs、DEmiRNAs和种子节点之间的关系，利用Cytoscape软件将上述种子映射到原始lncRNA-miRNA-mRNA网络中，并提取相关的DElncRNAs和DEmiRNAs构建网络。
3. 功能富集分析
免疫学评分、基质评分、ESTIMATE 评分和肿瘤纯度是通过 "ESTIMATE "程序确定的。为了量化 CRC 中每种细胞浸润的相对丰度，采用了 EPIC 算法：此外，还采用了六种独立算法：EPIC、MCP-COUNTER、XCELL、TIDE、CIBERSORT 和 QUANTISEQ 来计算癌症相关成纤维细胞（CAF）、上皮细胞和巨噬细胞的浸润水平。采用了肿瘤免疫功能障碍与排斥（TIDE）和免疫表观评分（IPS）算法来预测不同亚型 CRC 对免疫疗法的反应。
4. 潜在的治疗药物预测
从两个来源检索了癌细胞系的药物敏感性数据：癌症治疗反应门户网站（Cancer Therapeutics Response Portal，CTRP，https://portals.broadinstitute.org/ctrp/）和混合物中同时分析相对抑制数据集（Profiling Relative Inhibition Simultaneously in Mixtures dataset，PRISM https://depmap.org/portal/prism/）。排除了缺失值大于 20% 的化合物，并使用 K-nearest neighbor (k-NN) 估算法填补缺失的 AUC 值。为了预测临床样本中的药物反应，使用了 pRRophetic 软件包中的内置脊回归模型。预测的 AUC 值用于评估药物敏感性，AUC 值越低表明对治疗的敏感性越高。
5. 统计分析
所有生物信息学分析均使用 R 和 Graphpad Prism 进行（以均值 +/- SD 表示）。Student-t 检验和单因素方差分析分别用于比较两个或多个正态分布组之间的差异。对于非正态分布数据，分别采用 Wilcoxon 检验或 Kruskal-Wallis 检验来比较两组或多组之间的差异。连续变量之间的相关性采用斯皮尔曼相关分析法进行分析。分类变量之间的差异采用卡方检验（chi-square）和费雪精确检验（Fisher's exact tests）进行评估。使用 "pROC "R 软件包生成 ROC 曲线，使用 "CompareC "软件包比较不同变量的 C 指数。采用单变量 Cox 回归分析估算 CDRS 和 m6A 相关 PCD 基因的危险比 (HR)。使用多变量 Cox 回归分析确定了独立的预后因素。在所有分析中，双侧 P 值小于 0.05 被认为具有统计学意义。
研究结果
1. 确定并验证结直肠癌患者中与 N6-甲基腺苷相关的程序性细胞死亡基因
为了深入了解 m6A 相关 PCD 在 CRC 中的作用，首先分析了 m6A 调控因子在 CRC 和正常组织中的表达情况，为了探讨它们在 CRC 中的预后意义，进行了 Cox 回归分析和 Kaplan-Meier 生存分析。结果发现，ELAVL1、FMR1、HNRNPA2B1、LRPPRC、YTHDC2、YTHDF3 和 METTL14 是与总生存期显著相关的保护因子。相比之下，FTO、VIRMA 和 ZC3H13 被归类为危险因素（图 S1B）。此外，还通过斯皮尔曼相关性分析评估了这些 m6A 调控因子之间的相互调控。研究结果表明，在同一功能类别中，表达模式高度相关，书写者、擦除者和阅读者之间也存在显著相关性。这些结果强调了写入者、阅读者和擦除者调控因子之间错综复杂的交叉作用，这在 CRC 的发生和发展中起着至关重要的作用。接下来，从 TCGA-CRC 队列中提取了 m6A 基因和 PCD 基因的表达矩阵，并对其进行了斯皮尔曼相关性分析。与显著相关性（|R| ≥ 0.3 且 p 值<0.05）的 PCD 基因被定义为 m6A 相关 PCD 基因。共鉴定出 854 个 m6A 相关 PCD 基因（图 1A）。
2. 基于N6-甲基腺苷相关程序性细胞死亡基因构建风险模型
在TCGA-CRC队列中，共发现了276个差异表达基因（DEGs）165 个基因在 CRC 中上调，111 个基因下调（图 1B）。利用 TCGA-CRC 队列的 Cox 回归分析和 Kaplan-Meier 生存分析，确定了 14 个与 m6A 相关的 PCD 基因，这些基因与患者的总生存期（OS）有显著关联（图 1C）。选择了 10 个 m6A 相关 PCD 基因，以 TCGA-CRC 数据为训练数据集，通过 LASSO Cox 回归分析建立预测模型（图 1D-1F）。这些基因在其他验证数据集中也显示出与 m6A 调控因子的显著相关性（图 S2A）。GSE32323 中进一步验证 发现，除 TNFRSF10A、SFPQ 和 ITGA6 外，所有基因的表达趋势一致。此外，生存分析证实，在多个数据集中，除 SFPQ 和 MTM1 外，这些基因都与 CRC 预后有关。
接着用 qRT-PCR 验证了 INHBB 和 SNAI1 在细胞系和临床组织样本中的表达水平（图 1G 和 1H）。结果证实，INHBB 和 SNAI1 在 CRC 细胞系和新鲜组织样本中上调。
然后，根据 10 个选定的 m6A 相关程序性细胞死亡（PCD）基因的表达及其各自的回归系数，计算出每位患者的细胞死亡相关风险评分（CDRS）（图 1F）。CDRS 的计算公式如下：CDRS = 0.345 ∗ DAPK1 +0.179 ∗ INHBB +0.067 ∗SNAI1 - 0.028 ∗ TRAP1 - 0.034 ∗ ITGA6 - 0.046 ∗ GSKIP - 0.068 ∗ GLA - 0.09 ∗ MTM1- 0.111 ∗ SFPQ - 0.152 ∗TNFRSF10A （图 1F）。CDRS 与 14 个 m6A 相关 PCD 基因以及 23 个 m6A 基因之间的相关性。结果表明，CDRS 与 m6A 相关 PCD 基因和 m6A 基因呈显著正相关（图 S2G 和 S2H）。相反，CDRS 与其他 7 个 m6A 相关 PCD 基因、9 个 "阅读者 "和 5 个 "写作者 "呈显著负相关（图 S2H 和 S2I）。此外，还发现 CDRS 与各种临床特征（包括生存状况、性别、T 期、N 期、M 期和 AJCC 分期）显著相关（图 1I）。共识分子亚型（CMS）根据转录组图谱将 CRC 分成四个不同的亚型：CMS4 亚型患者的 CDRS 评分最高，其总生存率较差（图 1I）。
3. 作为结直肠癌独立风险因素的细胞死亡相关风险评分
根据 CDRS 的中位值将 CRC 患者分为高 CDRS 组和低 CDRS 组。在 TCGA-CRC 训练数据集和其他验证数据集中，高风险组患者的总生存期明显低于低风险组（图 2A）。为评估 CDRS 的独立预后价值，进行了单变量和多变量 Cox 分析（图 2B 和 2C）。单变量分析结果表明，高 CDRS 与总生存率降低之间存在显著相关性，所有数据集的 HR 均大于 1（图 2B）。在对年龄、性别、T期、N期、M期、AJCC分期和微卫星状态等现有临床特征进行调整后，多变量分析显示 CDRS 仍是一个具有统计学意义的独立预后因素（图 2C）。此外， AUC 和 c 指数分析表明 CDRS 在多个独立队列中始终表现出稳健的性能（图 2D 和 2E）。鉴于目前已有几种基于 m6A 或 CRC 个体死亡率模式的预后模型，进一步比较了 CDRS 与其他特征在预测 CRC 预后方面的表现。对 TCGA 和 GEO 数据集进行的单变量 Cox 回归、C-指数和 AUC 分析一致表明 CDRS 优于其他特征（图 2F、2G）。为了提高预测 CRC 患者总生存期的准确性，在 TCGA 队列中使用多变量 Cox 回归建立了一个列线图模型，以估计 1 年、3 年和 5 年的总生存期。该模型包括年龄、M、AJCC 分期和 CDRS。校准曲线证实了列线图在预测 1 年、3 年和 5 年总生存率方面的准确性。此外，决策曲线分析（DCA）表明，列线图模型优于本研究中使用的其他预测指标。
4. 与细胞死亡相关风险评分特征相关的生物通路
CDRS 在预测 CRC 预后方面的卓越能力激发了了解其潜在机制的好奇心。为了深入研究这些机制，采用了基因组变异分析（Gene Set Variation Analysis，GSVA）来探讨 CDRS 对所有数据集的生物通路的影响。分析表明，高 CDRS 组表现出基质和炎症相关通路的激活，包括但不限于上皮间质转化（EMT）、血管生成、顶端连接、炎症反应和伽马干扰素反应。相比之下，低 CDRS 组显示出与细胞生长和促癌机制相关的通路上调，包括 DNA 修复、G2M 检查点、PI3K/AKT/mTOR 信号转导和 Wnt/Beta-Catenin 信号转导（图 3A ）此外，相关性分析表明 CDRS 与 EMT1/2/3、癌相关成纤维细胞（CAFs）和泛成纤维细胞特征（PAN-F-TBR）等基质相关特征呈正相关，而与细胞生长相关通路呈负相关（图 3B）。虽然 CDRS 与 CD8 T 效应子和抗原处理机制（APM）特征没有明显的相关性，但它与免疫检查点阻断（ICB）抗性相关信号呈负相关（图 3B）。
5. 不同细胞死亡相关风险评分组的基因组特征
接着比较了 CDRS 高值组和低值组之间的体细胞突变和拷贝数改变，但结果显示两组之间没有显著差异。利用fish精确检验（阈值为 p <0.01），发现了 24 个不同的突变基因。虽然这些基因的突变率较低，但观察到低 CDRS 组的 APC 突变频率较高，而高 CDRS 组的 BRAF 突变频率较高（图 4A）。在 TCGA-CRC 和 GSE39582 数据集中，BRAF 突变组的 CDRS 较高，而 KRAS 和 TP53 的野生型组和突变组之间的 CDRS 没有明显差异（图 4B）。
6. 细胞死亡相关风险评分与免疫微环境的关系
鉴于 CDRS 与某些基质相关特征之间的重要关联，对 CDRS 与肿瘤微环境之间的关系进行了广泛的调查。首先，利用ESTIMATE算法计算了TCGA-CRC和meta-GEO数据集中的基质得分和免疫得分。研究结果表明，CDRS特征与基质评分、免疫评分和ESTIMATE评分（基质评分和免疫评分的整合）呈正相关，而与肿瘤纯度呈负相关（图5A）。接下来，采用 EPIC 方法量化了 TCGA-CRC 中的免疫细胞浸润情况。值得注意的是，基质微环境的关键成分，如癌症相关成纤维细胞（CAFs）和上皮细胞，与 CDRS 呈明显的正相关（图 S7A）。为了进一步评估免疫细胞浸润，使用了多种算法，包括 EPIC、MCP-COUNTER、XCELL、TIDE、CIBERSORT 和 QUANTISEQ。结果显示，CDRS 与 CAFs、上皮细胞和 M2 样巨噬细胞呈正相关（图 5B）。此外，与 m6A 相关的 PCD 基因 DAPK1、INHBB 和 SNAI1 与 CAFs、上皮细胞和 M2 样巨噬细胞呈正相关，而 GSKIO、ITGA6 和 TRAP1 则表现出相反的趋势（图 S7B）。已知转化生长因子 β（TGF-β）是 CAF 活化的重要调节因子。与 CDRS 低的患者相比，CDRS 高的患者表现出更高的 TGF-β 信号转导（图 5C ）。此外，还研究了免疫调节剂与 CDRS 值之间的相关性，发现 CDRS 主要与 TGFB1 呈正相关。此外，巨噬细胞相关趋化因子 CCL5 和 CX3CL1 在高 CDRS 组中高表达。使用 TIDE 方法，发现高 CDRS 组在 TCGA-CRC 和 meta-GEO 队列中都有较高的 TIDE、排除得分（仅在 meta-GEO 数据集中）和功能障碍得分（图 5D 和 5E）
7. 细胞死亡相关风险评分在预测药物敏感性中的作用
为了研究 CDRS 与药物敏感性之间的关系，利用癌症治疗反应门户网站（CRPT）和混合物中同时分析相对抑制（PRISM）数据集中来自数百个癌细胞系（CCL）的大规模药物敏感性和基因表达数据，建立了药物反应预测模型（图 6A）。利用基因表达谱，通过 "pRRophetic "软件包采用脊回归模型估算了每个临床样本中每种化合物的AUC值。结果发现，AUC 值与药物敏感性呈负相关。详细的工作流程如图 6A 所示。首先，在高 CDRS 组（最高十分位数）和低 CDRS 组（最低十分位数）之间进行了差异药物反应分析（|log2FC| > 0.10），以确定 AUC 值差异最显著的化合物（图 6A）。接下来，通过计算 AUC 值与 CDRS 之间的 Spearman 相关性来筛选化合物，选出 CTRP 和 PRISM 的 Spearman |R| 相关系数均大于 0.30 的化合物（图 6A）。最后，确定了 10 个源自 CTRP 的化合物和 11 个源自 PRISM 的化合物（图 6B 和 6C）。其中一种化合物达沙替尼是一种口服生物利用度高、前景广阔的治疗药物，用于治疗各种人类恶性肿瘤，不过它对 ABL、BCR-ABL 和 SFKs 具有相对特异性。33 虽然 CDRS 在 KRAS 突变型和野生型结直肠癌患者之间未显示出显著差异，但发现高 CDRS 组患者的达沙替尼估计 AUC 值较低（图 6B 和 6C）。同样，其他五种药物（romidepsin、panobinostat、vindesine、YM-155、棘霉素）在高 CDRS 组的估计 AUC 值也较低，这表明它们可能是高 CDRS CRC 患者的潜在治疗选择（图 6B 和 6C）。
