Methrix tutorial

#Load library
library(methrix)

#Genome of your preference to work with
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

library(BiocManager)

if(!requireNamespace("BSgenome.Hsapiens.UCSC.hg19")) {
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")
}
library(BSgenome.Hsapiens.UCSC.hg19) 

#Example bedgraph files
bdg_files <- list.files(
  path = system.file('extdata', package = 'methrix'),
  pattern = "*bedGraph\\.gz$",
  full.names = TRUE
)

print(basename(bdg_files))
#> [1] "C1.bedGraph.gz" "C2.bedGraph.gz" "N1.bedGraph.gz" "N2.bedGraph.gz"

#Generate some sample annotation table
sample_anno <- data.frame(
  row.names = gsub(
    pattern = "\\.bedGraph\\.gz$",
    replacement = "",
    x = basename(bdg_files)
  ),
  Condition = c("cancer", 'cancer', "normal", "normal"),
  Pair = c("pair1", "pair2", "pair1", "pair2"),
  stringsAsFactors = FALSE
)

print(sample_anno)
#>    Condition  Pair
#> C1    cancer pair1
#> C2    cancer pair2
#> N1    normal pair1
#> N2    normal pair2

#First extract genome wide CpGs from the desired reference genome
hg19_cpgs <- suppressWarnings(methrix::extract_CPGs(ref_genome = "BSgenome.Hsapiens.UCSC.hg19"))
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit
#> 
#> Attaching package: 'rtracklayer'
#> The following object is masked from 'package:BiocIO':
#> 
#>     FileForFormat
#> -Extracting CpGs
#> -Done. Extracted 28,217,448 CpGs from 25 contigs.

#Read the files 
meth <- methrix::read_bedgraphs(
  files = bdg_files,
  ref_cpgs = hg19_cpgs,
  chr_idx = 1,
  start_idx = 2,
  M_idx = 3,
  U_idx = 4,
  stranded = FALSE,
  zero_based = FALSE, 
  collapse_strands = FALSE, 
  coldata = sample_anno
)
#> ----------------------------
#> -Preset:        Custom
#> --Missing beta and coverage info. Estimating them from M and U values
#> -CpGs raw:      28,217,448 (total reference CpGs)
#> -CpGs retained: 28,217,448(reference CpGs from contigs of interest)
#> ----------------------------
#> -Processing:    C1.bedGraph.gz
#> --CpGs missing: 28,216,771 (from known reference CpGs)
#> -Processing:    C2.bedGraph.gz
#> --CpGs missing: 28,216,759 (from known reference CpGs)
#> -Processing:    N1.bedGraph.gz
#> --CpGs missing: 28,216,746 (from known reference CpGs)
#> -Processing:    N2.bedGraph.gz
#> --CpGs missing: 28,216,747 (from known reference CpGs)
#> -Finished in:  00:01:15 elapsed (51.3s cpu)

#Typing meth shows basic summary.
meth
#> An object of class methrix
#>    n_CpGs: 28,217,448
#> n_samples: 4
#>     is_h5: FALSE
#> Reference: hg19

methrix::methrix_report(meth = meth, output_dir = tempdir())

mpca <- methrix_pca(m = methrix_data, do_plot = FALSE)

#Plot PCA results
plot_pca(pca_res = mpca, show_labels = TRUE)

#Color code by an annotation
plot_pca(pca_res = mpca, m = methrix_data, col_anno = "Condition")

if(!require(bsseq)) {
BiocManager::install("bsseq")}

library(bsseq)
bs_seq <- methrix::methrix2bsseq(m = methrix_data)

bs_seq
#> An object of type 'BSseq' with
#>   743 methylation loci
#>   4 samples
#> has not been smoothed
#> All assays are in-memory

sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#> 
#> 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_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] bsseq_1.42.0                         MafDb.1Kgenomes.phase3.hs37d5_3.10.0
#>  [3] GenomicScores_2.18.0                 BSgenome.Hsapiens.UCSC.hg19_1.4.3   
#>  [5] BSgenome_1.74.0                      rtracklayer_1.66.0                  
#>  [7] BiocIO_1.16.0                        Biostrings_2.74.0                   
#>  [9] XVector_0.46.0                       methrix_1.20.0                      
#> [11] SummarizedExperiment_1.36.0          Biobase_2.66.0                      
#> [13] GenomicRanges_1.58.0                 GenomeInfoDb_1.42.0                 
#> [15] IRanges_2.40.0                       S4Vectors_0.44.0                    
#> [17] BiocGenerics_0.52.0                  MatrixGenerics_1.18.0               
#> [19] matrixStats_1.4.1                    data.table_1.16.2                   
#> 
#> loaded via a namespace (and not attached):
#>  [1] DBI_1.2.3                 bitops_1.0-9             
#>  [3] permute_0.9-7             rlang_1.1.4              
#>  [5] magrittr_2.0.3            compiler_4.4.1           
#>  [7] RSQLite_2.3.7             DelayedMatrixStats_1.28.0
#>  [9] png_0.1-8                 vctrs_0.6.5              
#> [11] pkgconfig_2.0.3           crayon_1.5.3             
#> [13] fastmap_1.2.0             dbplyr_2.5.0             
#> [15] labeling_0.4.3            utf8_1.2.4               
#> [17] Rsamtools_2.22.0          rmarkdown_2.28           
#> [19] UCSC.utils_1.2.0          bit_4.5.0                
#> [21] xfun_0.48                 zlibbioc_1.52.0          
#> [23] cachem_1.1.0              jsonlite_1.8.9           
#> [25] blob_1.2.4                highr_0.11               
#> [27] rhdf5filters_1.18.0       DelayedArray_0.32.0      
#> [29] Rhdf5lib_1.28.0           BiocParallel_1.40.0      
#> [31] parallel_4.4.1            R6_2.5.1                 
#> [33] bslib_0.8.0               RColorBrewer_1.1-3       
#> [35] limma_3.62.0              jquerylib_0.1.4          
#> [37] Rcpp_1.0.13               knitr_1.48               
#> [39] R.utils_2.12.3            Matrix_1.7-1             
#> [41] tidyselect_1.2.1          abind_1.4-8              
#> [43] yaml_2.3.10               codetools_0.2-20         
#> [45] curl_5.2.3                lattice_0.22-6           
#> [47] tibble_3.2.1              withr_3.0.2              
#> [49] KEGGREST_1.46.0           evaluate_1.0.1           
#> [51] BiocFileCache_2.14.0      pillar_1.9.0             
#> [53] BiocManager_1.30.25       filelock_1.0.3           
#> [55] generics_0.1.3            RCurl_1.98-1.16          
#> [57] BiocVersion_3.20.0        ggplot2_3.5.1            
#> [59] sparseMatrixStats_1.18.0  munsell_0.5.1            
#> [61] scales_1.3.0              gtools_3.9.5             
#> [63] glue_1.8.0                tools_4.4.1              
#> [65] AnnotationHub_3.14.0      locfit_1.5-9.10          
#> [67] GenomicAlignments_1.42.0  XML_3.99-0.17            
#> [69] rhdf5_2.50.0              grid_4.4.1               
#> [71] AnnotationDbi_1.68.0      colorspace_2.1-1         
#> [73] GenomeInfoDbData_1.2.13   HDF5Array_1.34.0         
#> [75] restfulr_0.0.15           cli_3.6.3                
#> [77] rappdirs_0.3.3            fansi_1.0.6              
#> [79] S4Arrays_1.6.0            dplyr_1.1.4              
#> [81] gtable_0.3.6              R.methodsS3_1.8.2        
#> [83] sass_0.4.9                digest_0.6.37            
#> [85] SparseArray_1.6.0         rjson_0.2.23             
#> [87] farver_2.1.2              memoise_2.0.1            
#> [89] htmltools_0.5.8.1         R.oo_1.26.0              
#> [91] lifecycle_1.0.4           httr_1.4.7               
#> [93] statmod_1.5.0             bit64_4.5.2