Methirx is failing on Bioconductor devel for a while now. This issue was a tricky one to track but should be fixed soon to avoid being kicked out from next release. Problem arises from changes BSgenome objects and the way they store metadata. methrix parses ref. build info from BSgenome object which is stored as an attribute.

Current implementation as on BC release and GitHub

#version: BSgenome.Hsapiens.UCSC.hg19_1.4.3
#Version: BSgenome_1.56.0
> ref_genome = BSgenome::getBSgenome(genome = 'hg19')

> names(attributes(x = ref_genome))
 [1] "pkgname"            "single_sequences"   "multiple_sequences"
 [4] "source_url"         "user_seqnames"      "injectSNPs_handler"
 [7] ".seqs_cache"        ".link_counts"       "nmask_per_seq"     
[10] "masks"              "organism"           "common_name"       
[13] "provider"           "provider_version"   "release_date"      
[16] "release_name"       "seqinfo"            "class"

>attributes(x = ref_genome)$provider_version
[1] "hg19"

In the newer versions (BC devel) this information is stored under metadata subslot (which results in NULL when one try to access)

#version: BSgenome.Hsapiens.UCSC.hg19_1.4.3
#Version: BSgenome_1.57.6

> names(attributes(x = ref_genome))
 [1] "pkgname"            "single_sequences"   "multiple_sequences"
 [4] "seqinfo"            "user_seqnames"      "injectSNPs_handler"
 [7] ".seqs_cache"        ".link_counts"       "metadata"          
[10] "class"

>attributes(x = ref_genome)$provider_version  #Results in NULL
NULL

> attributes(x = ref_genome)$metadata$genome #Workaround
[1] "hg19"

This line needs to be changed for the fix.

in documenatation: at STEP2

meth = methrix::read_bedgraphs(files = bdg_files[1:3], pipeline = "MethylcTools", collapse_starnds = TRUE, vect_batch_size = 3, ref_build = "Hs37d5", ref_cpgs = hs37d5_cpgs)

collapse_starnds --> collapse_strands

summary_data_tables<-list()
for(i in seq(granges_list)){

• summary_data_tables[[i]]<-get_region_summary(m = meth, regions = granges_list[[i]], type = "M", how = "mean")
  

• }
  -Checking for overlaps..
  -Summarizing by average
  Error in [.data.table(dat, , lapply(.SD, mean, na.rm = na_rm), by = yid, :
  Some items of .SDcols are not column names: [cells01-rep1, cells01-rep2, cells01-rep3, cells02-rep1, cells02-rep2, cells02-rep3, cells03-rep1, cells03-rep2, cells03-rep3, cells04-rep1, ...]

get_region_summary() will not work if sample names have underscore "_" as they are converted into "." in the usage of the function

methrix::plot_coverage(meth, type = "dens", size.lim = 1e+06)

warnings/errors:
The dataset is bigger than the size limit. A random subset of the object will be used that contains ~1e+06 observations.
Error in sample.int(length(x), size, replace, prob) :
cannot take a sample larger than the population when 'replace = FALSE'

further information:
meth object has 25679692 elements/CpGs,
plot funktion worked when size.lim was equal or larger than the methrix object and therefore did not need to select a random subset.

Here the code and the error message:

install_github( CompEpigen/methrix ,  ref = "devel_update")
library("methrix")

report.dir <- "/C010-projects/Pathway_27/SGBS_hydroxymethylation/WGoxBS/QC_all"

setwd(report.dir)
m <-readRDS("methrix_object_hs37d5.RDS")


#subset samples 
m_subset <-m[,m@ metadata$descriptive_stats$genome_stat$mean_cov>10]

#filter CpGs 
data_mask <- mask_methrix(m_subset , low_count = 2, high_quantile = 0.99)
m10masked <-coverage_filter(data_mask, cov_thr = 10, min_samples = 16)

#generate report and recalculate statistics
methrix_report(m10masked, output_dir=file.path(report.dir,"cov10" ), recal_stats = T)


**Error message:**
> _Step 1 of 5: Methylation/Coverage per chromosome  Error in merge.data.table(meth_stat, cov_stat, by = c("chr", "Sample_Name")) :    Elements listed in `by` must be valid column names in x and y In addition: Warning message: In merge.data.table(meth_stat, cov_stat, by = c("chr", "Sample_Name")) :   You are trying to join data.tables where 'x' and 'y' arguments are 0 columns data.table.
> --_

Please, fix the ordering in the chromosome-wise plots to make it natural (chr1,2,3,..,10,..,X,Y) instead of alpha-numeric (chr1,10,11,....,chr2,..).

Hey man !
Very cool plotly report !
I don't know if I am the only one to experience this, but the legend titles in some plot are cut out of the display windows, see below (upper right corner of the plots):

I am using Chromium Version 80.0.3987.87 (Build officiel) Built on Ubuntu , running on Ubuntu 16.04 (64 bits). It occurs on full screen mode half-splitted screen mode.
The issue is also affecting the exported PNG file.

This is no emergency.
If you wait enough maybe I will just fix it myself :) .
I guess adjusting plot's margin should solve this (check it first, and if setting new margins solves it you could just provide a fix on the right side margin setting a value on separated plots, or using nchar() on the given color variable name multiplied by a factor n that keeps the margin large enough to display the full legend title doesn't matter it's length).
Cheers !

Hi, I am attempting to use get_region_summary with a very large methrix object (almost 400 human samples) and a GRanges object with about 80,000 regions and am finding that it is constantly using all the RAM on my machine before being killed. Even on a workstation with almost 200 GB RAM, the problem persists. I guess by using the n_chunks parameter I can reduce the RAM being used, but I am wondering is there a rough guide for how many chunks to use depending on the number of samples and number of regions being evaluated?

Could not find BSgenome BSgenome.Hsapiens.UCSC.mm10
Found following BSgenome installtions. Provide the correct 'pkgname'

I've been having a great time using methrix over the past weeks. Thanks for creating such a nice package.

I've run into a bug when reading BedGraph files from Methyldackel using the pipeline = "MethyDackel preset and transforming the methrix file into a bsseq object using methrix2bsseq.

The resulting bsseq from methrix2bsseq will have non-integer M-Coverage values when retrieving them using:

meth <-  methrix::read_bedgraphs(n_threads = threads,
                                            files = bedgraph_fps,
                                            pipeline = "Methyldackel"
                                            ref_cpgs = mm10_cpgs,
                                            coldata = meta_data)
bsseq_obj <- methrix::methrix2bsseq(meth)                                
bsseq::getCoverage(bsseq_obj, type = M)

I've had a look at the code for methrix2bsseq and I think I've found the cause.

When using pipeline = "MethyDackel, read_bedgraph uses the Methyldackel calculated beta values, however these are rounded to the 2nd decimal place.
Example(4th column):

track type="bedGraph" description="SRR1182519.sorted CpG methylation levels"
1	25114	25115	66	2	1
1	25115	25116	100	3	0

methrix2bsseq calculates the M-Matrix of methylated reads by multiplying with the coverage, however using these rounded beta values will cause non-integer values to be inserted in the M-Coverage Matrix.

My workaround for now is setting the read_bedgraph fields manually since then methrix calculates exact beta values itself.
Then everything works out

methrix::read_bedgraphs(n_threads = threads,
                                           files = smoothed_bsseq_files,
                                           chr_idx = 1, start_idx = 2, end_idx = 3, M_idx =5, U_idx = 6,
                                           ref_cpgs = mm10_cpgs,
                                           coldata = meta_data,
                                            
  )

Maybe it would make sense having having methrix calculate the Beta values when reading from Methyldackel bedgraphs?

In the methrix_report in the image "Chromosome-wise depth of coverage, MEDIAN" it shows instead of median covarage, median chromosome-wise methylation.

meth <- methrix::remove_uncovered(m = meth)
Error in as.vector(x, "list") :
cannot coerce type 'closure' to vector of type 'list'
traceback)(
Error: unexpected ')' in "traceback)"
traceback()
No traceback available

To simplify the import of bedGraphs into UCSC Genome Browser it would be great to have the usual track lines (name, etc).

The loading report shows first a breakdown of CpGs by "raw", "filtered", "stranded". Downstream there is a "Processing" section, where also "CpGs missing" numbers are given after processing each bedGraph. It is unclear to which number from the first section (raw/filtered/stranded) the missing number should be related. Perhaps specifying this in parenthesis along with percentages could be useful.

Hello @tkik @MaxSchoenung @HeyLifeHD @lutsik ,

As a part of the package it would be nicer to have an example data for testing/demonstration. Bsseq has BS.chr22 which can be loaded with the command data(BS.chr22). In methrix I included a part of single cell dataset from bsmap protocol - data(mm9_bsmap). However @tkik suggested that this dataset is quite "boring" and we should use something more robust - and I agree.

Two options:

1. We convert the dataset bundled with bsseq to methrix object and use it.
2. We make our own custom dataset from well studied in-house data.

If you opt for 2, we should make sure that we have at-lease two samples while keeping the object size small. For example, bsseq data is quite small (3.8MB) with ~50K loci and 2 samples. Bioconductor has a package size limit of 5 MB and the core team is quite strict in this regard.

> data(BS.chr22)
> BS.chr22
An object of type 'BSseq' with
  494728 methylation loci
  2 samples
has not been smoothed
All assays are in-memory
> print(object.size(x = BS.chr22), units = "MB")
3.8 Mb

Suggestions welcome, and please post here if you have an example set in mind.

meth <- methrix::mask_methrix(m = meth, low_count=5, high_quantile = 0.99)
-Masked 5,498,862 CpGs due to low coverage.
Error in [<-(*tmp*, row_idx, value = NA) :
linear subassignment to a DelayedArray object 'x' (i.e. 'x[i] <-
value') is only supported when the subscript 'i' is a logical
DelayedArray object of the same dimensions as 'x' and 'value' an
ordinary vector of length 1)
traceback()
4: stop(wmsg(.linear_subassign_error_msg))
3: [<-(*tmp*, row_idx, value = NA)
2: [<-(*tmp*, row_idx, value = NA)
1: methrix::mask_methrix(m = meth, low_count = 5, high_quantile = 0.99)

Many bedGraph formats are not supported by IGV, would be great to have a "plug'n'play" export from methrix objects.

By default one gets a plot with both axes labeled as PC1

Some cosmetics: most graphs in the html report have no x and y axis labels.

After loading a methrix object from a RDS file ( readRDS() ), when I tried to remove SNPs this error message shown up (partially translated from French):

Error in h(simpleError(msg, call)) : 
  evaluation error of the argument 'x' when selecting a method for 'unique' function: evaluation error of the argument 'x' when selecting a method for 'unique' function: evaluation error of the argument 'x' when selecting a method for 'which' function : Sequence names M in 'ranges' not present in reference genome hs37d5.

For some reason the error for unique() is printed 2 times. The 3 errors are concatenated (without carriage returns). The last part of the error message sounds like an attempt to query M chromosome on the annotation MafDb.1Kgenomes.phase3.hs37d5.

Here is the complete code:

#Load library
library(methrix)
library(BSgenome.Hsapiens.UCSC.hg19)

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

#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
)

#First extract genome wide CpGs from the desired reference genome
hg19_cpgs <- methrix::extract_CPGs(ref_genome = "BSgenome.Hsapiens.UCSC.hg19")

meth <- readRDS(file = "methrix_object.RDS")

#Generate report
methrix::methrix_report(meth = meth, output_dir = "~/methrix_tuto/")

# #Remove uncovered locis
# meth = methrix::remove_uncovered(m = meth)

if(!require(MafDb.1Kgenomes.phase3.hs37d5)) {
  BiocManager::install("MafDb.1Kgenomes.phase3.hs37d5")}
if(!require(GenomicScores)) {
  BiocManager::install("GenomicScores")}

library(MafDb.1Kgenomes.phase3.hs37d5)
library(GenomicScores)

#Remove SNPs
meth_snps_filtered <- methrix::remove_snps(m = meth)

The error only happens for the last command.
I cannot reproduce the error when the methrix object is generated with read_bedgraphs(), instead of being loaded from a RDS file.

Let me know if you need anything else for the issue.
Best,

Thanks Anand, that is a great start!

I still think that the import of bed files should be a bit more structured. What I mean is that, no matter what exact format the bed has, the information which they provide is roughly the same (chromosome, start, end, strand, number of methylated, number of unmethylated, mean methylation, total coverage, anything else?). We could therefore first could try to map each input column to the set of predefined columns, and then write a common read-in procedure that does not depend on the format anymore. Back then I tried to do it this way in RnBeads, but not sure about my code anymore. Please do look/steal/ stuff from there (dataImport.R).

Hi,

I used HDF5 backend due to memory issues. I firstly tried to save a RDS file of the output from methrix::read_bedgraphs, but it gives me error that name = '/private/var/folders/0r/35wl9njn767_c2md3zy0h61w0000gn/T/Rtmpp14UbV/M_sink_1.h5', errno = 2, error message = 'No such file or directory'. I checked that there is the file named M_sink_1.h5 in the directory so I don't know what happened here. Also, I tried to remove the h5 file and rerun the code, it gives error that Error in DelayedArray::write_block(block = as.matrix(b$bdg[, .(beta)]), :
unused argument (sink = M_sink). Could you help me? :)

Thanks

methrix_report(m10masked, "/C010-projects/Pathway_27/SGBS_hydroxymethylation/WGoxBS/QC/methrix_reports_filtered2", recal_stats=T)

error message only occures if recal_stats=T:

Step 1 of 5: Methylation/Coverage per chromosome
Error in merge.data.table(meth_stat, cov_stat, by = c("chr", "Sample_Name")) :
Elements listed in by must be valid column names in x and y
In addition: Warning message:
In merge.data.table(meth_stat, cov_stat, by = c("chr", "Sample_Name")) :
You are trying to join data.tables where 'x' and 'y' arguments are 0 columns data.table.