---
title: "NCBI BLAST" 
author: "Author: FirstName LastName"
date: "Last update: `r format(Sys.time(), '%d %B, %Y')`" 
output:
  BiocStyle::html_document:
    toc_float: true
    code_folding: show
  BiocStyle::pdf_document: default
package: systemPipeR
vignette: |
  %\VignetteEncoding{UTF-8}
  %\VignetteIndexEntry{WF: BLAST Template}
  %\VignetteEngine{knitr::rmarkdown}
fontsize: 14pt
bibliography: bibtex.bib
editor_options: 
  chunk_output_type: console
---


```{css, echo=FALSE}
pre code {
white-space: pre !important;
overflow-x: scroll !important;
word-break: keep-all !important;
word-wrap: initial !important;
}
```

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
options(width=60, max.print=1000)
knitr::opts_chunk$set(
    eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
    cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")), 
    tidy.opts=list(width.cutoff=60), tidy=TRUE)
```

```{r setup, echo=FALSE, message=FALSE, warning=FALSE, eval=FALSE}
suppressPackageStartupMessages({
    library(systemPipeR)
})
```

# About the template
This section provides general description and how to use this BLAST workflow. 
In the actual analysis report, this section is usually **removed**. 

This BLAST workflow template is based on the 
[BLAST](https://blast.ncbi.nlm.nih.gov/Blast.cgi) based R package [rBLAST](https://github.com/mhahsler/rBLAST).

- The BLAST software can be downloaded from [NCBI](https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/).
  Please make sure it can be run from command-line.
- rBLAST can be installed with `install.packages('rBLAST', repos = 'https://mhahsler.r-universe.dev')`.

This workflow does:
1. Validate the BLAST installation
2. BLAST input fasta file with a reference genome. 
3. BLAST input with a certain database
4. BLAST sequence with general databases to find out the source organism(s).

All are written in R (`Linewise`) steps, but _BLAST+_ must be installed.

![](results/plotwf_spblast.jpg)

# Introduction

Users want to provide here background information about the design of their
BLAST project.

This report describes the analysis of a BLAST project studying ...

## Experimental design

Typically, users want to specify here all information relevant for the
analysis of their BLAST study. This includes detailed descriptions of
files, experimental design, reference genome, gene annotations,
etc.

# Workflow environment

To create a Workflow within _`systemPipeR`_, we can start by defining an empty
container and checking the directory structure:

```{r create_workflow, message=FALSE, eval=FALSE}
library(systemPipeR)
sal <- SPRproject()
sal
```

To load the entire workflow with default settings, run

```{r load_workflow_default, eval=FALSE}
sal <- importWF(sal, "SPblast.Rmd")
```


## Load packages 

```{r load_packages, eval=FALSE, spr=TRUE}
cat(crayon::blue$bold("To use this workflow, following R packages are expected:\n"))
cat(c("'rBLAST", "readr\n"), sep = "', '")
###pre-end
appendStep(sal) <- LineWise(
    code = {
        library(systemPipeR)
        library(rBLAST)
    }, 
    step_name = "load_packages"
)
```

## Test BLAST install

Molecules can be loaded or downloaded. This example dataset has 100 molecules.

```{r test_blast, eval=FALSE, spr=TRUE}
# Here, the dataset is downloaded. If you already have the data locally, change URL to local path.
appendStep(sal) <- LineWise(
    code = {
        # If you have a modular system, use following line
        moduleload("ncbi-blast") 
        # If not, comment out line above you need to install BLAST and configure the PATH. 
        blast_check <- tryCMD("blastn", silent = TRUE)
        if(blast_check  == "error") stop("Check your BLAST installation path.")
    }, 
    step_name = "test_blast",
    dependency = "load_packages"
)
```

## Load query sequence

Load query sequence from a `fasta` file.

In this template, an example fasta is provided, with 10 sequences from Arabidopsis, 
Cholera, Human, Mouse, and COVID-19, 2 for each. 

```{r load_query, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        query <- readDNAStringSet('data/example.fasta')
    }, 
    step_name = "load_query",
    dependency = "test_blast"
)
```


## BLAST against reference genome

In this step, we are trying to BLAST the query sequences to a reference genome 
and see if this genome contains the whole or part of the sequences.

In this example, a minimized `tair10` genome is used. In the real analysis, 
please replace it with a full genome `fasta` file. 

```{r build_genome_db, eval=FALSE, spr=TRUE}

appendStep(sal) <- LineWise(
    code = {
        reference <- 'data/tair10.fasta'
        # this command prepare BLAST-able database of genome
        makeblastdb(reference, dbtype='nucl')
    }, 
    step_name = "build_genome_db",
    dependency = "load_query"
)
```

Next BLAST is performed. Since there are only 2 Arabidopsis sequences in the 
example `fasta`. Only these two sequences are expected to return statistically 
meaningful BLAST results.

```{r blast_genome, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        bl_tair10 <- blast(db = reference, type='blastn')
        cl_tair10 <- predict(bl_tair10, query)
        readr::write_csv(cl_tair10, "results/blast_tair10.csv")
    }, 
    step_name = "blast_genome",
    dependency = "build_genome_db"
)
```


## BLAST existing databases

There are plenty of databases on [NCBI](https://ftp.ncbi.nlm.nih.gov/blast/db/) that one could 
download and run BLAST on. Once the databases are downloaded, unzip all files 
into one directory. We need to provide the path to the database. 

In this example, we want to know if COVID-19 is a beta coronavirus. Then, we can 
use some COVID sequence to BLAST all other existing beta coronavirus sequences and 
find the similarity. This resource is downloadable from NCBI. All downloaded 
`Betacoronavirus.XX.tar.gz` files are unzipped to `/srv/projects/db/ncbi/preformatted/20220131/`. 
Please change the path according to your project. Then, we can BLAST the last two sequence against
the database.


```{r blast_db, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        bl_covid <- blast(
            db = '/srv/projects/db/ncbi/preformatted/20220131/Betacoronavirus', 
            type='blastn'
        )
        cl_covid <- predict(bl_covid, query[9:10])
        readr::write_csv(cl_covid, "results/blast_covid.csv")
    }, 
    step_name = "blast_db",
    dependency = "load_query"
)
```

## BLAST to general databases
Sometimes we do not know the origin of a sequence, for example, a sequence comes 
from a contaminated sample, and we want to know the source. In such cases, we 
would need to BLAST the sequence to a more generic database. The most generic 
nucleotide BLAST database is the `nt` database. 

This database is extremely big and requires giant RAM and CPU cores to run. 
Please do not run the following example unless your system admin has provided 
you such store space and computational power. A better way for average the user 
is to use the website https://blast.ncbi.nlm.nih.gov/blast/Blast.cgi/ . 
The engine over there is optimized and can quickly
search for the species information. 


```{r blast_nt, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        bl_nt <- blast(db = '/srv/projects/db/ncbi/preformatted/20220131/nt', type='blastn')
        cl_nt <- predict(bl_nt, query[5])
        readr::write_csv(cl_nt, "results/blast_nt.csv")
    }, 
    step_name = "blast_nt",
    dependency = "load_query",
    run_step = "optional"
)
```


## Workflow session

```{r wf_session, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        sessionInfo()
    }, 
    step_name = "wf_session", 
    dependency = "blast_db")
```

# Manage the workflow

To run the workflow, use `runWF` function. It executes all the steps store in 
the workflow container. The execution will be on a single machine without 
submitting to a queuing system of a computer cluster. 

```{r runWF, eval=FALSE}
sal <- runWF(sal, run_step = "mandatory") # remove `run_step` to run all steps to include optional steps
```

- To use complex workflow control options, such as parallelization, subsetting samples, selecting steps, read the [documents](https://systempipe.org/sp/spr/sp_run/step_run/) on our website. 
- Explore [other details of the workflow object](https://systempipe.org/sp/spr/sp_run/sal_explore/).
- Create [logs and reports](https://systempipe.org/sp/spr/sp_run/step_reports/).
- [Visualize the workflow](https://systempipe.org/sp/spr/sp_run/step_vis/).

# About the workflow
## Tools used 

To check command-line tools used in this workflow, use `listCmdTools`, and use `listCmdModules`
to check if you have a modular system.

The following code will print out tools required in your custom SPR project in the report. 
In case you are running the workflow for the first time and do not have a project yet, or you 
just want to browser this workflow, following code displays the tools required by default.
```{r list_tools}
if(file.exists(file.path(".SPRproject", "SYSargsList.yml"))) {
    local({
        sal <- systemPipeR::SPRproject(resume = TRUE)
        systemPipeR::listCmdTools(sal)
        systemPipeR::listCmdModules(sal)
    })
} else {
    cat(crayon::blue$bold("Tools and modules required by this workflow are:\n"))
    cat(c("bowtie2/2.4.5", "samtools/1.14", "macs2"), sep = "\n")
}
```


## Session Info
This is the session information for rendering this report. To access the session information
of workflow running, check HTML report of `renderLogs`. 
```{r report_session_info, eval=TRUE}
sessionInfo()
```
