Use an existing database

You can download pretrained databases from NCBI at https://ftp.ncbi.nlm.nih.gov/blast/db/. Here we download the 16S rRNA database. To avoid multiple downloads of the same file, we use BiocFileCache in function blast_db_cache for download. The database compressed and needs to be expanded.

## download the 16S Microbial rRNA data base from NCBI
tgz_file <- blast_db_get("16S_ribosomal_RNA.tar.gz")
untar(tgz_file, exdir = "16S_rRNA_DB")

The extracted database consists of a folder with a set of database files.

list.files("./16S_rRNA_DB/")

Next, we open the database for querying using the blast() function which returns a BLAST database object.

bl <- blast(db = "./16S_rRNA_DB/16S_ribosomal_RNA")
bl

To demonstrate how to query the database, we read one sequence from an example FASTA file that is shipped with the package. Queries are performed using the predict() function. The result is a data.frame with one row per match. We will show the first 5 matches.

seq <- readRNAStringSet(system.file("examples/RNA_example.fasta",
    package = "rBLAST"
))[1]
seq

cl <- predict(bl, seq)
nrow(cl)
cl[1:5, ]

Additional arguments for BLAST can be passed on using the BLAST_args parameter for predict(). The output format can be specified using custom_format. In the following, we specify a custom format and that the sequences need to have 99% identity. See the BLAST Command Line Applications User Manual for details (https://www.ncbi.nlm.nih.gov/books/NBK279690/).

fmt <- paste(
    "qaccver saccver pident length mismatch gapopen qstart qend",
    "sstart send evalue bitscore qseq sseq"
)
cl <- predict(bl, seq,
    BLAST_args = "-perc_identity 99",
    custom_format = fmt
)
cl

Create a custom BLAST database

The makeblastdb utility can be used to create a BLAST database from a FASTA file with sequences. The package provides an R interface function of the same name. As an example, we will create a searchable database from a sequence FASTA file shipped with the package.

seq <- readRNAStringSet(system.file("examples/RNA_example.fasta",
    package = "rBLAST"
))
seq
#> RNAStringSet object of length 5:
#>     width seq                                               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, we write a FASTA file that will be used by blast to create the database.

writeXStringSet(seq, filepath = "seqs.fasta")

Next, we create a BLAST database from the file. We need to specify that the sequences contains RNA and thus nucleotides.

makeblastdb("seqs.fasta", db_name = "db/small", dbtype = "nucl")

Note that it is convenient to specify a folder and a name separated by a / to organize all index files in a folder.

We can now open the database and use it for queries.

db <- blast("db/small")
db

Check if it finds a 100 nucleotide long fragment from the first sequence in the training data.

fragment <- subseq(seq[1], start = 101, end = 200)
fragment

predict(db, fragment)

We see that the found sequence ID (ssequid) matches the query sequence ID (qsequid) and that it matches the correct region in the sequence (see sstart and send).

To permanently remove a database, the folder can be deleted.

unlink("seqs.fasta")
unlink("db", recursive = TRUE)