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().

bam_to_bigwig(
  bam_file,
  prefix = file.path(tempdir(), basename(bam_file)),
  min_unique_qual = FALSE,
  double_count = FALSE,
  overwrite = FALSE
)

Arguments

bam_file

A character(1) with the path to the input BAM file.

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.

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.

double_count

A logical(1) determining whether to count the overlapping ends of paired ends reads twice.

overwrite

A logical(1) specifying whether to overwrite the output file(s), if they exist already.

Value

A character() with the path to the output files(s).

Details

Note that this functionality is currently not supported on Windows by Megadepth.

Examples

## Install the latest version if necessary
install_megadepth(force = TRUE)
#> The latest megadepth version is 1.2.0
#> This is not an interactive session, therefore megadepth has been installed temporarily to 
#> /tmp/RtmpmgolT5/megadepth

## 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
}
#> GRanges object with 4 ranges and 1 metadata column:
#>         seqnames          ranges strand |     score
#>            <Rle>       <IRanges>  <Rle> | <numeric>
#>   [1]      chr10            0-10      * |      0.00
#>   [2]      chr10 8756697-8756762      * |     15.85
#>   [3]      chr10 4359156-4359188      * |      3.00
#>   [4] GL000219.1   168500-168620      * |      1.26
#>   -------
#>   seqinfo: 2 sequences from an unspecified genome; no seqlengths