getMappableRegions: Get mappable regions of a genome.

In fmicompbio/swissknife: Handy code shared in the FMI CompBio group

Get mappable regions of a genome.

Description

Given a k-mer length and the maximum number of allowed hits per k-mer, find all mappable regions in a genome.

Usage

getMappableRegions(
  genome,
  genomeIndex,
  kmerLength = 50,
  maxHits = 1,
  Ncpu = 2,
  quiet = TRUE
)

Arguments

genome	The genome sequence to work on. Either a BSgenome object, a character scalar with the name of an installed BSgenome or with a file path and name pointing to a fasta file with the genome sequence.
genomeIndex	character scalar with the path to the bowtie index and prefix to align against, in the form </path/to/index>/<prefix>, or the name of an installed Rbowtie index package created by the QuasR package for an installed BSgenome package.
kmerLength	numeric scalar specifying the k-mer length (width of overlapping windows in genome), usually set to the typical read length for which to get the mappable regions.
maxHits	numeric scalar specifying the maximum number of hits (matches) of a k-mer in the genome to be considered mappable.
Ncpu	numeric scalar specifying the number of CPU threads to use for alignments.
quiet	logical scalar indicating if progress information should be printed on the console.

Details

Sequences of all overlapping windows are extracted from the genome and aligned to the provided genome index using bowtie with parameters -f -v 0 -a -B 1 -m maxHits. If no more than maxHits hits are found, the window is defined mappable.

Value

A GRanges object with mappable regions. All plus-strand sequences in genome of length kmerLength with their start (leftmost) position overlapping the GRanges object do not generate more than maxHits hits when aligned to the genome.

Author(s)

Michael Stadler

See Also

bowtie in package Rbowtie used by getMappableRegions to align reads to the genome; bowtie_build in package Rbowtie for indexing a genome.

Examples

if (requireNamespace("Rbowtie", quietly = TRUE)) {
    library(Rbowtie)

    genomefile <- system.file("extdata", "getMappableRegions", "hg19sub.fa", package = "swissknife")
    indexdir <- tempfile()
    indexpre <- "index"
    indexname <- file.path(indexdir, indexpre)
    idx <- bowtie_build(genomefile, indexdir)

    mapgr <- getMappableRegions(genomefile, indexname, 50, quiet = FALSE)
    print(mapgr)
}

fmicompbio/swissknife documentation built on June 11, 2025, 4:17 p.m.