Wrapper functions to run Megadepth commands via
system2('megadepth', ...).
megadepth_cmd(...)
megadepth_shell(input = ".", ...)Arguments to be passed to system2('megadepth', ...), e.g.
annotation(path) is basically
megadepth_cmd(c('--annotation', path))
(i.e. run the command megadepth --annotation path).
A character(1) with the path to the input BAM, BigWig or
text file for Megadepth.
See base::system2() for the types of output you can generate.
A character() with the capture of the standard output stream
generated by Megadepth.
megadepth_cmd(): Run an arbitrary Megadepth command.
megadepth_shell(): Run an arbitrary Megadepth command.
megadepth_cmd() is based on blogdown::hugo_cmd() which is available at
https://github.com/rstudio/blogdown/blob/master/R/hugo.R.
megadepth_shell() is based on the shell_ls() example from cmdfun which is
available at https://snystrom.github.io/cmdfun/index.html.
## Install if necessary
install_megadepth()
#> The latest megadepth version is 1.2.0
#> This is not an interactive session, therefore megadepth has been installed temporarily to
#> /tmp/RtmpK16vxZ/megadepth
## Find version
## megadepth_shell() provides an interface more familiar to R users
megadepth_shell(version = TRUE)
#> [1] "megadepth 1.2.0"
## megadepth_cmd() requires using directly the command line syntax for
## Megadepth
megadepth_cmd("--version", stdout = TRUE)
#> [1] "megadepth 1.2.0"
## Compare the help files:
# megadepth_shell() captures the standard output and returns a character()
# megadepth_cmd() shows the standard output on the console
megadepth_shell("--help")
#> [1] "megadepth 1.2.0"
#> [2] ""
#> [3] "BAM and BigWig utility."
#> [4] ""
#> [5] "Usage:"
#> [6] " megadepth <bam|bw|-> [options]"
#> [7] ""
#> [8] "Options:"
#> [9] " -h --help Show this screen."
#> [10] " --version Show version."
#> [11] " --threads # of threads to do: BAM decompression OR compute sums over multiple BigWigs in parallel"
#> [12] " if the 2nd is intended then a TXT file listing the paths to the BigWigs to process in parallel"
#> [13] " should be passed in as the main input file instead of a single BigWig file (EXPERIMENTAL)."
#> [14] " --prefix String to use to prefix all output files."
#> [15] " --no-auc-stdout Force all AUC(s) to be written to <prefix>.auc.tsv rather than STDOUT"
#> [16] " --no-annotation-stdout Force summarized annotation regions to be written to <prefix>.annotation.tsv rather than STDOUT"
#> [17] " --no-coverage-stdout Force covered regions to be written to <prefix>.coverage.tsv rather than STDOUT"
#> [18] " --keep-order Output annotation coverage in the order chromosomes appear in the BAM/BigWig file"
#> [19] " The default is to output annotation coverage in the order chromosomes appear in the annotation BED file."
#> [20] " This is only applicable if --annotation is used for either BAM or BigWig input."
#> [21] ""
#> [22] "BigWig Input:"
#> [23] "Extract regions and their counts from a BigWig outputting BED format if a BigWig file is detected as input (exclusive of the other BAM modes):"
#> [24] " Extracts all reads from the passed in BigWig and output as BED format."
#> [25] " This will also report the AUC over the annotated regions to STDOUT."
#> [26] " If only the name of the BigWig file is passed in with no other args, it will *only* report total AUC to STDOUT."
#> [27] " --annotation <bed> Only output the regions in this BED applying the argument to --op to them."
#> [28] " --op <sum[default], mean, min, max> Statistic to run on the intervals provided by --annotation"
#> [29] " --sums-only Discard coordinates from output of summarized regions"
#> [30] " --distance (2200[default]) Number of base pairs between end of last annotation and start of new to consider in the same BigWig query window (a form of binning) for performance. This determines the number of times the BigWig index is queried."
#> [31] " --unsorted (off[default]) There's a performance improvement *if* BED file passed to --annotation is 1) sorted by sort -k1,1 -k2,2n (default is to assume sorted and check for unsorted positions, if unsorted positions are found, will fall back to slower version)"
#> [32] " --bwbuffer <1GB[default]> Size of buffer for reading BigWig files, critical to use a large value (~1GB) for remote BigWigs."
#> [33] " Default setting should be fine for most uses, but raise if very slow on a remote BigWig."
#> [34] ""
#> [35] ""
#> [36] "BAM Input:"
#> [37] "Extract basic junction information from the BAM, including co-occurrence"
#> [38] "If only the name of the BAM file is passed in with no other args, it will *only* report total AUC to STDOUT."
#> [39] " --fasta Path to the reference FASTA file if a CRAM file is passed as the input file (ignored otherwise)"
#> [40] " If not passed, references will be downloaded using the CRAM header."
#> [41] " --junctions Extract co-occurring jx coordinates, strand, and anchor length, per read"
#> [42] " writes to a TSV file <prefix>.jxs.tsv"
#> [43] " --all-junctions Extract all jx coordinates, strand, and anchor length, per read for any jx"
#> [44] " writes to a TSV file <prefix>.all_jxs.tsv"
#> [45] " --longreads Modifies certain buffer sizes to accommodate longer reads such as PB/Oxford."
#> [46] " --filter-in Integer bitmask, any bits of which alignments need to have to be kept (similar to samtools view -f)."
#> [47] " --filter-out Integer bitmask, any bits of which alignments need to have to be skipped (similar to samtools view -F)."
#> [48] " --add-chr-prefix Adds \"chr\" prefix to relevant chromosomes for BAMs w/o it, pass \"human\" or \"mouse\"."
#> [49] " Only works for human/mouse references (default: off)."
#> [50] ""
#> [51] "Non-reference summaries:"
#> [52] " --alts Print differing from ref per-base coverages"
#> [53] " Writes to a CSV file <prefix>.alts.tsv"
#> [54] " --include-softclip Print a record to the alts CSV for soft-clipped bases"
#> [55] " Writes total counts to a separate TSV file <prefix>.softclip.tsv"
#> [56] " --only-polya If --include-softclip, only print softclips which are mostly A's or T's"
#> [57] " --include-n Print mismatch records when mismatched read base is N"
#> [58] " --print-qual Print quality values for mismatched bases"
#> [59] " --delta Print POS field as +/- delta from previous"
#> [60] " --require-mdz Quit with error unless MD:Z field exists everywhere it's"
#> [61] " expected"
#> [62] " --head Print sequence names and lengths in SAM/BAM header"
#> [63] ""
#> [64] "Coverage and quantification:"
#> [65] " --coverage Print per-base coverage (slow but totally worth it)"
#> [66] " --auc Print per-base area-under-coverage, will generate it for the genome"
#> [67] " and for the annotation if --annotation is also passed in"
#> [68] " Defaults to STDOUT, unless other params are passed in as well, then"
#> [69] " if writes to a TSV file <prefix>.auc.tsv"
#> [70] " --bigwig Output coverage as BigWig file(s). Writes to <prefix>.bw"
#> [71] " (also <prefix>.unique.bw when --min-unique-qual is specified)."
#> [72] " Requires libBigWig."
#> [73] " --annotation <BED|window_size> Path to BED file containing list of regions to sum coverage over"
#> [74] " (tab-delimited: chrm,start,end). Or this can specify a contiguous region size in bp."
#> [75] " --op <sum[default], mean> Statistic to run on the intervals provided by --annotation"
#> [76] " --no-index If using --annotation, skip the use of the BAM index (BAI) for pulling out regions."
#> [77] " Setting this can be faster if doing windows across the whole genome."
#> [78] " This will be turned on automatically if a window size is passed to --annotation."
#> [79] " --min-unique-qual <int>"
#> [80] " Output second bigWig consisting built only from alignments"
#> [81] " with at least this mapping quality. --bigwig must be specified."
#> [82] " Also produces second set of annotation sums based on this coverage"
#> [83] " if --annotation is enabled"
#> [84] " --double-count Allow overlapping ends of PE read to count twice toward"
#> [85] " coverage"
#> [86] " --num-bases Report total sum of bases in alignments processed (that pass filters)"
#> [87] " --gzip Turns on gzipping of coverage output (no effect if --bigwig is passsed),"
#> [88] " this will also enable --no-coverage-stdout."
#> [89] ""
#> [90] "Other outputs:"
#> [91] " --read-ends Print counts of read starts/ends, if --min-unique-qual is set"
#> [92] " then only the alignments that pass that filter will be counted here"
#> [93] " Writes to 2 TSV files: <prefix>.starts.tsv, <prefix>.ends.tsv"
#> [94] " --frag-dist Print fragment length distribution across the genome"
#> [95] " Writes to a TSV file <prefix>.frags.tsv"
#> [96] " --echo-sam Print a SAM record for each aligned read"
#> [97] " --ends Report end coordinate for each read (useful for debugging)"
#> [98] " --test-polya Lower Poly-A filter minimums for testing (only useful for debugging/testing)"
#> [99] ""
#> [100] ""
megadepth_cmd("--help")