850K甲基化芯片数据的分析

作者是生信技能树组建的表观遗传学学习小组的小组长,前面已经发过一个:

学员分享-Chip-seq 实战分析流程

本文是看到生信技能树有个450K甲基化芯片数据处理传送门,我呢,恰好不久前用一个集成度很高的ChAMP包分析过850K的甲基化芯片数据。所以,就想着把自己的笔记整理下,可以和更多的小伙伴学习交流,还有个原因可能是因为这是四月份打算学生信时,接手的第一个任务,曲曲折折好几个月才跑通流程,遇到的坑也比较多,想记录下来。

我之前分析时是参考ChAMP包的源文档,非常详细的整个流程的介绍,但是,在笔记快整理完时突然发现作者的博客也写过一篇介绍的文章,博客里写的不像源文档很官方,这里面有很多作者很直白的解释和补充,还有作者一些很深刻的思考。看了之后发现自己对很多分析理解的还不是很深刻。所以如果想学甲基化芯片数据分析的小伙伴可以以官方源文档和作者的博客为主,这篇笔记仅仅作为额外的参考吧。

Illumina甲基化芯片目前仍是很多实验室做甲基化项目的首选,尤其是对于大样本研究而言,其性价比相当高;目前在临床上应用还是很广的。这种芯片的发展主要经历了27K、450K以及850K(27K,450K,850K指能测到的CpG甲基化位点),目前积累的数据主要是450K芯片的,之后850K可能会成为甲基化芯片的主流。楼主之前写过一篇450K芯片预处理的帖子,其中很详细介绍了这种芯片的基础知识以及流程图和代码,大家可以先看看。芯片的处理流程一般就是:数据读入——数据过滤——数据校正——下游分析

数据处理一种时基于GenomeStudio(illumina开发的软件),但是只对于小样本,另一种基于R的各种package,如lumi、minfi、wateRmelon、ChAMP等。

与测序相比,芯片的处理可能对计算资源的要求不算高,主要使用的工具就是R,但是R的使用比较耗内存,尤其是处理大批量数据的时候。

Step1: 基础知识的补充

在正式分析前,我结合作业先将有关甲基化和芯片的基础知识整理了一下。

Illumina 甲基化芯片的原理及探针的设计(I型探针和II型探针)

原理:简而言之,基于亚硫酸盐处理后的DNA序列杂交的信号探测。亚硫酸盐是甲基化探测的“金标准”,不管是芯片或者甲基化测序,都要先对DNA样品进行亚硫酸盐处理,使非甲基化的C变成U,而甲基化的C保持不变,从而在后续的测序或者杂交后区分出来。

450K和850K采用了两种探针Infinium Ⅰ 和Infinium Ⅱ对甲基化进行测定,Infinium I采用了两种bead(甲基化M和非甲基化U,如图显示),而II只有一种bead(即甲基化和非甲基化在一起),这也导致了它们在后续荧光探测的不同,450K采用了两种荧光探测信号(红光和绿光)(图1)。

甲基化概述:

DNA甲基化被认为是表观遗传调控的一种方式,如Cytosine methylation (5-mC)是研究最多的,被认为是哺乳动物中常见的甲基化方式, 最近有一些研究也发现了其他形式的甲基化,如2016年Nature上发表了一篇关于鼠的胚胎干细胞的m6A(N6-methyladenine)形式的甲基化。DAN甲基化被认为对基因表达,染色质重塑,细胞分化,疾病等都有重要影响(图2)。

甲基化的检测方法:

目前甲基化检测的方法可以概括为三种:芯片、测序、免疫沉淀。具体选择何种方法主要还是根据实验目的和实验室条件了。但目前来说,甲基化芯片技术从覆盖度,检测灵敏度和价格综合考虑,还是性价比相对高的(图3)。

关于甲基化芯片常见的Glossary:

CpG island: Defned as regions 500 bp, 55% GC and expected/observed CpG ratio of 0.65. 40% of gene promoters contain islands.

CpG shelves: ~4Kb from islands.

CpG shores: ~2Kb from islands, 75% of tissuespecifc differentially methylated regions found in shores. Methylation in shores shows higher correlation with gene expression than CpG islands.

Differentially methylated regions (DMR): Cell-, tissue-, and condition- specifc differences in methylation.

Enhancer: A short region of DNA that can activate transcription and is often regulated by methylation.

Hypermethylation: Most cytosines are methylated. Hypomethylation: Most cytosines do not have 5-mC. Euchromatin and active gene promoters are hypomethylated.

Beta value:通常的甲基化衡量方法被称为“Beta”值; 等于甲基化百分比,并定义为“Meth”除以“Meth + Unmeth”。

CGI: CpG island 即甲基化岛。

因为手头的数据是850K的甲基化数据,之前也只接触过ChAMP包,所以这里就以ChAMP包介绍850K甲基化数据分析。ChAMP包是一个集成度很高的包,它包括450K和EPIC(即通常所说的850K)两套分析流程,完整的包括了数据的载入,标准化,矫正,差异甲基化和富集分析等功能(图4)。

Step2:计算机资源的准备

作业1 安装好R软件及相应的包,下载R包的说明书,整理它们的官网链接。

R的使用真的很耗内存,我有28个样本(14个control, 14个case), 之前4G内存的电脑,本地分析总时半路电脑就卡死了。所以最好配置高一点,或者在服务器上下载安装R和Rstudio(这里最好安装Rstudio, 因为ChAMP包中有很多的GUI图形功能,Rstudio可以更好实现,或者含有X11功能的linux系统)。

软件的安装:

R和Rstudio 的本地安装很简单,直接到官网下载,只要注意安装时的路径不要有中文,Rstudio安装前要先安装R。

服务器版本的Rstudio安装好后,在网页地址栏输入访问地址:<服务器IP:8787,用户名和密码为Linux用户的用户名和密码。

具体安装方法可以参考生信宝典陈老师的一篇文章http://www.biotrainee.com/thread-1808-1-1.html。

下载R包:

下载ChAMP 包,官网给出了很详细的流程说明(https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html)。

source("https://bioconductor.org/biocLite.R")biocLite("ChAMP")

NOTE: ChAMP有很多依赖包,安装时,若报错有哪个包没有,就继续安装 biocLite("YourErrorPackage"),可能3-4次就可以安装成功。

导入ChAMP包并测试:

导入ChAMP包后,根据是450K的数据或者是850K的数据,导入测试数据集,走一下分析流程,检测包是否正常工作,更重要的是看该包的文档,理解每一步流程的意义。该包的文档很详细,建议大家看原文档,下面给出的啰啰嗦嗦的介绍基本上都来自官网的文档说明(https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html)。

library("ChAMP")#450K的数据导入:testDir=system.file("extdata",package="ChAMPdata")myLoad <- champ.load(testDir,arraytype="450K")#850K的数据data(EPICSimData)

Step 3: 数据读入

数据的读入这里可能坑最多,首先450K 和850K甲基化芯片的原始数据格式都是IDAT, 因为数组是用两种不同的颜色来测量的,所以每个样本都有两个文件,通常是扩展名Grn.idat和Red.idat。数据在载入时还需要一个SampleSheet.csv文件(图5)(也称做pd file), 这个文件很重要,它包含了样本的信息,可以对照测试数据的csv文件和自己的csv文件,对信息不全的地方进行补充。尤其要注意SampleGroup 这一列信息是否有,这一列信息代表你想比较的表型类型,比如癌和癌旁。另一个我遇到过的一个隐形坑在Sentrix_ID,这一列数因为数字串很长,在Excel中可能以科学计数法显示,然后本来是长数字串后两位不一样的数字串都变为一样的,在读入时就会报重复字符的错误,所以这里一定要核查下长数字串的信息,如果有错误,自己重新输入时以文档格式输入,或者前面加右单引‘。csv文件准备好后,将csv文件与所有样本的芯片数据(即IDAT文件)放在一个文件下,然后就可以正常读入了。

图5 Sample_Sheet.csv fiel

library("ChAMP")myLoad <- champ.load("F:/850K Methylation Chip/biotree_850K/methy_rawData",arraytype = "EPIC")save(myLoad,file="myLoad.rda")

champ.load()包含了 champ.import() 和champ.filter(),这里会自动过滤p值0.01; probes beadcount <3 in at least 5% of samples;NoCG;probes with SNPs; MultiHit; probes located on X,Y chromosome。

在读入数据之后,最好保存,后续重复读入时会加快速度。

Step 4: 质控和标准化

CpG overview:

质控前可以先看看CpG的分布,包括在染色体上的分布;CpG岛附近的 open sea, shelf,shore (参考图2,理解具体意思) ; UTR,TSS; I 型探针和II探针上的分布(图6),这个信息对后续DMP的分析有帮助。

CpG.GUI(arraytype="EPIC")

质控:

然后进行质控,有两种方式:champ.QC() 和 QC.GUI()。champ.QC会产生三种类型的图(点图,beta 分布图,聚类图)以pdf格式输出,QC.GUI产生5个图,多了一个I型、II型探针图和热图(图7)。所有的GUI功能都比较耗内存,且产生的是网页交互式的图片,每幅图的右上角给的都有保存按钮,要注意的是保存时文件名要加上.png的后缀(图7)。

#champ.QC()QC.GUI(arraytype="EPIC")

图7 QC Overview

标准化:

champ.norm 提供了四种方法:BMIQ, SWAN1, PBC2 and FunctionalNormliazation4。默认的方法是BMIQ, 且BMIQ对850K的标准化方法更好一点,所以这里我选择的是BMIQ的标准化方法,没有尝试其他的标准化方法。

myNorm <- champ.norm(arraytype="EPIC")QC.GUI(myNorm,arraytype="EPIC")save(myNorm,file="myNorm.rda")

SVD plot 和批次效应:

SVD(singular value decomposition) 这里用于评估数据集中变量的主要成分。这种成分可能确实是你感兴趣的生物因素,也可能是技术来源的一些变量成分(称为批次效应)(图8)。如果存在批次效应,就进行批次效应的矫正,矫正完之后可以再看看SVD plot。

champ.SVD()

图8 SVD Plot

Step 5: 差异甲基化分析(DMP & DMR & DMB)

差异分析是多数研究都要分析的,这里包括三种方法:DMP,DMR,DMB。DMP代表找出Differential Methylation Probe(差异化CpG位点),DMR代表找出Differential Methylation Region(差异化CpG区域),Block代表Differential Methylation Block(更大范围的差异化region区域)

简单来说,DMP是找出一个一个的差异甲基化CpG位点,DMR就是一个连续不断都比较长的差异片段,科学家们觉得,这样的连续差异片段,对于基因的影响会更加明显,只找这样的片段,可以使得计算生物学的打击精度更为准确,也可以让最终找出来的结论数据更少,便于实验人员筛选。另外一个类似的东西就是DMB,那个东西出现的原因是,有的科学家觉得,DMR这样的区域还不够显著,DNA上的甲基化出现变化,可能是绵延几千位点的!而且只会在基因以外的区域,但是这些基因以外的区域发生变化,却会导致基因的表达发生变化。你可以想象成,北京周边的河北在大炼钢铁,然后北京也跟着雾霾了,大概就是这意思。

DMP,DMR,DMB的结果都是基于的shiny的交互页面,左栏上方是 P-value 和 abs(logFC) ,可以选择想看的值,然后点submit, 右栏可以生成差异甲基化表,热图,feature&cgi, 左栏下方还有基因,CpG按钮,选择你想看的结果,submit, 右栏就会生成相应gene,CpG结果(图9)。

myDMP <- champ.DMP(arraytype="EPIC")save(myDMP,file="myDMP.rda")DMP.GUI()myDMR <- champ.DMR(arraytype = "EPIC",method="DMRcate",cores=1)save(myDMR,file="myDMR.rda")DMR.GUI(arraytype="EPIC")#myBlock <- champ.Block(arraytype = "EPIC")#Block.GUI(arraytype="EPIC",compare.group=c("PrEC_cells","LNCaP_cells"))

图9 DMP Overview

Step 6: 基因富集和网络分析(GSEA & EpiMod)

差异甲基化分析后,你可能想知道DMP,DMR中涉及到的基因是否可以富集到某个生物功能或通路,GSEA(Gene Set Enrichment Analysis)和EpiMod(Differential Methylated Interaction Hotspots)提供了可以寻找作用通路网络中的疾病关联小网络的功能 (图 10)。

myGSEA <- champ.GSEA(arraytype = "EPIC")save(myGSEA,file="myGSEA.rda")myEpiMod <- champ.EpiMod(arraytype="EPIC")save(myEpiMod,file="myEpiMod.rda")

Step 7: 拷贝数变异分析(CNA)

拷贝数变异,也就是有些基因片段被复制的此处过多或者过少,从而导致某些疾病。但是这个函数作者觉得有点粗糙,精度还不够。我试着跑了一下,时间超长(图11)。

myCNA <- champ.CNA(control = F,arraytype = "EPIC")save(myCNA,file=myCNA)

图11 Frequency Plot of Cancer Sample

小结:如果用ChAMP包对450K或850K甲基化数据进行分析时,一是最好有个配置高一点的电脑;二是初始数据导入时,注意csv文件的格式,且要和IDAT文件放在一个文件下;其余的流程很少会遇到bug, 但最关键的是理解每一步的意义,能够根据分析的结果挖掘出想要的东西。

ps: 这次作业提供的公共数据,有IDAT文件,也有个csv文件,但是这里的csv文件和我的csv文件差别很大,不是很明白这里的csv文件是什么,有什么作用。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-11-08

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏非典型技术宅

iOS传感器:实现一个随屏幕旋转的图片1. 加速计介绍2. 加速计的使用3. 获取加速计数据的两种方式4. 实现图片永远水平方向

1574
来自专栏Data Analysis & Viz

百年百图の中国(1900-1999):另类python爬虫和PIL拼图

标题有点长,也有点怪。前半部分文艺向,后半部分python技术向。目的就是用PIL库得到100张图的拼图(成果图见文末)。

742
来自专栏程序生活

Python爬虫系列(七)豆瓣图书排行榜(数据存入到数据库)

豆瓣用户每天都在对“读过”的书进行“很差”到“力荐”的评价,豆瓣根据每本书读过的人数 以及该书所得的评价等综合数据,通过算法分析产生了豆瓣图书250。 网址:豆...

3754
来自专栏无原型不设计

优秀原型设计欣赏:美食类App原型制作-Kitchen Stories

题材有Mockplus(摹客)团队提供,仅供参考学习。

2047
来自专栏机器人网

最常见的仪器仪表原理动图,工程师一看就能懂!

一、温度仪表原理 1.薄膜热电偶的结构 ? 2.固体膨胀式温度计 ? 3.热电偶补偿导线的外形图 ? 4.热电偶温度计 ? 5.热电阻的结构 ? 二、压力仪表原...

35711
来自专栏机器人网

为什么采用4~20mA的电流来传输模拟量?

大家可能会非常熟悉RS232,RS485,CAN等工业上常用的总线,他们都是传输数字信号的方式。那么,我们用什么方式来传输模拟信号呢?工业上普遍需要测量各类非电...

2628
来自专栏贾老师の博客

【笔记】读 JeffDean 分布式系统

1353
来自专栏落影的专栏

iOS音视频播放(Audio Unit播放音频+OpenGL ES绘制视频)

前言 相关文章: 使用VideoToolbox硬编码H.264 使用VideoToolbox硬解码H.264 使用AudioToolbox编码AAC 使...

5009
来自专栏zone7

用Python告诉你深圳房租有多高

最近各大一二线城市的房租都有上涨,究竟整体上涨到什么程度呢?我们也不得而知,于是乎 zone 为了一探究竟,便用 Python 爬取了房某下的深圳的租房数据,以...

1130
来自专栏龙行天下CSIEM

科学瞎想系列之六十四 双馈电机绕组故障诊断

双馈是大型风力发电的主流技术之一,目前已装机运行的并网型风力发电机组大多采用这一技术路线。通常双馈发电机绕组出现故障后很难在塔上维修,必须下塔。下塔!宝宝...

2855

扫码关注云+社区