王炸！15+高分文章教你如何只用公共数据库研究RNA修饰

m6A, m1A，选择性多聚腺苷酸化（APA）和腺苷到肌苷的RNA编辑（A-I），这四种常见的RNA腺苷修饰，主要是由“writer”酶介导的，构成免疫反应和肿瘤发生的表观遗传调控的关键机制。然而，这些“writer”在肿瘤微环境(TME)、药物敏感性、靶向治疗和免疫治疗中的相互作用和潜在作用仍不清楚。下面我们就来看一篇在结直肠癌里面的研究。这篇文章是：Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer（四种RNA修饰writers的交互定义了结直肠癌的肿瘤微环境和药物基因组景观）。发表于2021年的Molecular Cancer，IF是15.302，妥妥的高分文章啊。

本研究系统地描述了26个RNA修饰“writer”在结直肠癌(CRC)中的mRNA表达和基因改变，并评估了他们在8个数据集1697个CRC样本中的表达模式。使用无监督聚类方法将样本划分为两种RNA修饰“writer”的表达模式。随后，构建基于负责RNA修饰模式的差异表达基因(DEGs)的RNA修饰“writer”评分(WM_Score)模型，以量化个体肿瘤的RNA修饰相关亚型。此外，还对WM_Score与TME特征、共识分子亚型(CMSs)、临床特征、转录和转录后调控、药物反应和免疫治疗疗效进行了关联分析。本文的工作流程如下figure1 S所示

结果证实RNA修饰“writer”的多层改变与CRC患者生存和TME细胞浸润特征相关。发现WM_Score-high和-low两种不同的RNA修饰模式。WM_Score-high组与较差的患者总生存和抑制免疫细胞的浸润相关，如M2巨噬细胞、EMT活化和转移，而WM_Score-low组与生存优势、凋亡和细胞周期信号通路相关。WM_Score与CRC的转录和转录后事件调控高度相关。对于抗癌药物，WM_Score与靶向致癌相关通路如MAPK、EGFR、mTOR信号通路的药物呈高度负相关(药物敏感)，与靶向凋亡和细胞周期的药物呈正相关(耐药性)。重要的是，WM_Score与PD-L1阻断的疗效相关，这表明靶向这些“writers”的潜在药物的开发有助于免疫治疗的临床效益。（本文的图文总结如下figure7 S所示：）

下面我们就来看一下本文的详细结果是怎样呈现的

1. CRC中四种RNA修饰“writers”的遗传和转录改变

根据已发表数据，本研究共纳入26个RNA修饰“writers”，其中A-I修饰“writers”3个，m6A修饰“writers”7个，m1A修饰“writers”4个，APA修饰“writers”12个。

为了确定癌症中RNA修饰作者的基因改变，评估了26位writers的非沉默体细胞突变的患病率。在TCGA的LAML、PCPG和UVM三个癌症队列中，单个writer突变频率相对较低，而COAD队列中“writers”突变频率相对较高(图S2A)。在404个COAD样本中，119个(29.46%)存在RNA修饰“writers”突变(图1A)。然而，有这些“writers”突变的CRC患者总生存期比没有突变的患者短(图1B)，提示“writer”基因改变可能在CRC中起功能性作用。接下来，使用hallmark基因集进行了基因集变异分析（GSVA）的富集分析以比较“writers”突变组和非突变组。发现在突变组中更多的肿瘤标志物相关基因集被富集（图S2B）。 提示writer基因突变可能导致功能改变，影响结肠癌的生存预后。

然后，检测这些“writers”的体细胞拷贝数变化，发现CSTF1、CPSF1/4、ZC3H13和KIAA1429具有广泛的拷贝数变化(CNV)增益(图1C)。为了确定这些基因变异是否影响了CRC患者中RNA修饰“writer”的表达，比较了配对正常和CRC样本中调节因子mRNA的变化，发现大多数“writer”的表达在CRC中显著增加(图1D-G)。与正常结肠组织相比，CRC组织中CNV增加的RNA修饰“writers”(如CSTF1和CPSF1)更频繁(图1C和G)， 表明CNV可能是“writers” mRNA表达的调节因子 。然而，一些“writer”表达上调，但CNV丢失频率高。为探讨CNV值与肿瘤组织中mRNA表达的差异，选取8位CNV丢失率超过20%的研究者， 根据CNV值将大肠癌患者分为4组，包括CNV增加、CNV减少和CNV无显著改变 。然后，比较了这几组“writer”mRNA的表达情况(图S2C)。在这些“writer”中，CNV增加的患者比CNV减少的患者表达更高。与正常组织相比，METTL14、ADORB1、CPSF2、PABPN1在CNV减少组中表达显著下调或无显著变化。肿瘤的发生是一个复杂的过程，CNV的变化不能完全解释肿瘤与正常组织中“writer”表达的差异。虽然许多“writer”的表达变化可以用CNV来解释，但CNV是部分而不是唯一调控mRNA表达的因子。DNA甲基化和转录因子等其他特征，可以调节基因表达。

这一分析表明， 正常和CRC样本之间的基因景观和RNA修饰“writer”的表达具有高度的异质性，RNA修饰“writer”的表达失衡在CRC的发生和发展中具有潜在的作用 。

2. 与癌症特征和免疫浸润相关的RNA修饰“writers”的独特模式

为了全面了解参与CRC肿瘤发生的“writers”的表达模式，选择 8个包含临床信息的数据集(GSE41568、GSE39582、GSE13294、GSE14333、GSE18105、GSE20916、GSE21510、GSE37892) 中的1695个CRC样本进行进一步分析。单因素Cox回归显示，在GSE39582数据集中，26个RNA修饰“writers”中有10个与CRC预后相关(图S3A)。

为了探究writers之间的关系，计算了26个writers在CRC中的表达的成对相关性，发现正相关比负相关更频繁。发现不仅同一类别中RNA修饰“writers”的表达显著相关，而且不同类型的修饰“writers”之间也存在显著相关。ADARB2、ZC3H13和PABPN1的表达与其他“writers”呈负相关。而TRMT61B、TRMT6、KIAA1429、TRMT10C、CSTF2/3、CPSF4、CLP1的表达呈正相关(图2A)。因此，“writers”之间的交互可能对单个肿瘤之间产生不同的RNA修饰模式很重要。

接下来，根据26个RNA修饰“writers”的表达谱，应用共识聚类(Consensus Clustering)对RNA修饰模式定性不同的患者进行分类。经过非监督聚类，合并数据集中的727例CRC患者在Cluster_1中被识别，而其他968例患者在Cluster_2中被识别(图2B)。在RNA修饰模式的预后分析中，Cluster_2修饰模式的亚型中显示出特别显著的生存优势(图2C)。为了鉴定这些不同的RNA修饰模式的生物学意义，进行GSVA富集分析。 Cluster_1显著富集于基质和致癌激活通路，提示RNA修饰“writers”可能与肿瘤发生有关。Cluster_2在与增殖和凋亡相关的通路中富集(图2C) 。

大量研究证实了TME浸润免疫细胞与RNA修饰之间的关系。因此，本文试图研究“writers”在TME中的功能作用。为了比较RNA修饰模式之间免疫细胞成分的差异，使用 CIBERSORT方法 ，这是一种使用支持向量回归来确定肿瘤中免疫细胞类型的反卷积算法。该方法表明RNA修饰“writers”可能与TME细胞浸润有很强的相关性(图S3B)。分析两种RNA修饰簇在TME细胞浸润方面的差异。观察到Cluster_1中M2巨噬细胞、调节性细胞(Tregs)、T滤泡辅助细胞(Thf细胞)和T γδ细胞浸润率较高。在Cluster_2中，活化的树突状细胞、自然杀伤细胞和M1巨噬细胞的浸润增加(图2E)。总的来说， RNA修饰的Cluster_1在免疫抑制细胞中富集，如M2巨噬细胞和Treg，预示预后不良(图2D) 。 M2巨噬细胞在Cluster_1中显著富集 ，而M1巨噬细胞在Cluster_2中占优势(图2E-F)。同样，对巨噬细胞标记物的表达分析显示，M2巨噬细胞标记基因IL1R1、FIZ1、TGFB1、IL10、ARG1在Cluster_1中较Cluster_2中显著上调，而M1巨噬细胞标记基因IL12A、NOS2、IL23A、IL15RA显著下调(图2G)。这表明 RNA修饰模式影响了特定免疫细胞类型的浸润程度，但没有改变浸润免疫细胞的类型 。

3. 构建RNA修饰“writers”signature

为了进一步表征上述两种RNA修饰模式的功能作用，本文鉴定了463种RNA表型相关的DEGs并进行了富集分析。发现这些基因在生物过程中表现出富集，特别是与DNA复制、细胞周期和tRNA代谢过程相关的基因(图S4A)。它们也在信号通路中富集，特别是p53信号通路和IL-17信号通路(图S4B)。为了进一步验证这种差异调控，对与RNA修饰相关的463个基因进行了无监督聚类分析。该分析将CRC患者分为两个基因组亚型： 基因cluster_A和基因cluster_B 。与RNA修饰模式的聚类一致(图2B)，562例CRC患者中，205例聚集在与Cluster_1相关的基因cluster_A中，357例聚集在与Cluster_2相关的基因cluster_B中。此外，患者的基因Cluster_A的预后比基因Cluster_B的患者差(图S4C)，与Cluster_1的患者类似(图2D)。

考虑到RNA修饰的异质性和复杂性，基于这些表型相关基因构建了基于DEGs的评分模型，以量化个体CRC患者的RNA修饰模式；这个模型被称为 WM_Score (“writer”of RNA Modification_Score) 。发现，Cluster_1的WM_Score显著高于Cluster_2(图3B)。同样，基因cluster_A的WM_Score显著高于基因cluster_B(图3C)。为了评估WM_Score对TME的影响，比较WM_Score低和高组之间的免疫细胞浸润情况。发现WM_Score-high组M2巨噬细胞、Tregs、Tfh和Tγδ细胞的浸润更高，而WM_Score-low组活化的DCs和M1巨噬细胞的浸润更高(图S4D)。根据韦恩图和频率分布直方图对这三种不同的分类器进行重叠分析。如图S4所示，结果表明，这三种分类计算方法具有很高的符合度。

为了进一步评估WM_Score的临床相关性，将患者分为WM_Score-low和-high组。WM_Score-low的患者显示显著的生存益处(图3D)。WM_Score ROC曲线的AUC在总生存3、6和12个月时分别为0.64、066和0.64(图3H)。为了检验WM_Score是否可以作为一个独立的预后因素，使用患者的临床特征（年龄、性别和TNM状态）进行多因素Cox回归分析。发现WM_Score是评估患者预后的一个稳健且独立的预后生物标志物(图3F)。使用来自TCGA的562例CRC患者样本验证了WM_Score的可靠性。单变量Cox回归和和多变量Cox回归分析显示，WM_Scor-low组比WM_Score-high组有更好的总体生存优势(图3E,G)。这些结果提示 WM_Score可以反映RNA修饰模式并预测CRC患者的预后 。

4. CRC患者WM_Score的分子分型及临床特征

CRC可分为4种共识分子亚型(CMS1-4，consensus molecular subtypes) ，包括CMS1(微卫星不稳定)：具有免疫细胞浸润、高突变状态、不稳定的微卫星和强免疫激活特性；CMS2（典型）伴WNT和MYC信号的显著激活；CMS3（代谢异常）伴明显的代谢失调；CMS4（间充质）伴显著的TGF-β基质细胞的活化、侵袭和血管生成。在图2C中已经证明， Cluster_1 伴随EMT、TGF-β和VEGF信号通路的激活。此外， WM_Score与EMT评分呈正相关 (图S5A)。在TCGA-COAD/READ队列中，WM_Score-high组的EMT评分显著高于WM_score-low组(图S5B)。

为了检验WM_Score和CMS亚型之间的关系，分别比较了四个GEO数据集(GSE39582, GSE13294, GSE14333, GSE20916)和TCGA队列中不同CMS亚型的WM_Score。发现CMS亚型之间的WM_Scores存在显著差异，其中CMS4组得分最高(图4A-B)。WM_Score-high和-low组之间CMS亚型分布也有显著性差异。CMS4亚型多见于WM_Score-high组，而CMS1/2亚型多见于WM_Score-low组(图S5C-D)。进一步分析不同类型CMS的信号通路特征。CMS4中激活的信号通路主要是WNT、TGF-β、VEGF、EMT信号通路，而CMS1/2中激活的特征信号通路主要是与细胞周期和蛋白酶体相关的信号通路(图S 5E-F)。同样，在两个队列中，WM_Score-high组中CMS4相关信号通路的富集分数显著更高，而WM_Score-low组中CMS1/2相关信号通路的富集分数显著更高(图4C-D)。

CMS亚型反映在肿瘤进展和临床结果中。CMS4肿瘤倾向于在更晚期被诊断(图S5G)，并与较差的总生存率相关(图S5H)。本研究进一步证明了WM_Score在肿瘤分期中是不同的，并且在晚期的CRC中更高(图4E)，这表明包含WM_Score的参数参与了肿瘤进展。WM_Score与CMS4亚型的相关性以及WM_Score与CRC分期的相关性提示WM_Score可能通过反映肿瘤转移而与患者生存相关。转移性结直肠癌患者的WM_Score明显高于非转移性结直肠癌患者(图4F)。这些结果表明， WM_Score与CMS亚型密切相关，高WM_Score可能与EMT、TGF-β等介导肿瘤转移的信号通路的激活有关，提示预后较差 。

5. WM_Score参与转录和转录后调控

WM_Score是一个基于26个RNA修饰“writers”表达的评估模型，这些“writers”调控转录后修饰，包括RNA转运、定位、翻译等生物过程。为了进一步评估WM_Score在解释转录和转录后事件中的作用，关注RNA修饰“writers”相关过程(如APA、A-I编辑)。首先分析TCGA-COAD/READ队列中不同RNA修饰模式间miRNAs表达的差异。筛选出39个WM_Score-high和-low组间显著差异表达的miRNA，并对其靶基因信号通路进行富集分析。miRNA靶向基因与细胞周期、凋亡、EMT、PI3K-Akt等信号通路显著相关，且靶向基因在两组间存在差异表达。EMT信号通路中14个miRNA靶向的DEGs中有13个高表达(图5A)。相比之下，WM_Score-low组中表达较低的一组miRNA的靶基因在细胞周期和凋亡信号通路中表达上调并显著富集。这些分析表明，WM_Score与miRNA的表达和信号通路的调控有关。

为了探讨RNA修饰writers的功能作用，分析TCGA-COAD/READ队列中各基因的APA和A-I编辑事件，以观察转录后特征。鉴定不同RNA修饰模式之间APA差异的基因，并比较了与这些基因相关的生存，以确定3'UTR长度是否影响CRC患者的生存(图S6A)。大多数缩短APA事件的基因在WM_Score-low组中富集，并与较短的生存相关(图5B和S6A)。HEATR3和YBX2转录本具有统计学意义的缩短，这与CRC患者较差的生存率相关(图5C)。推测在WM_Score-high组中，由于HEATR3和YBX2的缩短，进而缩短了3'UTR, miRNA可能不能靶向相应的基因，从而激活基因表达，促进结CRC的起始和发展。

接下来鉴定159个TCGA-READ/COAD样本中WM_Score-high和-low组间A-to-I编辑差异的基因。A-to-I编辑率较低的基因在WM_Score-low组中富集，且与CRC患者生存期较差相关(图5D和S6B)。以NUP43和KCNE3基因为例，NUP43和KCNE3的A-to-I编辑与CRC患者更短的生存时间相关(图5E)。 这些基因在WM_Score-high和-low组之间的A-to-I编辑率的差异可能是由miRNA通过编辑3'UTR区域来调控的，从而影响CRC的发生发展 。

6. WM_Score的潜在治疗价值和预测PD-L1阻断剂免疫治疗

为了进一步了解WM_Score对药物反应的影响，评估WM_Score和癌细胞株对药物反应之间的关系。利用Spearman相关分析，在GDSC数据库（Genomics of Drug Sensitivity in Cancer）中发现42对WM_Score与药物敏感性显著相关的对(图6A)。进一步分析了这些药物所靶向基因的信号通路。发现敏感性与WM_Score-high相关的药物多针对MAPK、mTOR和VEGF信号通路。相比之下，敏感性与WM_Score-low相关的药物以凋亡和细胞周期信号通路为靶点(图6B)。 结果表明RNA修饰模式与药物敏感性相关。因此，WM_Score可能是建立适当治疗策略的潜在生物标志物 。

目前，肿瘤突变负荷（TMB）和PD-L1蛋白的表达等生物标志物的研究已经取得了很大进展。考虑到WM_Score似乎与肿瘤的免疫微环境有关(图S4D)，基于两个免疫治疗队列，检测WM_Score预测患者对 ICB（The immune checkpoint blockade，免疫检查点阻断） 治疗反应的能力 (图6C-I)。发现两个在 抗PD-L1 队列中（ IMvigor210队列 和膀胱癌队列），WM_Score-low患者表现出显著的临床获益，并显著延长了总生存期。IMvigor210队列中的348例患者对抗PD-L1受体阻滞剂表现出不同程度的反应，包括 完 全反应(CR) 、部分反应(PR)、疾病稳定(SD)和疾病进展(PD)，CR患者的WM_Score低于其他类型反应的患者(图6D)。在WM_Score-low和-high组之间进行的卡方检验也显示了WM_Score-low患者显著更好的治疗结果(图6E)。分析IMvigor210三个免疫亚型的WM_Score，包括“免疫炎症”、“免疫排除”和“免疫荒漠”(图6F)，观察到“免疫炎症”类型的WM_Score低于其他两组。此外，WM_Score-low组的TMB和新抗原负担明显高于WM_Score-high组(图6G-H)，这可能部分解释了WM_Score-low组的生存优势和ICB治疗的更大益处。

研究还发现，在WM_Score-high的肿瘤中，M2巨噬细胞的活化和基质细胞的TME明显升高，这些过程可能介导了肿瘤的免疫耐受。这一结果提示， WM_Score-high的肿瘤可能代表对免疫治疗有耐药性的“冷肿瘤” 。综上所述， 建立的RNA修饰评分模型的writers可能有助于CRC药物的选择和对抗PD-L1免疫治疗应答的预测。