专栏首页生信技能树甲基化芯片数据下载如何读入到R里面

甲基化芯片数据下载如何读入到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区域)

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-02-08

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • TCGA中GBM的RNA-seq和甲基化数据整合分析实践

    首先进入TCGA下载数据GBM的RNA-seq和甲基化数据,从下表可见GBM共有172套RNA-seq数据以及437套DNA甲基化数据,由于TCGA提供Infi...

    生信技能树
  • 甲基化信号值的差异分析也许不应该是看logFC

    但是差异分析大家还是首先limma,而limma这个包本来是针对log后的表达矩阵设计的,这样的话,如果我们的输入是甲基化信号矩阵,实际上出来的结果是有问题的。

    生信技能树
  • 甲基化相关的习题背景补充

    有学徒表示虽然看了我在B站免费分享的视频课程《甲基化芯片(450K或者850K)数据处理 》,详见:免费视频课程《甲基化芯片数据分析》,但是课程过于强调实操,很...

    生信技能树
  • 从 Anemometer BUG 到 FRM 文件的恢复

    最近深深体会到,目前的发展速度,数据库方面各种东西,原理层出不穷,一个礼拜不去看那些公众号去“滋养”,一下脑子,就发现新的概念不知道了。

    AustinDatabases
  • Lepus 天兔数据库监控

    Lepus是一套开源的数据库监控平台,目前已经支持MySQL、Oracle、SQLServer、MongoDB、Redis等数据库的基本监控和告警(MySQL已...

    小手冰凉
  • “我不碰比特币,我只做区块链”

    你可以在世界上几乎所有银行都听到这些话。错!讽刺的是,比特币区块链是金融机构最重要的区块链。每天,银行都会在比特币交易所和商户之间进行数百万美元的交易。一些银行...

    anthlu
  • 关于Python的JSON

    实际上json就是python字典的字符串表示,但是字典作为一个复杂对象是无法直接转换成定义它的代码的字符串,python有一个叫 simplejson的库可以...

    py3study
  • solidity智能合约开发工具Atom及其插件安装

    Atom简介 Atom代码编辑器支持Windows、Mac、Linux三大桌面平台,完全免费,并且已经在 GitHub 上开放了全部的源代码。 开发团队将Ato...

    程序新视界
  • 自学Python八 爬虫大坑之网页乱码

      Bug有时候破坏的你的兴致,阻挠了保持到现在的渴望。可是,自己又非常明白,它是一种激励,是注定要被你踩在脚下的垫脚石!

    叁金
  • 使用RIST或SRT进行实时云内容提取

    本篇是来自SMPTE 2019的演讲,演讲者是来自Net Insight的Doug Shelton和Mikael Wånggren,演讲题目是“Live Clo...

    用户1324186

扫码关注云+社区

领取腾讯云代金券