1简介

methcp是一种全基因组亚硫酸氢盐测序(WGBS)数据的差异甲基化区(DMR)检测方法。适用于广泛的实验设计。在本文档中，我们提供了两组比较和时间过程分析的例子。

methcp基于变化点检测识别DMRs，自然地分割基因组并提供区域级差异分析。我们把感兴趣的读者引向我们的报纸在这里

加载包。

库(bsseq)库(MethCP)
##警告:替换以前的导入'matrixStats::rowMedians'由## 'Biobase::rowMedians'加载'DSS'时
##警告:取代以前的导入'matrixStats::anyMissing'由## 'Biobase::anyMissing'时加载'DSS'

2双官能团的比较

在本节中，我们使用从GEO获得的具有登录号的拟南芥数据集中的CpG甲基化数据GSM954584．我们从每个样本中提取1号和2号染色体的子集，并在野生型植物和H2A之间进行差异分析。Z突变植物，我们称之为治疗而且控制在文件的其余部分。

2.1读取数据

我们使用完善的Bioconductor包bsseq加载存储原始数据。下面是如何使用读取原始计数的示例bsseq．

我们提供了一个辅助函数createBsseqObject当每个示例的数据存储在单独的文本文件中时，创建bsseq对象。

有关的更多操作bsseq对象，或创建bsseq对象自定义为您的文件格式，请参考其用户指南．

#数据集由6个样本组成。3个样品为H2A。Z突变#植株，3个样品为对照。sample_names <- c(paste0("control"， seq_len(3))， paste0("treatment"， seq_len(3)))) #获取文件路径和文件名的向量raw_files <- system。file("extdata"， paste0(sample_names， ".txt")， package = "MethCP") #加载数据bs_object <- createBsseqObject(files = raw_files, sample_names = sample_names, chr_col = 'Chr'， pos_col = 'Pos'， m_col = "M"， cov_col = 'Cov')

的bsseq对象。

bs_object
##一个“BSseq”类型的对象，其中## 82103甲基化位点## 6个样本##尚未平滑##所有分析都在内存中

其中一个示例的原始文件头。

Dt <- read。table(raw_files[1]， stringsAsFactors = FALSE, header = TRUE) head(dt)
## Chr Pos M Cov ## 1 chr1 109 3 3 ## 2 chr1 110 4 4 ## 3 chr1 115 3 3 ## 4 chr1 116 4 4 ## 5 chr1 161 0 1 ## 6 chr1 162 11 1

2．2计算统计数据

我们计算每个胞嘧啶统计使用两种不同的测试DSS而且methylKit使用函数calcLociStat．此函数返回一个MethCP对象，该对象用于未来的分割步骤。当数据集中有多个染色体时，我们允许并行计算。

#两个组的样本名称进行比较。它们应该是创建' bsseq '对象时提供的#示例名称的子集。group1 <- paste0("control"， seq_len(3)) group2 <- paste0("treatment"， seq_len(3)) #下面我们使用两个不同的# test ' DSS '和' methylKit '来计算每个胞嘧啶的统计量。用户可以为他们的#应用程序选择其中之一。# obj_DSS <- calcLociStat(bs_object, group1, group2, test = "DSS") obj_methylKit (bs_object, group1, group2, test = "methylKit")
obj_methylKit
## MethCP对象，2染色体，81936甲基化位点##测试:methylKit ## group1: control1 control2 control3 ## group2: treatment1 treatment2 treatment3 ##没有分段

如果用户希望将其预先计算的测试统计数据用于两组比较和时间过程数据以外的实验，我们将使用计算的统计数据并创建一个MethCP对象。

data <- data.frame(chr = rep("Chr01"， 5)， pos = c(2,5,9,10,18)，效果。size = c(1，-1, NA, 9, Inf)， pvals = c(0, 0.1, 0.9, NA, 0.02)) obj <- MethCPFromStat(data, test.name="myTest"， pvals. name="myTest")Field =" pvals"， effect.size. Field ="效果。”,seqnames大小。Field ="chr"， pos. Field ="pos")
##过滤NAs和无限值。##过滤前胞嘧啶计数:5。##过滤后胞嘧啶计数:3。
obj
## MethCP对象，1条染色体，3个甲基化位点## test: myTest ## group1:不适用## group2:不适用##未分割

2．3分割

segmentMethCP对MethCP对象。当数据集中有多个染色体时，我们允许并行计算。不同于calcLociStat函数中，我们没有对使用的核数施加任何限制。请参阅文档以调整分段中使用的参数。

# obj_DSS <- segmentMethCP(# obj_DSS, bs_object，区域。测试= "weighted-coverage") obj_methylKit <- segmentMethCP(obj_methylKit, bs_object，区域。Test = "fisher")
obj_methylKit
## MethCP对象与2染色体，81936甲基化位点##测试:methylKit ## group1: control1 control2 control3 ## group2: treatment1 treatment2 treatment3 ##已分割

2.4重要的地区

使用功能getSigRegion在一个MethCP对象获取DMRs的列表。

# region_DSS <- getSigRegion(obj_DSS) # head(region_DSS)
region_methylKit <- getSigRegion(obj_methylKit) head(region_methylKit)
## seqnames开始结束nC。有效nC均值。diff的意思。浸地区。Pval ## 1.2 chr1 26853 27083 10 10 -0.3126 6.8500 6.306067e-14 ## 1.5 chr1 30691 30768 11 11 -0.1776 8.3788 4.028997e- 1.13 chr1 47817 47966 16 16 -0.2058 7.1354 4.314515e-08 ## 1.21 chr1 65562 65735 10 10 -0.5742 5.2167 0.000000e+00 ## 1.38 chr1 92672 92891 10 10 0.4067 8.9667 0.000000e+00 ## 1.40 chr1 94712 95021 12 12 0.1304 7.9444 1.295822e-10

3.时间过程的例子

MethCP适用于各种各样的实验设计。我们将MethCP应用于从GEO获得的具有登录号的拟南芥种子萌发数据集https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE94712数据来自两个重复，分别是野生型植物(Col-0)和ros1 dml2 dml3 (rdd)三去甲基酶突变型植物的干种子和发芽种子，在吸胀4天(DAI)后的0-4天。对于基于胞嘧啶的统计，我们拟合甲基化比率的线性模型，并测试条件Col-0和条件rdd的时间系数之间的差异。

读取元数据。

Meta_file <- system。file("extdata"， "meta_data.txt"， package = "MethCP") meta <- read。table(meta_file, sep = "\t"， header = TRUE)
条件时间名时间复制采样名1 Col-0 DS 0 1 GSM2481288_mC_Col-0_DS_r1 ## 2 GSM2481289_mC_Col-0_DS_r2 ## 3 Col-0 0DAI 1 1 gsm2481290_mc_col - 0_dai_r1 ## 4 Col-0 0DAI 12 gsm2481291_mc_col - 0_dai_r2 ## 5 Col-0 1DAI 2 1 gsm2481292_mc_col - 0_dai_r1 ## 6 Col-0 1DAI 2 2 gsm2481293_mc_col - 0_dai_r1 ## 2 gsm2481293_mc_col - 0_dai_r1 ## 2

读取计数数据。

#获取文件路径的向量和文件名raw_files <- system。file("extdata"， paste0(meta$SampleName， ".tsv")， package = "MethCP") # read files bs_object <- createBsseqObject(files = raw_files, sample_names = meta$SampleName, chr_col = 1, pos_col = 2, m_col = 4, cov_col = 5, header = TRUE)

应用覆盖过滤器，以确保每个位点对每个条件的总覆盖率(在样本间的总和)大于3。

groups <- split(seq_len(nrow(meta))， meta$Condition) coverages <- as.data.frame(getCoverage(bs_object, type = "Cov")) filter <- rowsum (coverages[， meta$SampleName[groups[[1]]]] != 0) >= 3 & rowsum (coverages[， meta$SampleName[groups[[2]]]] != 0) >= 3 bs_object <- bs_object[filter，]

计算统计数据。元数据的数据帧将被传递给函数calcLociStatTimeCourse．注意，必须有命名的列条件，时间而且SampleName在数据框架中。

obj <- calcLociStatTimeCourse(bs_object, meta)
obj
组1:GSM2481288_mC_Col-0_DS_r1 GSM2481289_mC_Col-0_DS_r2 gsm2481291_mc_col - 0_dai_r1 gsm2481292_mc_col - 0_dai_r2 gsm2481292_mc_col - 0_dai_r1 gsm2481292_mc_col - 0_dai_r2 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r2 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r2 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1 gsm248129_mc_col - 0_dai_r1GSM2481300_mC_rdd_DS_r1 GSM2481301_mC_rdd_DS_r2 GSM2481302_mC_rdd_0DAI_r1 GSM2481303_mC_rdd_0DAI_r2 GSM2481304_mC_rdd_1DAI_r1 GSM2481305_mC_rdd_1DAI_r2 GSM2481306_mC_rdd_2DAI_r1 GSM2481307_mC_rdd_2DAI_r2 GSM2481308_mC_rdd_3DAI_r1 GSM2481309_mC_rdd_3DAI_r2 GSM2481310_mC_rdd_4DAI_r1 GSM2481311_mC_rdd_4DAI_r2 # #没有被分割

分割。

obj <- segmentMethCP(obj, bs_object, region. obj)Test = "stouffer")

去拿DMRs。

regions <- getSigRegion(obj)
头(地区)
## seqnames开始结束nC。有效nC均值。diff的意思。浸地区。Pval ## 31 1 238367 239326 69 69 -0.1042 5.7144 0.0012861236 ## 66 1 394175 395220 50 50 -0.2490 3.9292 0.0026789957 ## 70 1 405983 408116 52 52 -0.1079 6.3197 0.0082183067 ## 92 1 531037 532534 36 36 -0.1900 4.3507 0.0007397345 ## 187 1 962256 966028 145 145 -0.1049 5.9529 0.0029463070 ## 241 1 1258866 1259519 30 30 -0.1454 5.6153 0.0028433185

4会话信息

sessionInfo ()
## R版本4.2.1(2022-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.4 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]MethCP_1.11.0 bsseq_1.33.0 # # [3] SummarizedExperiment_1.27.2 Biobase_2.57.1 # # [5] MatrixGenerics_1.9.1 matrixStats_0.62.0 # # [7] GenomicRanges_1.49.1 GenomeInfoDb_1.33.5 # # [9] IRanges_2.31.2 S4Vectors_0.35.1 # # [11] BiocGenerics_0.43.1 BiocStyle_2.25.0 # # # #通过加载一个名称空间(而不是附加):## [11] permute_0.9-7 rhdf5filters_1.9.0 ## [13] fastseg_1.43.5 DNAcopy_1.71.0 ## [15] tidyselect_1.1.2 compiler_4.2.1 ## [17] cli_3.3.0 DelayedArray_0.23.1 ## [19] rtracklayer_1.57.0 bookdown_0.28 ## [21] sass_0.4.2 scales_1.2.1 ## [23] mvtnorm_1.1-3 string_1 .4.1 ## [25] digest_0.6.29 Rsamtools_2.13.4 ## [27] rmarkdown_2.15 ## [5] bslib_20.4 utf8_1.2.2 ## [27] rmarkdown_2.15 ##R.utils_2.12.0 # # [29] XVector_0.37.0 pkgconfig_2.0.3 # # [31] htmltools_0.5.3 sparseMatrixStats_1.9.0 # # [33] fastmap_1.1.0 bbmle_1.0.25 # # [35] limma_3.53.6 BSgenome_1.65.2 # # [37] rlang_1.0.4 DelayedMatrixStats_1.19.0 # # [39] generics_0.1.3 jquerylib_0.1.4 # # [41] BiocIO_1.7.1 jsonlite_1.8.0 # # [43] mclust_5.4.10 BiocParallel_1.31.12 # # [45] gtools_3.9.3 dplyr_1.0.9 # # [47] R.oo_1.25.0 rcurl_1.98 - 1.8 # # [49] magrittr_2.0.3 GenomeInfoDbData_1.2.8 # # [51] Matrix_1.4-1 fansi_1.0.3 # # [53]Rcpp_1.0.9 munsell_0.5.0 ## [55] Rhdf5lib_1.19.2 lifecycle_1.0.1 ## [59] yaml_2.3.5 MASS_7.3-58.1 ## [61] zlibbioc_1.43.0 rhdf5_2.41.1 ## [63] plyr_1.8.7 qvalue_2.29.0 ## [67] bdsmatrix_1.3 .3-6 crayon_1.5.1 ## [69] lattice_0.20-45 Biostrings_2.65.2 ## [71] splines_4.2.1 methylKit_1.23.0 ## [73] locfit_1. 40 ## [75] pillar_1.8.1 rjson_0.2.21 ## [77] reshape_1 .4.4 codetools_0.2-18 ## [79] glue_1.6.2## [81] evaluate_0.16 data.table_1.14.2 ## [83] BiocManager_1.30.18 vctrs_0.4.1 ## [85] purrr_0.3.4 gtable_0.3.0 ## [87] assertthat_0.2.1 cachem_1.0.6 ## [89] ggplot2_3.3.6 emdbook_1.3.12 ## [91] xfun_0.32 restfulr_0.0.15 ## [93] coda_0.19-4 tibble_3.1.8 ## [95] GenomicAlignments_1.33.1