---
title: "Exploring QuickBLAST"
author: "Vishvesh Karthik"
date: "`r Sys.Date()`"
urlcolor: blue
output:
  rmarkdown::html_vignette:
    highlight: espresso
    self_contained: true
    wrap: 60
    columns: 60
    toc: true
    toc_depth: 4
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Exploring QuickBLAST}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#> ",
  eval = arrow::arrow_with_parquet()
)
knitr::opts_chunk$set(error = TRUE)
options(try.outFile = stdout())
```

```{r, echo=FALSE, out.width="30%", fig.align="center"}
knitr::include_graphics("figures/logo.png")
```

# Introduction to QuickBLAST

R is widely used for data analysis, but running NCBI's standard BLAST tools within R has traditionally been slow. Because the NCBI C++ toolkit is massive and difficult to compile for R, existing packages are forced to run BLAST as an external subprocess, creating major read/write bottlenecks.

QuickBLAST solves this by building a direct bridge between R and the NCBI C++ toolkit via Rcpp. By bypassing traditional text-based formatting and transporting data directly into memory using Apache Arrow, QuickBLAST performs sequence comparisons exceptionally fast.

## Setup and Instance Management

QuickBLAST operates using "instances". An instance is an active C++ object pointer (QuickBLAST_XPtr) that holds your BLAST configuration open in memory.

```{r setup, echo=TRUE,include=TRUE, tidy=TRUE}
library(QuickBLAST)

QuickBLAST::isQuickBLASTLoaded()

# Create instances for different sequence types and programs.
# Note: seq_type = 1 (Protein), seq_type = 0 (Nucleotide)
blastp_inst <- QuickBLAST::CreateQuickBLASTInstance(seq_type = 1, strand = 0, program = "blastp", save_sequences = TRUE, options = "-evalue 100000")

blastn_inst <- QuickBLAST::CreateQuickBLASTInstance(seq_type = 0, strand = 0, program = "blastn", options = "-evalue 100000")

tblastx_inst <- QuickBLAST::CreateQuickBLASTInstance(seq_type = 0, strand = 0, program = "tblastx")

# Inspecting an instance reveals its underlying C++ pointer and attributes
blastn_inst
```

### Managing Memory

Because QuickBLAST creates objects in C++, we provide convenience functions to track, verify, and clean up instances to prevent memory leaks.

```{r inst_man, error = 0, echo=TRUE,include=TRUE, tidy=TRUE}
QuickBLAST::GetInstanceCount()

QuickBLAST::GetInstanceID(blastp_inst)

# Clean up memory by deleting specific instances
QuickBLAST::DeleteQuickBLASTInstance(QuickBLAST::GetQuickBLASTInstance(1))

QuickBLAST::DeleteQuickBLASTInstance(2)

QuickBLAST::GetInstanceCount()

# Attempting to fetch a deleted instance securely throws a C++ map error
# try(identical(try(QuickBLAST::GetQuickBLASTInstance(1), silent = TRUE, outFile = stdout()), blastn_inst), silent = TRUE, outFile = stdout())

# Recreate the blastn instance for later steps
blastn_inst <- QuickBLAST::CreateQuickBLASTInstance(seq_type = 0, strand = 0, program = "blastn")
```

## Raw Sequence Comparison

You can pass raw character strings directly to QuickBLAST. The engine evaluates the sequences in memory instantly. Notice that QuickBLAST includes automatic recovery: if a pointer is dead, it will attempt to reload it.

```{r blast2seqs, error = 0}
# (~)Successful TBLASTX execution : Garbage low information nucleotide sequences
try(QuickBLAST::BLAST2Seqs(tblastx_inst, "AAAAAAAAAAAATTTTTTTTTTTTGGGGGGGGGGGCCCCCCCCC", "TTTTTTTTTTTGGGGGGGGGGGG"), silent = TRUE, outFile = stdout())

# Expected empty results due to mismatched pointer type (Nucleotide vs Nucleotide with a blastp pointer)
try(QuickBLAST::BLAST2Seqs(blastp_inst, "AAAAAAAAAAAATTTTTTTTTTTTGGGGGGGGGGGCCCCCCCCC", "TTTTTTTTTTTGGGGGGGGGGGG"), silent = TRUE, outFile = stdout())

# (~)Successful Protein-Protein search : Same protein
QuickBLAST::BLAST2Seqs(blastp_inst, "MQILLVEDDNTLFQELKKELEQWDFNVAGIEDFGKVMDTFESFNPEIVILDVQLPKYDGFYWCRKMREVSNVPILFLSSRDNPMDQVMSMELGADDYMQKPFYTNVLIAKLQAIYRRVYEFTAEEKRTLTWQDAVVDLSKDSIQKGDDTIFLSKTEMIILEILITKKNQIVSRDTIITALWDDEAFVSDNTLTVNVNRLRKKLSEISMDSAIETKVGKGYMAHE", "MQILLVEDDNTLFQELKKELEQWDFNVAGIEDFGKVMDTFESFNPEIVILDVQLPKYDGFYWCRKMREVSNVPILFLSSRDNPMDQVMSMELGADDYMQKPFYTNVLIAKLQAIYRRVYEFTAEEKRTLTWQDAVVDLSKDSIQKGDDTIFLSKTEMIILEILITKKNQIVSRDTIITALWDDEAFVSDNTLTVNVNRLRKKLSEISMDSAIETKVGKGYMAHE")

# Successful NT-NT search
QuickBLAST::BLAST2Seqs(blastn_inst, "ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA", "ATGGAGGAGCCGCAGTCAGATCCTAGCCTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA")
```

## Remote NCBI Searching

If you do not want to compile databases locally, QuickBLAST can interface directly with NCBI's remote servers.

````markdown
```{r remoteblast, error = 0, eval = FALSE}
# Safely check if the machine has internet access before running
if (requireNamespace("curl", quietly = TRUE) && curl::has_internet()) {
  # Example of a remote timeout handling
  try(QuickBLAST::RemoteBLAST(blastn_inst, query_input = "AAAAAAAAAAAATTTTTTTTTTTTGGGGGGGGGGGCCCCCCCCC", database = "nr", input_type = 1, return_values = TRUE), silent = TRUE, outFile = stdout())

  # Successful Remote Protein Query against the PDB database
  QuickBLAST::RemoteBLAST(blastp_inst, query_input = "MQILLVEDDNTLFQELKKELEQWDFNVAGIEDFGKVMDTFESFNPEIVILDVQLPKYDGFYWCRKMREVSNVPILFLSSRDNPMDQVMSMELGADDYMQKPFYTNVLIAKLQAIYRRVYEFTAEEKRTLTWQDAVVDLSKDSIQKGDDTIFLSKTEMIILEILITKKNQIVSRDTIITALWDDEAFVSDNTLTVNVNRLRKKLSEISMDSAIETKVGKGYMAHE", database = "pdb", input_type = 1, return_values = TRUE)
} else {
  message("No internet connection available. Skipping remote BLAST examples.")
}
```
````

## File-to-File Comparisons

For large-scale genomics, sequences are often stored in FASTA files. QuickBLAST can compare two files directly via multi-threaded reading, without loading everything into the R global environment first.

```{r blast2files, error = 0}
safe_temp_dir <- normalizePath(tempdir(check = TRUE), mustWork = FALSE)

# Create temporary FASTA files
query_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)
subject_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)

# Human TP53 gene fragment (Query)
writeLines(c(
  ">Human_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), query_fasta)

# Slightly mutated version (Subject - e.g., an ortholog or variant)
writeLines(c(
  ">Variant_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCCTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), subject_fasta)

# Efficiently compare the files using the BLASTN instance
file_results <- QuickBLAST::BLAST2Files(
  blastn_inst,
  query = query_fasta,
  subject = subject_fasta,
  return_values = TRUE,
  num_threads = 1
)

print(file_results)
```

## Local Database Creation and Searching

For high-throughput homology searches against a fixed reference, compiling the sequences into an indexed BLAST database is much faster than raw file-to-file comparisons. QuickBLAST wraps the NCBI makeblastdb functionality directly into R.


```{r blastfile2db, error = 0}
safe_temp_dir <- normalizePath(tempdir(check = TRUE), mustWork = FALSE)
query_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)
subject_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)

# Human TP53 gene fragment (Query)
writeLines(c(
  ">Human_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), query_fasta)
# Mutant TP53 gene fragment (Subject)
writeLines(c(
  ">Variant_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCCTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), subject_fasta)

db_path <- file.path(safe_temp_dir, "My_Test_DB")
# db_path_sub <- file.path(safe_temp_dir, "My_Test_DB_sub")

# 1. Compile the local BLAST database
# db_type can be "nucl" or "prot"
# QuickBLAST::MakeBLASTDB(
#   ptr = blastn_inst,
#   input_file = query_fasta,
#   database_name = db_path,
#   parse_seqids = FALSE,
#   stdout_opt = "",
#   stderr_opt = FALSE
# )

try(QuickBLAST::MakeBLASTDB(
  ptr = blastn_inst,
  input_file = subject_fasta,
  database_name = db_path,
  parse_seqids = FALSE
), silent = TRUE, outFile = stdout())

# 2. Search a query file or query database against the custom local database
db_results <- try(QuickBLAST::BLAST2DBs(
  blastn_inst,
  query = query_fasta,
  subject = db_path,
  return_values = TRUE,
  num_threads = 2
), silent = TRUE, outFile = stdout())

print(db_results)
```

## High-Performance Output with Apache Arrow

By default, the return_values = TRUE flag seen above brings data directly into R as an Rcpp::List. However, when processing millions of alignments, memory can max out quickly.

To solve this, QuickBLAST writes results asynchronously using Apache Arrow formats (arrow::parquet, arrow::ipc and arrow::csv). This creates a highly compressed, blazingly fast columnar storage file, an inter-process-communicable file or a tabular CSV.


```{r outputformats, error = 0}
safe_temp_dir <- normalizePath(tempdir(check = TRUE), mustWork = FALSE)
query_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)
subject_fasta <- normalizePath(tempfile(fileext = ".fasta", tmpdir = safe_temp_dir), mustWork = FALSE)

# Human TP53 gene fragment (Query)
writeLines(c(
  ">Human_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), query_fasta)
# Mutant TP53 gene fragment (Subject)
writeLines(c(
  ">Variant_p53_fragment\nATGGAGGAGCCGCAGTCAGATCCTAGCCTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA"
), subject_fasta)

# Specify an output path
out_parquet <- normalizePath(tempfile(fileext = ".parquet"), mustWork = FALSE)
out_ipc <- normalizePath(tempfile(fileext = ".ipc"), mustWork = FALSE)
out_csv <- normalizePath(tempfile(fileext = ".csv"), mustWork = FALSE)

# Perform the search and write straight to an Arrow Parquet/IPC/CSV file
# We don't need it in RAM, we are writing to disk
QuickBLAST::BLAST2Files(
  blastn_inst,
  query = query_fasta,
  subject = subject_fasta,
  out_file = out_parquet,
  out_format = "parquet",
  return_values = FALSE,
  verbose = TRUE
)

QuickBLAST::BLAST2Files(
  blastn_inst,
  query = query_fasta,
  subject = subject_fasta,
  out_file = out_ipc,
  out_format = "ipc",
  return_values = FALSE
)

QuickBLAST::BLAST2Files(
  blastn_inst,
  query = query_fasta,
  subject = subject_fasta,
  out_file = out_csv,
  out_format = "csv",
  return_values = FALSE
)

# Load the hits back into an R data.frame natively via QuickBLAST
hits_df <- QuickBLAST::LoadBLASTHits(out_parquet, format = "parquet")
head(hits_df)

hits_df <- QuickBLAST::LoadBLASTHits(out_ipc, format = "ipc")
head(hits_df)

hits_df <- QuickBLAST::LoadBLASTHits(out_csv, format = "csv", sep = "\t", header = TRUE)
head(hits_df)
```

## Conclusion

QuickBLAST modernizes sequence alignment in R. By managing your instances appropriately, keeping data in native C++ structures, and offloading massive search results to Apache Arrow formats, you can run high-throughput bioinformatic pipelines natively without relying on slow Sys.Call() bash wrappers.

```{r sessionInfo}
sessionInfo()
```