学会利用Minfi包进行甲基化分析（一）
大家好，我是风间琉璃，这里是每周四的甲基化生信分析专栏。今天想跟大家介绍在进行甲基化分析中必不可少的R包——minfi包。可能大家没有使用过也应该听说过。那我们直接开始吧。
450K甲基化芯片简介
首先需要向大家介绍一点基础知识。
目前我们进行甲基化分析常规使用的是Illumina HumanMethylation450 BeadChip，也就是我们常说的450K芯片。送检的每个样本均在单独的阵列（包含红色和绿色两种通道）上进行测量，每一个阵列包含450,000个CpG位点。每一个位点具有两种不同的测量值：甲基化以及非甲基化测量值。
而这两种值是通过“Type I”或“Type II”中的一种方式进行测量。“Type I”用来只测量一种的颜色，而在这一个颜色通道种包括两个不同的探针来分别测量甲基化以及未甲基化值。而“Type II”只有一个探针，但包括双色通道来测量甲基化与未甲基化的值。我们需要注意的是在芯片种探针与CpG位点并非一一对应的，450K芯片一共有48万多个探针，而所包含的CpG位点差不多在45万个左右。这一点我们之后还会再提及。
实际上一个芯片（slide）包括12个阵列(array)，而每一个阵列能够分析一个样本，机器可以同时分析8张芯片，所有，一次性可以没有批次的分析96个样本。
M值和β值
接下来我们再介绍以下甲基化种的β值和M值。
首先在我们之前的分析中我们知道450K甲基化芯片能够对应一个CpG位点测出甲基化测量信号强度（M， methylated value）值和非甲基化信号强度（U，unmethylated value）。β值则等于M/(M + U + offset)，这里的offset是偏移量，以防止出现分母为0的情况。而M值则等于log2 (M/U)，也就是根据荧光信号进行log化。
在我们常规分析中，β值更加适合进行甲基化水平的定量，能够阐明生物学意义，任何等于或大于0.6的β值都被认为是完全甲基化的。
任何等于或小于0.2的β值被认为是完全未甲基化的。
β值在0.2和0.6之间被认为是部分甲基化的。
而M值更适合用于进行下游统计分析。介绍完理论知识之后，
Minfi包的使用
我们就介绍今天我们的主句Minfi包。首先下载和加载包
# if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")# BiocManager::install("minfi")library("minfi")library(minfiData)
其中的minfiData包存储着我们需要使用的数据。如果要了解minfi包，我们需要先了解这个包对数据的储存格式，先来看看例子。
#####数据类别RGsetEx## class: RGChannelSet ## dim: 622399 6 ## metadata(0):## assays(2): Green Red## rownames(622399): 10600313 10600322 ... 74810490 74810492## rowData names(0):## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02## 5723646053_R06C02## colData names(13): Sample_Name Sample_Well ... Basename filenames## Annotation## array: IlluminaHumanMethylation450k## annotation: ilmn12.hg19## RGsetEx: RGChannelSet, 622,399 featuresMsetEx <- preprocessRaw(RGsetEx)MsetEx## class: MethylSet ## dim: 485512 6 ## metadata(0):## assays(2): Meth Unmeth## rownames(485512): cg00050873 cg00212031 ... ch.22.47579720R## ch.22.48274842R## rowData names(0):## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02## 5723646053_R06C02## colData names(13): Sample_Name Sample_Well ... Basename filenames## Annotation## array: IlluminaHumanMethylation450k## annotation: ilmn12.hg19## Preprocessing## Method: Raw (no normalization or bg correction)## minfi version: 1.30.0## Manifest version: 0.4.0## MsetEx: MethylSet, 485,512 featuresGMsetEx <- mapToGenome(MsetEx)GMsetEx## class: GenomicMethylSet ## dim: 485512 6 ## metadata(0):## assays(2): Meth Unmeth## rownames(485512): cg13869341 cg14008030 ... cg08265308 cg14273923## rowData names(0):## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02## 5723646053_R06C02## colData names(13): Sample_Name Sample_Well ... Basename filenames## Annotation## array: IlluminaHumanMethylation450k## annotation: ilmn12.hg19## Preprocessing## Method: Raw (no normalization or bg correction)## minfi version: 1.30.0## Manifest version: 0.4.0## GMsetEx: GenomicMethylSet, 485,512 features####RsetEx<- ratioConvert(GMsetEx)RsetEx## class: GenomicRatioSet ## dim: 485512 6 ## metadata(0):## assays(2): Beta CN## rownames(485512): cg13869341 cg14008030 ... cg08265308 cg14273923## rowData names(0):## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02## 5723646053_R06C02## colData names(13): Sample_Name Sample_Well ... Basename filenames## Annotation## array: IlluminaHumanMethylation450k## annotation: ilmn12.hg19## Preprocessing## Method: Raw (no normalization or bg correction)## minfi version: 1.30.0## Manifest version: 0.4.0
我们这个一共有四个数据集，RGsetEx、MsetEx、GMsetEx、RsetEx，我们发现每一个的格式都不一样，第一个RGsetEx（RGChannelSet格式）、第二个MsetEx（MethylSet 格式）、第三个GMsetEx（GenomicMethylSet格式）、第四个（GenomicRatioSet格式）。我们发现minfi包具有有四种不同类型的数据格式。确实有点复杂，不过我们必须要记得，这对于我们之后的分析很有意义。
Minfi包的数据处理流程
首先我们先来看一张图。
我们可以看到，minfi包通过读取甲基化芯片的原始文件（IDAT文件）后，得到的是RGChannelSet格式，这很好理解，我们的探针读取本来就是通过红/蓝通道进行读取。然后通过preprocessRaw()命令得到了MethylSet 格式的文件，也就是我们示例中的MsetEx文件。到这儿一共出现了两条路，其实就是先进行映射到基因组得到GenomicMethylSet格式还是进行比值转换得到RatioSet格式的区别，最后得到GenomicRatioSet格式文件进行下游分析（既进行了映射又进行了比值转换）。
然后我们来介绍以下甲基化芯片的分析流程。
1.数据读取
2.数据质控
3.数据预处理
4.筛选探针
5.进一步分析
第一步：数据读取
#####reading data###baseDir <- system.file("extdata", package = "minfiData")baseDir#####展示我们读取的工作路径## [1] "D:/user/Documents/R/win-library/3.6/minfiData/extdata"list.files(baseDir)#####展示工作目录下的文件/文件夹## [1] "5723646052""5723646053""SampleSheet.csv"list.files(file.path(baseDir, "5723646052"))####展示baseDir下的5723646052文件夹中的文件## [1] "5723646052_R02C02_Grn.idat""5723646052_R02C02_Red.idat"## [3] "5723646052_R04C01_Grn.idat""5723646052_R04C01_Red.idat"## [5] "5723646052_R05C02_Grn.idat""5723646052_R05C02_Red.idat"
接下来我们需要构建targets文件，这对于我们读取数据很重要！！！！
######targets 文件targets <- read.metharray.sheet(baseDir)## [read.metharray.sheet] Found the following CSV files:## [1] "D:/user/Documents/R/win-library/3.6/minfiData/extdata/SampleSheet.csv"targets## Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID person age sex## 1 GroupA_3 H5 <NA> GroupA <NA> id3 83 M## 2 GroupA_2 D5 <NA> GroupA <NA> id2 58 F## 3 GroupB_3 C6 <NA> GroupB <NA> id3 83 M## 4 GroupB_1 F7 <NA> GroupB <NA> id1 75 F## 5 GroupA_1 G7 <NA> GroupA <NA> id1 75 F## 6 GroupB_2 H7 <NA> GroupB <NA> id2 58 F## status Array Slide## 1 normal R02C02 5723646052## 2 normal R04C01 5723646052## 3 cancer R05C02 5723646052## 4 cancer R04C02 5723646053## 5 normal R05C02 5723646053## 6 cancer R06C02 5723646053## Basename## 1 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646052/5723646052_R02C02## 2 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646052/5723646052_R04C01## 3 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646052/5723646052_R05C02## 4 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646053/5723646053_R04C02## 5 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646053/5723646053_R05C02## 6 D:/user/Documents/R/win-library/3.6/minfiData/extdata/5723646053/5723646053_R06C02####展示读取的方式# sub(baseDir, "", targets$Basename)
我们看一看我们电脑里面存储的Target.CSV长什么样。
其实我们可以很明显看到，Target文件的读取我们自己也可以做，不需要用它的read.metharray.sheet()命令，我们可以自己构建一个target文件，把对应的内容赋值粘贴进去就可以。当然自己构建的文件需要放进去分组信息、文件的路径（basename）、阵列（array）信息等。
####读取文件RGchannelSet文件格式，我们把它缩写成RGsetRGset <- read.metharray.exp(targets = targets)## Warning in readChar(con, nchars = n): 在non-UTF-8 MBCS语言环境里只能读取字节RGset## class: RGChannelSet ## dim: 622399 6 ## metadata(0):## assays(2): Green Red## rownames(622399): 10600313 10600322 ... 74810490 74810492## rowData names(0):## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02## 5723646053_R06C02## colData names(13): Sample_Name Sample_Well ... Basename filenames## Annotation## array: IlluminaHumanMethylation450k## annotation: ilmn12.hg19###展示分组信息pd <- pData(RGset)pd[,1:4]## DataFrame with 6 rows and 4 columns## Sample_Name Sample_Well Sample_Plate Sample_Group## <character> <character> <character> <character>## 5723646052_R02C02 GroupA_3 H5 NA GroupA## 5723646052_R04C01 GroupA_2 D5 NA GroupA## 5723646052_R05C02 GroupB_3 C6 NA GroupB## 5723646053_R04C02 GroupB_1 F7 NA GroupB## 5723646053_R05C02 GroupA_1 G7 NA GroupA## 5723646053_R06C02 GroupB_2 H7 NA GroupB
展示我们文件的注释信息
#####展示该甲基化文件的注释信息annotation(RGset)## array annotation ## "IlluminaHumanMethylation450k" "ilmn12.hg19"###quality controlmanifest <- getManifest(RGset)manifest####展示这个object的相关注释信息## IlluminaMethylationManifest object## Annotation## array: IlluminaHumanMethylation450k## Number of type I probes: 135476 ## Number of type II probes: 350036 ## Number of control probes: 850 ## Number of SNP type I probes: 25 ## Number of SNP type II probes: 40#####展示甲基化位点的探针信息head(getProbeInfo(manifest))## DataFrame with 6 rows and 8 columns## Name AddressA AddressB Color NextBase## <character> <character> <character> <character> <DNAStringSet>## 1 cg00050873 32735311 31717405 Red A## 2 cg00212031 29674443 38703326 Red T## 3 cg00213748 30703409 36767301 Red A## 4 cg00214611 69792329 46723459 Red A## 5 cg00455876 27653438 69732350 Red A## 6 cg01707559 45652402 64689504 Red A## ProbeSeqA ProbeSeqB nCpG## <DNAStringSet> <DNAStringSet> <integer>## 1 ACAAAAAAAC...ATAAACCCCA ACGAAAAAAC...ATAAACCCCG 2## 2 CCCAATTAAC...AAAACATACA CCCAATTAAC...AAAACGTACG 4## 3 TTTTAACACC...AAAAAAAACA TTTTAACGCC...AAAAAAAACG 3## 4 CTAACTTCCA...AACACAAACA CTAACTTCCG...AACGCGAACG 5## 5 AACTCTAAAC...AAAAAACTCA AACTCTAAAC...AAAAAACTCG 2## 6 ACAAATTAAA...ACAAAAAACA GCGAATTAAA...ACAAAAAACG 6####展示红色通道信息head(getRed(RGset))## 5723646052_R02C02 5723646052_R04C01 5723646052_R05C02## 10600313 1093 1160 681## 10600322 1968 2040 3571## 10600328 1646 2009 1297## 10600336 25841 31865 26460## 10600345 1671 1248 1124## 10600353 1057 2330 1129## 5723646053_R04C02 5723646053_R05C02 5723646053_R06C02## 10600313 1258 1518 1137## 10600322 2017 3175 1359## 10600328 2197 2206 1648## 10600336 22337 29089 26880## 10600345 944 1685 589## 10600353 1430 1580 1213
我们可以看到，RGset是IlluminaHumanMethylation450k进行分析的，注释比对用的基因组是ilmn12.hg19基因组。并且这个object中一共有135476 个type I型探针，350036个Type II型探针，850个control探针，SNP的type I和Type II型探针分别是24和40个。
好啦，今天就介绍到这里啦，这一期注释是对前期知识的储备和积累，所以显得比较啰嗦，能看到这里的同学为你们点个赞。下一次我将给大家分享接下来的四部分内容，展示如何使用一个包来处理整个甲基化芯片的分析流程。
好啦，我是风间琉璃，我们下期见~
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文