R/window.R

In cubar: Codon Usage Bias Analysis

#' slide window interval generator
#'
#' \code{slide} generates a data.table with start, center, and end columns
#' for a sliding window analysis.
#'
#' @param from integer, the start of the sequence
#' @param to integer, the end of the sequence
#' @param step integer, the step size
#' @param before integer, the number of values before the center of a window
#' @param after integer, the number of values after the center of a window
#' @return data.table with start, center, and end columns
#' @export
#' @examples
#' slide(1, 10, step = 2, before = 1, after = 1)
#'
slide <- function(from, to, step = 1, before = 0, after = 0){
    if(to - from + 1 < (1 + before + after)){
        stop('The window size is too large')
    }
    center <- seq(from + before, to - after, by = step)
    start <- center - before
    end <- center + after
    data.table::data.table(start = start, center = center, end = end)
}


#' sliding window of codons
#'
#' \code{slide_codon} generates a data.table with start, center, and end columns
#'   for a sliding window analysis of codons.
#'
#' @param seq DNAString, the sequence
#' @param step integer, the step size
#' @param before integer, the number of codons before the center of a window
#' @param after integer, the number of codons after the center of a window
#' @return data.table with start, center, and end columns
#' @export
#' @examples
#' x <- Biostrings::DNAString('ATCTACATAGCTACGTAGCTCGATGCTAGCATGCATCGTACGATCGTCGATCGTAG')
#' slide_codon(x, step = 3, before = 1, after = 1)
#'
slide_codon <- function(seq, step = 1, before = 0, after = 0){
    if(!inherits(seq, 'DNAString')){
        seq <- Biostrings::DNAString(seq)
    }
    slen <- (length(seq) %/% 3) * 3
    if(slen < 3 * (1 + before + after)){
        stop('The window size is too large')
    }
    slide(from = 1, to = slen, step = step * 3,
          before = before * 3, after = 2 + after * 3)
}


#' apply a cub index to a sliding window
#'
#' \code{slide_apply} applies a function to a sliding window of codons.
#'
#' @param seq DNAString, the sequence
#' @param .f function, the codon index calculation function to apply, for
#'   example, \code{get_enc}.
#' @param step integer, the step size in number of codons
#' @param before integer, the number of codons before the center of a window
#' @param after integer, the number of codons after the center of a window
#' @param ... additional arguments to pass to the function \code{.f}
#' @return data.table with start, center, end, and codon usage index columns
#' @importFrom data.table ':='
#' @export
#' @examples
#' slide_apply(yeast_cds[[1]], get_enc, step = 1, before = 10, after = 10)
#'
slide_apply <- function(seq, .f, step = 1, before = 0, after = 0, ...){
    index <- NULL
    if(!inherits(seq, 'DNAString')){
        seq <- Biostrings::DNAString(seq)
    }
    windt <- slide_codon(seq, step = step, before = before, after = after)
    winseq <- Biostrings::Views(seq, start = windt$start, end = windt$end)
    winseq <- Biostrings::DNAStringSet(winseq)
    windt[, index := .f(count_codons(winseq), ...)]
    windt[]
}


#' plot sliding window codon usage
#'
#' \code{slide_plot} visualizes codon usage in sliding window.
#'
#' @param windt data.table, the sliding window codon usage
#'   generated by \code{slide_apply}.
#' @param index_name character, the name of the index to display.
#' @return ggplot2 plot.
#' @export
#' @examples
#' sw <- slide_apply(yeast_cds[[1]], get_enc, step = 1, before = 10, after = 10)
#' slide_plot(sw)
slide_plot <- function(windt, index_name='Index'){
    center <- index <- NULL
    ggplot2::ggplot(windt, ggplot2::aes(x = ceiling(center/3), y = index)) +
        ggplot2::geom_line() +
        ggplot2::geom_point() +
        ggplot2::geom_smooth(method = 'loess', formula = y ~ x) +
        ggplot2::geom_hline(yintercept = stats::median(windt$index), linetype = 2, color = 'red') +
        ggplot2::labs(x = 'Codon position', y = index_name) +
        ggplot2::theme_classic(base_size = 12) +
        ggplot2::theme(axis.text = ggplot2::element_text(color = 'black'))
}