R/gsaregion.R

In missMethyl: Analysing Illumina HumanMethylation BeadChip Data

#' Generalised gene set testing for Illumina's methylation array data
#' 
#' Given a user specified list of gene sets to test, \code{gsaregion} tests
#' whether differentially methylated regions (DMRs) identified from Illumina's
#' Infinium HumanMethylation450 or MethylationEPIC array are enriched, taking
#' into account the differing number of probes per gene present on the array.
#' 
#' This function extends \code{goregion}, which only tests GO and KEGG pathways.
#' \code{gsaregion} can take a list of user specified gene sets and test whether
#' the significant DMRs are enriched in these pathways. This function takes a
#' \code{GRanges} object of DMR coordinates, maps them to CpG sites on the
#' array and then to Entrez Gene IDs, and tests for enrichment using Wallenius' 
#' noncentral hypergeometric test, taking into account the number of CpG sites 
#' per gene on the 450K/EPIC array.  If \code{prior.prob} is set to FALSE, then 
#' prior probabilities are not used and it is assumed that each gene is equally
#' likely to have a significant CpG site associated with it. 
#' 
#' The testing now also takes into account that some CpGs map to multiple genes. 
#' For a small number of gene families, this previously caused their associated 
#' GO categories/gene sets to be erroneously overrepresented and thus highly 
#' significant. If \code{fract.counts=FALSE} then CpGs are allowed to map to 
#' multiple genes (this is NOT recommended).
#' 
#' Genes associated with each CpG site are obtained from the annotation package
#' \code{IlluminaHumanMethylation450kanno.ilmn12.hg19} if the array type is
#' "450K". For the EPIC array, the annotation package
#' \code{IlluminaHumanMethylationEPICanno.ilm10b4.hg19} is used. To use a
#' different annotation package, please supply it using the \code{anno}
#' argument. 
#' 
#' In order to get a list which contains the mapped Entrez gene IDS,
#' please use the \code{getMappedEntrezIDs} function. The \code{topGSA} function 
#' can be used to display the top 20 most enriched pathways.
#' 
#' If you are interested in which genes overlap with the genes in the gene set, 
#' setting \code{sig.genes} to TRUE will output an additional column in the 
#' results data frame that contains all the significant differentially 
#' methylated gene symbols, comma separated. The default is FALSE.
#' 
#' @param regions \code{GRanges} Object of DMR coordinates to test for GO term
#' enrichment.
#' @param all.cpg Character vector of all CpG sites tested. Defaults to all CpG
#' sites on the array.
#' @param collection A list of user specified gene sets to test. Can also be a
#' single character vector gene set. Gene identifiers must be Entrez Gene IDs.
#' @param array.type The Illumina methylation array used. Options are "450K" or
#' "EPIC". Defaults to "450K".
#' @param plot.bias Logical, if true a plot showing the bias due to the
#' differing numbers of probes per gene will be displayed.
#' @param prior.prob Logical, if true will take into account the probability of
#' significant differentially methylation due to numbers of probes per gene. If
#' false, a hypergeometric test is performed ignoring any bias in the data.
#' @param anno Optional. A \code{DataFrame} object containing the complete
#' array annotation as generated by the \code{\link{minfi}}
#' \code{\link{getAnnotation}} function. Speeds up execution, if provided.
#' @param equiv.cpg Logical, if true then equivalent numbers of cpgs are used
#' for odds calculation rather than total number cpgs. Only used if 
#' \code{prior.prob=TRUE}.
#' @param fract.counts Logical, if true then fractional counting of cpgs is used
#' to account for cgps that map to multiple genes. Only used if 
#' \code{prior.prob=TRUE}.
#' @param genomic.features Character vector or scalar indicating whether the 
#' gene set enrichment analysis should be restricted to CpGs from specific 
#' genomic locations. Options are "ALL", "TSS200","TSS1500","Body","1stExon",
#' "3'UTR","5'UTR","ExonBnd"; and the user can select any combination. Defaults
#' to "ALL".
#' @param sig.genes Logical, if true then the significant differentially 
#' methylated genes that overlap with the gene set of interest is outputted 
#' as the final column in the results table. Default is FALSE.
#' 
#' @return A data frame with a row for each gene set and the following columns:
#' \item{N}{ number of genes in the gene set } \item{DE}{ number of genes that
#' are differentially methylated } \item{P.DE}{ p-value for over-representation
#' of the gene set } \item{FDR}{ False discovery rate, calculated using the
#' method of Benjamini and Hochberg (1995).  } \item{SigGenesInSet}{ Significant 
#' differentially methylated genes overlapping with the gene set of interest. }
#' @author Jovana Maksimovic
#' @seealso \code{\link{gometh},\link{goregion},\link{gsameth},
#' \link{getMappedEntrezIDs}}
#' 
#' @references Phipson, B., Maksimovic, J., and Oshlack, A. (2016). missMethyl:
#' an R package for analysing methylation data from Illuminas
#' HumanMethylation450 platform. \emph{Bioinformatics}, \bold{15};32(2),
#' 286--8. 
#' 
#' Geeleher, P., Hartnett, L., Egan, L. J., Golden, A., Ali, R. A. R.,
#' and Seoighe, C. (2013). Gene-set analysis is severely biased when applied to
#' genome-wide methylation data. \emph{Bioinformatics}, \bold{29}(15),
#' 1851--1857. 
#' 
#' Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A.
#' (2010). Gene ontology analysis for RNA-seq: accounting for selection bias.
#' \emph{Genome Biology}, 11, R14. 
#' 
#' Ritchie, M. E., Phipson, B., Wu, D., Hu, Y.,
#' Law, C. W., Shi, W., and Smyth, G. K. (2015). limma powers differential
#' expression analyses for RNA-sequencing and microarray studies. \emph{Nucleic
#' Acids Research}, gkv007. 
#' 
#' Benjamini, Y., and Hochberg, Y. (1995). Controlling
#' the false discovery rate: a practical and powerful approach to multiple
#' testing. \emph{Journal of the Royal Statistical Society Series}, B,
#' \bold{57}, 289-300.
#' @examples
#' 
#' \dontrun{ # to avoid timeout on Bioconductor build
#' library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
#' library(limma)
#' library(DMRcate)
#' library(ExperimentHub)
#' library(org.Hs.eg.db)
#' 
#' # Follow the example for the dmrcate function to get some EPIC data from
#' # ExperimentHub
#' eh <- ExperimentHub()
#' FlowSorted.Blood.EPIC <- eh[["EH1136"]]
#' tcell <- FlowSorted.Blood.EPIC[,colData(FlowSorted.Blood.EPIC)$CD4T==100 |
#'                                 colData(FlowSorted.Blood.EPIC)$CD8T==100]
#' detP <- detectionP(tcell)
#' remove <- apply(detP, 1, function (x) any(x > 0.01))
#' tcell <- tcell[!remove,]
#' tcell <- preprocessFunnorm(tcell)
#' #Subset to chr2 only
#' tcell <- tcell[seqnames(tcell) == "chr2",]
#' tcellms <- getM(tcell)
#' tcellms.noSNPs <- rmSNPandCH(tcellms, dist=2, mafcut=0.05)
#' tcell$Replicate[tcell$Replicate==""] <- tcell$Sample_Name[tcell$Replicate==""]
#' tcellms.noSNPs <- avearrays(tcellms.noSNPs, tcell$Replicate)
#' tcell <- tcell[,!duplicated(tcell$Replicate)]
#' tcell <- tcell[rownames(tcellms.noSNPs),]
#' colnames(tcellms.noSNPs) <- colnames(tcell)
#' assays(tcell)[["M"]] <- tcellms.noSNPs
#' assays(tcell)[["Beta"]] <- ilogit2(tcellms.noSNPs)
#' 
#' # Perform region analysis
#' type <- factor(tcell$CellType)
#' design <- model.matrix(~type) 
#' myannotation <- cpg.annotate("array", tcell, arraytype = "EPIC",
#'                              analysis.type="differential", design=design, 
#'                              coef=2)
#' # Run DMRCate with beta value cut-off filter of 0.1                              
#' dmrcoutput <- dmrcate(myannotation, lambda=1000, C=2, betacutoff = 0.1)
#' regions <- extractRanges(dmrcoutput)
#' length(regions)
#' 
#' ann <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
#' # All CpG sites tested (limited to chr 2)
#' allcpgs <- rownames(tcell)
#' # Use org.Hs.eg.db to extract a GO term
#' GOtoID <- suppressMessages(select(org.Hs.eg.db, keys=keys(org.Hs.eg.db), 
#'                                   columns=c("ENTREZID","GO"), 
#'                                   keytype="ENTREZID"))
#' keep.set1 <- GOtoID$GO %in% "GO:0010951"
#' set1 <- GOtoID$ENTREZID[keep.set1]
#' keep.set2 <- GOtoID$GO %in% "GO:0042742"
#' set2 <- GOtoID$ENTREZID[keep.set2]
#' keep.set3 <- GOtoID$GO %in% "GO:0031295"
#' set3 <- GOtoID$ENTREZID[keep.set3]
#' # Make the gene sets into a list
#' sets <- list(set1, set2, set3)
#' names(sets) <- c("GO:0010951","GO:0042742","GO:0031295")
#' 
#' # Testing with prior probabilities taken into account
#' # Plot of bias due to differing numbers of CpG sites per gene
#' gst <- gsaregion(regions = regions, all.cpg = allcpgs, collection = sets, 
#'                  array.type = "EPIC", plot.bias = TRUE, prior.prob = TRUE,
#'                  anno = ann)
#' topGSA(gst)
#' 
#' # Add significant genes in gene set to output
#' gst <- gsaregion(regions = regions, all.cpg = allcpgs, collection = sets, 
#'                  array.type = "EPIC", plot.bias = TRUE, prior.prob = TRUE,
#'                  anno = ann, sig.genes = TRUE)
#' topGSA(gst)
#' }
#' 
#' @export gsaregion
gsaregion <- function(regions, all.cpg=NULL, collection, 
                      array.type = c("450K","EPIC"), plot.bias=FALSE, 
                      prior.prob=TRUE, anno=NULL, equiv.cpg = TRUE, 
                      fract.counts=TRUE, 
                      genomic.features = c("ALL", "TSS200","TSS1500","Body",
                                           "1stExon","3'UTR","5'UTR","ExonBnd"),
                      sig.genes = FALSE)
    # Generalised version of goregion with user-specified gene sets 
    # Gene sets collections must be Entrez Gene ID
    # Takes into account probability of differential methylation based on
    # numbers of probes on array per gene
    # Jovana Maksimovic
    # 26 April 2019. Last updated 11 Setember 2020.
{
    
    array.type <- match.arg(toupper(array.type), c("450K","EPIC"))    
    genomic.features <- match.arg(genomic.features, c("ALL", "TSS200","TSS1500",
                                                      "Body", "1stExon","3'UTR",
                                                      "5'UTR","ExonBnd"), 
                                  several.ok = TRUE)
    
    if(length(genomic.features) > 1 & any(grepl("ALL", genomic.features))){
        message("All input CpGs are used for testing.") 
        genomic.features <- "ALL"      
    } 
    
    if(array.type == "450K" & any(grepl("ExonBnd", genomic.features))){
        stop("'ExonBnd' is not an annotated feature on 450K arrays,
           please remove it from your genomic.feature parameter
           specification.") 
    }
    
    if(is.null(anno)){
        if(array.type=="450K"){
            anno <- minfi::getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19::IlluminaHumanMethylation450kanno.ilmn12.hg19)
        } else {
            anno <- minfi::getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19::IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
        }
    }
    
    if(!is.null(all.cpg)){
        anno <- anno[all.cpg,]
    }
    
    cpgs <- GenomicRanges::GRanges(seqnames = anno$chr, 
                                   ranges = IRanges::IRanges(start = anno$pos, 
                                                             end = anno$pos),
                                   strand = anno$strand,
                                   name = anno$Name)
    
    overlaps <- GenomicRanges::findOverlaps(cpgs,regions)
    sig.cpg <- cpgs$name[from(overlaps)]
    
    result <- gsameth(sig.cpg=sig.cpg, all.cpg=all.cpg, collection=collection, 
                      array.type=array.type, plot.bias=plot.bias, 
                      prior.prob=prior.prob, anno=anno, equiv.cpg=equiv.cpg,
                      fract.counts=fract.counts,
                      genomic.features = genomic.features,
                      sig.genes = sig.genes)
    result
}