甲基化分析系列之CHAMP包的使用（一）
大家好，我是风间琉璃，这里是每周四的甲基化生信分析专栏。前面有小伙伴在评论使用CHAMP包对甲基化芯片进行分析，没错，那今天它来了！CHAMP的功能作为集成式一站式分析甲基化的包，其作用非常强大，并具有minfi包的一部分功能。在进行分析的同时以图片的方式提供分析的结果。首先我们先来看一张流程图，详细的提供了我们使用CHAMP包分析甲基化芯片的详细流程。
流程图中的不同的颜色代表甲基化芯片分析的不同阶段。
蓝色部分代表着甲基化数据的准备阶段，比如加载数据（loading）、标准化（normalization）、质控（Quality control）等等。
红色部分是分析的总流程，包括找出甲基化差异位点（Differentially Methylated Positions, DMPs）、找出差异甲基化区域（Differentially Methylated Regions, DMRs）、拷贝数变异分析（champ.CNA）、甲基化通路富集分析（champ.GSEA）等等。
黄色部分则是可以进一步提供图形用户界面（Graphical User Interface, GUI）功能，这正是CHAMP包分析的强大之处，我们之后还会进一步介绍。
另外灰色的实线代表着我们流程中必须进行的步骤，而灰色的虚线则不是必要步骤。而灰色实线边缘具有绿色的阴影部分是CHAMP开发者推荐的主要分析流程。
来，根据这张流程图我们再复习一下甲基化芯片的处理流程：
1.数据导入（champ.import）,包括甲基化位点筛选（champ.filter）和缺失值的插补（champ.impute）
2.数据质控（champ.QC）
3.数据标准化（champ.norm）
4.奇异值分析了解批次效应（champ.SVD）
5.批次校正（champ.runCombat）和细胞类型异质性校正（champ.refbase）
6.甲基化差异位点分析（champ.DMP）和差异区域分析（champ.DMR）
7.甲基化位点/区域的通路富集分析（champ.GSEA）
安装包
那我们复习一遍之后，我们就直接开始吧！！首先我们先安装包和加载包。
#####安装包，并加载#if (!requireNamespace("BiocManager", quietly = TRUE))#   install.packages("BiocManager")# BiocManager::install("ChAMP")library("ChAMP")
CHAMP作为一站式集成式甲基化芯片分析包，下载过程中小伙伴们需要一点耐心。因为每次运行这个命令可能都会提示“xx package”没有安装的情况。这时候我们只需要根据它的指令，一步步安装它的提示缺少的包。过程可能会比较长。安装好之后，我们加载CHAMP包的时候，如果出现以下提示，证明我们加载成功。不然还是需要把提示缺乏的包一个一个安装上。
第一步：数据导入和筛选
#####首先找到我们测试数据的位置testDir=system.file("extdata",package="ChAMPdata")#展示位置testDir## [1] "D:/user/Documents/R/win-library/4.0/ChAMPdata/extdata"####champ.load加载数据，8个样本：4个肺癌，4个正常myLoad <- champ.load(testDir,arraytype="450K")## Failed CpG Fraction.## C1 0.0013429122## C2 0.0022162171## C3 0.0003563249## C4 0.0002842360## T1 0.0003831007## T2 0.0011946152## T3 0.0014953286## T4 0.0015447610#####full pipeline全部流程 一次性跑完，但会报错。#champ.process(directory = testDir)
当然上面还提供了一站式全流程分析的命令-champ.process，但是很容易报错，所以并不建议大家使用。
这一步需要花一会儿的时间，因为这里CHAMP包一步步给我们提示我们它的每一步做了什么。因为给出的提示很多，我截取了最重要的部分。
这一步骤主要包括2个模块。
模块一
section 1:读入表型信息文件（pdata）。
section 2：读入IDAT文件（包括计算P值、M值矩阵等等）。
section 3：注释探针成为甲基化位点。
模块二
section 1:查看输入的数据是否正确 。
section 2: 进行筛选，（包括排除具有高比例失败的探针的样本、筛查P值、筛查SNP等等） 。
这些我们在minfi包需要自己计算的过程，CHAMP已经自动给我们做了，不愧是一站式集成式的分析神器! 当然，champ.load如果只能根据CHAMP内置的筛选标准来进行导入数据，那就太不智能，所以我们可以在命令中设置自己的筛选标准，比如：
#myLoad <- champ.load(testDir,arraytype="450K",filterDetP=TRUE,# autoimpute=TRUE,filterBeads=TRUE,filterSNPs=TRUE,# filterMultiHit=TRUE,filterXY=TRUE)####查看数据集的pdata，分组信息等myLoad$pd## Sample_Name Sample_Plate Sample_Group Pool_ID Project Sample_Well Slide## 1 C1 NA C NA NA E09 7990895118## 2 C2 NA C NA NA G09 7990895118## 3 C3 NA C NA NA E02 9247377086## 4 C4 NA C NA NA F02 9247377086## 5 T1 NA T NA NA B09 7766130112## 6 T2 NA T NA NA C09 7766130112## 7 T3 NA T NA NA E08 7990895118## 8 T4 NA T NA NA C09 7990895118## Array## 1 R03C02## 2 R05C02## 3 R01C01## 4 R02C01## 5 R06C01## 6 R01C02## 7 R01C01## 8 R01C02还有另外一种方式进行。##我们换另外一种方式进行导入数据# myImport <- champ.import(testDir)# myLoad <- champ.filter()
这里向大家提出一个问题，如果我们只有β值和M值得矩阵，能否进行相关得分析呢？（不知道β和M值得同学请查阅上一期：《手把手教你甲基化生信分析—甲基化minfi包的使用（一）》）。
为什么会有这样得问题呢？因为我们在GEO数据库进行数据下载得时候，有些GSE提供得并不是IDAT格式得文件，而是β值或者M值得矩阵。这时候就需要想一想别的方法进行数据导入了。CHAMP提供了一个很好得方法。
###如果只有beta值矩阵,我们需要提供β矩阵以及表型数据#我们先看一看提供数据的格式是怎么样的class(myLoad$beta)## [1] "matrix" "array"head(myLoad$beta)## C1 C2 C3 C4 T1 T2## cg00000957 0.79274273 0.81628392 0.87269772 0.84761398 0.87328478 0.82818422## cg00001349 0.67334505 0.62750070 0.67697354 0.70975670 0.46656322 0.76167976## cg00001583 0.07940828 0.07251949 0.37411396 0.15932854 0.44702168 0.14756981## cg00002028 0.02977963 0.02320829 0.04002213 0.02730524 0.03504087 0.09662921## cg00002719 0.03607880 0.04530982 0.07608838 0.02276532 0.23375170 0.34626219## cg00002837 0.21105214 0.27023541 0.34336824 0.43784615 0.13593510 0.41579813## T3 T4## cg00000957 0.69923664 0.75353389## cg00001349 0.24434991 0.45570105## cg00001583 0.39635270 0.39144462## cg00002028 0.02083333 0.01201803## cg00002719 0.07642563 0.30255255## cg00002837 0.18168243 0.24598131##matrix格式的矩阵，行是cg位点，列是样本class(myLoad$pd)## [1] "data.frame"head(myLoad$pd)## Sample_Name Sample_Plate Sample_Group Pool_ID Project Sample_Well Slide## 1 C1 NA C NA NA E09 7990895118## 2 C2 NA C NA NA G09 7990895118## 3 C3 NA C NA NA E02 9247377086## 4 C4 NA C NA NA F02 9247377086## 5 T1 NA T NA NA B09 7766130112## 6 T2 NA T NA NA C09 7766130112## Array## 1 R03C02## 2 R05C02## 3 R01C01## 4 R02C01## 5 R06C01## 6 R01C02#数据框格式，样本名、样本分组、Slide等信息##使用champ.filter,同样包含了过滤等步骤，可以自定义#myLoad=champ.filter(beta = myLoad$beta,pd = myLoad$pd)###如果是M值呢？#myLoad=champ.filter(M=myLoad$M,pd = myLoad$pd)
那么不论是β/M值矩阵又或者IDAT原始文件，我们都可以通过CHAMP包进行一站式读取，是不是很方便呢？
好啦，总结一下，我们这一期详细讲解了CHAMP的分析流程，以及如何进行数据导入。我们下期将详细讲解整个CHAMP的分析流程。
我是风间琉璃，我们下期见。
回复“风间琉璃甲基化03”，得到今天代码。
欢迎大家关注解螺旋生信频道-挑圈联靠公号~

继续阅读
阅读原文