Get genome

DESCRIPTION

@@ -20,7 +20,7 @@ Encoding: UTF-8
 License: GPL-3 + file LICENSE
 biocViews: Genetics, Epigenetics, ATACSeq, RNASeq, SingleCell
 Depends:
-   R (>= 4.1.0),
+   R (>= 4.4.0),
    SummarizedExperiment,
    RaggedExperiment,
    HDF5Array

NAMESPACE

@@ -23,6 +23,7 @@ export(getChrs)
 export(getCorMatrix)
 export(getDenoisedCorMatrix)
 export(getDomainInflections)
+export(getGenome)
 export(getGlobalMeans)
 export(getMatrixBlocks)
 export(getOpenSeas)

R/getABSignal.R

@@ -21,7 +21,7 @@
 #' #Modified from https://www.biostars.org/p/225520/
 #' random_genomic_int <- data.frame(chr = rep("chr14", 100))
 #' random_genomic_int$start <- apply(random_genomic_int, 1, function(x) {
-#'   round(runif(1, 0, getSeqLengths(chr = x)[[1]]), 0)
+#'   round(runif(1, 0, getSeqLengths(getGenome("hg19"), chr = x)[[1]]), 0)
 #' })
 
 #' random_genomic_int$end <- random_genomic_int$start + runif(1, 1, 1000)

R/getBinMatrix.R

@@ -26,9 +26,10 @@
 #' 
 #' #Generate random genomic intervals of 1-1000 bp on chr1-22
 #' #Modified from https://www.biostars.org/p/225520/
+#' genome.gr <- getGenome("hg19")
 #' random_genomic_int <- data.frame(chr = rep("chr14", 100))
-#' random_genomic_int$start <- apply(random_genomic_int, 1, function(x) { 
-#'   round(runif(1, 0, getSeqLengths(chr = x)[[1]]), 0)
+#' random_genomic_int$start <- apply(random_genomic_int, 1, function(x) {
+#'   round(runif(1, 0, getSeqLengths(genome.gr, chr = x)[[1]]), 0)
 #' })
 #' random_genomic_int$end <- random_genomic_int$start + runif(1, 1, 1000)
 #' random_genomic_int$strand <- "*"
@@ -60,21 +61,7 @@ getBinMatrix <- function(x, genloc, chr = "chr1", chr.start = 0,
   }
 
   #which genome do we have
-  genome <- match.arg(genome)
-
-  if (is.null(chr.end)) {
-    if (genome %in% c("hg19", "hg38", "mm9", "mm10")) {
-      chr.end <- switch(genome,
-                        hg19 = getSeqLengths(genome = "hg19", chr = chr),
-                        hg38 = getSeqLengths(genome = "hg38", chr = chr),
-                        mm9 = getSeqLengths(genome = "mm9", chr = chr),
-                        mm10 = getSeqLengths(genome = "mm10", chr = chr))
-    }
-    else {
-      message("Don't know what to do with ", genome)
-      stop("If you'd like to use an unsupported genome, specify chr.end to an appropriate value...")
-    }
-  }
+  chr.end <- chr.end %||% getSeqLengths(getGenome(genome), chr = chr)
 
   start <- seq(chr.start, chr.end, res) #Build the possible bin ranges given resolution
   end <- c(start[-1], chr.end) - 1L #Set the end ranges for desired resolution

R/getCorMatrix.R

@@ -20,7 +20,7 @@
 #' # Modified from https://www.biostars.org/p/225520/
 #' random_genomic_int <- data.frame(chr = rep("chr14", 100))
 #' random_genomic_int$start <- apply(random_genomic_int, 1, function(x) { 
-#'   round(runif(1, 0, getSeqLengths(chr = x)[[1]]), 0)
+#'   round(runif(1, 0, getSeqLengths(getGenome("hg19"), chr = x)[[1]]), 0)
 #' })
 #' random_genomic_int$end <- random_genomic_int$start + runif(1, 1, 1000)
 #' random_genomic_int$strand <- "*"

R/plotAB.R

@@ -33,7 +33,7 @@
 #' # Modified from https://www.biostars.org/p/225520/
 #' random_genomic_int <- data.frame(chr = rep("chr14", 100))
 #' random_genomic_int$start <- apply(random_genomic_int, 1, function(x) { 
-#'   round(runif(1, 0, getSeqLengths(chr = x)[[1]]), 0)
+#'   round(runif(1, 0, getSeqLengths(getGenome("hg19"), chr = x)[[1]]), 0)
 #' })
 #' random_genomic_int$end <- random_genomic_int$start + runif(1, 1, 1000)
 #' random_genomic_int$strand <- "*"

R/utils.R

@@ -168,39 +168,60 @@ removeEmptyBoots <- function(obj) {
   Filter(Negate(anyNA), obj)
 }
 
-#' Get the seqlengths of a chromosome
+#' Get a GRanges object from bundled compartmap genomes
+#'
+#' @param genome The desired genome to use ("hg19", "hg38", "mm9", "mm10")
+#' @param type The type of data - full genome or open sea regions
+#'
+#' @return Granges of the genome
+#'
+#' @examples
+#' hg19 <- getGenome(genome = "hg19")
+#'
+#' @export
+getGenome <- function(
+  genome = c("hg19", "hg38", "mm9", "mm10"),
+  type = "genome"
+) {
+  genome.name <- match.arg(genome) |> tryCatch(error = function(e) {
+    e <- gsub("'arg'", "'genome'", e)
+    msg <- paste0(e, "Only human and mouse genomes are supported for the time being.")
+    stop(msg)
+  })
+  gr <- switch(type,
+    genome = paste0(genome.name, ".gr"),
+    openseas = paste0("openSeas.", genome.name)
+  )
+  return(get(gr))
+}
+
+#' Get the seqlengths of a chromosome from a given genome's GRanges
 #'
 #' The goal for this function is to eliminate the need to lug around
 #' large packages when we only want seqlengths for things.
 #'
-#' @param genome The desired genome to use ("hg19", "hg38", "mm9", "mm10")
+#' @param genome.gr A GRanges object of the genome (from `getGenome()`)
 #' @param chr What chromosome to extract the seqlengths of
 #'
 #' @return The seqlengths of a specific chromosome
 #'
 #' @importFrom GenomeInfoDb seqlengths seqlevels
-#' @importFrom utils data
 #' @import GenomicRanges
 #'
 #' @examples
-#' hg19.chr14.seqlengths <- getSeqLengths(genome = "hg19", chr = "chr14")
+#' hg19.chr14.seqlengths <- getSeqLengths(getGenome('hg19'), chr = "chr14")
 #'
 #' @export
-getSeqLengths <- function(
-  genome = c("hg19", "hg38", "mm9", "mm10"),
-  chr = "chr14"
-) {
-  # eventually we should support arbitrary genomes
-  genome.name <- match.arg(genome)
-  # check if the genome used exists in what is currently supported, stopping if not
-  if (!genome.name %in% c("hg19", "hg38", "mm9", "mm10")) stop("Only human and mouse are supported for the time being.")
-  # import
-  genome.gr <- get(paste0(genome.name, ".gr"))
-
-  # make sure that the chromosome specified exists in the seqlevels
-  if (!chr %in% seqlevels(genome.gr)) stop("Desired chromosome is not found in the seqlevels of ", genome.name)
-  # get the seqlengths
+getSeqLengths <- function(genome.gr, chr = "chr14") {
   sl <- seqlengths(genome.gr)[chr]
+  if (is.na(sl)) {
+    genome.build <- unique(genome(genome.gr))
+    msg <- paste(
+      chr, "not found in seqlevels of", genome.build,
+      "- check that the 'genome' and 'chr' arguments are correct"
+    )
+    stop(msg)
+  }
   return(sl)
 }
 
@@ -463,21 +484,15 @@ filterOpenSea <- function(
   genome = c("hg19", "hg38", "mm10", "mm9"),
   other = NULL
 ) {
-  # get the desired open sea loci given the genome
-  genome.name <- match.arg(genome)
-  if (is.null(other)) {
-    openseas.genome <- get(paste0("openSeas.", genome.name))
-  } else {
-    # check if it's a GRanges flavored object
-    if (!is(other, "GRanges")) stop("The 'other' input needs to be a GRanges of open sea regions")
-    openseas.genome <- other
-  }
+  # get the desired open sea loci given the genome GRanges
+  openseas.genome <- other %||% getGenome(genome, type = "openseas")
+  stopifnot("The 'other' input needs to be a GRanges of open sea regions" = is(openseas.genome, "GRanges"))
+
   # Subset by overlaps
   message("Filtering to open sea CpG loci...")
   # subset to just CpG loci if CpH or rs probes still exist
   obj <- obj[grep("cg", rownames(obj)), ]
-  obj.openseas <- subsetByOverlaps(obj, openseas.genome)
-  return(obj.openseas)
+  subsetByOverlaps(obj, openseas.genome)
 }
 
 #' Gather open sea CpG from a GRanges of CpG islands