纯生信分析套路 胃肠道癌m6A研究
胃肠癌PI3K/Akt/mTOR信号通路中m6A RNA修饰调节的研究
研究概述| 腺苷(m6A) N6位点的甲基化是哺乳动物蛋白编码mRNA中最常见的RNA修饰，是一种具有多种重要生物学功能的可逆修饰。m6A在胃肠道(GI)癌中的潜在修饰机制尚不清楚。本文对胃肠道中已知的9种m6A writer、erasers和reader蛋白进行了全面的分子层面分析。
数据概览| TCGA数据库 + GEO数据库
分析方向|用cBioportal进行基因改变及通路分析。STRING分析m6A调控因子及其相关通路蛋白网络，MEGA7构建系统发生树，SRAMP预测m6A修饰位点。m6ASNP分析m6A相关SNP。此外实验部分使用m6A-seq、real-time PCR和phospho - mapk阵列验证m6A在其通路的调控。
下面将以具体方法和数据结果阐述：
到目前为止，在7000多个哺乳动物基因的转录本中，超过12000个m6A位点已经通过甲基化RNA免疫沉淀测序(MeRIP-Seq)进行了鉴定。一些关键研究表明，在胶质母细胞瘤、肺癌、肝癌、乳腺癌等肿瘤中，核心基因的改变影响着肿瘤的发生、癌细胞的增殖、肿瘤微环境和预后。然而，在大多数人类癌症中，这些m6A调控因子的生物学功能和关键靶基因仍不清楚。
近年来，SNP，表观修饰遗传等作为研究对象的胃肠道肿瘤的重要调控因子逐渐被发现，但是m6A在恶性肿瘤中的生物学功能和机制的研究进展还很少。所以，本文研究了五种GI亚型(食道、胃、结直肠、肝脏和胰腺)与m6A writer、reader和eraser蛋白之间的相关性。
首先，从TCGA数据库(http:// cancergenome.nih.gov)获取1920/10037(19.1%) GI癌症样本，共有9,628个样本用于通路改变分析。并同时获取基因表达数据和蛋白富集数据（http://www cbioportal.org/）。SNP数据从千人基因组中获取美国,东亚洲和欧洲人的SNP。
基于TCGA数据采用Reverse Phase Protein Array (RPPA) 法收集蛋白组数据，用cBioportal对27种人类癌症进行了分析。富集蛋白采用Log2比值的对数比值和查询事件结果p值<0.05来确定。另外使用DAVID推断靶基因和KEGG途径，p value和富集数表示该通路水平是否显著。
在cBioportal 和Prism6中计算TCGA数据中基因共表达和DNA甲基化的相关性使用皮尔森相关系数，为了明确目标基因的系统发育特征，从GenBank中征集了66条全基因组序列，构建了系统发生树。Neighbour-Joining分别采用MEGA7中的NJ算法，随机添加分类单元和树分重连接分支交换。
SRAMP作为m6A位点预测工具，基于三种类型的随机森林分类器，即二进制(核苷酸序列的位置二进制编码)，KNN (k近邻)和光谱(核苷酸对光谱)最终通过加权求和，确定其结果。另外注意假阳性样本纠正，避免显著的偏倚，保持1:10的正阴性比例。对于具有相同基因座的多个转录本，只有m6A数量最多的最长转录本被保留。
----------------
根据上述方法步骤，首先可获得GI cancer样本的writer、erasers和reader分布情况，如下图一所示。比较了m6A相关基因的表达模式，对1,782例癌症样本和148例正常对照样本进行热图分析，发现m6A调控因子在胃肠道中表达显著升高。进一步比较五种类型胃肠道肿瘤与正常样本的差异(食道、胃、结直肠、肝脏和胰腺)(图一C)显示，与正常对照组相比，除胰腺癌外，胃肠道患者的阅读和写作水平大多上调。cBioPortal进一步分析了m6A调节因子对癌症遗传因子的影响，发现约11.9%的临床样本(204/1710)存在m6A调控因子的基因扩增，近4.6%的样本(79/1710)存在基因组突变。其中，m6a结合蛋白(reader)， YTHDF1, YTHDF3，和HNRNPA2B1拷贝数扩增频率较高。
图一，m6A writer、erasers和reader在GI cancer中的全面分子谱分析
接下来，获得临床相关的验证，对TCGA数据库中五种亚型的GI癌症m6A调控因子进行OS分析(如下图二)。Kaplan-Meier分析显示，m6A调控因子的差异表达与OS显著相关。在食管癌和结直肠癌中，METTL3的高表达与较差的生存率显著相关(p<0.05)。相比之下，在胃癌、肝癌和胰腺癌中，METTL3低表达与较短的生存期显著相关。在结肠直肠癌中，高表达的ALKBH5与低OS显著相关，而在胃癌和胰腺癌中，它与较好的生存期相关。综上所述，这些结果表明m6A调节因子可能对胃肠道肿瘤的生存有重要影响。
图二，m6A甲基化调节因子的改变与GI患者总体生存率的关系
然后进行GI 和其他癌症中蛋白表达的m6A通路分析，可见下图三，其中热图展示了41种通路途径，磷脂酰肌醇-3-激酶(PI3K)/Akt，哺乳动物雷帕霉素靶蛋白(mTOR)，叉头盒转录因子(FoxO) O类，催乳素，雌激素，胰岛素和ErbB信号通路，通常存在于27种主要类型的人类癌症中。与m6A修饰在GI中呈正相关的5个通路(PI3K/Akt, mTOR, FoxO，丝裂原活化蛋白激酶[MAPK]和p53)与KEGG癌症通路重叠。
图三，与m6A失控调控相关的信号通路
接下来，为了评估肿瘤中PI3K/Akt、FoxO、mTOR、MAPK和p53信号通路的改变，作者通过cBioPortal对TCGA中27种不同类型肿瘤的患者(n= 9628)体细胞突变和DNA拷贝数变化(扩增、浅缺失和深缺失)进行了肿瘤基因组学研究。对于每一种癌症类型，如果一个或多个基因参与某一通路，则至少有一种改变的靶基因被认为在该通路中发生了改变。如下图四所示，27种癌症中，PI3K/Akt通路的改变频率最高(75%)。该通路中消化道肿瘤亚型的改变频率分别为食管癌(95%)、胃癌(83%)、结直肠癌(77%)、肝癌(82%)、胰腺癌(80%)。mTOR通路在胃肠道中也有较高的改变率(75%)，在五种癌症亚型中分别为95%、74%、62%、86%和70%。并用热图分析其改变频率。PIK3CA、PIK3R1、AKT1、PTEN、BCL2和SYK是PI3K/AKT通路中最常见的突变基因。
图四，选择的信号通路和m6A调节因子之间的关系
之后，为了了解这5种通路在肿瘤和正常样本中是否受到DNA甲基化的调控，我们使用DiseaseMeth version 2.0基于TCGA Human- Methylation450的数据集进行了进一步的分析。系统图树显示，AKT1和YTHDF1 (reader)之间有很强的关系，表明它们可能是一个共同的祖先或相同的进化起源点。此外，进行了m6A位点预测。特殊DRACH (D = A, G或U;R = A或G;可以将H = A, C或U)和GAC consensus motif识别为m6A位点，通过SRAMP的两种随机森林预测模式进行计算m6A位点普遍存在于PI3K/Akt的组成分子中(AKT1/2/3, PI3KCA，PTEN, PI3KR1/2, BCL2, BCL2L1, IRS1, GSK3B, MYC)和mTOR通路(mTOR, AKT1S1, TSC1，RICTOR, RPTOR, RPS6KB1, PRKCA)，以及来自FoxO, MAPK, p53途径分子。
综上所述，这些结果表明m6A在调控PI3K/Akt/mTOR等肿瘤通路中发挥着基础性作用。
-------------------
至此生信研究结束，而作者为了更有说服力的结果支持生信发现，进行了敲除实验。下面是实验方法部分：
为了证实m6A在人类癌症中的作用，采用与免疫沉淀实验不同输入FASTQ文件的MeT-DB v2.0生物信息学工具。U2SO细胞(GSE48037)和Hela细胞(GSE46705) 、HEK293T细胞(GSE55572)和A549细胞(GSE76367)的METTL3敲低后表达谱高通量测序数据。通过JBrowse轨迹筛选面板筛选感兴趣的靶基因。在Origin9内重建并计算了目标基因外显子峰区输出。
m6ASNP使用来自1000基因组计划数据库的遗传变异文件，预测m6A变异位点。根据文献，共有55,548个人类m6A位点。人m6A相关的高置信度为13,703，中置信度为54,222，低置信度为243,880。为避免偏差，使用5种软件将各变量的有害程度从0 ~ 5排序，得分越高，表明危害可能性越大。
在4℃条件下孵育2 h，得到10μg的亲和纯化抗m6A多克隆抗体(突触系统)放进450ml的免疫沉淀(IP)缓冲液(150 mMNaCl, 0.1% NP-40 (v/v)， 10mM Tris-HCl, pH 7.4)和300 U mL−1 RNase抑制剂(Promega)进行m6A-seq coverage分。并进行哺乳动物细胞培养，转染对照siRNA，48 h后，Q-PCR检测METTL3、AKT1、mTOR、PIK3CA、PTEN mRNA水平。之后进行m6A甲基化定量，测量mRNA半衰期。
然后，让我们看看作者的验证结果。
为了研究m6A对PI3K/Akt/mTOR通路的直接影响，作者使用MeT-DB V2.0用于检测基于GEO数据库的基因组表达谱。由于m6A的存在可诱导真核生物的mRNA降解，因此我们重点研究了基于GEO数据库的m6Awriter下调的上调基因。METTL3在U2SO细胞(GSE48037)、Hela细胞(GSE46705)、HEK293T细胞(GSE55572)和A549细胞(GSE76367)导致PI3K/表达相对升高Akt/mTOR信号通路蛋白:AKT1, AKT2, PTEN, PIK3CA, BCL2如下图五所示。在AKT1中鉴定了具有dbSNP注释的SNPs (rs185142105、rs189600949、rs5027535、rs113244357、rs192598160和rs111541557)突变体，这些突变体在相同位点或不同位点干扰m6A修饰。
图五，m6A对PI3K/Akt/mTOR通路分子的调控
为了在转录组水平上鉴定和定位m6A位点，作者从人胃癌细胞株(TMK1)中纯化RNA进行m6A-seq。经鉴定，mTOR、AKT1、PIK3CA和PTEN具有特定的m6ARNA甲基化序列元件。在mTOR和PIK3CA中，3 ' UTR和相邻段在m6A峰中富集最多，如下图六所示。AKT1 m6A主要修饰靶密码子区域。
图六， PI3K/Akt/mTOR通路RNA修饰的m6A-seq捕获
最后，作者通过甲基转移酶(METTL3-METTL14复合物)的活性来验证m6A的解除是否会导致GI癌症中PI3K/Akt/mTOR信号通路的改变抑制s -腺苷同型半胱氨酸(SAH)或METTL3 siRNA，并抑制去甲基化酶的活性(FTO)与甲氯灭酸(MA)活性。PTEN、AKT1 PIK3CA,和mTOR的信使RNA和蛋白质水平显著上调，在SAH治疗GI细胞下,其中包括胃癌细胞TMK1,结直肠癌细胞系HCT116、食道癌细胞关冲,胰腺癌细胞PANC-1肝癌细胞系Huh-7， TMK1和HCT116细胞的敲除METTL3情况。此外，甲基转移酶或去甲基酶诱导的m6A修饰的代偿反馈可以尝试用4个基因间动态变化解释。为了检验m6A是否调节PI3K/Akt/mTOR信号通路组分的降解，在SAH或MA处理48小时后，使用Act D处理TMK1和HCT116。在TMK1细胞中，信号密度分析显示磷酸化在SAH和MA治疗后，AKT1、GSK3 inhibitor / akt、p38 inhibitor和mTOR分别明显被抑制或激活，如下图七所示。
图七，m6A对PI3K/Akt/mTOR通路调控的验证
到这里，本文研究就已经结束了。有分析意向（http://gaptechsxr.mikecrm.com/1vdMmqy）生信人WX公众号
文章从公共数据入手，着手研究甲基化对于GI患者的总OS的影响，并实验证明使用药物可抑制甲基化产生。首次证明m6A修饰通过调节GI中关键激酶，直接调节PI3K/Akt和mTOR信号通路的活性。有助于肿瘤的早期发现和治疗策略的制定。
所以，如果你有类似的想法，不如现在就尝试一下，说不定，就有新的发现。
