#' 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.
#' @param filter_in A `integer(1)` used to filter in read alignments. See
#' https://github.com/ChristopherWilks/megadepth#bamcram-processing and
#' https://samtools.github.io/hts-specs/SAMv1.pdf for further documentation on
#' how to apply this parameter.
#' @param filter_out A `integer(1)` used to filter out read alignments. See
#' https://github.com/ChristopherWilks/megadepth#bamcram-processing and
#' https://samtools.github.io/hts-specs/SAMv1.pdf for further documentation on
#' how to apply this parameter.
#'
#' @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,
filter_in = 65535,
filter_out = 260,
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,
"filter-in" = filter_in,
"filter-out" = filter_out
)
)
## 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.