Getting Started with BiocRefgetStore

BiocRefgetStore provides a BSgenome-compatible interface to reference genomes backed by GA4GH refget stores. Instead of managing FASTA files, you connect to a refget store (local or remote) and access sequences by digest.

This tutorial uses the 2023 Human Pangenome Reference -- a remote refget store containing 47 haplotype-resolved assemblies hosted on S3. Sequences are downloaded on-demand and cached locally, so you don't need to download the entire dataset upfront.

Installation

BiocRefgetStore depends on gtars (>= 0.10.0), the Rust-backed sequence-store backend. gtars is not on CRAN or Bioconductor, so install it first, then install BiocRefgetStore.

1. Install gtars (prerequisite)

On Linux x86_64, the easiest route is the prebuilt R binary attached to the gtars release -- no Rust toolchain required:

install.packages(c("data.table", "BiocManager"))
BiocManager::install(c("BiocGenerics", "IRanges"))

# Install the prebuilt gtars R binary directly from the release asset:
install.packages(
  "https://github.com/databio/gtars/releases/download/gtars-r-v0.10.0/gtars_0.10.0_R_x86_64-pc-linux-gnu.tar.gz",
  repos = NULL
)

macOS (Apple Silicon / arm64). A prebuilt binary is also published for Apple Silicon Macs:

install.packages(c("data.table", "BiocManager"))
BiocManager::install(c("BiocGenerics", "IRanges"))

install.packages(
  "https://github.com/databio/gtars/releases/download/gtars-r-v0.10.0/gtars_0.10.0_R_arm64-apple-darwin.tgz",
  repos = NULL
)

Prebuilt binaries are built for R 4.4 and won't load under a different R minor version. Check the releases page (tag gtars-r-v0.10.0) for the exact asset names if a URL 404s.

Intel macOS / other platforms / other R versions. No prebuilt binary is published for these, so build from source. gtars lives in a Rust monorepo and its R package depends on sibling crates by relative path, so remotes::install_github() does not work -- you need a full clone plus a Rust toolchain (cargo/rustc):

# git clone https://github.com/databio/gtars.git   # full monorepo, not a subdir
install.packages("path/to/gtars/gtars-r", repos = NULL, type = "source")

2. Install BiocRefgetStore

remotes::install_github("databio/BiocRefgetStore")
# Or from a local clone:
# install.packages("path/to/BiocRefgetStore", repos = NULL, type = "source")

Connect to a remote pangenome store

Load a genome from the Human Pangenome Reference store. Nothing is downloaded upfront: opening fetches only the store/collection metadata. Sequence data is then pulled on demand when you ask for it.

How remote retrieval works: two modes

Against a remote store, sequence access uses one of two on-demand modes depending on whether you ask for a region or a whole sequence:

• Region reads (lean byte-range). getSeq(genome, name, start, end), extractRegions(), and the GRanges/vectorized forms pull only the bytes covering the requested region, over an HTTP byte-range request. The whole sequence never lands on disk, and these reads persist nothing by default -- each call re-fetches just the bytes it needs. This is the right default for sparse, random access: you can query a genome-scale collection while transferring only the handful of bases you touch.

• Whole-sequence reads (streaming). Coordinate-less access -- genome[["chr"]] or getSeq(genome, name) with no start/end -- streams the entire sequence on demand with constant peak memory, then returns it. Reach for this only when you genuinely need a full chromosome; for random-access extraction, prefer region reads so you don't transfer bases you won't use.

Local stores (RefgetGenome.from_fasta() or an on-disk store) serve both modes straight from disk/memory, so the same code works unchanged.

library(BiocRefgetStore)

# 2023 Human Pangenome Reference (47 haplotype-resolved assemblies)
pangenome_url <- "https://refgenie.s3.us-east-1.amazonaws.com/pangenome_refget_store"

# One assembly from the pangenome (481 sequences)
genome <- RefgetGenome.from_remote(
  cache_path = "~/.cache/refget/pangenome",
  remote_url = pangenome_url,
  digest = "0qveCdMlbF_kYn6XWb7YBy-FtRZ6gSAL"
)
genome
#> RefgetGenome with 481 sequences
#>   collection_digest: 0qveCdMlbF_kYn6XWb7YBy-FtRZ6gSAL
#>   seqnames: JAHEPG010000017.1, JAHEPG010000014.1, JAHEPG010000009.1, JAHEPG010000029.1, JAHEPG010000005.1 ... (476 more)

Use the names reported by seqnames(genome) for all access below; this vignette uses the first one, JAHEPG010000017.1, as the example.

The data-access chunks below are kept eval=FALSE only because they require network access (and the recent gtars noted above); they are not run when this vignette is knit. Run them in a live session to fetch real sequence from the remote store.

Basic sequence access

Extract a full sequence or a region by coordinates:

# Full sequence (returns a DNAString)
seq <- genome[["JAHEPG010000017.1"]]

# Region by coordinates (1-based, inclusive)
region <- getSeq(genome, "JAHEPG010000017.1", start = 1000, end = 2000)

# Force character output
region_chr <- getSeq(genome, "JAHEPG010000017.1", start = 1000, end = 2000,
                     as.character = TRUE)

Negative-strand extraction applies the reverse complement:

rc <- getSeq(genome, "JAHEPG010000017.1", start = 1000, end = 2000, strand = "-")

Multiple regions at once

Pass vectors of names, starts, and ends:

seqs <- getSeq(
  genome,
  names = c("JAHEPG010000017.1", "JAHEPG010000014.1", "JAHEPG010000009.1"),
  start = c(100, 200, 300),
  end   = c(199, 299, 399)
)

If you have a GRanges object, pass it directly:

library(GenomicRanges)
gr <- GRanges(c("JAHEPG010000017.1:100-199:+", "JAHEPG010000014.1:200-299:-"))
seqs <- getSeq(genome, gr)

Bulk extraction from a data.frame

extractRegions accepts a data.frame with chrom, start, end columns:

regions <- data.frame(
  chrom = c("JAHEPG010000017.1", "JAHEPG010000017.1", "JAHEPG010000014.1"),
  start = c(100, 5000, 200),
  end   = c(199, 5099, 299)
)
seqs <- extractRegions(genome, regions, as.character = TRUE)

Export sequences to FASTA

Write extracted regions or full sequences to a FASTA file:

# Regions to FASTA
extractToFasta(genome, regions, "extracted_regions.fa")

# Specific sequences
exportChromosomes(genome, c("JAHEPG010000017.1", "JAHEPG010000014.1"), "subset.fa")

# All sequences
exportChromosomes(genome, output_path = "full_assembly.fa")

Genome metadata

Inspect sequences and their properties:

seqnames(genome)       # sequence names
seqlengths(genome)     # named integer vector of lengths
seqinfo(genome)        # full Seqinfo object
length(genome)         # number of sequences
collection_digest(genome)   # seqcol digest
coordinate_system(genome)   # sorted_name_length_pairs digest
sequence_digests(genome)    # per-sequence SHA512t24u digests

Working with local FASTA files

You can also create a genome directly from a local FASTA file. This builds an in-memory refget store, computes sequence digests, and creates a Seqinfo object automatically:

genome <- RefgetGenome.from_fasta("genome.fa")
genome

Persistent on-disk store

For large genomes you access repeatedly, use an on-disk store so sequences are indexed once and reused across sessions:

# First time: create the store from FASTA
store <- gtars::refget_store_on_disk("~/.local/share/refget/hg38")
result <- gtars::add_fasta(store, "hg38.fa")
# Save the digest: result$digest

# Later: reload without re-parsing
genome <- RefgetGenome.from_directory(
  "~/.local/share/refget/hg38",
  digest = "saved_digest_string"
)

The same workflow in Python

The remote-only workflow shown here mirrors the Python refget package. For a side-by-side reference -- open by URL + cache, list sequences, fetch by collection + name, batch-extract from a BED file, and export to FASTA -- see the Python demo at data_loaders/demo_remote_store.py in the refget repository.

Next steps

See the Reference vignette for complete documentation of every function and method in the package.