学会利用Minfi包进行甲基化分析（二）
大家好，我是风间琉璃。上一周我们在《手把手教你甲基化生信分析——甲基化minfi包的使用（一）》介绍了甲基化芯片450k的工作原理，minfi包所包含的数据结构类型，我们怎么读取甲基化的原始数据IDAT文件，如何自定义自己的target文件。并且介绍了我们甲基化芯片的分析流程。
1.数据读取
2.数据质控
3.数据预处理
4.筛选探针
5.进一步分析
今天我们就继续介绍下面四步的分析流程。
第二步：数据质控
########quality control########制作QC报告qcReport(RGset, sampNames=colnames(RGset), sampGroups=pData(RGset)$rep, pdf= "qcReport.pdf")## png ## 2#####conversion probe plotcontrolStripPlot(RGset)
###QCMset=preprocessRaw(RGset)QC=getQC(Mset)####展示样本的QC图plotQC(QC)
#####展示密度图densityPlot(Mset, sampGroups =pd$rep)
densityBeanPlot(Mset, sampGroups = pd$rep)
### Type plot and mds plotplotBetasByType(Mset[,1], main=colnames(Mset)[1])
mdsPlot(Mset, numPositions = 1000, sampNames=colnames(Mset), sampGroups=pData(Mset)$rep)
 可以看到我们的样本整体质量还是很不错的。
第三步：数据的预处理
对于数据的预处理，我们分为两步，一步是筛选不合格的探针和不合格的样本。第二步是对样本进行标准化。一共有六种方法将将我们读取的RGChannelSet格式文件标准化并转化为MethylSet格式的文件。我们首先介绍第一步：筛选不合格的探针和样本。我们可以通过detectionP命令对每一个样本的每一个探针检测其P值。在这里我们默认P值大于0.01是不可信的。所以我们要剔除它。这一步其实包括两个内容，通过探针得到的P值筛选合格的探针，以及通过P值筛选合格的样本。
我们这里的标准是：
（1）某个样本呢所有的探针Pvalue平均值大于0.01，那么我们认为这个样本不合格；
（2）如果某一个探针，在任何一个样本出现Pvalue值大于0.01，那么这个探针不合格。
所以使用如下代码。
#######filtering###sample filter by detectionP#######得到Pvalue矩阵detP <- detectionP(RGset)##根据列计算，如果Pvalue均值小于0.01，则为Truekeep_sample = apply(detP,2,mean)<0.01##类似于数据框格式，可以行提取合格探针，列提取合格样本RGsetEx_cutSam <- RGset[,keep_sample]###find failed probe by detectionP##得到Pvalue矩阵detP <- detectionP(RGset)##根据行计算，如果所有的Pvalue均小于0.01，则为Truekeep_probe = apply(detP,1,function(x) {all (x< 0.01)})##类似于数据框格式，可以行提取合格探针，列提取合格样本RGset_cutprob = RGset[keep_probe,]
进行探针筛选之后，我们就需要进行数据预处理和标准化了，这一步，我们通过以下六个方法，把RGChannelSet格式文件转化为MethylSet格式的文件。
## [preprocessFunnorm] Background and dye bias correction with noob## [preprocessFunnorm] Mapping to genome## [preprocessFunnorm] Quantile extraction## [preprocessFunnorm] Normalization## [preprocessQuantile] Mapping to genome.## [preprocessQuantile] Fixing outliers.## [preprocessQuantile] Quantile normalizing.
进行数据的标准化之后，我们就需要对探针进行筛选也就是我们的第四步。
第四步：筛选探针
这时候我们使用之前的preprocessRaw()命令得到Mset.raw进行下面的计算，当然这里我们同样可以选择其他方法预处理后得到的MethylSet格式格式文件进行分析。在筛选探针的过程中，我们主要筛选排除具有SNP位点、存在在性染色体（X/Y）上的探针，当然还有其他的探针，我们之后再说。这里还需要强调一点，我们目前得到的是MethylSet格式的文件，但这个文件还没有映射到基因组中，但排除SNP以及X/Y染色体上的探针都是需要我们知道每一个探针在基因组上的位置信息才能够进行下面的步骤，所以我们需要再熟悉minfi包的包含的数据格式。
#####map to genes#####GMSet=mapToGenome(Mset.raw)##得到GenomicMethylSet格式文件###得到注释文件annotation <- getAnnotation(GMSet)#head(annotation) ###可以展示annotation文件#####remove the SNPs#####得到SNP的注释信息snps = getSnpInfo(GMSet)##展示snps文件head(snps)## DataFrame with 6 rows and 6 columns## Probe_rs Probe_maf CpG_rs CpG_maf SBE_rs SBE_maf## <character> <numeric> <character> <numeric> <character> <numeric>## cg13869341 NA NA NA NA NA NA## cg14008030 NA NA NA NA NA NA## cg12045430 NA NA NA NA NA NA## cg20826792 NA NA NA NA NA NA## cg00381604 NA NA NA NA NA NA## cg20253340 NA NA NA NA NA NA#####去除三种SNP位点探针GRSet = dropLociWithSnps(GMSet,snps = c("CpG","SBE","Probe"))
这里我们可以看到一共有三种SNP注释信息，分别是Probe_rs、CpG_rs、SBE_rs。CpG_rs 指的是这个CpG位点同时也是一个snp位点，SBE_rs指的是CpG位点下游的第一个碱基是snp位点，Probe_rs 指的是探针区域覆盖到的snp位点。这三类我们都需要去除。去除SNP之后，我们接下来需要提出性染色体上的探针。
####remove the sex chromosome#######得到性染色体上的探针sex_probe = rownames(annotation) [annotation$chr %in% c("chrX","chrY")]###返回逻辑值True/Falsekeep = !(featureNames(GRSet) %in% sex_probe)#####类似于数据框提取行的方式进行筛选探针GRSet = GRSet[keep,]
这样，探针筛选就结束了，接下来就到了大家都喜闻乐见的一步，差异分析找出差异甲基化位点。
第五步：差异分析
其实往往大家觉得最想做得一步是最简单的，最主要的还是在前面的数据预处理阶段。接下来我们看看差异分析需要用到的代码。
#######差异化位点#########得到β值beta = getBeta(GRSet)M=getM(GRSet)######得到分组信息group = pd$Sample_Group#####dmp differentially methylated positions的简称dmp = dmpFinder(beta,pheno = group,type = "categorical")# pheno参数用于指定样本的分组情况。#type参数用于指定分组类型：case/control不同的实验处理就是categorical；分组是连续型变量###就是continuous。###展示dmphead(dmp)## intercept f pval qval## cg01976913 0.65724404 6782.310 1.303074e-07 0.01721145## cg09496161 0.79161242 5208.394 2.208961e-07 0.01721145## cg00128718 0.80535201 4352.225 3.162741e-07 0.01721145## cg20663364 0.15947351 3923.227 3.891588e-07 0.01721145## cg16485975 0.18167801 3454.334 5.018627e-07 0.01721145## cg26729197 0.07409764 3176.724 5.933096e-07 0.01721145
我们发现一共四个指标，除了第一列的cg开头的甲基化位点信息，还包括intercept、f、pval、qval，其中的qvalue相当于adj.Pvalue，所以我们就可以找到所有adj.Pvalue<0.05的甲基化位点啦。另外这尔再加一个彩蛋，怎么计算差异甲基化区域，differentially methylated region(dmr)。
看代码
## [bumphunterEngine] Using a single core (backend: doSEQ, version: 1.5.1).## [bumphunterEngine] Computing coefficients.## [bumphunterEngine] Finding regions.## Warning in regionFinder(x = beta, chr = chr, pos = pos, cluster = cluster, : NAs## found and removed. ind changed.## [bumphunterEngine] Found 31133 bumps.
好啦，我们使用minfi包进行分析的450k甲基化芯片的流程就到这里啦，回复“风间琉璃甲基化02”得今天的代码。接下来我将给大家在介绍功能更加强大的甲基化处理分析包——CHAMP包。大家可以期待一下~我们下期见。
欢迎大家关注解螺旋生信频道-挑圈联靠公号~

继续阅读
阅读原文