R/bam_to_bigwig.R

In LieberInstitute/megadepth: megadepth: BigWig and BAM related utilities

#' Convert a BAM file to a BigWig
#'
#' Given an input BAM file, convert this to the BigWig format which can then be
#' used in `get_coverage()`.
#'
#' Note that this functionality is currently not supported on Windows by
#' Megadepth.
#'
#' @param bam_file A `character(1)` with the path to the input BAM file.
#' @param prefix A `character(1)` specifying the output file prefix. This
#'   function creates a BigWig file called `{prefix}.all.bw`. By default, the
#'   prefix is the BAM file name and the file is created in the `tempdir()` and
#'   will be deleted after you close your R session.
#' @param min_unique_qual A `integer(1)` specifying a mapping quality threshold
#'   and only bases above this will be used to generate the BigWig. If set to
#'   `FALSE` this argument is not used by Megadepth. Otherwise it will generate
#'   files `{prefix}.unique.bw` and `{prefix}.unique.tsv`.
#' @param double_count A `logical(1)` determining whether to count the
#'   overlapping ends of paired ends reads twice.
#' @param overwrite A `logical(1)` specifying whether to overwrite the output
#' file(s), if they exist already.
#'
#' @return A `character()` with the path to the output files(s).
#' @export
#'
#' @examples
#'
#' ## Install the latest version if necessary
#' install_megadepth(force = TRUE)
#'
#' ## Find the example BAM file
#' example_bam <- system.file("tests", "test.bam",
#'     package = "megadepth", mustWork = TRUE
#' )
#'
#' ## Create the BigWig file
#' ## Currently Megadepth does not support this on Windows
#' if (!xfun::is_windows()) {
#'     example_bw <- bam_to_bigwig(example_bam, overwrite = TRUE)
#'
#'     ## Path to the output file(s) generated by bam_to_bigwig()
#'     example_bw
#'
#'     ## Use the all.bw file in get_coverage(), first find an annotation file
#'     annotation_file <- system.file("tests", "testbw2.bed",
#'         package = "megadepth", mustWork = TRUE
#'     )
#'
#'     ## Compute the coverage
#'     bw_cov <- get_coverage(
#'         example_bw["all.bw"],
#'         op = "mean",
#'         annotation = annotation_file
#'     )
#'     bw_cov
#' }
bam_to_bigwig <- function(
        bam_file,
        prefix = file.path(tempdir(), basename(bam_file)),
        min_unique_qual = FALSE,
        double_count = FALSE,
        overwrite = FALSE) {
    if (xfun::is_windows()) {
        warning(
            "Megadepth currently does not support creating bigWig files on ",
            "Windows. This function will still try to run the code, as a ",
            "newer version of Megadepth might support this feature. Please ",
            "let us know if code works for you on Windows! Thanks!",
            call. = FALSE
        )
    }

    ## Hm.. the unique.tsv file is empty with the current example data
    ## even with min_unique_qual = 0
    expected_ext <- c("all.bw", "unique.bw", "unique.tsv")
    if (!overwrite) prefix_exists(prefix, expected_ext)

    megadepth_shell2(
        bam_file,
        list(
            "prefix" = prefix,
            "bigwig" = TRUE,
            "min-unique-qual" = min_unique_qual,
            "double-count" = double_count
        )
    )

    ## Find the prefix files
    new_files <- prefix_files(prefix)
    ## List only those that match the expected file names
    new_files <- new_files[names(new_files) %in% expected_ext]

    return(new_files)
}