0．1简介

蜂鸟是一个使用全基因组亚硫酸氢盐测序(WGBS)或还原亚硫酸氢盐测序(RRBS)实验数据识别病例组和对照组之间差异甲基化区域(DMRs)的包。

的蜂鸟包使用贝叶斯隐马尔可夫模型(HMM)来检测DMRs。它每次适合一条染色体的贝叶斯HMM。的最终输出蜂鸟为在给定染色体中检测到的起始和结束位置的DMRs, DMRs的方向(高或低)，以及这些DMRs中的cpg的数量。

0.2功能

的蜂鸟包包含以下三个功能:

1.hummingbirdEM:该函数读取输入数据，设置初始值，执行贝叶斯HMM的期望最大化(EM)算法，并推断甲基化状态的最佳序列。

它以三个参数作为输入:对照组数据、案例组数据和容器大小。对这个函数的调用如下所示:hummingbirdEM(experimentInfoControl, experimentInfoCase, binSize = 40)，其中experimentInfoControl和experimentInfoCase分别包含控制组和用例组的输入数据。该输入数据包括整个染色体的每个CpG位置的甲基化读取数和未甲基化读取数。它们的格式需要是summarizeexperiment对象的格式。在第3节(样本数据集)中，我们提供了关于如何将数据组织到summarizeexperiment对象中的详细信息。第三个参数binSize是用户想要的bin大小。我们默认的bin大小是40个碱基对。料仓尺寸越小，DMR边界预测越准确。更大的容器大小导致更快的计算时间。默认的bin大小是通过平衡这两个因素来选择的。更多关于统计模型和如何选择一个好的箱子大小的详细信息可以在Ji(2019)中找到。

2.hummingbirdPostAdjustment:该函数通常在执行后执行hummingbirdEM．它允许研究人员对DMR提出三个额外的要求:1)DMR的最小长度，2)DMR中cpg的最小数量，3)DMR中任何两个相邻cpg之间的最大距离(以碱基对为单位)。

的hummingbirdPostAdjustment函数有六个参数。对这个函数的调用如下所示:hummingbirdPostAdjustment(experimentInfoControl, experimentInfoCase, emInfo, minCpGs = 10, minLength = 500, maxGap = 300)，其中experimentInfoControl和experimentInfoCase采用与函数相同的输入数据hummingbirdEM；emInfo是运行函数的结果hummingbirdEM；minCpGs, minLength, maxGap是前面提到的三个额外需求。默认值分别为10、500和300。

3.hummingbirdGraph:生成用户指定区域的观测和预测图。它通常在执行后调用hummingbirdEM而且hummingbirdPostAdjustment功能。

这个函数hummingbirdGraph需要五个参数。对这个函数的调用如下所示:hummingbirdGraph(experimentInfoControl, experimentInfoCase, postAdjInfoDMRs, coord1, coord2)，其中experimentInfoControl和experimentInfoCase为输入数据，如hummingbirdEM．postAdjInfoDMRs是从函数的结果中检测到的DMRs的读取hummingbirdPostAdjustment坐标1和坐标2是绘图的开始和结束基因组位置。这个函数的执行产生两个图形，我们称之为观察图和预测图。观察图显示病例组和对照组的平均甲基化率。预测图显示了按仓预测，其中“0”表示预测的正常仓;“1”表示预测的高甲基化bin;“-1”表示预测的低甲基化仓。

0.3样本数据集

该软件包提供了一个名为“exampleHummingbird”的示例数据集作为示例。

具体来说，它是Chen Z. et al(2017)描述的大后代综合征(LOS)研究中29号染色体的部分数据。本研究中WGBS实验的原始FASTQ文件可在基因表达综合(GEO)数据库中公开获取，登录号为。GSE93775。

在本节中，我们将使用这个示例数据来演示如何以正确的格式组织数据蜂鸟包中。我们的包需要R版本4.0和Rcpp包。

图书馆(genome ranges)图书馆(summarizeexperiment)图书馆(蜂鸟)数据(例如蜂鸟)

首先，我们分别使用“abnormUM”、“abnormM”、“normM”和“normUM”来表示四个矩阵，它们包含异常组的未甲基化读取数、异常组的甲基化读取数、正常组的甲基化读取数和正常组的未甲基化读取数。对于这四个矩阵中的每一个，每行都是一个CpG位置，每列都是一个生物复制(例如，一个病人、一只老鼠等)。在LOS研究中，异常组有4头牛，正常组也有4头牛。因此，这四个矩阵每个都包含四列。我们要求这四种基质只包含在相同基因组位置上共有的cpg。在分析之前，不被所有生物重复共享的cpg被移除。下面显示了范数矩阵的前6行。

头(normM)

# # [1] [2] [3] [4] # # [1] 8 7 12 10 # # [2] 4 4 2 4 # # [3] 0 1 0 4 # # [4] 2 2 0 2 # # [5] 1 1 1 1 # # [6] 8 0 0 7

我们的贝叶斯HMM对每个治疗组的生物重复数量没有最低要求。病例组(或异常组)和对照组(或正常组)可以有一个或多个重复。两组的重复次数可能不相等。

其次，我们使用向量pos来包含这些cpg在上述四个矩阵“abnormUM”、“abnormM”、“normM”和“normUM”中的基因组位置。

头(pos)

# # # # (1) (1) 271 # # [2] 331 # # [3] 363 # # [4] 386 # # (5) 418 # # [6] 464

使用蜂鸟包中，我们需要把四个矩阵和向量pos放在两个summarizeexperiment对象中，一个用于案例组，一个用于对照组。可以这样做:

pos <- pos[，1] assaysControl <- list(normM = normM, normUM = normUM) assaysCase <- list(abnormM = abnormM, abnormUM = abnormUM) exampleSEControl <- summarizeexperiment (assaysControl, rowRanges = GPos("chr29"， pos)) exampleSECase <- summarizeexperiment (assaysCase, rowRanges = GPos("chr29"， pos))

示例esecontrol和exampleSECase已经准备好供蜂鸟包中。

要在summarizeexperiment对象中显示数据，可以执行以下操作。

CpG的位置是:

rowRanges (exampleSEControl)

## UnstitchedGPos对象，4746个位置，0元数据列:## seqnames pos strand ##    ## [1] chr29 271 * ## [2] chr29 331 * ## [3] chr29 363 * ## [4] chr29 386 * ## [5] chr29 418 * ## ... ... ... ...## [4742] chr29 399795 * ## [4743] chr29 399802 * ## [4744] chr29 399833 * ## [4745] chr29 399864 * ## [4746] chr29 399987 * ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths

rowRanges (exampleSECase)

## UnstitchedGPos对象，4746个位置，0元数据列:## seqnames pos strand ##    ## [1] chr29 271 * ## [2] chr29 331 * ## [3] chr29 363 * ## [4] chr29 386 * ## [5] chr29 418 * ## ... ... ... ...## [4742] chr29 399795 * ## [4743] chr29 399802 * ## [4744] chr29 399833 * ## [4745] chr29 399864 * ## [4746] chr29 399987 * ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths

正常组甲基化和非甲基化读计数数据的矩阵如下:

头(化验(exampleSEControl)[[“normM”]])

# # [1] [2] [3] [4] # # [1] 8 7 12 10 # # [2] 4 4 2 4 # # [3] 0 1 0 4 # # [4] 2 2 0 2 # # [5] 1 1 1 1 # # [6] 8 0 0 7

头(化验(exampleSEControl)[[“normUM”]])

# # [1] [2] [3] [4] # # 7 [1] 4 2 4 # # [2] 12 11 11 10 # # [3] 10 10 8 7 # # 13 # # [4] 8 11 10 [5] 7 11 6 17 # # 8 9 7 8 [6]

异常组甲基化和非甲基化读计数数据的矩阵如下:

头(化验(exampleSECase)[[“abnormM”]])

# # [1] [2] [3] [4] # # [1] 10 7 10 13 # # [2] 6 2 6 8 # # 3 0 3 0 # # [3] [4] 0 1 1 0 # # [5] 1 1 2 2 # # 6 4 8 7 [6]

头(化验(exampleSECase)[[“abnormUM”]])

# # [1] [2] [3] [4] # # [1] 6 3 3 3 # # [2] 9 5 6 12 # # [3] 12 11 8 20 # # [4] 8 13 12 15 # # [5] 10 12 12 19 # # 8 7 6 14 [6]

0.4例子

本节使用上述示例数据集来展示如何使用我们的蜂鸟包来推断甲基化状态。首先，我们需要加载蜂鸟包和exampleHummingbird数据集。exampleHummingbird数据集包含summarizeexperiment对象、exampleSEControl和exampleSECase，这些对象已经准备好供蜂鸟包中。

库(蜂鸟)数据(exampleHummingbird)

emInfo <- hummingbirdEM(experimentInfoControl = exampleSEControl, experimentInfoCase = exampleSECase, binSize = 40)

##读取输入…##容器大小:40。##总行数:4746，正常组总重复数:4，异常组总重复数:4##处理输入完成…##计算初始值…##初始值计算…## EM开始…##迭代:1 ##迭代:2 ##迭代:3 ##迭代:4 ##迭代:5 ##迭代:6 ##迭代:7 ##迭代:8 ##迭代:9 ##迭代:10 ##迭代:11 ##迭代:13 ##迭代:14 ##迭代:16 ##迭代:17 ##迭代:18 ##迭代:19 ##迭代:20 ##迭代:21 #迭代:22 ##状态计算……EM经过22次迭代后收敛。##保存输出… ## ****** Program ended. ******

emInfo

## seqnames range对象，包含3296个范围和4个元数据列:## seqnames ranges strand | distance norm abnorm ##    |    [1] chr29 311-350 * | 40 0.672414 0.711864 ## [2] chr29 351-390 * | 40 0.156250 0.104348 ## [5] chr29 431-470 * | 40 0.333333 0.421875 ## ... ... ... ... . ... ... ...## [3292] chr29 399671- 399771 * | 40 0.250000 0.294118 ## [3294] chr29 399791-399830 * | 40 0.700000 0.743902 ## [3295] chr29 399831-399870 * | 40 0.448980 0.416667 ## [3296] chr29 399951-399990 * | 120 0.108696 0.268293 ##方向## <整数> ## [1]1 ## [2]1 ## [4]0 ## [5]0 ## # ... ...## [3292] 0 ## [3293] 0 ## [3294] 0 ## [3295] 0 ## [3296] 0 ## ------- ## seqinfo:来自未指定基因组的1个序列;没有seqlengths

emInfo是一个GenomicRanges对象，包含每个bin的起始和结束位置、当前bin与前面bin之间的距离、正常和异常基团的平均甲基化率以及预测的甲基化变化方向(“0”表示预测的正常bin;“1”表示预测的高甲基化仓;“-1”表示预测的低甲基化仓)。

hummingbirdPostAdjustment调整emInfo，使每个检测到的DMR具有用户定义的最小长度、最小cpg数量以及每个DMR中相邻cpg之间的最大间隙。如果用户没有定义，默认值为minLength=500, minCpGs=10, maxGap=300。

postAdjInfo <- hummingbirdPostAdjustment(experimentInfoControl = exampleSEControl, experimentInfoCase = exampleSECase, emInfo = emInfo, minCpGs = 10, minLength = 100, maxGap = 300)

##读取输入…最小cpg: 10，最小长度:100，最大差距:300。##工作地点差价调整开始…工作岗位调整完成…##输出DMRs…##共有3个DMRs。显示前3个。##区域:开始，结束，长度，方向，CpGs ## 1: 98391, 98590,200， - 1,10 ## 2: 107991, 108350, 360， - 1,12 ## 3: 110551, 110870,320， - 1,10 ##保存输出…## ******程序结束。＊＊＊＊＊＊

postAdjInfo dmr美元

## seqnames范围链|长度方向CpGs ##    |    ## [1] chr29 98391-98590 * | 200 -110 ## [2] chr29 107991-108350 * b| 360 -1 12 ## [3] chr29 110551-110870 * | 320 -110 ## ------- ## seqinfo: 1个序列来自一个未指定的基因组;没有seqlengths

postAdjInfo是两个GenomicRanges对象的列表，DMRs和obsPostAdj。具体来说，DMRs包含基于用户定义的参数(minLength、minCpGs和maxGap)的检测区域。它包含了细化后的DMRs，包括起始基因组位置、结束基因组位置、区域长度、预测甲基化变化方向(“0”表示无显著变化，“1”表示预测超甲基化，“-1”表示预测低甲基化)以及cpg的数量。obsPostAdj对象包含每个CpG位点的甲基化状态。

最后，我们使用蜂鸟图对用户定义的基因组区域进行可视化观察和预测。在观察图中，横轴表示基因组位置;纵轴分别显示每个箱子正常组和异常组的样本平均甲基化率。预测图显示了每个箱子中异常组和正常组之间的样本平均差值。数字(“0”，“1”，“-1”)表示预测。

接下来的两张图(前者是观察图，后者是预测图)将上述输出中的第二个DMR可视化。

蜂鸟图(experimentInfoControl = exampleSEControl, experimentInfoCase = exampleSECase, postAdjInfoDMRs = postAdjInfo$DMRs, coord1 = 107991, coord2 = 108350)

0．5引用

如果你使用蜂鸟包装，请引用以下论文，包括统计模型和拟合算法:

• Ji(2019)用于检测差异甲基化区域的贝叶斯隐马尔可夫模型。生物识别技术75(2):663 673。

0.6参考

LOS研究的真实数据来自以下论文:

• (2017)与DNA甲基化组外突变不耦合的基因的整体调控异常是先天性过度生长综合征的特征。科学报告7,12667。

0.7会话信息

本文分析的内容包括:

sessionInfo ()

## R版本4.2.1(2022-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。utf - 8 LC_IDENTIFICATION = C附加基本包:# # # # # # [1]stats4统计图形grDevices跑龙套数据集方法# # # # # #[8]基地其他附加包:# # [1]hummingbird_1.8.0 SummarizedExperiment_1.28.0 # # [3] Biobase_2.58.0 MatrixGenerics_1.10.0 # # [5] matrixStats_0.62.0 GenomicRanges_1.50.0 # # [7] GenomeInfoDb_1.34.0 IRanges_2.32.0 # # [9] S4Vectors_0.36.0 BiocGenerics_0.44.0 # # [11] BiocStyle_2.26.0 # # # #通过加载一个名称空间(而不是附加):## [1] Rcpp_1.0.9 highr_0.9 bslib_0.4.0 ## [4] compiler_4.2.1 BiocManager_1.30.19 jquerylib_0.1.4 ## [7] XVector_0.38.0 bitops_1.0-7 tools_4.2.1 ## [10] jsonlite_1.8.3 evaluate_0.17 rlang_1.0.6 ## [13] Matrix_1.5-1 DelayedArray_0.24.0 cli_3.4.1 ## [19] magick_2.7.3 yaml_2.3.6 xfun_0.34 ## [22] fastmap_1.1.0 GenomeInfoDbData_1.2.9 string_1 .4.1 ## [28] R6_2.5.1 rmarkdown_2.17 bookdown_0.29 ## [31]RCurl_1.98-1.9 cachem_1.0.6 . htmltools_0.5.3 stringi_1.7.8