rBLAST: R Interface for the Basic Local Alignment Search Tool

Michael Hahsler and Anurag Nagar

library("rBLAST")
library("Biostrings")

Overview

rBLAST connects Bioconductor sequence containers to the NCBI BLAST+ command-line programs. A typical local search has three steps:

  1. Create a database from FASTA sequences with makeblastdb().
  2. Open the database with blast().
  3. Search it with predict() using a DNAStringSet, RNAStringSet, or AAStringSet query.

predict() returns a data frame with one row per BLAST hit. rBLAST can also open downloaded NCBI databases and submit remote searches to NCBI.

System requirements

BLAST+ is external software and must be installed separately. See the package installation instructions and the BLAST+ installation guide.

For the nucleotide workflow below, R must find all three programs:

required_tools <- c("blastn", "makeblastdb", "blastdbcmd")
Sys.which(required_tools)
#>      blastn makeblastdb  blastdbcmd 
#>          ""          ""          ""

Note: A complete local BLAST+ installation was not available when this vignette was built. The code is shown, but output from BLAST-dependent chunks is omitted.

An empty path means that the corresponding executable is unavailable. If BLAST+ is installed outside the normal search path, update PATH before using rBLAST:

Sys.setenv(
    PATH = paste(
        Sys.getenv("PATH"),
        "/path/to/ncbi-blast/bin",
        sep = .Platform$path.sep
    )
)

Use blast_help("blastn") or run system2("blastn", "-version") to inspect the installed command-line program.

Quick start: build and query a local database

The following example is self-contained. It uses five RNA sequences bundled with rBLAST and writes all generated files to a unique temporary directory.

queries <- readRNAStringSet(system.file(
    "examples/RNA_example.fasta",
    package = "rBLAST"
))
queries
#> RNAStringSet object of length 5:
#>     width sequence                                          names               
#> [1]  1481 AGAGUUUGAUCCUGGCUCAGAAC...GGUGAAGUCGUAACAAGGUAACC 1675 AB015560.1 d...
#> [2]  1404 GCUGGCGGCAGGCCUAACACAUG...CACGGUAAGGUCAGCGACUGGGG 4399 D14432.1 Rho...
#> [3]  1426 GGAAUGCUNAACACAUGCAAGUC...AACAAGGUAGCCGUAGGGGAACC 4403 X72908.1 Ros...
#> [4]  1362 GCUGGCGGAAUGCUUAACACAUG...UACCUUAGGUGUCUAGGCUAACC 4404 AF173825.1 A...
#> [5]  1458 AGAGUUUGAUUAUGGCUCAGAGC...UGAAGUCGUAACAAGGUAACCGU 4411 Y07647.2 Dre...

First, write the sequences to FASTA and create a nucleotide database. RNA and DNA sequences both use dbtype = "nucl".

work_dir <- tempfile("rBLAST-vignette-")
dir.create(work_dir)
fasta <- file.path(work_dir, "sequences.fasta")
db_path <- file.path(work_dir, "database")

writeXStringSet(queries, fasta)
makeblastdb(
    fasta,
    db_name = db_path,
    dbtype = "nucl",
    verbose = FALSE
)

db_path is the common prefix of the database files, not the name of one individual index file. Open that prefix with blast():

db <- blast(db_path, type = "blastn")
print(db, info = FALSE)

Now search for a 100-nucleotide fragment taken from the first database sequence:

fragment <- subseq(queries[1], start = 101, end = 200)
hits <- predict(db, fragment)
hits

The most commonly used result columns are:

Column Meaning
qseqid, sseqid Query and database-sequence identifiers
pident Percentage of identical positions in the alignment
length Alignment length
qstart, qend Alignment coordinates in the query
sstart, send Alignment coordinates in the database sequence
evalue Expected number of chance matches; smaller is stronger
bitscore Normalized alignment score; larger is stronger

In this example, sseqid identifies the sequence from which the fragment was taken, and sstart and send locate the fragment within that sequence.

Multiple queries and custom output

An XStringSet can contain one or many query sequences. BLAST options are passed in command-line form through BLAST_args. Output fields are selected with space-separated BLAST field specifiers in custom_format:

selected_hits <- predict(
    db,
    queries[1:2],
    BLAST_args = "-perc_identity 99",
    custom_format = paste(
        "qseqid sseqid pident length",
        "qstart qend sstart send evalue bitscore"
    )
)
selected_hits

A successful search without matches returns a zero-row data frame. A failed BLAST command instead produces an R error containing the executable’s exit status and the query and output paths.

Threads and parallel searches

Use BLAST’s -num_threads option to parallelize one local search:

hits <- predict(
    db,
    queries,
    BLAST_args = "-num_threads 8"
)

This differs from running several predict() calls concurrently. Each concurrent call creates its own query and output files and starts a separate BLAST process. Account for both levels of concurrency when choosing worker and thread counts; for example, 8 workers with 8 BLAST threads can request up to 64 CPU threads.

Temporary files and troubleshooting

By default, rBLAST removes a prediction’s temporary FASTA and output files. Retain them for one diagnostic run with:

hits <- predict(
    db,
    queries[1],
    verbose = TRUE,
    keep_tmp = TRUE
)

list.files(tempdir())

verbose = TRUE prints the generated filenames and BLAST command. keep_tmp = TRUE leaves the files in tempdir() until they are manually removed or the R session ends. Repeated or parallel searches can retain large amounts of data, potentially exhausting temporary storage or a user quota. Use keep_tmp = FALSE for routine work.

Common failures can be diagnosed as follows:

On systems with a small /tmp, select a larger temporary filesystem before starting R:

TMPDIR=/path/with/space Rscript analysis.R

Prebuilt NCBI databases

blast_db_get() downloads and caches archives from the NCBI BLAST database directory. Because these archives can be large and require network access, this example is not executed while building the vignette:

archive <- blast_db_get("16S_ribosomal_RNA.tar.gz")
database_dir <- "/path/to/16S_rRNA_DB"
dir.create(database_dir, recursive = TRUE, showWarnings = FALSE)
untar(archive, exdir = database_dir)

ncbi_db <- blast(file.path(database_dir, "16S_ribosomal_RNA"))
hits <- predict(ncbi_db, queries[1])

Large databases may be divided into numbered archives. Download and extract every part into the same directory. See the BLAST command-line manual for database layouts and search options.

Remote searches

A BLAST object can refer to an NCBI database without downloading it:

remote_db <- blast("nt", remote = TRUE, type = "blastn")
remote_hits <- predict(remote_db, queries[1])

Remote searches use a shared NCBI service, may be queued, and are usually much slower than local searches. They also require network access. Prefer a local database for large or repeated analyses.

Session information

sessionInfo()
#> R version 4.6.1 Patched (2026-06-24 r90190)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.7.8
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] rBLAST_1.9.2         Biostrings_2.81.9    Seqinfo_1.3.2       
#> [4] XVector_0.53.0       IRanges_2.47.5       S4Vectors_0.51.9    
#> [7] BiocGenerics_0.59.12 generics_0.1.4      
#> 
#> loaded via a namespace (and not attached):
#>  [1] crayon_1.5.3    cli_3.6.6       knitr_1.52      rlang_1.3.0    
#>  [5] xfun_0.60       otel_0.2.0      jsonlite_2.0.0  htmltools_0.5.9
#>  [9] sass_0.4.10     rmarkdown_2.32  evaluate_1.0.5  jquerylib_0.1.4
#> [13] fastmap_1.2.0   yaml_2.3.12     lifecycle_1.0.5 compiler_4.6.1 
#> [17] digest_0.6.39   R6_2.6.1        bslib_0.12.0    tools_4.6.1    
#> [21] cachem_1.1.0