minfi 1.36.0
的minfi包提供了分析Illumina甲基化阵列的工具,特别是450k和EPIC(也称为850k)阵列。我们对旧的27k数组有部分支持。
在这个包中处理的任务包括预处理,QC评估,感兴趣的甲基化位点的识别和绘图功能。分析这些类型的阵列是我们和其他小组正在进行的研究。
这个包的输入数据是IDAT文件,在标准化之前表示两个不同的颜色通道。这是最完整的数据类型,因为IDAT文件包括控制探针上的测量值。可以将Genome Studio文件与此包中包含的数据结构一起使用,但只有一些功能可用,因为Genome Studio输出不包含控制探针信息。此外,通常Genome Studio的输出是使用在Genome Studio中实现的方法进行标准化的,这些方法通常被认为是较差的。
MINFI包包含由不同的非重叠作者在多个手稿中描述的方法。这使得引用包有点困难,所以这里有一个指南。
getSex ()
,质量控制使用getQC ()
,分位数归一化使用preprocessQuantile ()
,使用bumphunter ()
块查找使用blockFinder ()
.convertArray ()
而且combineArrays ()
,扩展NOOB工作在单样本模式以及使用estimateCellCounts ()
使用来自450k数组的参考数据来估计EPIC数据的单元类型组成。preprocessQuantile ()
,这将是体贴的也引用(Touleimat and Tost 2012),因为本文描述的方法本质上与preprocessQuantile ()
.bumphunter ()
,你也应该引用原始的bump hunter出版物(Jaffe et al. 2012).preprocessSWAN ()
请引用(Maksimovic, Gordon, and Oshlack 2012).preprocessNoob ()
,请引用(Triche et al. 2013).preprocessFunnorm ()
,请引用(Fortin et al. 2014).中默认的preprocessFunnorm ()
就是做新手背景校正。如使用,请同时注明(Triche et al. 2013).隔间()
而且extractAB ()
,请引用(Fortin and Hansen 2015).如果你使用Bibtex,你可以得到引用在这种格式通过
toBibtex(引用(“minfi”))
文献往往有点不具体。DNA甲基化微阵列的术语。
对于450k微阵列,每个样品都在一个阵列上测量,在两个不同的颜色通道(红色和绿色)上。每个阵列测量大约45万个CpG位置。每个CpG都与两个测量相关联:甲基化测量和“非”甲基化测量。这两个值可以用两种方法之一来测量:使用“类型I”设计或“类型II设计”。使用I型设计测量的cpg使用单一颜色测量,在同一颜色通道中使用两个不同的探针提供甲基化和非甲基化测量。使用II型设计测量的cpg使用单个探针测量,两种不同的颜色提供甲基化和非甲基化的测量值。实际上,这意味着在这个数组上有不探针和CpG位置之间的一对一对应关系。因此,我们试图精确地描述这一点,当我们提到单碱基基因组位点时,我们指的是“位点”(或“CpG”),我们将其与“探针”区分开来。上一代27k甲基化阵列仅使用Type I设计,而EPIC阵列同时使用Type I和Type II。
样本之间DNA甲基化的差异既可以是在单个CpG上,称为差异甲基化位置(DMP),也可以是在区域水平上,称为差异甲基化区域(DMR)。
物理上,每个样本都是在一个“阵列”上测量的。对于450k的设计,在一个物理“幻灯片”上有12个数组(以6 × 2的网格组织)。幻灯片被组织成最多包含8个幻灯片(96个数组)的“板”。EPIC阵列每个幻灯片有8个阵列,每个板有64个阵列。
该文档具有以下依赖项
库(minfi)库(minfiData)
MINFI包被设计为非常灵活的方法开发人员。2021欧洲杯体育投注开户这种灵活性对用户来说是有代价的;他们需要了解一些不同的数据类:
RGChannelSet
: IDAT文件的原始数据;这些数据是在探针(而不是CpG位点)水平上组织的。该数据有两个通道:红色和绿色。MethylSet
:由CpG位点水平组织的数据,但未映射到基因组。该数据有两个通道:Meth(甲基化)和Unmeth(非甲基化)。RatioSet
:由CpG位点水平组织的数据,但未映射到基因组。数据至少具有两个通道中的一个:Beta和/或M (Beta logratio)。它可以选择包含CN通道(拷贝号)。GenomicMethylSet
:像MethylSet
,但映射到一个基因组。GenomicRatioSet
:像RatioSet
,但映射到基因组。类层次结构如下所示
为了更清楚地说明这一点,让我们看一下示例数据RGsetEx
从minfiData包中。这是我们在类层次结构中移动时的特性和类的数量(代码未运行):
RGsetEx
##类:RGChannelSet ## dim: 622399 6 ##元数据(0):## assays(2):绿红## 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…注释:IlluminaHumanMethylation450k注释:ilmn12.hg19
## RGsetEx: RGChannelSet, 622,399 features MsetEx <- preprocessRaw(RGsetEx) ## MsetEx: MethylSet, 485,512 features GMsetEx <- mapToGenome(MsetEx) ## GMsetEx: GenomicMethylSet, 485,512 features
注意特性的数量是如何变化的。在RGChannelSet
一个特征是一个探针(它不同于CpG轨迹,请参阅术语部分)。在MethylSet
每个特征现在都是一个甲基化位点,它的特征更少,因为一些位点是用两个探针测量的。最后,GenomicMethylSet
尺寸和MethylSet
,但原则上它可以更小,以防你使用的注释说一些探针不能映射到基因组(在这种情况下是hg19)。
最后我们可以转换成aRatioSet
通过ratioConvert ()
.这两个函数ratioConvert ()
而且mapToGenome ()
如上面的类层次结构图所示。许多预处理函数在图中经历了几个步骤,例如,如果函数需要知道探针的基因组位置(一些预处理函数以不同的方式处理在性染色体上测量的探针)。
此包支持对IDAT文件的分析,其中包含汇总的珠子信息。
根据我们的经验,大多数实验室使用“样本表”CSV文件来描述实验的布局。这是基于Illumina提供的示例表文件。我们的管道假设存在这样一个文件,但是如果它不可用的话,使用例如Excel创建这样一个文件相对容易。
我们使用了一个包含6个样本的示例数据集,分布在两张幻灯片上。首先,我们获得IDAT文件的系统路径;这需要一点,因为数据来自已安装的包
baseDir <- system。file("extdata", package = "minfiData") list.files(baseDir)
## [1] "5723646052" "5723646053" "SampleSheet.csv"
这显示了450k数据的典型布局:每个“幻灯片”(包含12个数组,参见术语)存储在一个单独的目录中,具有一个数字名称。顶层目录包含示例表文件。在幻灯片目录中,我们可以找到IDAT文件(可能还有一些JPG图像或其他文件):
list.files(文件。路径(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”
每个数组的文件都有另一个数字,由一个红色和一个Grn(绿色)IDAT文件组成。注意,对于这个示例数据,每张幻灯片只包含3个数组,而不是12个。这样做是因为文件大小的限制,因为我们只需要6个数组来说明包的功能。
首先我们读一下样品单。我们提供了一个方便的函数来读取这个文件read.metharray.sheet ()
.这个函数有几个吸引人的附加功能。让我们看看输出
目标<- read.metharray.sheet(baseDir)
# # (read.metharray。找到以下CSV文件:
## [1] "/home/biocbuild/bb -3.12-bioc/R/library/minfiData/extdata/SampleSheet.csv"
目标
# # Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID人年龄性别# # 1 GroupA_3 H5 < NA > GroupA < NA > id3 83 # # 2 GroupA_2 D5 < NA > GroupA F < NA > id2 58 # # 3 GroupB_3 C6 < NA > GroupB < NA > id3 83 # # 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 F < NA > id2 58 # #状态数组幻灯片# # 1正常R02C02 5723646052 # # 2正常R04C01 5723646052 # # 3癌症R05C02 5723646052 # # 4正常R05C02 R04C02 5723646053 # # 51 /home/biocbuild/ bbbs -3.12-bioc/R/library/minfiData/extdata/5723646052/5723646052_R02C02 ## 2/ home/biocbuild/ bbbs -3.12-bioc/R/library/minfiData/extdata/5723646052/5723646052_R04C01 ## 3/ home/biocbuild/ bbbs -3.12-bioc/R/library/minfiData/extdata/5723646052/5723646052_R05C02 ## 4 /home/biocbuild/ bbbs -3.12-bioc/R/library/minfiData/extdata/5723646052/5723646052 _r04c02 ##6 /home/biocbuild/ bbs3.12 -bioc/R/library/minfiData/extdata/5723646053/5723646053_R06C02 ## 6 /home/biocbuild/ bbs3.12 -bioc/R/library/minfiData/extdata/5723646053/5723646053_R06C02
首先输出:这是adata.frame
.它包含一列Basename
描述了对应于示例的IDAT文件的位置,以及两列数组
而且幻灯片
.在Illumina提供的样例表中,命名了这两列Sentrix_Position
而且Sentrix_ID
,但我们重新命名它们。下面将提供关于此函数使用的更多细节。的Basename
列往往太大,不便于显示,这里是相对简化baseDir
:
sub(baseDir, "",目标$Basename)
## [1] "/ 5723646052/5723646052_r02c02 " "/ 5723646052/5723646052_r04c01 " ## [3] "/ 5723646052/5723646052_r05c02 " " ## [5] "/ 5723646053/5723646053_r05c02 "
(这只是为了展示)。
用这个data.frame
,从数据中很容易读出
RGset <- read.metharray。Exp(目标=目标)
## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串#### readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串##
让我们看看相关的pheno数据,这实际上只是上面的targets对象中包含的信息。
RGset
##类:RGChannelSet ## dim: 622399 6 ##元数据(0):## assays(2):绿红## 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…注释:IlluminaHumanMethylation450k注释:ilmn12.hg19
pd <- pData(RGset) pd[,1:4]
##数据帧,6行4列## Sample_Name Sample_Well Sample_Plate Sample_Group ## <字符> <字符> <字符> <字符> ## 5723646052_R04C01 GroupA_2 H5 NA groupa# # 5723646052_R04C01 GroupA_2 D5 NA groupa# # 5723646052_r04c02 GroupB_3 C6 NA groupb# # 5723646053_R04C02 GroupB_1 F7 NA groupb# # 5723646053_R06C02 GroupB_2 H7 NA groupb# # 5723646053_R06C02 GroupB_2 H7 NA groupb# # 5723646053_R06C02 GroupB_2 H7 NA groupb# # 5723646053_R06C02 GroupB_2
的read.metharray.exp ()
函数还可以读取整个目录或目录树(使用递归
设置为真正的
)通过只使用函数的实参基地
而且目标=零
,就像
RGset2 <- read.metharray.exp(文件。路径(baseDir,“5723646052”))
## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串##
RGset3 <- read.metharray。exp(baseDir, recursive = TRUE)
## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串#### readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串## readChar警告(con, nchars = n):嵌入空的截断字符串## readChar(con, nchars = n):嵌入空的截断字符串##
表中唯一重要的列data.frame
用于目标
的参数read.metharray.exp ()
函数是列名Basename
.通常,这样的对象也有命名的列数组
,幻灯片
,和(可选)板
.
我们使用了基于Illumina提供的样例表数据文件构建的表数据文件。这是一个带有头的CSV文件。在这种情况下,我们假设表型数据从第一行开始(数据)
(或者没有头文件)。
使用该函数手动读取示例表也很容易read.csv ()
.在这里,我们知道要跳过文件的前7行。
Targets2 <- read.csv(file. csv)targets2 . path(baseDir, "SampleSheet.csv"), stringsAsFactors = FALSE, skip = 7
# # Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID Sentrix_ID # # 1 GroupA_3 H5 NA GroupA NA 5723646052 # # 2 GroupA_2 D5 NA GroupA NA 5723646052 # # 3 GroupB_3 C6 NA GroupB NA 5723646052 # # 4 GroupB_1 F7 NA GroupB NA 5723646052 # # 5 GroupA_1 G7 NA GroupA NA 5723646053 # # 6 GroupB_2 H7 NA GroupB NA 5723646053 # # Sentrix_Position人年龄性别地位# # 1正常R02C02 id3 83 # # 2 R04C01 id2 58 # # 3 F正常R05C02 id3 83癌症# # 4 R04C02 id1 75 F # # 5 R05C02 id1 75 F正常R06C02 id2 58 F癌症
我们现在需要填充Basename
列。关于可能的方法如下
targets2$Basename <- file。path(baseDir, targets2$Sentrix_ID, paste0(targets2$Sentrix_ID, targets2$Sentrix_Position))
最后,MINFI包含一个基于文件的解析器:read.metharray ()
.返回对象表示样本的红色和绿色通道测量值。我们从包中得到的一个有用的函数Biobase是结合()
这将组合(“添加”)两组样本。这允许用户手动建立一个RGChannelSet
.
对于甲基化阵列,我们有两种类型的注释包:包含阵列设计的“manifest”包和包含甲基化位点在基因组上的位置、它们映射到哪些基因组特征以及它们是否可能重叠任何已知snp的“annotation”包。
您可以看到哪些包正在被使用
注释(RGsetEx)
## IlluminaHumanMethylation450k“ilmn12.hg19”
本讨论面向希望了解MINFI内部原理的包开发人员或用户。2021欧洲杯体育投注开户
一组450k的数据文件最初将被读入RGChannelSet
,表示原始强度为两个矩阵:一个是绿色通道,一个是红色通道。这是一个非常类似于ExpressionSet
或者一个NChannelSet
.的RGChannelSet
是,在一起了吗IlluminaMethylationManifest
对象,预处理为MethylSet
.的IlluminaMethylationManifest
对象包含阵列设计,并描述探针和颜色通道如何配对在一起以测量特定CpG的甲基化水平。对象还包含关于控制探测(也称为QC探测)的信息。的MethylSet
包含标准化数据,基本上由两个矩阵组成,其中包含每个CpG的甲基化和非甲基化证据。只有RGChannelSet
包含关于控制探测的信息。
方法分析Affymetrix表达式数组的范例与前一段中描述的过程非常相似affy包(一个AffyBatch
被预处理成ExpressionSet
使用存储在CDF环境(包)中的数组设计信息)。
MINFI中的预处理是由一系列名为preprocessXX
.不同的函数有不同的输出类。
目前,我们有
preprocessRaw
:不处理。preprocessIllumina
: Illumina预处理,由Genome Studio执行(由我们进行逆向工程)。preprocessSWAN
: SWAN归一化,在(Maksimovic, Gordon, and Oshlack 2012).preprocessQuantile
:分位数归一化(适用于DNA甲基化阵列),在(Touleimat and Tost 2012,阿尔耶等人(2014))preprocessNoob
: Noob预处理,详见(Triche et al. 2013).preprocessFunnorm
:函数归一化,如(Fortin et al. 2014).补救办法:讨论文学
Aryee, Martin J, Andrew E Jaffe, Hector Corrada Bravo, Christine Ladd-Acosta, Andrew P Feinberg, Kasper D Hansen和Rafael A Irizarry. 2014。“Minfi:用于Infinium DNA甲基化微阵列分析的灵活而全面的Bioconductor包。”生物信息学30(10): 1363 - 9。https://doi.org/10.1093/bioinformatics/btu049.
让-菲利普·福廷,卡斯珀·D·汉森,2015。“利用表观遗传数据中的长程相关性重建Hi-C揭示的A/B区室。”基因组生物学16:180。https://doi.org/10.1186/s13059-015-0741-y.
Fortin, Jean-Philippe, Aurélie Labbe, Mathieu Lemire, Brent W Zanke, Thomas J Hudson, Elana J Frtig, Celia MT Greenwood, Kasper D Hansen. 2014。“450k甲基化阵列数据的功能规范化提高了大型癌症研究的复制能力。”基因组生物学15(11): 503。https://doi.org/10.1186/s13059-014-0503-2.
福廷,让-菲利普,小蒂莫西·特里切,卡斯珀·D·汉森,2017。Illumina人类甲基化epic阵列与Minfi的预处理、标准化和集成生物信息学33(4): 558 - 60。https://doi.org/10.1093/bioinformatics/btw691.
贾菲,安德鲁·E,彼得·村上,李华金,杰弗里·T·莱克,M·丹尼尔·法林,安德鲁·P·范伯格,拉斐尔·A·伊瑞扎里,2012。“在表观遗传流行病学研究中寻找差异甲基化区域。”国际流行病学杂志41(1): 200 - 209。https://doi.org/10.1093/ije/dyr238.
Maksimovic, Jovana, Lavinia Gordon和Alicia Oshlack, 2012。“SWAN: Illumina Infinium HumanMethylation450串珠芯片阵列归一化子集分位数。”基因组生物学13 (6): R44机身内部。https://doi.org/10.1186/gb-2012-13-6-r44.
Touleimat, Nizar, Jörg Tost, 2012。“使用子集分位数归一化进行精确DNA甲基化估计的Infinium()人类甲基化450K珠子芯片数据处理的完整管道。”表观基因组学4(3): 325 - 41。https://doi.org/10.2217/epi.12.21.
特里谢,蒂莫西J,丹尼尔J魏森伯格,大卫范登伯格,彼得W莱尔德,金伯利D齐格蒙德。2013。Illumina Infinium DNA甲基化串珠阵列的低水平加工核酸研究41 (7): e90。https://doi.org/10.1093/nar/gkt090.