甲基化芯片数据下载如何读入到R里面

数据是一切的开始，前面我们介绍了一些背景知识，主要是理解什么是DNA甲基化，为什么要检测它，以及芯片和测序两个方向的DNA甲基化检测技术。具体介绍在：甲基化的一些基础知识，也了解了甲基化芯片的一般分析流程 。既然要开始甲基化芯片数据挖掘实战，那么首先要有数据咯！需要区别的是甲基化芯片样本的idat原始文件，以及甲基化信号值矩阵。前面我们介绍了如何在GEO里面下载甲基化数据，拿到的数据文件必须要导入到R里面才能分析，现在我们就讲一下不同数据如何导入R里面。

首先你需要成功下载哦。如果你相信作者对他自己的甲基化芯片数据处理，就可以直接使用其 _series_matrix.txt.gz 存储的甲基化信号矩阵。如果你不相信作者，就下载他上传的idat芯片原始数据，然后自己走minfi或者champ流程，自己拿甲基化信号矩阵走下游分析。

网速好就可以使用GEOquery可以直接下载甲基化信号值矩阵

如果你网速非常好（比如海外用户），使用GEOquery可以直接下载甲基化信号值矩阵，取决于你是否相信作者对芯片原始数据的处理。代码很简单，如下：

require(GEOquery)    
require(Biobase)
GSE80559 <- getGEO("GSE80559")
beta.m <- exprs(GSE80559[[1]])
# 这样就拿到了甲基化信号值矩阵，而且在R环境里面啦。

再次强调，这个方法适用于数据集的研究者处理好了idat芯片原始数据，而且处理的格式符合要求哈。大概率上，你还是得自己去下载idat芯片原始数据走minfi流程的。

其实就是使用了这个数据集存放在GEO里面的 _series_matrix.txt.gz 文件而已，这个文件直接读入到R即可，没什么好说的了。所以你可以找朋友帮你下载好 _series_matrix.txt.gz 文件，存放在当前目录，使用getGEO指定当前目录，这样的话，这个getGEO函数就会读取你下载好 _series_matrix.txt.gz 文件，而不是重新下载！

GSE75679 <- getGEO("GSE75679",destdir = './')
beta.m <- exprs(GSE75679[[1]])

再次强调，你在中国大陆，基本上不可能下载成功这个 _series_matrix.txt.gz 文件，我对中国大陆严峻的科研环境是深恶痛绝！读取本地文件的log日志如下：

> GSE75679 <- getGEO("GSE75679",destdir = './',AnnotGPL = F)
Found 1 file(s)
GSE75679_series_matrix.txt.gz
Using locally cached version: .//GSE75679_series_matrix.txt.gz
|=================================================================| 100%  231 MB
Parsed with column specification:
cols(
  .default = col_double(),
  ID_REF = col_character()
)
See spec(...) for full column specifications.
|=================================================================| 100%  231 MB
File stored at: 
.//GPL13534.soft

这个时候，你关注的数据集的甲基化信号值矩阵，就被加载到R里面啦。后面我们再介绍后续处理。

其实 ChAMP也可以直接导入这个 _series_matrix.txt.gz 文件代表的甲基化信号值矩阵哦！

然后如果下载了芯片的idat原始文件

可以使用minfi包的read.metharray.exp函数读取，你前面下载的该数据集的RAW.tar 里面的各个样本的idat文件，就被批量加载到R里面啦。（必须注意的的是， 下面代码里面GSE68777/idat文件夹里面有idat文件哦！！！）

library("minfi")
rgSet <- read.metharray.exp("GSE68777/idat")
rgSet 
save(rgSet,file = 'GSE68777_minfi_rgSet.Rdata')
# 你需要学习 rgSet 所表示的对象 class: RGChannelSet ，才能进行后续处理

也可以使用The Chip Analysis Methylation Pipeline，但是不得不说，每次安装 ChAMP 都得脱一层皮，它的依赖包实在是太多了。其中一个ChAMPdata_2.18.0.tar.gz就是680M文件。首先可以读取R包自带的芯片的idat原始文件，如下

# BiocManager::install("ChAMP",ask = F,update = F)
library("ChAMP")
testDir=system.file("extdata",package="ChAMPdata")
myLoad <- champ.load(testDir,arraytype="450K")

这个数据集就在ChAMPdata_2.18.0.tar.gz就是680M文件，很可怕！champ.load函数作用于 testDir 这个文件夹，我们就去看看这个 testDir 文件夹里面的内容，发现除了idat芯片原始数据之外，还有一个csv文件。

如果你要读取自己下载好的idat芯片数据，合理的思维迁移，你就应该明白， 你也需要制作自己的csv文件啦。当然了，这个csv文件的规则， 可以去看champ包说明书，也可以自己揣测。

那么，我们自己的GSE68777/idat文件呢，其实也不小，如下

7.7M Feb  7 23:21 GSM1681154_5958091019_R03C02_Grn.idat
7.7M Feb  7 23:21 GSM1681154_5958091019_R03C02_Red.idat
7.7M Feb  7 23:21 GSM1681155_5935446005_R05C01_Grn.idat
7.7M Feb  7 23:21 GSM1681155_5935446005_R05C01_Red.idat
7.7M Feb  7 23:21 GSM1681156_5958091020_R01C01_Grn.idat
7.7M Feb  7 23:21 GSM1681156_5958091020_R01C01_Red.idat
7.7M Feb  7 23:21 GSM1681157_5958091020_R03C02_Grn.idat
7.7M Feb  7 23:21 GSM1681157_5958091020_R03C02_Red.idat
7.7M Feb  7 23:21 GSM1681158_5935403032_R05C01_Grn.idat
7.7M Feb  7 23:21 GSM1681158_5935403032_R05C01_Red.idat
7.7M Feb  7 23:21 GSM1681159_5958091019_R04C02_Grn.idat
7.7M Feb  7 23:21 GSM1681159_5958091019_R04C02_Red.idat
7.7M Feb  7 23:21 GSM1681160_5935446005_R06C01_Grn.idat
7.7M Feb  7 23:21 GSM1681160_5935446005_R06C01_Red.idat
7.7M Feb  7 23:21 GSM1681161_5958091020_R02C01_Grn.idat
7.7M Feb  7 23:21 GSM1681161_5958091020_R02C01_Red.idat
7.7M Feb  7 23:21 GSM1681162_5958091020_R04C02_Grn.idat
7.7M Feb  7 23:21 GSM1681162_5958091020_R04C02_Red.idat 

可以看到文件名是有规律的。然后样本名对应的是不同的分组，需要自己仔细看该数据集的文献啦。

总之，你需要耗费至少半个小时去理解如何制作自己的csv文件，以及理解你想要挖掘的数据，然后才有可能使用champ读取那些idat挖掘咯。

champ能够自动化处理你所有的数据~

Generating beta Matrix
  Generating M Matrix
  Generating intensity Matrix
  Calculating Detect P value
  Counting Beads

后面我们再介绍后续处理。

如果是TCGA数据库的甲基化芯片数据

通常呢，tcga数据库的样本数量很大，而idat芯片原始文件太大，所以一般就直接下载甲基化信号矩阵即可。

通常我建议大家在UCSC的XENA浏览器下载。

Level 3 DNA methylation data of GC samples evaluated on the Illumina Infinium HumanMethylation450 platform (450K array) which assesses 482,421 CpG sites throughout the genome were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/)

These data consist of pre-processed data via TCGA pipelines in the form of β values, which are a ratio between methylated probe intensities and total probe intensities. Probe-level data were condensed to a summary beta value for each gene using the Methylation450_single_value function in TCGA-Assembler, which calculates the aver- age methylation value for all CpG sites associated with a gene.

甲基化信号矩阵文件非常大，如果全部的tcga的1万多个样本，文件是34G, 通常大家没必要做pan-cancer研究啦，但是下载其中一个癌症也是不小，几个G的文件，读取到R里面到没有问题。大家可能更关心的是这个甲基化信号矩阵如何被minfi或者champ读取成为对象。因为你不想重复造轮子，想使用minfi或者champ大量的质控函数，统计可视化函数，就必须把你的数据搞成为minfi或者champ的对象！

数据文件导入R之后呢？

其实目前甲基化信号值矩阵差异分析的R包非常多，比如 IMA, minfi, methylumi, RnBeads and wateRmelon，以及我们一直强推的ChAMP！不过我没有那么多时间去一一解读啦，我相信minfi或者champ就够用了。你的目的其实是做完甲基化信号值矩阵有3个层次的差异分析：DMP，DMR，DMB：

• DMP代表找出Differential Methylation Probe（差异化CpG位点）
• DMR代表找出Differential Methylation Region（差异化CpG区域）
• Block代表Differential Methylation Block（更大范围的差异化region区域）