R/plotPeaks.R

In plbaldoni/mixHMM: Differential Peak Calling via Hidden Markov Model With Mixture of Negative Binomial Distributions

#' Create a plot of ChIP read counts with called peaks and posterior probabilities
#'
#' `plotPeaks()` imports the output from `mixHMM()` and plots the peak
#' calls of a given genomic region.
#'
#' @param output An output from `mixHMM()`
#' @param ranges A vector with two elements that specifies the genomic windows to be plotted.
#' @param ChIP M*N matrix of ChIP read counts, where M is the number of windows in the analyzed genome and N is the number of conditions*replicates
#' @param group vector of length N with condition labels (e.g. c(1,1,2,2,3,3) for three conditions and two replicates each)
#' @param offset M*N matrix of offsets (default is offset = matrix(0,nrow=M,ncol=N))
#' @param peaks An optional vector of length M with zeros (background), ones (differential peaks), and twos (consensus peaks) representing the enriched windows. If not provided,
#' `plotPeaks()` will use the Viterbi sequence from `output`.
#'
#' @return A ggplot
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#' @references \url{https://github.com/plbaldoni/mixHMM}
#'
#' @importFrom RColorBrewer brewer.pal
#' @importFrom tidyr gather
#' @importFrom dplyr as.tbl
#' @importFrom plyr mapvalues
#' @importFrom ggpubr ggarrange
#' @importFrom scales comma
#' @importFrom ZIMHMM createOffset
#' @importFrom RColorBrewer brewer.pal
#' @importFrom SummarizedExperiment assay
#'
#' @examples
#' data(H3K36me3)
#' ChIP = SummarizedExperiment::assay(H3K36me3)
#' offset = ZIMHMM::createOffset(ChIP,method = 'loess',span = 1)
#' group = c(1,1,1,2,2,2,3,3,3)
#' B = 2^(length(unique(group)))-2
#' \dontrun{output = mixHMMConstr(ChIP,Control=NULL,offset,group,B,control = controlPeaks(epsilon.em = rep(1e-6,4)))}
#' \dontrun{plotPeaks(output = test,ranges = c(which.min(abs(start(rowRanges(H3K36me3))-17800000)),
#' which.min(abs(end(rowRanges(H3K36me3))-18150000))),ChIP = ChIP,group = group,offset = offset)}
#'
#' @export
#'
plotPeaks = function(output,ranges,ChIP,group,offset = NULL,peaks = NULL){
    if(length(ranges)!=2 | ranges[2]<ranges[1]){stop('Argument ranges must be a vector of two integers (x,y) with y>x')}
    ranges = ranges[1]:ranges[2]

    color = RColorBrewer::brewer.pal(3,'Set1')[2]

    if(is.null(offset)){offset = matrix(0,nrow=nrow(ChIP),ncol=ncol(ChIP))}

    DT <- as.data.table(ChIP/exp(offset))
    for(i in unique(group)){
        DT[,paste0('Group',i) := rowSums(.SD),.SDcols = which(group==i)]
    }
    DT <- DT[,paste0('Group',unique(group)),with=F]

    datamelt <- data.table::melt(DT,measure.vars=seq_len(ncol(DT)),value.name = 'Counts',variable.name = 'Groups')
    setattr(datamelt$Groups,"levels",paste0('Group ',1:ncol(DT)))
    datamelt[,Window := rep(seq_len(nrow(DT)),ncol(DT))]

    datapeak = data.table::data.table(Window = seq_len(nrow(DT)),Groups = unique(datamelt$Groups)[1])
    if(is.null(peaks)){
        datapeak[,Peaks := output$Viterbi]
    } else{
        datapeak[,Peaks := peaks]
    }
    datapeak$Peaks = plyr::mapvalues(datapeak$Peaks,from=0:2,to=c(NA,'Differential',NA))

    ### Figure 1: Read Counts and Peak Calls ###
    maxy = max(datamelt[(Window%in%ranges),Counts])*1.1
    fig.ChIP = ggplot2::ggplot(data=datamelt[(Window%in%ranges),],ggplot2::aes(x=Window,y=Counts))+
        ggplot2::geom_line()+
        ggplot2::facet_grid(rows=ggplot2::vars(Groups))+
        ggplot2::geom_segment(inherit.aes=FALSE,data = datapeak[(Window%in%ranges) & !is.na(Peaks),],ggplot2::aes(x=Window,xend=Window+1,y=maxy,yend=maxy),size=2,color=color)+
        theme_bw()+xlab('Genomic Window')+ylab('Adjusted Read Counts')+
        theme(axis.title.x = element_blank(),axis.line.x = element_blank(),axis.ticks.x = element_blank(),axis.text.x = element_blank())

    ### Figure 2: Post. Probabilities ###
    PostProb = data.table::data.table(Window=1:nrow(DT),output$Prob,Label='PP')
    PostProb <- tidyr::gather(dplyr::as.tbl(PostProb),Component,PP,PostProb1:PostProb3)
    PostProb$Component = plyr::mapvalues(PostProb$Component,from=paste0('PostProb',1:3),to=c(NA,'Differential',NA))

    fig.Prob = ggplot2::ggplot(data=PostProb[PostProb$Window%in%ranges & !is.na(PostProb$Component),],ggplot2::aes(x=Window,y=PP))+
        ggplot2::facet_grid(rows=ggplot2::vars(Label))+
        ggplot2::geom_area(position='identity',alpha=0.5,color=color,fill=color)+
        ggplot2::scale_y_continuous(limits = c(0,1),breaks=c(0,0.5,1))+
        ggplot2::scale_x_continuous(limits = range(PostProb[PostProb$Window%in%ranges,'Window']),labels = scales::comma)+
        ggplot2::xlab(paste0('Genomic Window'))+ylab('Post. Prob.')+
        ggplot2::theme_bw()+
        ggplot2::theme(legend.position = 'bottom',legend.direction = 'horizontal',strip.text.y = element_text(colour = alpha('grey',0.0)))+
        ggplot2::theme(legend.position="none")

    fig = ggpubr::ggarrange(fig.ChIP,fig.Prob,ncol=1,nrow=2,heights = c(0.8,0.2),common.legend = T,legend='bottom')

    return(fig)
}