Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
example: GenomicFiles MAP on GRangesList
library(GenomicFiles)
register(SerialParam())
# use 4 BigWigs from
# ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeCshlLongRnaSeq/
files <- list.files("~/data/bigwig/",full=TRUE)
ranges0 <- GRanges("chr1",IRanges(11:20 * 1e6, width=1e5))
f <- factor(rep(1:2, each=length(ranges0)/2))
ranges <- split(ranges0, f)
MAP <- function(range, file) {
rle <- import(file, which = range, as = "Rle", format = "bw")[range]
mean(rle)
}
REDUCE <- function(mapped) {
Reduce(cbind, mapped)
}
(res <- reduceByRange(ranges, files, MAP, REDUCE))
(resOut <- do.call(rbind, res))
sessionInfo()
> register(SerialParam())
> # use 4 BigWigs from
> # ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeCshlLongRnaSeq/
> files <- list.files("~/data/bigwig/",full=TRUE)
> ranges0 <- GRanges("chr1",IRanges(11:20 * 1e6, width=1e5))
> f <- factor(rep(1:2, each=length(ranges0)/2))
> ranges <- split(ranges0, f)
> MAP <- function(range, file) {
+ rle <- import(file, which = range, as = "Rle", format = "bw")[range]
+ mean(rle)
+ }
> REDUCE <- function(mapped) {
+ Reduce(cbind, mapped)
+ }
> (res <- reduceByRange(ranges, files, MAP, REDUCE))
$`1`
init
chr1 0.08394 0.08899 4.58022 3.49198
chr1 0.04923 0.04763 10.90287 8.61417
chr1 0.03631 0.04303 0.00780 0.00727
chr1 0.03408 0.04532 2.16463 1.69084
chr1 0.01378 0.00825 0.01869 0.00581
$`2`
init
chr1 0.13007 0.07382 3.83418 3.69978
chr1 0.91577 0.59110 50.07868 71.86047
chr1 0.01693 0.04954 0.18872 0.17385
chr1 0.00671 0.00805 0.20328 0.17623
chr1 1.12673 0.94558 0.05567 0.04120
> (resOut <- do.call(rbind, res))
init
chr1 0.08394 0.08899 4.58022 3.49198
chr1 0.04923 0.04763 10.90287 8.61417
chr1 0.03631 0.04303 0.00780 0.00727
chr1 0.03408 0.04532 2.16463 1.69084
chr1 0.01378 0.00825 0.01869 0.00581
chr1 0.13007 0.07382 3.83418 3.69978
chr1 0.91577 0.59110 50.07868 71.86047
chr1 0.01693 0.04954 0.18872 0.17385
chr1 0.00671 0.00805 0.20328 0.17623
chr1 1.12673 0.94558 0.05567 0.04120
> sessionInfo()
R Under development (unstable) (2014-09-25 r66681)
Platform: x86_64-apple-darwin12.5.0 (64-bit)
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] stats4 parallel stats graphics grDevices datasets utils
[8] methods base
other attached packages:
[1] GenomicFiles_1.1.19 BiocParallel_0.99.22 futile.logger_1.3.7
[4] rtracklayer_1.25.16 Rsamtools_1.17.34 Biostrings_2.33.14
[7] XVector_0.5.8 GenomicRanges_1.17.42 GenomeInfoDb_1.1.23
[10] IRanges_1.99.28 S4Vectors_0.2.4 BiocGenerics_0.11.5
[13] RUnit_0.4.27 devtools_1.6 slidify_0.4.5
[16] knitr_1.6 BiocInstaller_1.15.5
loaded via a namespace (and not attached):
[1] base64enc_0.1-2 BatchJobs_1.4 BBmisc_1.7
[4] bitops_1.0-6 brew_1.0-6 checkmate_1.4
[7] codetools_0.2-9 DBI_0.3.1 digest_0.6.4
[10] evaluate_0.5.5 fail_1.2 foreach_1.4.2
[13] formatR_1.0 futile.options_1.0.0 GenomicAlignments_1.1.30
[16] iterators_1.0.7 lambda.r_1.1.6 markdown_0.7.4
[19] RCurl_1.95-4.3 RSQLite_0.11.4 sendmailR_1.2-1
[22] stringr_0.6.2 tools_3.2.0 whisker_0.3-2
[25] XML_3.98-1.1 yaml_2.1.13 zlibbioc_1.11.1
>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment