.expand <- function(granges, size){
if(any(is.na(seqlengths(granges)))) stop("seqlengths must be specified")
if(identical(length(size), 1L)){
size <- rep(size, length(granges))
}
starts <- pmax(start(granges)-size, 1L)
ends <- pmin(end(granges)+size, seqlengths(granges)[chromosome(granges)])
d1 <- start(granges)-starts
d2 <- ends-end(granges)
i <- which(d1 > 0 & d2 > 0 & d1 != d2)
if(length(i) > 0){
dd <- pmin(d1[i], d2[i])
starts[i] <- start(granges)[i]-dd
ends[i] <- end(granges)[i]+dd
}
start(granges) <- starts
end(granges) <- ends
granges
}
#' Expand a genomic interval
#'
#' Expand a genomic interval by a user-specified number of basepairs.
#'
#' @examples
#' library(GenomeInfoDb)
#' library(BSgenome.Hsapiens.UCSC.hg19)
#' library(GenomicRanges)
#' si <- keepSeqlevels(seqinfo(BSgenome.Hsapiens.UCSC.hg19), "chr1")
#' gr <- GRanges("chr1", IRanges(15e3, 20e3), seqinfo=si)
#' gr2 <- expandGRanges(gr, 5e3)
#' ## will not go beyond the seqlengths
#' gr3 <- GRanges("chr1", IRanges(seqlengths(gr)-5000, seqlengths(gr)-2000),
#' seqinfo=si)
#' gr4 <- expandGRanges(gr3, 10e3)
#' end(gr4) == seqlengths(gr)
#'
#'
#' @export
#' @param granges a \code{GRanges} object
#'
#' @param size a length-one numeric vector specifying the number of
#' basepairs to subtract from the start and add to the ends of the
#' intervals in object \code{granges}
expandGRanges <- function(granges, size) .expand(granges, size)
#' Find the fraction of genomic interval overlapping with another
#' genomic interval.
#'
#'
#' Compute the fraction of a query genomic interval that overlaps with another
#' set of genomic intervals. We divide the total width of the intersection with
#' the width of the query.
#'
#' @examples
#' library(GenomicRanges)
#' autosomes <- function() paste0("chr", 1:22)
#' del <- GRanges(c("chr1", "chr2"), IRanges(c(1, 1), c(10, 10)))
#' filters <- GRanges(c("chr1", "chr2"), IRanges(c(5, 8), c(20, 10000)))
#' test <- intOverWidth(del, filters)
#' identical(test, c(0.6, 0.3))
#'
#' @return A numeric vector of the same length as the query
#' \code{GRanges}
#'
#' @export
#' @param query a \code{GRanges} object
#' @param subject a \code{GRanges} object
intOverWidth <- function(query, subject){
subject <- reduce(subject)
p <- rep(0, length(query))
hits <- findOverlaps(query, subject)
if(length(hits) == 0) return(p)
## vectorized intersect
i <- queryHits(hits)
j <- subjectHits(hits)
intersection <- width(pintersect(query[i], subject[j]))
numer <- sapply(split(intersection, i), sum)
denom <- width(query)[as.integer(names(numer))]
p[as.integer(names(numer))] <- numer/denom
p
}
#' Uncouple linked bins
#'
#' Two genomic intervals that are linked are represented by a
#' \code{GRanges} object where a single elemnt provides the genomic
#' interval of a bin and a genomic interval of a bin that it is linked
#' to. The latter is included in the metadata (mcols) variable
#' 'linked.to'. This function represents every genomic interval as a
#' single element. Hence, a \code{GRanges} of linked intervals with
#' length one would have length-two after uncoupling the bins.
#'
#' @examples
#' linked <- GRanges("chr1", IRanges(1,5))
#' linked$linked.to <- setNames(GRanges("chr1", IRanges(10, 20)), "b")
#' names(linked) <- "a"
#' uncouple(linked)
#'
#' @return a \code{GRanges} object
#' @export
#'
#' @param linked_bins a \code{GRanges} representation of two genomic
#' intervals that are linked.
uncouple <- function(linked_bins){
gr <- c(granges(linked_bins), granges(linked_bins$linked.to))
## order by rearrangement
L <- seq(1, length(linked_bins)*2, 2)
R <- seq(2, length(linked_bins)*2, 2)
gr <- gr[c(L, R)]
if(!is.null(names(linked_bins))){
names(gr)[L] <- names(linked_bins)
}
if(!is.null(names(linked_bins$linked.to))){
names(gr)[R] <- names(linked_bins$linked.to)
}
gr <- gr[!duplicated(gr)]
gr
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.