早！今天小编和大家分析一篇23年4月发表在Clin Epigenetics. （IF: 7.259）杂志的文章《Comprehensive analyses of molecular features, prognostic values, and regulatory functionalities of m6A-modified long non-coding RNAs in lung adenocarcinoma 》。作者通过检测总RNA中m6A的修饰水平和探索m6A调节因子的表达水平来评估和证实m6A修饰与LUAD之间的关系。利用M6A修饰的LncRNA表位转录组微阵列分析m6A修饰的lncRNA的特征，并筛选出LUAD中具有差异甲基化水平的lncRNA。结合癌症基因组图谱(TCGA)数据库的分析，开发和验证了6个m6A调节的lncRNA(6-MRlncRNA)风险评分来预测LUAD患者的总生存期(OS)。最后预测了差异m6A修饰的lncRNAs在LUAD中的竞争性内源性调节作用。
背景
根据2020年GLOBOCAN的统计数据，肺癌是世界上第二常见的肿瘤类型，死亡率居世界首位。尽管在过去三十年中，抗肿瘤治疗有了很大的改进，但区域性和广泛播散的肺癌患者的5年生存率分别只有27%和4%。肺腺癌（LUAD）是最普遍的肺癌组织学亚型，约占病例的40%，且近年来发病率不断增加。为了提高治疗效果，延长患者生存期，需要新的有效生物标志物用于LUAD的早期检测、预后和监测，以及探索该疾病的分子机制。
在共有基序DRm6ACH(D=G/A/U，R=G/A，H=A/U/C)内发现了人类mRNA中最普遍的转录后修饰N6-甲基腺苷(m6A)，并且在长的内部外显子中、终止密码子周围以及在mRNA的3’非翻译区(3’UTR)中高度富集。这种细胞修饰是受三种类型的RNAm6A调节剂(包括甲基转移酶(METTL3，METTL14，METTL16，VIRMA，RBM15，RBM15B，WTAP，ZC3H13和CBLL1)，脱甲基酶(FTO和ALKBH5)和结合功能蛋白质)调节的高度动态状态的一部分(YTHDF1/2/3，YTHDC1/2，HNRNPC/A2B1，RBMX，IGF2BP1/2/3，FMR1，EIF3A，ABCF1，ELAVL1，G3BP1/2，ZNF217和LRPPRC)。这些调节器类型分别介导m6A修饰的发生、去除和功能。一些研究发现，RNA的异常丰度m6修饰，可能由m的突变或异常表达介导6调节因子是某些肿瘤中细胞分化和增殖中的一种新的表观遗传基因调控机制。
长非编码RNA（lncRNA）被定义为一种长度超过200个核苷酸且蛋白质编码能力有限或不存在的RNA分子，在物种之间保守性差，并且在各种组织和细胞中表现出动态和特异性表达模式。由于发现lncRNA在各种人体器官中的丰度通常远低于mRNA，因此最初认为它们是一种转录副产物。然而，研究人员后来证明，lncRNA可以显著影响肿瘤细胞中的细胞增殖和分化，因为它们在某些癌基因的转录和转录后调节中具有多种作用。目前，在许多lncRNA中发现了m6A修饰，并且可以通过影响自身稳定性，核积累，RNA蛋白相互作用和结合microRNA(miRNA)的亲和力来支持其调节肿瘤相关基因的表达水平。然而，在LUAD中m6A修饰的lncRNA谱仍然不清楚。
方法
1.临床组织样本采集和细胞培养
2.核糖核酸提取和逆转录-qPCR
3.总m6ARNA的量化
4.蛋白质印迹
5.LUAD表达谱的获取及靶基因的表达分析
6.m6A调节因子的相关性分析和蛋白质-蛋白质交互网络
7.M6A修饰的LncRNA表观转录组微阵列芯片数据分析
8.MeRIP-qPCR
9.M6A调控的LncRNA预后特征的构建与验证
10.对预后特征的独立性评估和分层分析
11.M6A诱导的竞争性内源性RNA(CeRNA)网络的构建
研究结果
LUAD与低水平的总m6A修饰RNA有关
为了探讨m6A修饰在LUAD中的潜在作用，首先检测了总RNA样品中的m6A水平。发现与邻近的非肿瘤肺组织相比，大多数早期LUAD肿瘤组织表现出RNAm6A甲基化水平降低(N=20，P<0.01)(图1A)。与正常肺上皮细胞BEAS-2B相比，6种LUAD细胞系RNA的m6A修饰丰度显著降低(P<0.05)；在H838，PC9和A549细胞系中尤其明显(图1B)。这些结果表明，LUAD组织和细胞中RNAm6A的修饰水平降低。
M6A调节因子在LUAD组织和细胞系中的表达
通过TCGA和GTEx数据库的联合评估30个已报道的m6A调节因子的基因表达水平，这些调节因子的转录表达水平在LUAD肿瘤组织中全部紊乱(除METTL14外，所有P<0.05)，但所有差异均不显着(|log2FC|<1)(图2A)。RT-qPCR检测表明，LUAD组METTL3和FTO降低(P<0.01)，IGF2BP3升高(P<0.01)。METTL14，ALKBH5和IGF2BP1没有显示出显著的统计学差异(图2B)。总的来说，在8对配对的临床LUAD样本中，6种调节因子的无序表达模式显示出与数据库分析结果相似的趋势。接下来，蛋白质印迹分析显示，LUAD肿瘤组织中六种调节因子的蛋白质水平均显着增加(图2C)，这与具有BEAS-2B或HBE细胞作为正常对照的LUAD细胞系中的结果一致(图2D)。只有IGF2BP3在转录本和蛋白质水平上表现出一致的差异，这两者在LUAD中都升高。这些结果提示，m6A调节因子的无序表达模式可能参与了LUAD的发生发展。
为了进一步探索30个m6A调节因子的潜在关系，分别用Pearson系数和PPI网络分析研究了它们的RNA表达关系和蛋白质相互作用。除了METTL16之外，大多数m6A调节剂在肺组织的转录水平上与至少一种其他调节剂具有相关性(|r|>0.3，P<0.001)，并且最强的正负共表达是METTL3-RBM15B(r=-0.76，P<0.001)和YTHDF1-RBM15B(r=0.82，P<0.001)(图2E)。此外，大多数m6A调节剂彼此之间具有频繁的蛋白质相互作用，特别是WTAP，METTL14，ALKBH5和HNRNPC(图2F)。
LUAD组织中m6A修饰的LncRNA微阵列数据
实验表明，m6A调节因子在LUAD中异常表达，推测这些调节因子调控的lncRNA也可能发生特征性改变。在12,496个lncRNA特异性探针中，微阵列杂交并鉴定了总共2846个m6A修饰的lncRNA转录物，其来源分布主要是外显子重叠(48%)，基因间(24%)，天然反义(17%)和双向(13%)(图3A)。2142m6A修饰的lncRNA转录物(≤3000bp)的长度分布大多为500至1000bp(29%)，而其他长度范围没有显着比例差异，均为10%至17%(图3B)。M6A修饰的lncRNA转录物定位于所有染色体，染色体1和Y分别具有最高(13.5%)和最低(0.5%)丰度(图3C)。此外，来自正义和反义链的转录本数量没有明显差异(图3C)。
总的来说，123个高甲基化和20个低甲基化的lncRNA转录物在LUAD肿瘤组织中具有不同的m6A甲基化水平(|log2FC|>0.585，P<0.05)(图3D-E)。143个差异修饰的lncRNA转录物的表达分析表明，32个上调，48个下调(|log2FC|>1，P<0.05)(图3F-G)。
综合分析LUAD组织中的m6A修饰谱和表达谱
利用UpSet图将143个差异修饰的lncRNA转录物根据修饰和表达的改变分为六类(图4A)。我们发现“hyper–notdifferent”是最多的类别(N=54)，其次是“hyper–down”(N=42)，“hyper–up”(N=27)，“hypo–notdifferent”(N=9)，“hypo–down”(N=6)和“hypo–up”(N=5)。分别在正常组织、肿瘤组织以及正常组织-肿瘤组织中进行了143个lncRNA转录物的表达水平和m6A水平之间的Spearman相关分析。结果表明，在LUAD中，m6A甲基化水平较高的m6A修饰的lncRNA可能表达水平较低(图4B)。随后，选择4个差异m6A修饰的lncRNAs（SOCAR，PCAT19，SNHG8和COLCA1）进行RT-qPCR验证，以评估微阵列数据结果的准确性。如图4C所示，临床组织样品显示在LUAD中SNHG8的表达增强，SOCAR和PCAT19的表达降低，COLCA1的表达未改变。此外，MeRIP-qPCR测定验证了SOCAR和SNHG8转录物被高甲基化，并且PCAT19和COLCA1转录物在LUAD肿瘤组织中被低甲基化(图4D)。SNHG8和SOCAR的实验验证结果与微阵列数据的结果一致，表明微阵列数据具有准确性，其分析具有参考意义。这些结果提示lncRNA表达模式的紊乱与lncRNA上m6A修饰的异常改变有关。
M6A调控的LncRNA预后特征的构建与评价
对TCGA-LUAD中m6A修饰的lncRNA和30m6A调节剂进行共表达分析，将215个lncRNA定义为m6A调节的lncRNA(图5A)。然后筛选出13个与单变量Cox回归分析设定的TCGA-LUAD患者的OS相关的m6A调节的lncRNAs(P<0.01)(图5B)。图5C显示了13m6A调节的lncRNA和m6A调节因子的共表达。使用LASSOCox分析对13个m6A调节的lncRNA进行分析，(图6A，B)。获得了一个包含六个lncRNAs的m6A调控lncRNA预后特征，以及每个lncRNAs的系数，并计算出了风险评分。(图6C)。然后，采用RT-qPCR方法对LUAD临床组织标本进行检测，如图6D所示，LUAD肿瘤组织中OS危险因素AL590666.2(HR=1.232,95%CI:1.104-1.375)的表达增加(P<0.05)，而保护因子LINC01936的表达(HR=0.721,95%CI:0.567-0.917)降低(P<0.05)。
为了评估预后风险相关特征的表现，TCGA-LUAD患者根据风险评分的中值被分配到低风险和高风险亚组。KM生存曲线显示高风险亚组的生存概率较低(训练组P<0.001;验证组P=0.005)(图6E，H)。ROC曲线显示，1年、2年和3年OS的AUC值均大于0.66(图6F，I)。此外，风险评分分布图和散点图显示，高风险亚组的生存时间短于低风险亚组(图6G，J)。聚类热图显示了两个亚组之间6个预后相关的m6A调节的lncRNA的表达差异(图6G，J)。
特征的独立性评估和基于临床特征的分层分析
单变量Cox回归分析显示，风险评分可以预测不良OS以及AJCC阶段(图7A，B)。多变量分析同时证实了风险评分和AJCC分期是独立预后因素(图7A，B)。
采用分层分析来确认6-MRlncRNA预后特征的预后价值。将TCGA-LUAD患者根据年龄（≥65岁和65岁<）、性别（女性和男性）和肿瘤分期（I-II和III-IV）分为不同组，并根据中位风险评分将每组分为高低风险亚组。KM生存曲线显示，在临床特征各层中，高风险亚组的OS均较差（P<0.05；图C-H），表明风险评分可以准确区分不同临床特征人群的预后。
建造m6A诱导的CeRNA网络
在根据RNALocate记录的143个不同m6A修饰的lncRNAs的43个细胞定位中，40个lncRNAs在细胞外泌体中被检测到。4个lncRNAs只存在于细胞质中，12个只存在于细胞核中，12个同时存在于细胞质和细胞核中（图8A）。用miRcode确定了201个miRNAs与16个细胞质lncRNAs相互作用。使用三个数据库预测了miRNAs靶向的1269个mRNAs。通过TCGA和GTEx的联合分析，共筛选出4702个LUAD中差异表达的mRNAs（|log2FC|>1，P<0.05）。最终获得11个m6A调节因子、7个差异m6A修饰的lncRNAs、30个miRNAs和110个靶向mRNAs的ceRNA网络（图8B，C）。并且筛选出11个连接节点数大于或等于3的mRNAs（图8D），这些mRNAs可能被LUAD中m6A差异修饰的lncRNA诱导的ceRNA网络所调控并在LUAD肿瘤组织和对照组织中具有显著差异（|log2FC|>1，P<0.001）（图8D）。
