Search
Question: Simplifying a RaggedExperiment using qreduceAssay and weighted means.
2
10 weeks ago by
Levi Waldron690
CUNY Graduate School of Public Health and Health Policy, New York, NY
Levi Waldron690 wrote:

The qreduceAssay() example from the RaggedExperiment vignette demonstrates reducing multiple ranges using a weightedmean function, but there are not actually any ranges requiring weighted means:

library(GenomicRanges)
library(RaggedExperiment)
sample1 <- GRanges(
c(A = "chr1:1-10:-", B = "chr1:8-14:+", C = "chr2:15-18:+"),
score = 3:5)
sample2 <- GRanges(
c(D = "chr1:1-10:-", E = "chr2:11-18:+"),
score = 1:2)
colDat <- DataFrame(id = 1:2)
ragexp <- RaggedExperiment(sample1 = sample1,
sample2 = sample2,
colData = colDat)
query <- GRanges(c("chr1:1-14:-", "chr2:11-18:+"))

weightedmean <- function(scores, ranges, qranges)
{
isects <- pintersect(ranges, qranges)
sum(scores * width(isects)) / sum(width(isects))
}

(res <- qreduceAssay(ragexp, query, simplifyReduce = weightedmean))

##              sample1 sample2
## chr1:1-14:-        3       1
## chr2:11-18:+       5       2


Here is a simpler RaggedExperiment that contains two ranges falling within one of the query ranges:

ragexp2 <- RaggedExperiment(sample1=GRanges(
c(A = "chr1:1-10:-", B = "chr1:11-14:-"),
score = 3:4))


This throws an error:

try(qreduceAssay(ragexp2, query, simplifyReduce = weightedmean))

Error in .pintersect_GRangesList_GRanges(x, y, drop.nohit.ranges = drop.nohit.ranges,  :
'y' must have the length of 'x' or length 1


But what I was hoping to have returned in this example was:

matrix(c((3 * 10 + 4 * 4) / 14, 0), ncol=1, dimnames=dimnames(res[, 1, drop=FALSE]))

##               sample1
## chr1:1-14:-  3.285714
## chr2:11-18:+ 0.000000


Here is a gist containing this code.

> sessionInfo()

R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.5

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] RaggedExperiment_1.5.2 GenomicRanges_1.33.6   GenomeInfoDb_1.17.1    IRanges_2.15.14
[5] S4Vectors_0.19.17      BiocGenerics_0.27.1

loaded via a namespace (and not attached):
[1] lattice_0.20-35             matrixStats_0.53.1          bitops_1.0-6
[4] grid_3.5.0                  zlibbioc_1.27.0             XVector_0.21.3
[7] Matrix_1.2-14               BiocParallel_1.15.7         tools_3.5.0
[10] Biobase_2.41.1              RCurl_1.95-4.10             DelayedArray_0.7.19
[13] yaml_2.1.19                 compiler_3.5.0              SummarizedExperiment_1.11.5
[16] GenomeInfoDbData_1.1.0
>

modified 10 weeks ago by Martin Morgan ♦♦ 22k • written 10 weeks ago by Levi Waldron690
1
10 weeks ago by
Martin Morgan ♦♦ 22k
United States
Martin Morgan ♦♦ 22k wrote:

The weightedmean function on the help page is

weightedmean <- function(scores, ranges, qranges)
## weighted average score per query range
sum(scores * width(ranges)) / sum(width(ranges))

(should the denominator be width(qranges)?) and the results seem to be correct for both your examples, adjusted so that ?

> (res <- qreduceAssay(ragexp, query, simplifyReduce = weightedmean))
sample1 sample2
chr1:1-14:-        3       1
chr2:11-18:+       5       2
> try(qreduceAssay(ragexp2, query, simplifyReduce = weightedmean))
sample1
chr1:1-14:-  3.285714
chr2:11-18:+       NA

and for a modified version of the help page example...

> sample1 <- GRanges(
+     c(A = "chr1:1-10:-", B = "chr1:8-14:-", C = "chr2:15-18:+"),
+     score = 3:5)
> sample2 <- GRanges(
+     c(D = "chr1:1-10:-", E = "chr2:11-18:+"),
+     score = 1:2)
> query <- GRanges(c("chr1:1-14:-", "chr2:11-18:+"))
> colDat <- DataFrame(id = 1:2)
> ragexp3 <- RaggedExperiment(sample1 = sample1,
+                            sample2 = sample2,
+                            colData = colDat)
> qreduceAssay(ragexp3, query, simplifyReduce = weightedmean)
sample1 sample2
chr1:1-14:-  3.411765       1
chr2:11-18:+ 5.000000       2

Hey all, How does this commit look?

https://github.com/Bioconductor/RaggedExperiment/commit/9fdcfbf0779d3508af7ae0afc288f946758b60d6

Thanks, Marcel