[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",
representation(m="integer"),
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="")
callNextMethod()
})
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",
prototype=prototype(elementType="Tuple"))
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.
Martin
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