Simplifying a RaggedExperiment using qreduceAssay and weighted means.
1
2
Entering edit mode
@levi-waldron-3429
Last seen 6 weeks ago
CUNY Graduate School of Public Health a…

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   
>
raggedexperiment genomicranges • 704 views
ADD COMMENT
1
Entering edit mode
@martin-morgan-1513
Last seen 3 days ago
United States

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 COMMENT
0
Entering edit mode
ADD REPLY

Login before adding your answer.

Traffic: 210 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6