结合WGCNA的ceRNA分析

以往研究表明miRNA能通过与mRNA上的miRNA应答元件（MRE）结合来调节靶基因的表达。竞争性内源RNA（competing endogenous RNAs），简称ceRNA，是从功能上定义的一类RNA的总称。ceRNA具有miRNA结合位点，能够竞争性地结合miRNA，抑制miRNA对靶基因的调控。ceRNA包括mRNA、lncRNA、circRNA等。 通过进行ceRNA分析，我们能从一个更为宏观的角度来解释转录体如何构建基因表达调控网络，从而进一步挖掘基因在其中的调控机制。目前，已经证明ceRNA原理与肝癌、胃癌、乳腺癌、结肠癌、胰腺癌及膀胱癌等多种类型的癌症发展有关。

ceRNA-miRNA-mRNA之间有两个重要特征：

① 表达关系的相关性，在ceRNA体系中，由于调控机制的存在ceRNA-miRNA及mRNA-miRNA之间分别为表达负相关，ceRNA-mRNA之间为表达正相关。表达联动越强，越可能有显著的调控关系；

② 序列结合关系，miRNA通过与mRNA上MRE结合调节mRNA的表达，而ceRNA通过与miRNA结合减少miRNA与mRNA的结合，因而ceRNA-miRNA及mRNA-miRNA之间均存在序列结合关系。 通过这两个特征，我们可以用生物信息的方法进行ceRNA预测。 WGCNA其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module)，可以较好地分析基因之间表达相关性。而序列结合关系可通过ceRNA数据库或Miranda等工具预测。 接下来，我们将通过分享一篇文章展示WGCNA与ceRNA结合的生信分析思路及技巧。这篇文章PMID为31164492。发表在期刊Aging上（具有5.179的影响因子）。

文章思路介绍

这篇文章于2018年11月发表在THE LANCET（影响因子59.102）

1、数据获取

研究用到了来自GTEx的407个正常血液样本及来自TCGA的151个急性骨髓性白血病样本的RNA表达数据。 TCGA是目前最大的癌症基因信息数据库，收录了包括罕见癌型在内的30多种癌症。其数据的全面也体现在数据类型丰富上，包括临床数据、基因表达数据、miRNA表达数据、拷贝数变异、DNA甲基化、SNP 等。

TCGA数据一般通过美国国家癌症研究所（GDC，NCI Genomic Data Commons）的基因数据共享门户（DataPortal）、数据存档库（Legacy Archive）访问。GDC Data Portal的基因数据已经统一比对到GRCh38 (hg38)参考基因组。而GDC LegacyArchive基因数据的参考基因组版本为GRCh37 (hg19)和GRCh36 (hg18)。

TCGA数据根据处理级别（raw、normalized、integrated）和访问级别（controlled、open access）提供分级访问，其中完成度较高的3、4级数据不需要任何身份验证和授权就可以访问，可以满足绝大部分的研究需要。

TCGA中存储了大量癌症患者各个部位取得的癌样本与癌旁样本，然而对于白血病、脑、卵巢等特殊部位组织，通常不会对癌旁进行采集。GTEx数据库存储了大量来自健康人类捐献者的组织样本的表达数据，包含大约714个捐献者的接近11688个RNA-seq样本，共涉及53个组织、12个脑分区以及2个细胞系。可以在缺乏癌旁数据时用于对照。GTEx数据可以通过访问https://gtexportal.org/获取（更具体的数据下载步骤见文尾教程）。

2、基因差异表达分析与差异基因的功能富集

获得表达数据后，进一步对407个正常血液样本及151个急性骨髓性白血病样本的mRNA表达数据进行RNA表达差异分析，作者使用R包edgeR将表达数据标准化为TMM值并进行差异分析，共发现2667个显著上调mRNA及2456个显著下调mRNA。除了edgeR外，差异表达的常用R包还有DESeq2（同样适用于测序数据）及Limma（适用于芯片数据）。

图1.mRNA差异表达火山图

为了进一步分析上下调基因的生理功能，作者对差异基因进一步进行了GO富集分析，并用GSEA工具进行了Pathway分析。

图2.上调差异基因的GO富集结果

图3.上调差异基因的GSEAPathway分析

3、WGCNA分析基因共表达模块

接下来，作者基于WGCNA对方差较大的mRNA进行了加权共表达分析，鉴定了19个mRNA共表达模块。

图4.mRNA的WGCNA模块的聚类树状图与TOM热图

并进一步通过各个模块中的基因表达谱分析模块与性状（AML、正常）的相关性，最终挑出了与性状高度相关的两个mRNA模块（cyan、turquoise）。

图5.mRNA的WGCNA模块与性状相关性热图

并同样对方差较大的lncRNA进行了WGCNA加权共表达分析，鉴定了8个共表达的lncRNA模块。通过与性状进行相关性分析，挑选出了与性状高度相关的turquoise lncRNA模块。

4、ceRNA分析与mRNA挑选

对于挑选出的与性状高度相关的lncRNA共表达模块（turquoise）中的lncRNA（2662），作者通过miRcode数据库预测其可能作为ceRNA结合的miRNA，并与TCGAmiRNA表达数据中表达量最高的400个miRNA取交集，共获得47个目标miRNA。 接下来，作者通过starBase、miRDB、miRTarBase、Targetscan等数据库获取了这47个miRNA可能的靶向mRNA，并与WGCNA中与性状高度相关的两个mRNA模块中的mRNA、显著差异表达的mRNA取交集，最终获取了111个上调mRNA和9个下调mRNA。

图6.靶mRNA、WGCNA-cyan-turquoise mRNA以及显着上下调mRNA的集合分析

5、Cox生存分析

接下来，对于交集获得的120个mRNA，作者进一步在TCGA患者中进行了单元Cox生存分析，筛选了22个显著与生存相关的基因（Pvalue<0.05，NCBI ID<15000），然后对这22个基因进行多元Cox分析，构建了能够预测患者3年预后的基因Panel模型，共包括HOXA9、INSR、KRIT1、MYB、SPRY2、UBE2V1、WEE1、ZNF711八个基因。结合这8个基因的多元回归结果计算了风险指数，根据风险指数将TCGA患者分为高风险组与低风险组，并绘制生存曲线与ROC曲线验证这8个基因的模型对预后的预测准度。

图7.构建基因panel的风险指数生存曲线

图8.构建基因panel的风险指数ROC曲线

6、最终lncRNA-miRNA-mRNAceRNA网络

对于挑选出的基因panel中的8个mRNA，作者抓取上文中对其有靶向调控作用的miRNA，再抓取能够作为ceRNA结合这些miRNA的lncRNA，将lncRNA与差异表达的lncRNA取交集，最终构建了108个lncRNA、10个miRNA以及8个mRNA的ceRNA网络。

图9.最终获得的ceRNA网络

相关生信方案介绍

与这篇文章的思路相比，我们的WGCNA分析流程更加清晰和深入，在ceRNA预测上不仅考虑了序列结合也考虑了不同种基因间（mRNA-miRNA-lncRNA）的共表达关系。在miRNA靶基因预测上，不仅可以像文章中利用数据库查找，也可以利用工具进行独立的预测。基本流程如下：

1.    数据、分组和临床信息检索和下载

2.    分析其中lncRNAs、mRNAs、miRNAs表达的情况，对差异表达明显RNA进行分析；

3.    差异表达mRNAs的GO分析、GO富集分析和KEGG分析，KEGG富集分析；于miRNA靶向预测，分析lncRNA-miRNA、miRNA-mRNA靶基因分析，并构建靶向预测的RNA调节网络，推测关键miRNAs；

4.    基于WGCNA，构建基因共表达网络，确定包含差异基因最多的模块,对该模块进行关键基因分析, 共表达找到的关键基因与4中的关键基因取交集，

5.    将4-5的结果关联，构建针对关键mRNA、miRNAs，lncRNA的ceRNAs网络，挑选对于生存分析最具显著性意义关键miRNAs，并绘制网络图；

6.    验证筛选的关键miRNAs在癌中的生存率分析。从分期、年龄、性别，miRNA的表达分析其与预后的关系；

7.    利用cox回归、lasso等算法筛选能够预测癌症病人生存预后的基因集，计算病人的风险指数，构建预测模型；

8.    风险指数模型联合临床信息进行分层的生存预后分析。

此外根据数据与疾病特征的不同，具体流程会进行合理的调整。

相关数据库教程

整个分析流程主要涉及两类数据库：①获取样本及表达组数据的TCGA、GTEx；②miRcode、starBase、miRDB等ceRNA相关的数据库。下文我们将对其中较为介绍，并着重讲解GTEx数据的获取。

1、GTEx

GTEx全称Genotype-TissueExpression Project（基因组-组织表达项目），由美国国立卫生研究院（NIH）于2010年9月推出，旨在为研究人员提供研究不同组织间基因表达情况以及基因遗传变化如何导致常见疾病的一手数据。

GTEx样本来源于生前健康的捐献者的尸体，截止今日已更新到V7版，包含大约714个捐献者的接近11688个RNA-seq样本，共涉及53个组织、12个脑分区以及2个细胞系。还包含了48个组织的表达数量性状基因座分析结果（eQTL，研究与单个基因mRNA表达量相关的DNA突变）。此外还提供了各个样本的组织学切片图像。GTEx数据通常通过GTEx portal（https://gtexportal.org/）访问。

图10.GTEx portal首页

①获取基因表达数据

在首页点击菜单栏“Datasets”-“Data Download”即可进入数据下载界面。

图10.GTEx portal数据下载界面

其中用途最广泛的为基因表达数据（上图中红框圈出的496M文件），此外，因为下载矩阵中样本编号没有注释，为了筛选需要的组织部位，同时需要下载样本注释文件（上图中红框圈出的7.9M的文件）。

表达数据下载后为gct格式，为标准的以基因为行，样本为列表达矩阵。其中第一列为基因的Ensembl编号，而第二列为注释后的标准基因名。第一行为固定注释，第二行为基因数和样本数。进行后续分析前可以通过R或其他工具删除第一列以及第一、二行。

样本注释表格下载后为txt格式，第一列“SAMPID”为与表达矩阵对应的样本编号，另一个关键列为“SMTS”为样本对应的组织部位，可以通过“SMTS”列筛选所需部位的样本，再通过R或者其他工具提取表达矩阵中对应样本的表达数据。

②探索对基因表达量产生影响的基因突变（eQTL）

eQTL 分析是单个SNP和单个基因表达的相关性分析，可以为基因在人体生理和疾病中的作用机制，基因间的相互作用研究提供思路。如果对某个基因感兴趣，可以通过GTEx查找跟这个基因表达量相关的SNP。

在eQTL portal首页点击菜单栏“QTLs & Browsers”-“Single-Tissue eQTLs ”，进入基因eQTL列表界面。

图11.GTEx portal基因eQTL列表界面

在上图①处可以选择组织部位，②处可以输入感兴趣的基因。我们以在所有组织中以“KRAS”基因为例检索eQTL数据，获得了一个表格，如下图。

图12.eQTL检索结果，KRAS的eQTLs列表

点击左上角的CSV可以下载检索获得的eGene-eQTL表格，通过注释表格中的SNP信息，或筛选组织部位可以进行进一步的深入分析。

③获取样本的组织切片图像

通过GTEx portal还可以获取样本的组织切片图像，可以作为同部位患病组织切片的正常对照。

在eQTL portal首页点击菜单栏“Sample Data”-“Histology Viewer ”，进入样本切片图像访问界面。

图13.GTEx portal样本切片图像访问界面

点击“Hide search options”可以展开详细的搜索界面，可以根据捐赠者的年龄、性别、死亡原因及组织部位等进行筛选，符合筛选条件的样本会显示在搜索界面下方的表格中（附带各种样本信息）。点击表格中的样本即可显示对应的样本切片的显微图像。点击图片上的“Aperio Image”即可下载对应源图像。

图14.GTEx portal获取样本切片显微图像

2、starBase

starBase（http://starbase.sysu.edu.cn/）由中山大学团队研发，Starbase是做lncRNA/circRNA/microRNA等研究常用的强大数据库，在谷歌学术上显示获得了2500个以上的引用。目前已更新到3.0版本。包含23个物质的大量RNA相互作用数据，并可以进行基因的功能注释、差异表达分析、生存分析和共表达分析。

图15.starBase首页

首页点击“miRNA-Target”可以根据miRNA查找对应的靶向mRNA、lncRNA、circRNA等数据。左边边栏可以基于物种、miRNA、靶基因、CLIP-seq严格等级、涉及癌症类别等进行数据筛选。页面主题表格列出了符合条件的miRNA、靶基因信息及关系对来源。上方快速检索栏显示了目前的检索条件，可以点击黄色字体部分进行修改。表格右上角点击Search可进行自定义检索，点击“EXCEL”或“TXT”可进行下载。（其他功能界面类似故不再截图）

图16.starBase中的miRNA-mRNA数据界面

点击“Degradome-RNA”访问degradome-seq数据，查找miRNA介导的剪切降解片段深度测序数据，从中筛选miRNA作用的靶基因。

点击“RNA-RNA”访问RNA-RNA相互作用组高通量测序数据鉴定的RNA-RNA相互作用网络，涉及LIGR-Seq, PARIS, SPLASH等技术数据。

点击“ceRNA-Network”访问通过分析CLIP-seq数据得到的的ceRNA网络数据。

点击“RBP-Target”访问通过CLIP-seq数据获得的RNA结合蛋白与RNA结合的靶向关系对。

点击“RBP-Motif”访问通过CLIP-seq数据获得的各种RNA结合蛋白结合的特征蛋白序列（Motif）。

点击“RBP-Disease”可以查找RBPs与44种组织和366种疾病中基因的体细胞突变之间的潜在相关性，包括癌症和罕见疾病。

点击“Pathway”可以查找miRNA、RNA结合蛋白、RNA及ceRNA基因相关的KEGG通路。

点击“Pan-Cancer”可以在线进行差异表达分析、生存曲线分析、共表达分析等等。

图17.starBase中的癌症数据分析界面

3、miRcode

miRcode（http://www.mircode.org/index.php）由哥德堡大学生物医学研究所研发，提供全转录组规模（包含lncRNA）的人类miRNA靶基因预测数据，基于hg19版本的基因组, 转录本注释来源于Gencode v11版本，并使用TargetScan6工具。

图18.miRcode首页

在miRcode检索界面可以直接输入miRNA或者靶基因名字进行检索，并可以基于靶基因类型、靶基因位置进行筛选。也可以点击“Download”直接访问数据下载界面。

图19.miRcode下载界面

可以根据需要下载高置信度的预测结果，或者中等置信度的预测结果，数据表格较大，需要注意打开的方式。

4、miRDB

miRDB（http://mirdb.org/）时miRNA靶基因预测以及功能注释的在线数据库。miRDB中的靶基因预测基于MirTarget。miRDB包括了人类、小鼠、大鼠、狗和鸡五种物种的miRNA靶基因预测数据，用户也可以自己提供数据进行预测。

图18.miRDB首页

miRDB可以通过首页输入miRNA或者靶基因名字进行检索，也可以点击“Data Download”下载所有预测关系对，再根据自己的需要进行筛选。