R/bam_to_junctions.R

In megadepth: megadepth: BigWig and BAM related utilities

#' Extract junctions from a BAM file
#'
#' Given a BAM file, extract junction information including co-ordinates,
#' strand, anchor length for each junction read. For details on the format of
#' the output TSV file, check
#' <https://github.com/ChristopherWilks/megadepth#junctions>.
#'
#' @inheritParams bam_to_bigwig
#'
#' @param prefix A `character(1)` specifying the output file prefix. This
#'   function creates a file called `prefix.jxs.tsv`. 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 all_junctions A `logical(1)` indicating whether to obtain all
#'   junctions.
#' @param junctions A `logical(1)` indicating whether to obtain co-occurring jx
#'   coordinates.
#' @param long_reads A `logical(1)` indicating whether to increase the buffer
#'   size to accommodate for long-read RNA-sequencing.
#'
#' @return A `character(1)` with the path to the output junction tsv file.
#'
#' @export
#'
#' @examples
#'
#' ## Install if necessary
#' install_megadepth()
#'
#' ## Find the example BAM file
#' example_bam <- system.file("tests", "test.bam",
#'     package = "megadepth", mustWork = TRUE
#' )
#'
#' ## Run bam_to_junctions()
#' example_jxs <- bam_to_junctions(example_bam, overwrite = TRUE)
#'
#' ## Path to the output file generated by bam_to_junctions()
#' example_jxs
bam_to_junctions <- function(bam_file,
    prefix = file.path(tempdir(), basename(bam_file)),
    all_junctions = TRUE,
    junctions = FALSE,
    long_reads = FALSE,
    overwrite = FALSE) {
    expected_ext <- c("all_jxs.tsv", "jxs.tsv")
    if (!overwrite) prefix_exists(prefix, expected_ext)

    megadepth_shell2(
        bam_file,
        list(
            "prefix" = prefix,
            "junctions" = junctions,
            "all-junctions" = all_junctions,
            "long-reads" = long_reads
        )
    )

    ## 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)
}