Merge branch 'dev' into dev-nathan-2

Author: nleroy917

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: GenomicDistributions
-Version: 1.1.10
-Date: 2021-08-28
+Version: 1.3.2
+Date: 2022-01-01
 Title: GenomicDistributions: fast analysis of genomic intervals with Bioconductor
 Description: If you have a set of genomic ranges, this package can help you with
  visualization and comparison. It produces several kinds of plots, for example:
@@ -38,6 +38,7 @@ Imports:
     Biostrings,
     plyr,
     dplyr,
+    scales,
     GenomeInfoDb
 Suggests:
     AnnotationFilter,

NAMESPACE

@@ -16,10 +16,11 @@ export(calcFeatureDistRefTSS)
 export(calcGCContent)
 export(calcGCContentRef)
 export(calcNearestGenes)
+export(calcNearestNeighbors)
 export(calcNeighborDist)
-export(calcOpenSignal)
 export(calcPartitions)
 export(calcPartitionsRef)
+export(calcSummarySignal)
 export(calcWidth)
 export(dtToGr)
 export(genomePartitionList)
@@ -39,9 +40,9 @@ export(plotExpectedPartitions)
 export(plotFeatureDist)
 export(plotGCContent)
 export(plotNeighborDist)
-export(plotOpenSignal)
 export(plotPartitions)
 export(plotQTHist)
+export(plotSummarySignal)
 export(retrieveFile)
 import(dplyr)
 import(ggplot2)

NEWS

@@ -2,6 +2,12 @@
 All notable changes to this project will be documented in this file. Here we
 will document changes to major new releases only (not point releases).
 
+## [1.3.2] -- 2020-01-01
+  - Bioconductor released new version
+  - Cell specificity plots are now more genric signal summary plots - the calc function is now "calcSummarySignal", plot function "plotSummarySignal" 
+  - Default plotting of neighbor distance is now on linear scale
+  - Corrected plotting of partition overlap for multiple region sets - now sums to 100 by group
+
 ## [1.1.2] -- 2020-07-07
   - Added functions to calculate and plot dinucleotide frequencies
 

R/chrom-plots.R

@@ -174,10 +174,9 @@ calcChromBins = function(query, bins) {
 #' @param binCount Number of bins to divide the chromosomes into
 #' @return A data.table showing the distribution of regions across bins of the
 #' reference genome.
-#' @export
 #' @examples 
 #' ChromBins = calcChromBinsRef(vistaEnhancers, "hg19")
-calcChromBinsRef = function(query, refAssembly, binCount=10000) {
+calcChromBinsRefSlow = function(query, refAssembly, binCount=3000) {
     .validateInputs(list(refAssembly="character", 
                            query=c("GRanges","GRangesList")))
     # Bin the genome
@@ -189,13 +188,64 @@ calcChromBinsRef = function(query, refAssembly, binCount=10000) {
     return(calcChromBins(query, genomeBins))
 }
 
+
+#' Returns the distribution of query over a reference assembly
+
+#' Given a query set of elements (a GRanges object) and a reference assembly
+#' (*e.g. 'hg38'), this will aggregate and count the distribution of the query
+#' elements across bins of the reference genome. This is a helper function to
+#' create features for common genomes. It is a wrapper of
+#' \code{calcChromBins}, which is more general.
+
+#' @param query A GenomicRanges or GenomicRangesList object with query regions
+#' @param refAssembly A character vector that will be used to grab chromosome
+#'     sizes with \code{getChromSizes}
+#' @param binCount Number of bins to divide the chromosomes into
+#' @return A data.table showing the distribution of regions across bins of the
+#' reference genome.
+#' @export
+#' @examples 
+#' ChromBins = calcChromBinsRef(vistaEnhancers, "hg19")
+calcChromBinsRef = function(query, refAssembly, binCount=3000) {
+   .validateInputs(list(refAssembly="character",
+                           query=c("GRanges","GRangesList")))
+    if (is(query, "GRangesList"))  {
+        # Recurse over each GRanges object
+        x = lapply(query, calcChromBinsRef, refAssembly)
+        # To accommodate multiple regions, we'll need to introduce a new 'name'
+        # column to distinguish them.
+        nameList = names(query)
+       if(is.null(nameList)) {
+            nameList = seq_along(query) # Fallback to sequential numbers
+        }
+        # Append names
+        xb = rbindlist(x)
+        xb$name = rep(nameList, vapply(x, nrow, integer(1)))
+        return(xb)
+    }        
+   # Bin the genome
+    chromSizes = getChromSizes(refAssembly)
+    binnedDT = binChroms(binCount, chromSizes)
+    queryDT = grToDt(query)
+    setnames(binnedDT, "idCol", "chr")
+    queryDT[, midpoint:=start + (end-start)]
+    # Here I use a non-equi join to get the overlaps
+    res = binnedDT[queryDT, .(chr, regionID=ubinID, withinGroupID=x.binID, start=x.start, end=x.end), 
+                    on=.(chr, start<=midpoint, end>=midpoint), nomatch=0L][, list(.N), by=list(chr, start, end, regionID, withinGroupID)][order(regionID),]
+    res[, chr:=factor(chr, levels=unique(res$chr))]
+    return(res)
+}
+
+
+
 #' Plot distribution over chromosomes
 #' 
 #' Plots result from \code{genomicDistribution} calculation
 #' @param genomeAggregate The output from the genomicDistribution function
 #' @param binCount Number of bins (should match the call to
 #'     \code{genomicDistribution})
 #' @param plotTitle Title for plot.
+#' @param ylim Limit of y-axes. Default "max" sets limit to N of biggest bin.
 #' @return A ggplot object showing the distribution of the query 
 #'     regions over bins of
 #' the reference genome.
@@ -206,7 +256,7 @@ calcChromBinsRef = function(query, refAssembly, binCount=10000) {
 #' ChromBins = plotChromBins(agg)
 #' 
 plotChromBins = function(genomeAggregate, binCount=10000, 
-                           plotTitle="Distribution over chromosomes") {
+                           plotTitle="Distribution over chromosomes", ylim="max") {
     .validateInputs(list(genomeAggregate=c("data.table","data.frame")))
 
     if ("name" %in% names(genomeAggregate)){
@@ -231,8 +281,13 @@ plotChromBins = function(genomeAggregate, binCount=10000,
         theme(panel.spacing=unit(0, "lines")) + # Reduce whitespace
         theme(strip.text.y=element_text(size=12, angle=0)) + # Rotate labels
         geom_hline(yintercept=0, color="#EEEEEE") + # Light chrom lines
-        scale_y_continuous(breaks=c(max(genomeAggregate$N)), 
-                            limits=c(0, max(genomeAggregate$N))) +
+        {if (ylim == "max") {
+            scale_y_continuous(breaks = c(max(genomeAggregate$N)),
+                               limits = c(0, max(genomeAggregate$N)))
+        } else {
+            scale_y_continuous(breaks = ylim,
+                               limits = c(0, ylim))
+        }} +
     scale_x_continuous(breaks=c(0, binCount), labels=c("Start", "End")) +
     theme(plot.title=element_text(hjust=0.5)) + # Center title
     ggtitle(plotTitle) +

R/feature-plots.R

@@ -17,7 +17,7 @@
 # @return A vector of genomic distances for each query region relative to its 
 #         closest feature.
 calcFeatureDistBioc = function(query, features) {
-    .validateInputs(list(query=x("GRangesList","GRanges")))
+    .validateInputs(list(query=c("GRangesList","GRanges")))
     if (is(query, "GRangesList")) {
         # Recurse over each GRanges object
         x = lapply(query, calcFeatureDist, features)
@@ -337,6 +337,7 @@ cutDists = function(dists, divisions=NULL, nbins=50,
     labels = labelCuts(sort(divisions), collapse=" to ", infBins=infBins)
     cuts = cut(dists, divisions, labels)
     df = as.data.frame(table(cuts))
+    setDT(df)
     return(df)
 }
 

R/neighbor-distances.R

@@ -2,21 +2,27 @@
 
 
 #' Group regions from the same chromosome together and
-#' calculate the distances between neighboring regions. 
+#' calculate the distances of a region to its upstream and
+#' downstream neighboring regions.  
 #' Distances are then lumped into a numeric vector. 
 #'
 #' @param query A GRanges or GRangesList object.
-#'
+#' @param correctRef A string indicating the reference genome
+#' to use if distances are corrected for the number of 
+#' regions in a regionSet. 
+#' 
 #' @return A numeric vector or list with different vectors containing the
-#'  distances within neighboring regions.
+#'  distances of regions to their upstream/downstream neighbors.
 #' @export
 #' @examples 
 #' dist = calcNeighborDist(vistaEnhancers)
-calcNeighborDist =  function(query) {
-    .validateInputs(list(query=c("GRanges","GRangesList")))
+calcNeighborDist =  function(query, correctRef="None") {
+    .validateInputs(list(query=c("GRanges","GRangesList"),
+                         correctRef=c("character")))
     # lapply if a GRangeslist is provided
     if (is(query, "GRangesList")) {
-        dist = lapply(query, calcNeighborDist)
+        dist = lapply(query,
+                      function(x){calcNeighborDist(x, correctRef = correctRef)})
         namelist = names(query)
         if (is.null(namelist)) {
             newnames = seq_along(query)
@@ -30,10 +36,20 @@ calcNeighborDist =  function(query) {
     querydts = splitDataTable(querydt, "chr")
     distanceVectors = lapply(querydts, neighbordt)
     d = as.vector(unlist(distanceVectors))
-    # remove overlaps for log10 transformation
-    dcvec = d[!(d == "0")] 
-    dcvec = log10(dcvec)
-    return(dcvec)
+    # remove overlaps
+    dcvec = d[!(d == "0")]
+    # Correct for number of regions
+    if (!correctRef=="None") {
+        chromSizes = getChromSizes(correctRef)
+        genomelen = sum(chromSizes)
+        meanWidth = mean(calcWidth(query))
+        expectedDist = genomelen/nrow(querydt) - meanWidth
+        correctedDist = log10(dcvec/expectedDist)
+        return(correctedDist)
+    # If we just want to look at the raw neighbor distances
+    } else {
+        return(dcvec)
+    }
 }
 
 #' Internal helper function to calculate distance 
@@ -46,61 +62,133 @@ neighbordt = function(querydt)  {
     # there should be at least 2 regions for each chr
     if (nrow(querydt) > 1) {
         endVect = abs(querydt[, diff(end)])
-        regionWidth = querydt[, (end-start)]
+        regionWidth = querydt[, (end-start+1)]
         distancesVector = endVect - regionWidth[-1]
         # neg values represent overlaps between neighbor regions, set those to 0
         distancesVector[which(distancesVector < 0)] = 0 
         return(distancesVector)
   }
 }
-  
 
-#' Plot the distances between neighboring regions.The distance in the 
-#' x axis is log10 transformed for ease of comparison between 
-#' different regionsets and to account for outliers. 
+
+#' Group regions from the same chromosome together and
+#' compute the distance of a region to its nearest neighbor. 
+#' Distances are then lumped into a numeric vector. 
+#'
+#' @param query A GRanges or GRangesList object.
+#' @param correctRef A string indicating the reference genome
+#' to use if Nearest neighbor distances are corrected for the 
+#' number of regions in a regionSet. 
+#'
+#' @return A numeric vector or list of vectors containing the
+#'  distance of regions to their nearest neighbors.
+#' @export
+#' @examples 
+#' Nneighbors = calcNearestNeighbors(vistaEnhancers)
+calcNearestNeighbors = function(query, correctRef="None") {
+    .validateInputs(list(query=c("GRanges","GRangesList"),
+                         correctRef=c("character")))
+    # lapply if a GRangeslist is provided
+    if (is(query, "GRangesList")) {
+        dist = lapply(query,
+                      function(x){calcNearestNeighbors(x, correctRef = correctRef)})
+        namelist = names(query)
+        if (is.null(namelist)) {
+            newnames = seq_along(query)
+            namelist = newnames
+            # Append names
+            names(dist) = namelist 
+        }
+        return(dist)
+    }
+    # Calculate nearest neighbors in a vectorized manner
+    dist = calcNeighborDist(query)
+    upstream = dist[-length(dist)]
+    downstream = dist[-1]
+    dt = data.table(i=upstream, j=downstream)
+    pairmins = dt[, pmin(i, j)]
+    # First and last distances are default nearest neighbors
+    nNeighbors = c(dist[1], pairmins, dist[length(dist)])
+    # Correct for number of regions
+    if (!correctRef=="None") {
+        chromSizes = getChromSizes(correctRef)
+        genomelen = sum(chromSizes)
+        meanWidth = mean(calcWidth(query))
+        expectedDist = genomelen/length(query) - meanWidth
+        correctedDist = log10(nNeighbors/expectedDist)
+        return(correctedDist)
+    } else {
+        return(nNeighbors)
+    }
+}
+
+#' Plot the distances from regions to their upstream/downstream neighbors
+#' or nearest neighbors. Distances can be passed as either raw bp or
+#' corrected for the number of regions (log10(obs/exp)), but this has
+#' to be specified in the function parameters. 
 #' 
-#' @param dcvec A numeric vector or list with vectors containing distances 
-#' between neighbor regions. Produced by \code{calcNeighborDist}
+#' @param dcvec A numeric vector or list of vectors containing distances 
+#' to upstream/downstream neighboring regions or to nearest neighbors. 
+#' Produced by \code{calcNeighborDist} or \code{calcNearestNeighbors}
+#' @param correctedDist A logical indicating if the plot axis should
+#' be adjusted to show distances corrected for the number of regions
+#' in a regionset.
+#' @param Nneighbors A logical indicating whether legend should be adjusted
+#' if Nearest neighbors are being plotted. Default legend shows distances
+#' to upstream/downstream neighbors.  
 #'
 #' @return A ggplot density object showing the distribution of
-#' log10 transformed distances.  
+#' raw or corrected distances.  
 #' @export
 #' @examples
 #' numVector = rnorm(400, mean=5, sd=0.1)
 #' d = plotNeighborDist(numVector)
-plotNeighborDist = function(dcvec) {
+plotNeighborDist = function(dcvec, correctedDist=FALSE,
+                            Nneighbors=FALSE) {
     .validateInputs(list(dcvec=c("numeric","list")))
-    # if input is list, conver it to a data frame with 
+    # if input is list, convert it to a data frame with 
     # value and region set name, if input is vector - make a single
     # columns data.frame
-    if (is(dcvec, "list")){
-      nameList = names(dcvec)
-      vectorLengths = unlist(lapply(dcvec, length))
-      distReshaped = data.frame(value = unlist(dcvec),
-                              regionSet = rep(nameList, vectorLengths))
-    } else {
-      distReshaped = data.frame(value = dcvec)
-    }
-  
     if (is(dcvec, "list")) {
-        g = ggplot2::ggplot(distReshaped, aes(x=value, 
-                                              fill=regionSet, 
+        nameList = names(dcvec)
+        vectorLengths = unlist(lapply(dcvec, length))
+        distReshaped = data.frame(value = unlist(dcvec),
+                                  regionSet = rep(nameList, vectorLengths))
+        g = ggplot2::ggplot(distReshaped, aes(x=value,
+                                              fill=regionSet,
                                               colour=regionSet)) +
-          geom_density(alpha=0.5) +
-          theme_classic() +
-          theme(legend.position = "bottom")
+          geom_density(alpha=0.4)
     } else {
+        distReshaped = data.frame(value = dcvec)
         g = ggplot2::ggplot(distReshaped, aes(x=value)) +
-          geom_density(alpha=0.4) + 
-          theme_classic()
+          geom_density() 
+    }
+    if (correctedDist==TRUE) {
+        g = g + 
+          xlab(expression(log[10](over(Obs, Exp)))) +
+          geom_vline(xintercept = 0, linetype="dashed") +
+          ggtitle("Corrected neighboring regions distance distribution") 
+    } else {
+        g = g + 
+          xlab(expression("bp distance")) +
+          scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
+                        labels = scales::trans_format("log10", 
+                                                  scales::math_format(10^.x))) +
+          ggtitle("Neighboring regions distance distribution") 
     }
     g = g + 
-        xlab(expression(log[10]*("bp distance"))) +
-        xlim(0, 10) +
-        theme(aspect.ratio=1) +
-        theme_blank_facet_label() +
-        ggtitle("Neighboring regions distance distribution") +
-        theme(plot.title = element_text(hjust=0.5)) 
+      theme_classic() +
+      theme(aspect.ratio=1,
+            plot.title = element_text(hjust=0.5),
+            legend.position = "bottom") +
+      theme_blank_facet_label()
+  
+  # Adjust legend if plotting nearest neighbors 
+    if (Nneighbors==TRUE){
+        g = g + 
+          labs(fill="regionSet Nneighbors", 
+               colour="regionSet Nneighbors")
+    }
     return(g)
 }