R/pgxFreqplot.R

In progenetix/pgxRpi: R wrapper for Progenetix

#' Plot CNV frequency data
#'
#' Thie function plots the frequency of deletions and duplications
#'
#' @param data CNV frequency object returned by the `pgxLoader` or `segtoFreq` functions.
#' @param chrom A vector specifying which chromosomes to plot. If NULL, the plot will cover the entire genome. 
#' If specified, the frequencies are plotted with one panel for each chromosome. Default is NULL.
#' @param layout Number of columns and rows in plot. Only used in plot by chromosome. Default is c(1,1).
#' @param filters Index or string value indicating which filter to plot. The length of filters
#' is limited to one if the parameter `circos` is FALSE. Default is the first filter.
#' @param circos A logical value indicating whether to return a circos plot. If TRUE, it returns a circos plot 
#' that can display and compare multiple filters. Default is FALSE.
#' @param highlight Indices of genomic bins to be highlighted in red.
#' @param assembly A string specifying the genome assembly version to apply to CNV frequency plotting. 
#' Allowed options are "hg19" and "hg38". Default is "hg38".
#' @return The binned CNV frequency plot
#' @importFrom grDevices dev.cur dev.new devAskNewPage
#' @importFrom graphics abline axis mtext par polygon rect text title
#' @importFrom utils data
#' @importFrom methods is
#' @export
#' @examples
#' ## load necessary data (this step can be skipped in real implementation)
#' data("hg38_cytoband")
#' ## get frequency data
#' freq <- pgxLoader(type="cnv_frequency", output ='pgxfreq', filters="NCIT:C3512")
#' ## visualize
#' pgxFreqplot(freq)

pgxFreqplot <- function(data,chrom=NULL,layout=c(1,1),filters=NULL,circos = FALSE,highlight=NULL,assembly = 'hg38'){
    if (circos){
        return(cplotpgxFreq(data,filters,highlight,assembly))
    }

    if (length(filters) > 1){
        stop("The length of filters exceeds limit")
    }

    type <- ifelse(is.null(chrom),"genome","bychrom")
    switch(type,
            genome = genomeFreq(data=data,id=filters,highlight=highlight,assembly = assembly),
            bychrom = chromosomeFreq(data=data,chrom=chrom,layout=layout,id=filters,highlight=highlight,assembly = assembly))
}


cplotpgxFreq  <- function(data,filters,highlight,assembly){
    if (is.null(filters)) filters <- 1
    data_list <- list()
    # input is pgxfreq
    if (is(data, "CompressedGRangesList")){
        meta <- S4Vectors::mcols(data)
        id_name <- rownames(meta[filters,])
        for (i in id_name){
          ind_data <- data[[i]]
          data_list[[i]] <- data.frame(chr=paste0('chr',as.character(GenomicRanges::seqnames(ind_data))),
                                       start=GenomicRanges::start(ind_data),
                                       end=GenomicRanges::end(ind_data),
                                       gain=S4Vectors::mcols(ind_data)$gain_frequency,
                                       loss=S4Vectors::mcols(ind_data)$loss_frequency)
        } 
    } else if (is(data, "RangedSummarizedExperiment")){
        meta <- SummarizedExperiment::colData(data)
        id_name <- rownames(meta[filters,])
        range.info <- unique(GenomicRanges::granges(SummarizedExperiment::rowRanges(data)))
        for (i in id_name){
            gain.freq <- SummarizedExperiment::assay(data)[which(SummarizedExperiment::rowData(data)$type == "DUP"),i]
            loss.freq <- SummarizedExperiment::assay(data)[which(SummarizedExperiment::rowData(data)$type == "DEL"),i]
            data_list[[i]] <-  data.frame(chr=paste0('chr',as.character(GenomicRanges::seqnames(range.info))),
                                             start=GenomicRanges::start(range.info),
                                             end=GenomicRanges::end(range.info),
                                             gain=gain.freq,
                                             loss=-loss.freq)
        }
    } else{
        stop("\n The input is invalid \n")
    }
 

  # Identify highlight location
    if (!is.null(highlight)){
        h_chr <- data_list[[1]][['chr']][highlight]
        h_start <- data_list[[1]][['start']][highlight]
    }


  # set text size and color
    cex_text <- ifelse(length(data_list) == 1, 0.8, 0.6)
    col.gain <-'#FFC633'
    col.loss <- '#33A0FF'
    bar.col.gain <- col.gain
    bar.col.loss <- col.loss

      # plot
    for (i in seq_len(length(id_name))){
        id <- id_name[i]
        data_cir <- data_list[[id]]
        # Initialize
        if (i == 1){
            circlize::circos.clear()
            circlize::circos.par(start.degree=90, gap.degree=c(rep(1,23),10))
            circlize::circos.initializeWithIdeogram(species = assembly)
        }
        circlize::circos.genomicTrack(data_cir, ylim = c(-100, 100),
                                      numeric.column = c("gain", "loss"),
                                      panel.fun = function(region, value, ...){
                                          value1 <- unlist(value[1])
                                          value2 <- unlist(value[2])
                                          pos <- unlist((region[2]+region[1])/2)
                                          # check if highlight exists
                                          if (!is.null(highlight)){
                                          # check if the plotting chr is highlight chr
                                              idx <- as.numeric(which(h_chr == circlize::CELL_META$sector.index))
                                          if (length(idx) != 0){
                                              bar.col.gain <- rep(col.gain,length(value1))
                                              bar.col.gain[which(unlist(region[1]) %in% h_start[idx])] <- 'red'
                                              bar.col.loss <- rep(col.loss,length(value1))
                                              bar.col.loss[which(unlist(region[1]) %in% h_start[idx])] <- 'red'
                                          }
                                        }
                                        circlize::circos.barplot(value=value1,pos=pos, border = bar.col.gain, col = bar.col.gain)
                                        circlize::circos.barplot(value=value2,pos=pos, border = bar.col.loss, col = bar.col.loss)
                                        y <- seq(-100,100,by=20)
                                        circlize::circos.segments(0,y,circlize::CELL_META$xlim[2], y,col=ifelse(y>0,col.gain,ifelse(y == 0,'#909090',col.loss)))
                                      })
        if (i <= 2 & length(id_name) <= 2){text(0,-0.1*(i-1),id,cex = cex_text)}
    }

}