Search
Question: Simplifying a RaggedExperiment using qreduceAssay and weighted means.
2
gravatar for Levi Waldron
6 days ago by
Levi Waldron600
CUNY Graduate School of Public Health and Health Policy, New York, NY
Levi Waldron600 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   
>
ADD COMMENTlink modified 5 days ago by Martin Morgan ♦♦ 22k • written 6 days ago by Levi Waldron600
1
gravatar for Martin Morgan
5 days 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

 

ADD COMMENTlink written 5 days ago by Martin Morgan ♦♦ 22k

Hey all, How does this commit look?

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

Thanks, Marcel

ADD REPLYlink modified 5 days ago • written 5 days ago by Marcel Ramos ♦♦ 210
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 295 users visited in the last hour