R/genemodel2.R

In vjcitn/TFutils: TFutils

reqNS = function(pname) if (!requireNamespace(pname)) {
  stop(paste("please install", pname, "to use this function"))
}

#' use EnsDb to generate an exon-level model of genes identified by symbol
#' @param sym a character() vector of gene symbols
#' @param resource should be or inherit from EnsDb, 
#' answering exons(), with AnnotationFilter::SymbolFilter 
#' as filter parameter
#' @param columnsKept character vector used as `columns` param in exons()
#' @param \dots passed to exons()
#' @note There are many approaches available to acquiring 'gene models'
#' in Bioconductor; this one emphasizes the use of the exons method
#' for Ensembl annotation.
#' @return data.frame instance with exons in rows
#' @examples
#' if (requireNamespace("EnsDb.Hsapiens.v75")) {
#'  orm = genemodelDF("ORMDL3", EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75)
#'  dim(orm)
#' }
#' head(orm)
#' @export
genemodelDF = function(sym, resource, 
      columnsKept=c("gene_id", "tx_id"), ...) {
  if (is.character(resource)) reqNS(resource)
  reqNS("AnnotationFilter")
  reqNS("GenomicFeatures")
  if (!is(resource, "EnsDb")) warning("resource not of class EnsDb")
  exs = GenomicFeatures::exons(resource, filter=AnnotationFilter::SymbolFilter(sym), 
         columns=columnsKept, ...)
  names(exs) = NULL
  as(exs, "data.frame")
}

#' create a GeneRegionTrack instance for selected symbols
#' @param sym character vector of gene symbols, should be neighboring genes
#' @param id_elem vector of names of columns generated by genemodelDF to be
#' used to label transcripts
#' @param resource should be or inherit from EnsDb, 
#' answering exons(), with AnnotationFilter::SymbolFilter 
#' as filter parameter
#' @note This function helps to display the locations of TF binding sites
#' in the context of complex gene models.  A complication is that we
#' have nice visualization of quantitative affinity predictions for TFs
#' in the vignette, based on ggplot2, but it is not clear how to use
#' that specific code to work with Gviz.
#' @param \dots passed to genemodelDF
#' @return instance of Gviz GeneRegionTrack
#' @examples
#' if (requireNamespace("EnsDb.Hsapiens.v75") &
#'     requireNamespace("Gviz")) {
#'  orm = genemodForGviz("ORMDL3", resource= EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75)
#'  orm
#'  Gviz::plotTracks(orm, showId=TRUE) # change id_elem for shorter id string
#' }
#' @export
genemodForGviz = function(sym="ORMDL3", id_elem=c("symbol", "tx_id"),
        resource = EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75, ...) {
 reqNS("Gviz")
 exdf = genemodelDF(sym=sym, resource=resource, ...)
 cind = which(names(exdf)=="seqnames")
 if (length(cind)==1) names(exdf)[cind] = "chromosome"
 Gviz::GeneRegionTrack(exdf, transcript=exdf$tx_id, 
     symbol = apply(exdf[,id_elem,drop=FALSE],1, paste, collapse=" "), ...)
}