#' @title Assert chromosome exists
#'
#' @description Checks whether a chromosome exists in the target bigwig. Throws
#' error if it does not.
#' @param bigwig_file a string pointing to a bigwig file
#' @param chromosome chromosome name
#' @rdname assert_chromosome_exists
#' @return silent if true, else error
#' @export
assert_chromosome_exists <- function(chromosome, bigwig_file) {
bw_seqnames <- GenomicRanges::seqnames(
rtracklayer::seqinfo(rtracklayer::BigWigFile(bigwig_file))
)
pass <- all(chromosome %in% bw_seqnames)
if (!pass) {
not_incl <- setdiff(chromosome, bw_seqnames)
stop(paste("chromosome",
paste(not_incl, collapse = ", "),
"does not exist in", bigwig_file))
}
}
#' Class profile_views
#'
#' Class \code{profile_views} holds views of ranges on a dataset.
#'
#' @slot query_ranges a \code{\link[GenomicRanges]{GRanges-class}} that holds
#' all the query ranges
#' @slot profiles a \code{\link[IRanges]{ViewsList-class}} object that holds the data
#' views
#' @slot idx a lookup index used by acessors for the class
#'
#' @name profile_views-class
#' @rdname profile_views-class
#' @importClassesFrom GenomicRanges GRanges
#' @importClassesFrom IRanges ViewsList
#' @exportClass profile_views
methods::setClass("profile_views",
slots = c(query_ranges = "GRanges",
profiles = "ViewsList",
idx = "matrix")
)
#' @title extract regions from file
#'
#' @description Extracts regions from a bigwig over a set of ranges from a
#' \code{\link[GenomicRanges]{GRanges-class}} object and build a
#' \code{\link{profile_views}} object
#'
#' @inheritParams assert_chromosome_exists
#' @param query_ranges A \code{\link[GenomicRanges]{GRanges-class}} object
#'
#' @importClassesFrom GenomicRanges GRanges
#'
#' @return A list of run length encodings
#'
#' @name profile_views
#' @rdname profile_views
#'
#' @export
profile_views <- function(bigwig_file, query_ranges) {
# Restrict to active seqnames
GenomeInfoDb::seqlevels(query_ranges) <- GenomeInfoDb::seqlevelsInUse(query_ranges)
# Get unique chromosome names in bins
chrom <-
unique(S4Vectors::runValue(GenomicRanges::seqnames(query_ranges)))
# Check that all query chromosomes are present in bigwig
assert_chromosome_exists(chromosome = chrom, bigwig_file)
# Split queries into list
names(query_ranges) <- seq_along(query_ranges)
query_list <- split(IRanges::IRanges(GenomicRanges::start(query_ranges),
GenomicRanges::end(query_ranges),
names = names(query_ranges)),
f = GenomicRanges::seqnames(query_ranges))
# Import bigwig, only chromosomes to be queried
imported_bw <- rtracklayer::import.bw(
con = bigwig_file, which = query_ranges,
as = "RleList")[names(query_list)]
# Iterate over chromosomes
tx_views <- IRanges::Views(
subject = imported_bw,
query_list)
# Extract the originating row and final location in views object
origin <- as.integer(unlist(lapply(query_list, names)))
list_idx <- rep(seq_along(query_list), unlist(lapply(query_list, length)))
row_idx <- unlist(lapply(query_list, seq_along))
lookup_idx <- cbind(origin = origin, list_idx = list_idx, row_idx = row_idx)
rownames(lookup_idx) <- NULL
lookup_idx <- lookup_idx[order(lookup_idx[, 1]), ]
return(methods::new(Class = "profile_views",
query_ranges = query_ranges,
profiles = tx_views,
idx = lookup_idx
))
}
#' @title subsetting
#'
#' @description Convenience subsetting for views object. Retrieves view based on index of
#' query range.
#'
#' @param x \link{profile_views-class} object
#' @param i subset entries of views object
#' @param j unsupported
#' @param drop unsupported
#' @param ... just for compatibility with S4 method
#'
#' @aliases [
#' @return A views object
#' @export
methods::setMethod("[", c(x = "profile_views", i = "numeric", j = "missing", drop = "ANY"),
function(x, i, j, ..., drop=TRUE){
i <- as.integer(i)
return(x@profiles[[x@idx[i, 2]]][x@idx[i, 3]])
}
)
#' @title show
#'
#' @description Display object
#'
#' @inheritParams methods::show
#' @importMethodsFrom methods show
#'
#' @export
methods::setMethod("show", c(object = "profile_views"),
function(object) {
print(paste("profile_views object with", length(object@query_ranges),
"views"))
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.