[Bioc-devel] Help in designing class based on SummarizedExperiment

Martin Morgan mtmorgan at fhcrc.org
Sat Feb 15 02:55:24 CET 2014

Hi Peter --

I'm not really sure that I can provide good advice here; it feels a bit like a 
'round peg square hole' fit with SummarizedExperiment, since your 'pos' is not 
characterized by genomic ranges, and columns of 'assay' does not correspond to 
separate samples.

One possibility is to create a Tuple class

     .Tuple <- setClass("Tuple",
         contains="SimpleList", )

     Tuple <- function(seqnames, pos, count)
         ## note constructor: base class(es) as first and unnamed arg
         .Tuple(SimpleList(seqnames=seqnames, pos=pos, count=count), m=ncol(pos))

     setMethod(show, "Tuple", function(object) {
         cat(length(object$seqnames), " x ", object at m, "-", sep="")

     tuple <- Tuple(test_data$seqnames,
                    matrix(unlist(test_data[2:4]), ncol=3),
                    matrix(unlist(test_data[5:12]), ncol=8))

You could implement a validity method to enforce correct list membership and 
element dimensions, and a subset method to allow for row-wise subsetting across 
elements. You could store the seqnames as an Rle, and the matrix as a sparse 
matrix (Matrix::Matrix(..., sparse=TRUE)) or as a DataFrame with columns as 
Rle's (but many of your operations will presumably be matrix-like, so you'll end 
up coercing to a (dense) matrix anyway). I guess your overall data is a 
TupleList, something like

     .TupleList <- setClass("TupleList", contains="CompressedList",

The pass-by-reference approach used in SummarizedExperiment is not exported from 
GenomicRanges, so it's not easy to get the pass-by-reference / copy-on-change 
semantics; worth noting that the copying cost is incurred when the data is 
_modified_, not just accessed, and it's often the case that the data is not 
updated after creation. A not-obvious paradigm when updating several slots is to 
call the initialize method with the original object, maybe like

     setMethod("[", c("Tuple", "ANY", "missing"),
         function(x, i, j, ..., drop=FALSE)
         if (!missing(drop))
             warning("'drop' ignored when subsetting ", sQuote(class(x)))
         initialize(x, SimpleList(seqnames=x$seqnames[i],
             pos=x$pos[i,,drop=FALSE], count=x$count[i,,drop=FALSE]))

This is more memory-efficient than updating slots one at a time, and presumes 
that you have avoided writing an initialize method that breaks this behavior.

The assays slot of SummarizedExperiment can be anything that has matrix-like 
semantics, so for instance you could save a sparse matrix

     SummarizedExperiment(SimpleList(counts=Matrix(m, sparse=TRUE)))

There is not a matrix representation based on Rle's, though (maybe there's a 
sparse representation in the Matrix package that does a good job of representing 
your data anyway?).

Hope that helps, and is not too misleading. Perhaps others on the list will have 
additional ideas.


On 02/09/2014 06:20 PM, Peter Hickey wrote:
> Hi all,
> Apologies up front for the rather long post.
> I'm designing a class to store what I call co-methylation m-tuples. These are based on a very simple tab-delimited file format.
> For example, here are 1-tuples (m = 1):
> chr		pos1	M	U
> chr1	57691	0	1
> chr1	59276	1	0
> chr1	60408	1	0
> chr1	63495	1	0
> chr1	63568	2	0
> chr1	63627	3	0
> 2-tuples (m = 2):
> chr		pos1	pos2	MM	MU	UM	UU
> chr1	567438	567570	0	0	0	2
> chr1	567501	567549	0	0	0	35
> chr1	567549	567558	0	1	0	139
> 3-tuples (m = 3):
> chr		pos1	pos2	pos3	MMM	MMU	MUM	MUU	UMM	UMU	UUM	UUU
> chr1	13644	13823	13828	1		0		0		0		0		0		0		0
> chr1	14741	14747 	14773	1		0		0		0		0		0		0		0
> etc.
> 1-tuples are basically the standard input to an analysis of BS-seq data.
> I think of these files as being comprised of 3 parts: the 'chr' column (chr), the 'pos' matrix (pos1, pos2, pos3) and the 'counts' matrix (MMM, MMU, MUM, MUU, UMM, UMU, UUM, UUU), when m = 3. For a given value of 'm' there is one 'chr' column, m 'pos' columns and 2^m 'counts' columns.
> I want to implement a class for these objects as I'm writing a package for the analysis of this type of data. I'd like a GRanges-type object storing the genomic information and a matrix-like object storing the counts. After tinkering around for a while, and doing some reading of the code in packages such as GenomicRanges and bsseq, I decided to extend the SummarizedExperiment class.  I now have a prototype but I have some questions and would appreciate feedback on some of my design choices before I translate my existing functions to work with this class of object.
> Here is the code for the prototype:
> #####################################################################
> library(GenomicRanges)
> setClass("CoMeth", contains = "SummarizedExperiment")
> CoMeth <- function(seqnames, pos, counts, m, methylation_type, sample_name, strand = "*", seqlengths = NULL, seqinfo = NULL){
>    # Argument checks, etc. go here #
>    gr <- GRanges(seqnames = seqnames, ranges = IRanges(start = pos[[1]], end = pos[[length(pos)]]), strand = strand, seqlengths = seqlengths, seqinfo = seqinfo) # The width of each element is defined by the first and last 'pos', e.g. for 3-tuples it is defined by pos1 and pos3.
>    # Need to store the "extra" positions if m > 2. Each additional position is stored as a separate assay
>    if (m > 2){
>      extra_pos <- lapply(seq(2, m - 1, 1), function(i, pos){
>        pos[[i]]
>      }, pos = pos)
>      names(extra_pos) <- names(pos)[2:(m-1)]
>    } else {
>      extra_pos <- NULL
>    }
>    assays <- SimpleList(c(counts, extra_pos))
>    colData <- DataFrame(sample_name = sample_name, m = m, methylation_type = paste0(sort(methylation_type), collapse = '/'))
>    cometh <- SummarizedExperiment(assays = assays, rowData = gr, colData = colData)
>    cometh <- as(cometh, "CoMeth")
>    return(cometh)
> }
> And here's some example data:
> # A function that roughly imitates the output of a call to scan() to read in BS-seq m-tuple data
> # m is the size of the m-tuples
> # n is the number of m-tuples
> # z is the proportion of each column of 'counts' that is zero
> make_test_data <- function(m, n, z){
>    seqnames <- list(seqnames = rep('chr1', n))
>    pos <- lapply(1:m, function(x, n){matrix(seq(from = 1 + x - 1, to = n + x - 1, by = 1), ncol = 1)}, n = n) # Need these to be matrices rather than vectors
>    names(pos) <- paste0('pos', 1:m)
>    # A rough hack to simulate counts where a proportion (z) are 0 and the rest are sampled from Poisson(lambda). Small values of lambda will inflate the zero-count.
>    counts <- mapply(FUN = function(i, z, n, lambda){
>      nz <- floor(n * (1 - z))
>      matrix(sample(c(rpois(nz, lambda), rep(0, n - nz))), ncol = 1)
>      }, i = 1:(2 ^ m), z = z, n = n, lambda = 4, SIMPLIFY = FALSE) # Need these to be matrices rather than vectors
>    names(counts) <- sort(do.call(paste0, expand.grid(lapply(seq_len(m), function(x){c('M', 'U')}))))
>    return(c(seqnames, pos, counts))
> }
> m <- 3 # An example using 3-tuples
> n <- 1000 # A typical value for 3-tuples from a methylC-seq experiment is n = 17,000,000
> z <- c(0.2, 0.6, 0.6, 0.7, 0.6, 0.8, 0.8, 0.7) # Typical proportions of each column of 'counts' that are zero when using 3-tuples for a methylC-seq experiment
> test_data <- make_test_data(n = n, m = m, z = z)
> cometh <- CoMeth(seqnames = test_data[['seqnames']], pos = test_data[grepl('pos', names(test_data))], counts = test_data[grepl('[MU]', names(test_data))], m = m, methylation_type = 'CG', sample_name = 'test_data')
> sessionInfo()
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-apple-darwin10.8.0 (64-bit)
> locale:
> [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets  methods   base
> other attached packages:
> [1] GenomicRanges_1.14.4 XVector_0.2.0        IRanges_1.20.6       BiocGenerics_0.8.0
> loaded via a namespace (and not attached):
> [1] stats4_3.0.2 tools_3.0.2
> #####################################################################
> Questions
> 	1. How can I move the 'extra_pos' columns from the assay slot but keep the copy-on-change behaviour? From a design perspective, I think it would make more sense for the 'extra_pos' columns, i.e. ('pos2') for 3-tuples and ('pos2', 'pos3') for 4-tuples etc., to be in their own slot rather than in the assays slot, after all, they aren't assays but rather are additional genomic co-ordinates. The 'extra_pos' fields are fixed (at least until I start subsetting or combining multiple CoMeth objects). My understanding of the the SummarizedExperiment class is that the assays slot is a reference class to avoid excessive copying when changing other slots of a SummarizedExperiment object. So if the 'extra_pos' columns were stored outside of the assays slot then these would have to be copied when any changes are made to the other slots of a CoMeth object, correct? Is there a way to avoid this, i.e. so that these 'extra_pos' columns are stored separately from the assays slot but with the
>   copy-on-change behaviour of the assays slot?
> 	2. Is the correct to compute something based on the 'counts' data via the assay() accessor? For example, I might want a helper function getCounts(cometh) that does the equivalent of sapply(X = 1:(2^m), function(i, cometh){assay(cometh, i)}, cometh = cometh). Similarly, I might want to compute the coverage of an m-tuple, which would be the equivalent of rowSums(getCounts(cometh)). Is this the correct way to do this sort of thing?
> 	3. How do I measure the size of a SummarizedExperiment/CoMeth object? For example, with the test data, print(object.size(cometh), units = "auto") < print(object.size(assays(cometh)), units = "auto"), so it seems that the size of the assays slot isn't counted by object.size().
> 	4. Is it possible to store an Rle-type object in the assays slot of a SummarizedExperiment? 20-80% of the entries in each column of 'counts' are zero and there are often runs of zeros. So I thought that perhaps an Rle representation (column-wise) might be more (memory) efficient. But I can't seem to get an Rle object in the assays slot (I tried via DataFrame); is it even possible?
> 	5. Are there matrix-like objects with Rle columns? I found this thread started by Kasper Hansen (https://stat.ethz.ch/pipermail/bioconductor/2012-June/046473.html) discussing the idea of matrix-like object where the columns are Rle's; I could imagine using such an object for a CoMeth object containing multiple samples, i.e. MMM is a matrix-like object with ncol = # of samples, MMU is matrix-like object with ncol = # of samples, etc. Was anything like this ever implemented? My reading of the previous thread was to use a DataFrame but the "matrix API", e.g. rowSums, doesn't work with DataFrames (and see (4) as to whether it's even possible to store such objects in the assays slot).
> Many thanks for your help in answering these questions. Any other suggestions on the design of the CoMeth class are appreciated.
> Thanks,
> Pete
> --------------------------------
> Peter Hickey,
> PhD Student/Research Assistant,
> Bioinformatics Division,
> Walter and Eliza Hall Institute of Medical Research,
> 1G Royal Parade, Parkville, Vic 3052, Australia.
> Ph: +613 9345 2324
> hickey at wehi.edu.au
> http://www.wehi.edu.au
> ______________________________________________________________________
> The information in this email is confidential and inte...{{dropped:16}}

More information about the Bioc-devel mailing list