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

R语⾔实例操作分析GEO数据库甲基化芯⽚⼩伙伴们，上次为⼤家解读了⼀篇GEO甲基化芯⽚相关的SCI⽂献(Aberrantly methylated-diferentially expressed genes and pathways in colorectal cancer)，今天，⼩编打算带领⼤家⽤R软件实例操作分析GEO甲基化芯⽚。作为⽬前最⼤的芯⽚数据库，GEO数据库提供给我们了海量的数据，但是，错综复杂的数据交织在⼀起，如何选择数据是摆在我们⾯前最重要的问题，读完今天这篇⽂章，我相信⼤家都能学会GEO甲基化芯⽚的分析。下⾯，就和⼤家⼀起跑⼀遍R，希望⼤家喜欢这篇⽂章！⾸先是GEO甲基化芯⽚的下载和预处理。进⼊GEO数据库主页，也可以通过NCBI官⽹的GEODataSets或GEO Profiles进⼊（进⼊NCBI数据库后下拉搜索框）。主页搜索框输⼊关键词“Methylation”，点击search，出来两⾏英语，⼀般选择第⼀⾏的数字，点进去。页⾯跳转到GEO DataSets，这和NCBI中直接进⼊是⼀样的，直接输⼊需要检索的肿瘤，或者如下图点击Advanced后，在⾼级检索中分别输⼊“Colorectal cancer”和“Methylation”，根据⾃⼰的需要选择合适的甲基化芯⽚。如果⼤家知道芯⽚的GSE号，也可以直接根据GSE号来检索。例如，本帖根据研究内容选择了“GSE29490”这张芯⽚。点击芯⽚标题，则可弹出这张芯⽚的全部注释信息。下拉该页⾯，可以看到该芯⽚的探针平台信息，样本信息，以及矩阵⽂件(名字Series MatrixFile(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芯⽚来⼀打！百味芝⼠百味芝⼠祝您猪年⼤吉