genomicFilesFromBeds: create a GenomicFiles from the results of sesamizeToBeds (or...

In trichelab/h5testR: Tests for HDF5-backed DNAme, RNAseq, and MultiAssayExperiment objects

Description

see sesamizeToBeds for an easy way to go from IDATs to tabixed BED files.

Usage

genomicFilesFromBeds(beds, rr = NULL, cd = NULL, verbose = TRUE)

Arguments

beds	a character vector of files or data.frame with the column $bed
rr	rowRanges for GenomicFiles object (default: union BED ranges)
cd	colData for GenomicFiles object (default: skeletal info)
verbose	squawk? (TRUE)

Details

TODO: pull down the seqinfo for a genome assembly and add it to rowRanges

Value

a GenomicFiles object

Examples

# median methylation across a range
mm <- function(range, file, ...) { 
  res <- rtracklayer::import(file, which=range, ...)$score
  if (!is.null(res)) res <- median(res, na.rm=TRUE)
  if (is.null(res)) is.na(res) <- TRUE 
  return(res)
}

# EPIC data, hg38 mappings, parallel
if (require("minfiDataEPIC")) {

  baseDir <- system.file("extdata", package = "minfiDataEPIC")
  sheet <- minfi::read.metharray.sheet(baseDir) 

  # don't panic, feedback from sesamizeToBeds will appear eventually 
  beds <- sesamizeToBeds(sheet, refversion="hg38", BPPARAM=MulticoreParam(3))

  # package up the resulting beds 
  gf <- genomicFilesFromBeds(beds)
  show(gf)

  # get one probe's worth of data
  cg21870274 <- GRanges("chr1", IRanges(69591), "*", name="cg21870274")
  names(cg21870274) <- "cg21870274"
  genome(cg21870274) <- "hg38"
  unlist(reduceByRange(ranges=cg21870274, files=files(gf), MAP=mm)[[1]])

  # summarize across some CpG islands
  data("HMM_CpG_islands.hg38")
  set.seed(1234)
  someIndices <- sample(seq_along(HMM_CpG_islands.hg38), 20)
  cpgi <- HMM_CpG_islands.hg38[someIndices]
  se <- reduceByRange(ranges=cpgi, files=files(gf), MAP=mm, summarize=TRUE)
  show(se)
  head(assays(se)$data)

  # information extracted from each BED
  colData(gf)

}

trichelab/h5testR documentation built on July 12, 2020, 5:18 p.m.