Illumina HumanMethylation450 BeadChip (甲基化450k芯片) 预处理初探

旭日早升

声名：本文旨在提供一个初步的450K芯片处理需要考虑的问题，它并不完美，希望大家可以更新完善，初学者可以借鉴。

一、Illumina HumanMethylation450 BeadChip （甲基化450k芯片）简介
Illumina最早的甲基化芯片是27K（K代表1000，表示大概可以测到的CpG位点数）的数据，后来增加到了450K（主流的甲基化芯片），而目前illumina已经出了新一代产品EPIC（我习惯称之为850K），但是技术核心在450K已经成熟了，所以分析流程这里以450K为主（也是目前数据库主流的甲基化芯片数据）。
450K芯片的芯片图形及其原理可以以下图（来源于illumina.com）展示

1、芯片：一张芯片包括12个array（如图显示），也就是一张芯片可以做12个sample，一台机子一次可以跑8张芯片，也就是一共96个sample，每个样本可以测到超过450，000个CpG位点的甲基化信息（大概人所有的1%，但是覆盖了多数CpG岛和启动子区），芯片本身包含一些控制探针可以做质控。
2、原理：简而言之，基于亚硫酸盐处理后的DNA序列杂交的信号探测。亚硫酸盐是甲基化探测的“金标准”，不管是芯片或者甲基化测序，都要先对DNA样品进行亚硫酸盐处理，使非甲基化的C变成U，而甲基化的C保持不变，从而在后续的测序或者杂交后区分出来。450K采用了两种探针对甲基化进行测定，Infinium I采用了两种bead（甲基化M和非甲基化U，如图显示），而II只有一种bead（即甲基化和非甲基化在一起），这也导致了它们在后续荧光探测的不同，450K采用了两种荧光探测信号（红光和绿光）。
更多信息可以在官网查看 http://www.illumina.com/products ... _beadchip_kits.html


二、分析需要考虑的问题
1、背景校正
2、红光和绿光的校正
3、控制芯片的使用（illumina450K本身有一些控制芯片，可以用来做质控，如亚硫酸盐处理效率）
4、探针类型（I型和II型）的校正（不同探针类型产生的数据不同）
这个问题我们之前关注很多，这里附上两篇文献供大家参考，最终我们选择BMIQ的方法（基于ebayes的原理将II型探针的甲基化水平拉伸到I型水平，如下图显示）来做矫正。
文献1： http://bioinformatics.oxfordjournals.org/content/29/2/189.short
文献2： http://www.tandfonline.com/doi/abs/10.4161/epi.24008
图片来源于第1篇文献

5、位置的校正（芯片上的不同位置产生的数据可能会有偏差）
6、批次的校正（不同的批次做的数据会有偏差）
7、探针序列本身是否可靠（有些探针本身位于repeat区或者包含snp等就会影响杂交及最后的结果，应该去除，附上一片参考文献，里边有list可以用来去除不好的探针）
文献：http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-51


三、数据处理
1、GenomeStudio
GS是illumina开发的软件，基于图形界面的操作处理，适合于没有R及编程基础的人使用，但是使用GS需要权限，本人试过挺难破解的，如果有人可以破解，可以提供给论坛小伙伴使用。
2、基于R和bioconductor的pipeline
bioconductor里开发了很多package供大家使用，如果你会R，那么处理这个甲基化芯片的数据将变得简单。
可以处理450K芯片的package有lumi、minfi、wateRmelon、ChAMP等，没有哪一种就特别好，大家都在不断改进，所以只要你知道大概的流程和需要注意的问题，那么你也可以自己写代码处理，只是package可以帮你省很多事情。下边我会附上我的处理流程图和详解，还有我的代码，我的代码是基于minfi的，我再次强调，代码只是手段，你可以用其他的package（例如我的同事很多使用ChAMP），也许会更好。
流程图：

流程图详解：我们一般下载或者iscan后的原始数据格式为Idat，首先可以得到每个sample的每个probe的p值和bead数，根据p值和bead数可以进行样本和探针的过滤，过滤之后需要用BMIQ的方法进行I和II型探针的校正，矫正之后去掉那些包含snp之类的不好的探针，最后对数据做batch的校正。校正之后的数据就是预处理后的数据了，可以用于更下游的分析，如差异甲基化和甲基化与表达的关联分析等。
我的代码： https://github.com/wkl1990/Illumina450K-methylation，再次强调，我的代码并不完美，如果有bug，请告诉我修改，构建你自己的代码，让自己敲起来。

最后，感谢坛主大神的邀请和鼓励，他的学习和分享精神一直是我的榜样，希望在这里我可以书写和学习更多的知识，结识更多的朋友。

作者：生信技能树-旭日早升
转载请说明来源：生信技能树-biotrainee.com

Jimmy

Chips were scanned using an Illumina HiScan on a two-color channel to detect Cy3 labeled probes on the green channel and Cy5 labeled probe on the red channel.


Illumina GenomeStudio Software 2011.1 was used to read the array output and conduct background normalization.


illumina的软件自动就把IDAT文件变成了每个样本3列信息的txt文件啦。当然啦，minfi包的一个函数read.450k.exp也可以直接读取 .IDAT files
The signalA, signalB and probe intensity were exported for autosomal probes and read into R.


M values were generated using the Bioconductor (http://www.bioconductor.org; Fred Hutchinson Cancer Research Center, Seattle, WA, USA) methylumi package, M = log2(intensity m + 1/intensity u + 1) since this value has been shown to be valid for statistical analyses [42].
Following correction for chip to chip color bias using the Bioconductor lumi package [43] and probe type correction using subset-quantile within array normalization (SWAN) [44], M values were converted to ß values using the equation ß = (2M/(2M + 1)). 它这里给的公式太复杂了，这里有一个简单的，https://www.bioconductor.org/hel ... /minfi_BioC2014.pdf就是(average β = signal B/(signal A + signal B + 100))而已

The ß value is a number ranging from 0 to 1 that is directly proportional to percentage DNAm; thus to ease interpretation, we have reported results as ß values.


The microarray data used in this article was submitted to the NCBI GEO under accession number [GSE:42409].


Probes with a detection p value >0.01 in any sample, probes with no ß value in any sample, all rs and ch probes, all sex chromosome and non-specific probes were removed prior to analyses.


The level of DNAm for 428,216 probes in our sample dataset was intersected with the expanded annotation for further analyses.

Jimmy

给你两篇文献：
2015
Analysis pipelines and packages for Infinium HumanMethylation450 BeadChip (450k) data
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4304832/

2013
A data-driven approach to preprocessing Illumina 450K methylation array data
https://bmcgenomics.biomedcentra ... 86/1471-2164-14-293
还推荐看看这个；https://github.com/elevenfly/Min ... kflow_1_17_2015.rmd

kevinerin

感谢楼主,450k的应用很广泛！！收益挺多，之前自己的分析对校正并没有太在意。
我使用的是bioconductor的RnBeads, 好处是使用说明非常的详细，还有R的脚本，地址http://rnbeads.mpi-inf.mpg.de/。除了找到Different methylation regions外，RnBeads可以查找任意区间的甲基化水平，不过注释的区间长度有限。
问楼主个问题，如果从GEO上下载数据来处理的话，经常没有idat文件而只是个处理后的txt文件，这时候用哪个package处理，有什么需要注意的？

旭日早升

kevinerin 发表于 2016-9-25 10:42
感谢楼主,450k的应用很广泛！！收益挺多，之前自己的分析对校正并没有太在意。
我使用的是bioconductor的Rn ...

RnBeads之前也经常看到，应该很不错。确实GEO很多下的文件都不是idat，这时候我一般就注意有没有给p值和nbead值，如果有就可以做过滤，一般会给出p值，实在没有就看看他给出的txt是什么值，例如是beta值还是信号值，确定之后一般就不会用固定的package的，就根据需要来做。

kevinerin

谢谢楼主的回复。自己还只是比葫芦画瓢，要向您请教。在芯片数据预处理和质控方面，方法是有一些，不过我不知道评判标准是什么，如果质控看着不怎么好，如何舍弃些芯片的数据。

旭日早升

kevinerin 发表于 2016-9-27 02:02
谢谢楼主的回复。自己还只是比葫芦画瓢，要向您请教。在芯片数据预处理和质控方面，方法是有一些，不过我不 ...

大家互相交流学习，我也不敢说自己有多了解。illumina会提供一些control的探针，例如亚硫酸盐的探针，可以用来评估，之前看到的标准就是用测到的control探针的上下四分位数作为上下边界。另外根据p值过滤，例如超过1%的探针都是p>0.01，这样的sample可能就要去掉。还可以再看下PCA和MDS或者聚类之类的图，如果某些样本和其他样本差距很大，也最好去掉。

BUNO19

楼主您好，看了您的分析，收益颇多！
还想向您请教，如果只是β值的数据、探针的信息，应该如何分析？

旭日早升

BUNO19 发表于 2017-7-6 13:09
楼主您好，看了您的分析，收益颇多！
还想向您请教，如果只是β值的数据、探针的信息，应该如何分析？ ...

如果有beta值和探针信息，就可以做下游分析了，例如差异甲基化位点分析等等

bioant_2017

总结很丰富，mark了