集识别差异甲基化和拷贝数变异于一体的R包--ChAMP

摘要：ChAMP是一个功能异常强大的R包，包括了从甲基化芯片原始数据预处理、标准化到差异的识别等全面的功能。本文以450K数据为实例进行全方面的演示，非常详细的展示了每一步数据处理的过程及结果。

简介：ChAMP是一个针对HumanMethylation450 BeadChip和EPIC甲基化数据的R包，可以实现数据预处理、质量控制、数据标准化、校正批次效应、识别差异甲基化位点/差异甲基化区域、识别拷贝数变异区域等功能。值得注意的是使用相同的450K数据做差异甲基化分析和拷贝数变异分析可以在很大程度上降低不同平台数据对分析结果的影响。ChAMP包适用于各个平台，需要的R版本为R-3.2及以上，一般200个样本运行程序时需要8G内存。本文以450K数据为例介绍ChAMP的基本用法。

安装

>source('http://bioconductor.org/biocLite.R')

>biocLite(c('minfi','DNAcopy','impute','marray',

'limma','preprocessCore','RPMM','sva',

'IlluminaHumanMethylation450kmanifest',

'wateRmelon','isva','quadprog',

'bumphunter','doParallel','qvalue',

'RefFreeEWAS','GenomicRanges',

'plyr','ChAMP'))

>Library(ChAMP)

使用

1、数据准备

（1）level1阶段的450数据，格式为.idat

（2）样本说明文件，格式为.csv。文件中必须包括样本名、样本类别（用于识别差异甲基化和拷贝数变异）、Sentrix_ID和Sentrix_Position，其余为可选。文件格式如下图所示，[Header]为文件说明（可以删除），[Data]为数据部分。

*注：样本文件和说明文件必须在同一个目录下。

2、载入数据，数据预处理

>myLoad=champ.load(directory='/input/dir',

resultsDir='/output/dir/',methValue = 'B',

filterXY = TRUE, QCimages = TRUE,

filterDetP = TRUE,detPcut = 0.01,

filterBeads=TRUE,beadCutoff=0.05,

filterNoCG=FALSE,filterSNPs=TRUE,

filterMultiHit=TRUE,arraytype='450K'))

#参数含义：

# methValue：选择输出beta值或M值

# filterXY：是否过滤X和Y染色体

# QCimages：是否保留质控图片

# filterDetP：是否对detection P value进行筛选，dePcut指定P value阈值

# filterBeads：根据beadCutoff确定是否去除在至少n%的样本中beadcount

# filterNoCG：是否过滤非CG位点

# filterSNPs：是否过滤SNP

# filterMultiHit:是否过滤匹配到多个位置上的探针

# arraytype：选择输入文件类型，可选项为450K或EPIC

常用筛选条件：

（1）去除在一个或多个样本中P>0.01的探针

（2）去除在至少5%的样本中，beadcount

（3）去除包含SNP的探针

（4）去除映射到多个位置上的探针

（5）去除X、Y染色体上的探针

*注：当样本量较大时，载入数据较慢，在执行完champ.load函数后可以将其保存，方便下次直接使用。

>save(myLoad,file='/output/dir/currentStudyloadedData.RData')

>load('currentStudyloadedData.RData')

#输出结果：

（1）failedSample.txt文件（如下图）中给出每个样本中不符合条件的探针比例，如果比例超过5%，则考虑删除该样本，重新执行champ.load()函数。

（2）raw_densityPlot.pdf文件展示不同组别样本未标准化的beta值的密度分布。

（3）raw_mdsPlot.pdf展示利用前1000个样本间差异最大的位点绘制的multidimensionalscaling图，可以反映样本的相似性。

（4）raw_SampleCluster.jpg样本的层次聚类图（当样本数多于65个时不展示聚类图）

3、数据标准化，校正type-2 bias

ChAMP包中数据标准化的可选方法包括BMIQ(Teschendorff,2013), SWAN (Maksimovic, 2012), PBC (Dedeurwaerder, 2011)或NONE（不进行校正），默认使用BMIQ

>myNorm=champ.norm(resultsDir = '/output/dir')

#输出结果：

（1）标准化的beta值密度分布图

（2）标准化后的multidimensionalscaling图

（3）标准化后的样本层次聚类图

4、用奇异值分解（SVD）方法识别变异组分（包括生物学因素或者技术变异）

SVD不对数据做任何处理，只是给出各变异组分的显著性，输出结果是一个热图，颜色越深表示p值越小，即技术因素造成的变异较大，如果经SVD分析发现技术因素造成的变异较大，则应使用Combat等方法校正批次效应。

>champ.SVD(resultsDir = '/output/dir')

#输出结果：

SVDsummary.pdf

5、校正批次效应

当SVD分析显示技术因素造成的变异较大时，需要根据样本信息中SentrixID进行校正,校正之后再使用champ.SVD()函数检测变异组分，确定校正效果。

>batchNorm=champ.runCombat()

>champ.SVD(resultsDir = '/output/dir',beta= batchNorm$beta)

6、识别甲基化可变位点（Methylation Variable Positions，MVPs）

根据样本信息中的Sample_Group筛选组间差异甲基化位点，如果文件中有多个组，则默认筛选前两个组间的差异甲基化位点。

>limma=champ.MVP(resultsDir = '/output/dir',beta.norm = myNorm$beta)

*注：如果没有用champ.runCombat()校正，则beta.norm = batchNorm$beta，或者不指定该参数，程序会自动选择最后一次出现的标准化的甲基化数据。

#输出结果：

（1）MVP_ALL_CvsT_BHadjust.txt，全部位点的结果，包括原始P值、校正后的P值、map到的gene、位点在基因上的位置、位点相对于island的位置等信息（如下图）。

（2）MVP_0.05_CvsT_BHadjust_3882.bed，校正P值后组件显著差异的位点结果，包括位点名、染色体和染色体上的位置（如下图）。

7、识别差异甲基化区域（DMRs）

方法1：Bumphunter，不依赖于之前的输出结果

>bump=champ.DMR(method='Bumphunter',arraytype='450K')

#输出结果：

（1）Bumphunter方法识别出的差异甲基化位点

（2）Bumphunter方法识别出的差异甲基化区域

方法2：ProbeLasso，依赖于上一步输出的差异甲基化位点结果limma

>lasso=champ.DMR(resultsFile=limma,method='ProbeLasso',arraytype='450K')

#输出结果：

（1）ProbeLasso方法识别到的差异甲基化探针

（2）ProbeLasso方法识别到的差异甲基化区域

8、基于环状二元分割方法识别拷贝数变异区域

>CNA=champ.CNA(resultsDir = '/output/dir',controlGroup='C')

*注：参数中controlGroup的值应该与样本文件Sample_Group一致，否则会报错，默认情况下是“Control”。

#输出结果：

（1）实验组（T）中各染色体上拷贝数扩增和拷贝数缺失的样本比例

（2）实验组每个样本（T1/T2/T3/T4）中拷贝数扩增、拷贝数缺失在各染色体上的分布情况

（3）实验组（T）各样本的拷贝数扩增、缺失详细结果，下图中每一列分别表示区域ID、染色体、区域起始位置、区域终止位置、区域中包含的CG位点个数及该区域上拷贝数变异的度量值seg.mean，一般认为seg.mean≥0.3为拷贝数扩增，seg.mean≤-0.3为拷贝数缺失（有的研究中也以0.2为阈值）。

*注：以上流程可以用一个函数实现，但是由于程序运行过程可能出现一些错误，另外不利于具体参数的设置，所以建议分部执行。

>champ.process(directory ='/input/dir')

9、获得基因的拷贝数变异谱

由于ChAMP包得到的每个实验组样本单独的拷贝数区域结果，为了方便后续分析，可以使用CNTools和cghMCR将ChAMP识别到的拷贝数变异区域合并，比对到参考基因，以获得所有基因的拷贝数变异谱，执行此过程前，需要先把多个独立的拷贝数变异结果文件合并成一个文件，本文中合并后的文件名为segAll.txt。此外，该过程还需要一个参考基因文件，格式如下图：每一列分别表示染色体、起始位置、终止位置、Entrez ID和gene symbol（列的顺序不能更改）。

>library(CNTools)

>library(cghMCR)

>sampleData

>refseq

>set.seed(1234)

>sample.names

>number.of.samples

>convertedData

>rdseg

>cn.profile

>write.table(cn.profile,'can_profile.txt',row.names=F,col.names=T,quote=F)

#结果文件

每一行表示一个基因，每一列表示一个样本，数值为每个基因在各个样本上的seg.mean。

参考文献：

Morris TJ, Butcher L, Feber A, Teschendor A, Yuan Tian, Chakravarthy A, Beck S.The ChAMP Package.April 18, 2016.

Tiffany J. Morris1, Lee M.Butcher, Andrew Feber, Andrew E, et al.ChAMP: 450k Chip Analysis Methylation Pipeline.Dec 12, 2013.