R语言实例操作分析GEO数据库甲基化芯片

本文首发于“百味科研芝士”微信公众号，转载请注明：百味科研芝士，Focus科研人的百味需求。

小伙伴们，上次为大家解读了一篇GEO甲基化芯片相关的SCI文献(Aberrantly methylated-diferentially expressed genes and pathways in colorectal cancer)，今天，小编打算带领大家用R软件实例操作分析GEO甲基化芯片。作为目前最大的芯片数据库，GEO数据库提供给我们了海量的数据，但是，错综复杂的数据交织在一起，如何选择数据是摆在我们面前最重要的问题，读完今天这篇文章，我相信大家都能学会GEO甲基化芯片的分析。下面，就和大家一起跑一遍R，希望大家喜欢这篇文章！

首先是GEO甲基化芯片的下载和预处理。进入GEO数据库主页(网址https://www.ncbi.nlm.nih.gov/geo/)，可以直接打开链接进入，也可以通过NCBI官网的GEO DataSets或GEO Profiles进入（进入NCBI数据库后下来搜索框）。

主页搜索框输入关键词“Methylation”，点击search，出来两行英语，一般选择第一行的数字，点进去。

页面跳转到GEO DataSets，这和NCBI中直接进入是一样的，直接输入需要检索的肿瘤，或者如下图点击Advanced后，在高级检索中分别输入“Colorectal cancer”和“Methylation”，根据自己的需要选择合适的甲基化芯片。如果大家知道芯片的GSE号，也可以直接根据GSE号来检索。

例如，本帖根据研究内容选择了“GSE29490”这张芯片。

点击芯片标题，则可弹出这张芯片的全部注释信息

下拉该页面，可以看到该芯片的探针平台信息，样本信息，以及矩阵文件(名字Series Matrix File(s)的TXT格式)和原始文件(TAR格式的文件)。一般选择下载矩阵文件，如果下载原始文件，需要我们自行整理文矩阵文件，还是比较麻烦的！

点击Series Matrix File(s)后，选择文件的路径点击保存。

对下载好的矩阵文件解压，使用EXCEL表格打开，如下图，其中感叹号开头的是注释文件，将其删除。

将注释文件删除后，把EXCEL里面的矩阵粘贴在txt文档里，命名为M.txt。

同时在EXCEL中建立两列以sample,group分组的表格，对样本进行分组，C为肿瘤组，T为病例组。将表格粘贴在txt文档里面，以“group.txt”命名。

将M.txt和group.txt放在同一文件夹里，文件夹名就叫M吧，文件的准备已告一段落。

接下来，就可以打开Rstdio了(和R软件运行一样，本质也是R软件，只是界面不同)，做分析之前，需要安装甲基化芯片相关的包，这个过程一般比较慢，大概2h。而通过以下的代码就可以实现安装。对于这些包的安装和功能，我们可以参考Bioconductor网站。

分析GEO甲基化芯片需要安装如下包：

加载安装包，设置工作目录(注意R软件中工作目录需要用”\\”或“/”，不可用“\”来设置工作目录)，工作目录直接设置M文件夹的路径，然后读取txt文件。

对数据进行标准化处理并输出结果

我们可以看出，标准化处理之前的箱线图的中位值未处于同一水平线上。

而标准化处理后的箱线图中所有样本均处于同一水平线上，使各种实验条件下的测量可以相互比较，消除测量间的非实验差异。这样所有的样本就具有可比性了。

对芯片进行质量控制(QC)，这里需要输出两幅质量控制的图：DensityBean图，MDS图。

DensityBean图:图中可以看出，峰值主要出现在0附近，说明这张芯片的甲基化水平较低。

MDS图：选出样本中1000个变异最大的位点，观察肿瘤组和对照组样本分布情况，在我们实验中，对于这种分布明显偏离组内其他样本的样本要予以剔除，因为这样的样本对实验结果影响较大。

甲基化位点的差异分析

dmpFinder函数，其参数设置如下，对于这些函数的设置在R软件中输入“?dmpFinder”,运行代码即跳转到参数设置页面，非常方便。

打开文件夹，即可看见dmpDiff命名的甲基化位点的表格，而我们挑选的差异甲基化位点的值则是根据q-val＜0.05来判定的。

甲基化差异位点制作热图：这里是根据M值(甲基化的率)来做差异甲基化位点的热图。

甲基化差异区域分析

运行代码后，得到名为dmrs的表格。

我们在做甲基化差异区域注释时候需要五列数据：chr，start，end，Ref，Alt。其中Ref，Alt这两列数据缺失，需要手动用0补充。

补充好了之后，就可以对甲基化芯片差异区域进行注释了。

甲基化区域的注释，我们采用wANNOVAR网站(http://wannovar.wglab.org/)，进入主页后需要输入邮箱(机构邮箱)，这里小编也是在网上随便机构邮箱，这个网站有点不地道，不支持个人邮箱！当然，只是为了使用这个网站，不需要邮箱来接受消息。Sample identifer栏随便输入英文名称即可。将上述准备好的5列数据（chr，start，end，Ref，Alt）粘贴在第三个检索框里。下拉页面，在Input Fomat栏里选择ANNOVAR，然后点击Submit提交数据。

大概几分钟的时间，注释结束，便可下载我们刚刚注释的甲基化区域了，其中第一行为外显子区域结果，第二行为整个基因上面的注释。

打开文件，我们可以看到第6列是甲基化位于基因上的位置，第7列是甲基化区域所在的基因，第8行是对甲基化基因的注释。到这里，我们针对一张GEO甲基化分析已完成。

当我们得到这些甲基化的基因后，可以对这些基因进行GO、KEGG富集分析，生存分析等，或者去联合基因表达的芯片进行分析。作为最大的芯片数据库，因GEO芯片来自用户的上传，我们需要对芯片的质量做质控，并且做标准化的处理，以尽可能消除对因实验条件不同造成的实验误差，这一点在TCGA数据库中则不需要。不需要花费经费，只需要大家动脑动手就能发SCI，哈哈，这样的GEO芯片来一打！