R/waldUMI4C.R

In Pasquali-lab/UMI4Cats: UMI4Cats: Processing, analysis and visualization of UMI-4C chromatin contact data

#' DESeq2 Wald test for differential contacts
#' 
#' Using a \code{UMI4C} object, infers the differences between conditions 
#' specified in \code{design} using a Wald Test from \code{DESeq2} package. 
#' @param umi4c UMI4C object as generated by \code{makeUMI4C} or the
#' \code{UMI4C} constructor.
#' @param query_regions \code{GRanges} object containing the coordinates of the 
#' genomic regions you want to use to perform the analysis in specific genomic 
#' intervals. Default: NULL.
#' @param subset If \code{query_regions} are provided, how to subset the UMI4C
#' object:  "sum" for summing raw UMIs in fragments overlapping 
#' \code{query_regions} (default) or "overlap" for selecting overlapping fragments.
#' @param design A \code{formula} or \code{matrix}. The formula expresses how
#' the counts for each fragment end depend on the variables in \code{colData}.
#' See  \code{\link[DESeq2]{DESeqDataSet}}.
#' @param normalized Logical indicating if the function should return normalized
#'  or raw UMI counts.  Default: TRUE.
#' @param padj_method The method to use for adjusting p-values, see
#' \code{\link[stats]{p.adjust}}.  Default: fdr.
#' @param padj_threshold Numeric indicating the adjusted p-value threshold to
#' use to define significant differential contacts.  Default: 0.05.
#' @return \code{UMI4C} object with the DESeq2 Wald Test results, which can be
#' accessed using \code{\link{resultsUMI4C}}.
#' @export
#' @examples 
#' data("ex_ciita_umi4c")
#' 
#' umi_dif <- waldUMI4C(ex_ciita_umi4c)
waldUMI4C <- function(umi4c,
                      query_regions=NULL,
                      subset="sum", #overlap
                      design=~condition,
                      normalized = TRUE,
                      padj_method = "fdr",
                      padj_threshold = 0.05) {
  umi4c_ori <- umi4c
  ## Obtain only fragments in query_regions, if available
  if(!is.null(query_regions)) {
    if (subset=="overlap") umi4c <- subsetByOverlaps(umi4c, query_regions)   # Option 1: by overlaps
    else if (subset=="sum") umi4c <- combineUMI4C(umi4c, query_regions)   # Option 2: sum
  }
  
  ## Convert to dds
  dds <- UMI4C2dds(umi4c=umi4c,
                   design=design)
  
  ## Run DESeq2
  dds <- DESeq2::DESeq(dds, test="Wald")
  
  ## Convert back to UMI4C
  umi4c_final <- dds2UMI4C(umi4c=umi4c_ori,
                           dds=dds,
                           normalized=normalized,
                           padj_method=padj_method,
                           padj_threshold=padj_threshold)
  
  return(umi4c_final)
}