---
title: "VAR-Seq Workflow Template"
author: "First Last"
date: last-modified
sidebar: tutorials
bibliography: bibtex.bib
---

<p align="right">
  <a href="https://raw.githubusercontent.com/tgirke/GEN242/main/tutorials/spvarseq/spvarseq_index.qmd">
    <img src="https://img.shields.io/badge/Download-.qmd-blue?style=for-the-badge&logoColor=white" alt="Download qmd">
  </a>
</p>

<!--
- Render Rmd from command-line
Rscript -e "rmarkdown::render('spvarseq.Rmd', c('BiocStyle::html_document'), clean=F); knitr::knit('spvarseq.Rmd', tangle=TRUE)"
- Render qmd from command line:
quarto render spvarseq_index.qmd
Render qmd from R:
quarto::quarto_render("spvarseq_index.qmd")
-->


```{r setup_dir}
#| echo: false
#| include: false
#| message: false
#| warning: false
unlink(".SPRproject/", recursive = TRUE)

```


```{r setup_libraries}
#| echo: false
#| message: false
#| warning: false
suppressPackageStartupMessages({
    library(systemPipeR)
})

```


## Introduction

### Overview

This workflow template is for analyzing VAR-Seq data. It is provided by
[systemPipeRdata](https://bioconductor.org/packages/devel/data/experiment/html/systemPipeRdata.html), 
a companion package to [systemPipeR](https://www.bioconductor.org/packages/devel/bioc/html/systemPipeR.html) [@H_Backman2016-bt].
Similar to other `systemPipeR` workflow templates, a single command generates
the necessary working environment. This includes the expected directory
structure for executing `systemPipeR` workflows and parameter files for running
command-line (CL) software utilized in specific analysis steps. For learning
and testing purposes, a small sample (toy) data set is also included (mainly
FASTQ and reference genome files). This enables users to seamlessly run the
numerous analysis steps of this workflow from start to finish without the
requirement of providing custom data. After testing the workflow, users have
the flexibility to employ the template as is with their own data or modify it
to suit their specific needs. For more comprehensive information on designing
and executing workflows, users want to refer to the main vignettes of
[systemPipeR](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html)
and
[systemPipeRdata](https://www.bioconductor.org/packages/devel/data/experiment/vignettes/systemPipeRdata/inst/doc/systemPipeRdata.html). 

The `Rmd` file (`systemPipeVARseq.Rmd`) associated with this vignette serves a dual purpose. It acts
both as a template for executing the workflow and as a template for generating
a reproducible scientific analysis report. Thus, users want to customize the text
(and/or code) of this vignette to describe their experimental design and
analysis results. This typically involves deleting the instructions how to work
with this workflow, and customizing the text describing experimental designs,
other metadata and analysis results.

### Experimental design

Typically, the user wants to describe here the sources and versions of the
reference genome sequence along with the corresponding annotations. The standard
directory structure of `systemPipeR` (see [here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#3_Directory_structure)), 
expects the input data in a subdirectory named `data`
and all results will be written to a separate `results` directory. The Rmd source file 
for executing the workflow and rendering its report (here `systemPipeVARseq.Rmd`) is 
expected to be located in the parent directory. 

The test (toy) data set used by this template ([SRP010938](http://www.ncbi.nlm.nih.gov/sra/?term=SRP010938))
contains 18 paired-end (PE) read sets from *Arabidposis thaliana*
[@Howard2013-fq]. To minimize processing time during testing, each FASTQ
file has been reduced to 90,000-100,000 randomly sampled PE reads that
map to the first 100,000 nucleotides of each chromosome of the *A.
thaliana* genome. The corresponding reference genome sequence (FASTA) and
its GFF annotation files have been reduced to the same genome regions. This way the entire 
test sample data set is less than 200MB in storage space. A PE read set has been
chosen here for flexibility, because it can be used for testing both types
of analysis routines requiring either SE (single end) reads or PE reads.

To use their own VAR-Seq and reference genome data, users want to move or link the
data to the designated `data` directory and execute the workflow from the parent directory 
using their customized `Rmd` file. Beginning with this template, users should delete the provided test
data and move or link their custom data to the designated locations.
Alternatively, users can create an environment skeleton (named `new` [here](https://www.bioconductor.org/packages/devel/data/experiment/vignettes/systemPipeRdata/inst/doc/new.html)) or
build one from scratch. To perform a VAR-Seq analysis with new FASTQ files
from the same reference genome, users only need to provide the FASTQ files and
an experimental design file called 'targets' file that outlines the experimental 
design. The structure and utility of targets files is described in `systemPipeR's`
main vignette [here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#4_The_targets_file).

### Workflow steps

The default analysis steps included in this VAR-Seq workflow template are listed below. Users 
can modify the existing steps, add new ones or remove steps as needed.

__Default analysis steps__

1. Read preprocessing
    + Quality filtering (trimming)
    + FASTQ quality report
2. Alignments: `BWA` (or any other DNA aligner)
3. Alignment stats 
4. Variant calling 
    + with `GATK`, substeps 1-10
    + with `BCFtools`, 1 substep
5. Filter variants
    + Filter variants called by `GATK` 
    + Filter variants called by `BCFtools` 
6. Annotate filtered variants
7. Combine annotation results among samples
8. Calculate Summary statistics of variants
9. Plot Venn diagram of variants
10. Plot variants of interests

### Setup of workflow environment

<font color="red">NOTE: this section</font> describes how to set up the proper 
environment (directory structure) for running `systemPipeR` workflows in the GEN242
class. This routine applies to all workflows. After mastering this task the workflow 
run instructions <font color="red">can be deleted</font> since they are not expected 
to be included in a final HTML/PDF report of a workflow.

1. If a remote system or cluster is used, then users need to log in to the
   remote system first. The following applies to an HPC cluster (_e.g._ HPCC
   cluster). 

   A terminal application or RStudio server via onDemand needs to be used to log in to a user's cluster account. Next, one
   can open an interactive session on a computer node with `srun`. More details about
   argument settings for `srun` are available in this [HPCC
   manual](http://hpcc.ucr.edu/manuals_linux-cluster_jobs.html#partitions) or
   the HPCC section of this website
   [here](../linux/index.qmd#job-submission-with-sbatch).
   Next, load the R version required for running the workflow with `module load`. Sometimes it may be necessary to
   first unload an active software version before loading another version, _e.g._ `module unload R`.

From command-line

```{bash shell_01}
#| eval: false
srun --x11 --partition=gen242 --account=gen242 --mem=20gb --cpus-per-task 8 --ntasks 1 --time 20:00:00 --pty bash -l
module unload R; module load R/4.4.2

```

2. Load a workflow template with the `genWorkenvir` function. This can be done 
from the command-line or from within R. However, only one of the two options needs to be used.

The environment for this VAR-Seq workflow is auto-generated below with the 
`genWorkenvir` function (selected under `workflow="varseq"`). It is fully populated 
with a small test data set, including FASTQ files, reference genome and annotation data. The name of the
resulting workflow directory can be specified under the `mydirname` argument.
The default `NULL` uses the name of the chosen workflow. An error is issued if
a directory of the same name and path exists already. After this, the user’s R
session needs to be directed into the resulting `varseq` directory (here with
`setwd`). 

From command-line

```{bash shell_02}
#| eval: false
Rscript -e "systemPipeRdata::genWorkenvir(workflow='varseq')"
cd varseq

```

From R

```{r gen_workflow_envir}
#| eval: false
library(systemPipeRdata)
genWorkenvir(workflow = "varseq", mydirname = "chipseq")
setwd("varseq")

```

3. If the user wishes to use another `Rmd` file than the template instance 
provided by the `genWorkenvir` function, then it can be copied or downloaded 
into the root directory of the workflow environment (_e.g._ with `cp`, `download.file` 
or `wget`). For the first introduction to VAR-Seq analysis in GEN242, we will use the 
`Rmd` file obtained by `genWorkenvir`. Note, the `Rmd` source file of this tutorial page
is linked on the top.  

```{bash shell_03}
#| eval: false
# From command-line 
wget <*.Rmd> -O <*.Rmd> 
# From R
download.file("<*.Rmd>", "*.Rmd>")

```

4. Now one can open from the root directory of the workflow the corresponding R Markdown script (_e.g._ systemPipeVARseq.Rmd) using an R IDE, such as _nvim-r_, _ESS_ or RStudio. 
   Subsequently, the workflow can be run as outlined below. 

#### Import custom functions

Custom functions for the challenge projects can be imported with the source
command from a local R script (here [challengeProject_Fct.R](https://raw.githubusercontent.com/tgirke/GEN242/main/content/en/tutorials/spchipseq/challengeProject_Fct.R)). Skip this step if such a script is not available. Alternatively, these 
functions can be loaded from a custom R package. 

```{r load_custom_fct2}
#| eval: false
#| message: false
source("challengeProject_Fct.R")

```

#### Input data: `targets` file

The `targets` file defines the input files (e.g. FASTQ or BAM) and sample
comparisons used in a data analysis workflow. It can also store any number of
additional descriptive information for each sample. The following shows the first
four lines of the `targets` file used in this workflow template.


```{r load_targets_file}
targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4,-c(5,6)]

```

To work with custom data, users need to generate a _`targets`_ file containing
the paths to their own FASTQ files. [Here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#4_The_targets_file) is a detailed description of the structure and 
utility of `targets` files.


### Please Note

This template provides some common steps for a `VARseq` workflow. One can add, remove,
modify workflow steps by operating on the `SYSargsList` workflow object.
For more details of all the features and utilities, please consult the [main tutorial](
https://girke.bioinformatics.ucr.edu/GEN242/tutorials/systempiper/systempiper/).

In this template, code chunks with the option `spr = TRUE'` will be added to the workflow. 
Other R code chunks without this option will be ignored. The option `eval = FALSE` 
can be ignored when imported and build the workflow object. Please be aware of this 
possibility. 

The template can provide more than one alternative for each step, such as 
different mapping methods, that will receive the `mandatory` or `optional` flag. 
One can run just the `mandatory` steps, `ALL`, or `optional` steps when running 
the workflow.

Also, each one of the steps can be run on compute clusters (`compute` option)
or on the current session, here called `management` session. 
For the demonstration of this template, a `management` session will be chosen. 


## Quick start {#quick-start}

After a workflow environment has been created with the above `genWorkenvir` 
function call and the corresponding R session directed into the resulting directory (here `varseq`), 
the `SPRproject` function is used to initialize a new workflow project instance. The latter 
creates an empty `SAL` workflow container (below `sal`) and at the same time a
linked project log directory (default name `.SPRproject`) that acts as a
flat-file database of a workflow. Additional details about this process and 
the SAL workflow control class are provided in `systemPipeR's` main vignette 
[here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#11_Workflow_control_class)
and [here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#5_Detailed_tutorial).

Next, the `importWF` function imports all the workflow steps outlined in the
source Rmd file of this vignette (here `systemPipeVARseq.Rmd`) into the `SAL` workflow container.
An overview of the workflow steps and their status information can be returned
at any stage of the loading or run process by typing `sal`.

```{r project_varseq}
#| eval: false
library(systemPipeR)
sal <- SPRproject() 
sal <- importWF(sal, file_path = "systemPipeVARseq.Rmd", verbose = FALSE)
sal

```

After loading the workflow  into `sal`, it can be executed from start to finish
(or partially) with the `runWF` command. Running the workflow will only be
possible if all dependent CL software is installed on a user's system. Their
names and availability on a system can be listed with `listCmdTools(sal,
check_path=TRUE)`. For more information about the `runWF` command, refer to the
help file and the corresponding section in the main vignette
[here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#61_Overview). 

Running workflows in parallel mode on computer clusters is a straightforward
process in `systemPipeR`. Users can simply append the resource parameters (such
as the number of CPUs) for a cluster run to the `sal` object after importing
the workflow steps with `importWF` using the `addResources` function. More
information about parallelization can be found in the corresponding section at
the end of this vignette [here](#paralellization) and in the main vignette
[here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#63_Parallel_evaluation). 

```{r run_varseq}
#| eval: false
sal <- runWF(sal)

```
Workflows can be visualized as topology graphs using the `plotWF` function.

```{r plot_varseq}
#| eval: false
plotWF(sal)

```

![Figure 1: Toplogy graph of RNA-Seq workflow.](results/plotwf_varseq.png){fig-align="center" width=50%}

Scientific and technical reports can be generated with the `renderReport` and
`renderLogs` functions, respectively. Scientific reports can also be generated
with the `render` function of the `rmarkdown` package. The latter option with
`rmarkdown::render` is often more flexible and preferred for most users, since
it provides the advantage that any modifications to the Rmd file are instantly
reflected in the HTML report, eliminating the necessity to update the `sal`
object. The technical reports are based on log information that `systemPipeR` 
collects during workflow runs.

```{r report_varseq}
#| eval: false
# Scientific report
sal <- renderReport(sal)
rmarkdown::render("systemPipeVARseq.Rmd", clean = TRUE, output_format = "BiocStyle::html_document") 

# Technical (log) report
sal <- renderLogs(sal)

```

The `statusWF` function returns a status summary for each step in a `SAL` workflow instance.

```{r status_varseq}
#| eval: false
statusWF(sal)

```

## Workflow steps

The data analysis steps of this workflow are defined by the following workflow code chunks.
They can be loaded into `SAL` interactively, by executing the code of each step in the 
R console, or all at once with the `importWF` function used under the Quick start section. 
R and CL workflow steps are declared in the code chunks of `Rmd` files with the 
`LineWise` and `SYSargsList` functions, respectively, and then added to the `SAL` workflow 
container with `appendStep<-`. Their syntax and usage is described 
[here](https://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#52_Constructing_workflows). 

### Required packages and resources

The first step loads the `systemPipeR` package.

```{r load_SPR}
#| message: false
#| eval: false
#| spr: true
# Some toy data used here do not work well in the VAR-Seq workflow since 
# they take too much time to process. To save processing time,
# the toy data set is reduced to the first 13 samples.
# Please REMOVE the next two lines in your real analysis
cat(crayon::red$bold("Some samples in targets are removed for test workflow. Please change the template to disable this in your real analysis.\n"))
writeLines(readLines("targetsPE.txt")[1:13], "targetsPE.txt")

cat(crayon::blue$bold("To use this workflow, following R packages are expected:\n"))
cat(c("'GenomicFeatures", "VariantAnnotation", "GenomicFeatures", "ggbio", "ggplot2'\n"), sep = "', '")
###pre-end
appendStep(sal) <- LineWise(
    code = {
        library(systemPipeR)
        }, 
    step_name = "load_SPR"
)

```

### FASTQ quality report

The following `seeFastq` and `seeFastqPlot` functions generate and plot a series of
useful quality statistics for a set of FASTQ files, including per cycle quality
box plots, base proportions, base-level quality trends, relative k-mer
diversity, length, and occurrence distribution of reads, number of reads above
quality cutoffs and mean quality distribution. The results can be exported to
different graphics formats, such as a PNG file, here named `fastqReport.png`.
Detailed information about the usage and visual components in the quality plots
can be found in the corresponding help file (see `?seeFastq` or `?seeFastqPlot`).

The following step creates a quality report prior to trimming the fastq reads.
In a later step another quality report is generated after the reads are trimmed. 
It is recommended to inspect the first report in order to decide whether read 
trimming is required at all. 

Please note that initial targets files are being used here. In this case, 
it has been added to the first step, and later, we used the function `getColumn` 
to extract input file names in a named vector. 

```{r fastq_report_pre}
#| eval: false
#| message: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        targets <- read.delim("targetsPE.txt", comment.char = "#")
        updateColumn(sal, step = "load_SPR", position = "targetsWF") <- targets
        fq_files <- getColumn(sal, "load_SPR", "targetsWF", column = 1)
        fqlist <- seeFastq(fastq = fq_files, batchsize = 10000, klength = 8)
        png("./results/fastqReport.png", height = 162, width = 288 * length(fqlist))
        seeFastqPlot(fqlist)
        dev.off()
    },
    step_name = "fastq_report_pre",
    dependency = "load_SPR"
)

```

![Figure 2: FASTQ quality report for 18 samples.](results/fastqReport.png)

### Read preprocessing

#### Read trimming with Trimmomatic

Next, we need to populate the object created with the first step in the
workflow. Here, an example of how to perform this task using parameters template
files for trimming FASTQ files with [Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) 
software [@Bolger2014-yr].
For this step, the `SYSargsList` function has been used to build the command-line
and append to `sal` object. For more details of all the features and utilities, 
please consult the [main vignette](http://www.bioconductor.org/packages/release/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html).

If GATK (default) is used for variant calling, any type of fastq trimming is 
strongly depreciated. GATK have internal function to handle low quality posistions.

```{r trimmomatic}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "trimmomatic",
    targets = "targetsPE.txt",
    wf_file = "trimmomatic/trimmomatic-pe.cwl",
    input_file = "trimmomatic/trimmomatic-pe.yml",
    dir_path = "param/cwl",
    inputvars = c(
        FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"
    ),
    dependency = c("fastq_report_pre"),
    run_step = "optional"
)

```

#### Preprocessing with _`preprocessReads`_ function

The function _`preprocessReads`_ allows to apply predefined or custom
read preprocessing functions to all FASTQ files referenced in a
_`SYSargsList`_ container, such as quality filtering or adaptor trimming
routines. Internally, _`preprocessReads`_ uses the _`FastqStreamer`_ function from
the _`ShortRead`_ package to stream through large FASTQ files in a
memory-efficient manner. The following example performs adaptor trimming with
the _`trimLRPatterns`_ function from the _`Biostrings`_ package.

Here, we are appending this step at the _`SYSargsList`_ object created previously. 
All the parameters are defined on the _`preprocessReads/preprocessReads-pe.yml`_ file.

```{r preprocessing}
#| message: false
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "preprocessing",
    targets = "targetsPE.txt", dir = TRUE,
    wf_file = "preprocessReads/preprocessReads-pe.cwl",
    input_file = "preprocessReads/preprocessReads-pe.yml",
    dir_path = "param/cwl",
    inputvars = c(
        FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"
    ),
    dependency = c("fastq_report_pre"),
    run_step = "optional"
)

```

After the trimming step, the `outfiles` files can be used to generate the new 
targets files containing the paths to the trimmed FASTQ files. The new targets 
information can be used for the next workflow step instance, _e.g._ running the 
NGS alignments with the trimmed FASTQ files.

The following example shows how one can design a custom read _'preprocessReads'_
function using utilities provided by the _`ShortRead`_ package, and then run it
in batch mode with the _'preprocessReads'_ function. For here, it is possible to 
replace the function used on the `preprocessing` step and modify the `sal` object. 
Because it is a custom function, it is necessary to save the part in the R object, 
and internally the `preprocessReads.doc.R` is loading the function. If the R object
is saved with a different name (here `"param/customFCT.RData"`), please replace 
that accordingly in the `preprocessReads.doc.R`.

Please, note that this step is not added to the workflow, here just for demonstration.

First, we defined the function in the workflow:

```{r custom_preprocessing_function}
#| eval: false
appendStep(sal) <- LineWise(
    code = {
        filterFct <- function(fq, cutoff = 20, Nexceptions = 0) {
            qcount <- rowSums(as(quality(fq), "matrix") <= cutoff, na.rm = TRUE)
            # Retains reads where Phred scores are >= cutoff with N exceptions
            fq[qcount <= Nexceptions]
        }
        save(list = ls(), file = "param/customFCT.RData")
    },
    step_name = "custom_preprocessing_function",
    dependency = "preprocessing"
)

```

After, we can edit the input parameter:

```{r editing_preprocessing}
#| message: false
#| eval: false
yamlinput(sal, "preprocessing")$Fct
yamlinput(sal, "preprocessing", "Fct") <- "'filterFct(fq, cutoff=20, Nexceptions=0)'"
yamlinput(sal, "preprocessing")$Fct ## check the new function
cmdlist(sal, "preprocessing", targets = 1) ## check if the command line was updated with success

```

## FASTQ quality after Trimming

This is the post-trimming fastq quality report. If the trimming step is included, it is 
recommended to add this step to compare trimming of fastq before and after.

```{r fastq_report_pos}
#| eval: false
#| message: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        fq_files <- getColumn(sal, "preprocessing", "outfiles", column = 1) ## get outfiles path
        fqlist <- seeFastq(fastq = fq_files, batchsize = 10000, klength = 8)
        png("./results/fastqReport_pos.png", height = 18, width = 4 * length(fqlist))
        seeFastqPlot(fqlist)
        dev.off()
    },
    step_name = "fastq_report_pos",
    dependency = "trimmomatic",
    run_step = "optional"
)

```

### Read mapping with `BWA-MEM` 

The NGS reads of this project are aligned against the reference genome
sequence using the highly variant tolerant short read aligner `BWA-MEM`
[@Li2013-oy; @Li2009-oc]. The parameter settings of the aligner are
defined in the `param/cwl/gatk/bwa-pe.cwl`.

This test code uses untrimmed fastq files since the demo data is minimal and 
limited. However, it is best to test with `FASTQ quality report` function provided 
above to verify your real data first.

#### Build index and dictionary files for BWA and GATK

Build the index and dictionary files for BWA and GATK to run. 

```{r bwa_index}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "bwa_index",
    dir = FALSE, targets = NULL,
    wf_file = "gatk/workflow_bwa-index.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = "load_SPR"
)

```

Create reference `fasta` dictionary.

```{r fasta_index}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "fasta_index",
    dir = FALSE, targets = NULL,
    wf_file = "gatk/workflow_fasta_dict.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = "bwa_index"
)

```

Create dictionary index.

```{r faidx_index}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "faidx_index",
    dir = FALSE, targets = NULL,
    wf_file = "gatk/workflow_fasta_faidx.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = "fasta_index"
)

```

#### Mapping reads with BWA

```{r bwa_alignment}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "bwa_alignment",
    targets = "targetsPE.txt",
    wf_file = "gatk/workflow_bwa-pe.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(
        FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"
    ),
    dependency = c("faidx_index")
)

```

### Read and alignment stats

The following provides an overview of the number of reads in each sample
and how many of them aligned to the reference.

```{r align_stats}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        bampaths <- getColumn(sal, step = "bwa_alignment", "outfiles", column = "samtools_sort_bam")
        fqpaths <- getColumn(sal, step = "bwa_alignment", "targetsWF", column = "FileName1")
        read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE)
        write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE, sep = "\t")
    },
    step_name = "align_stats",
    dependency = "bwa_alignment",
    run_step = "optional"
)

```

### Create symbolic links for viewing BAM files in IGV

The `symLink2bam` function creates symbolic links to view the BAM alignment files in a
genome browser such as IGV. The corresponding URLs are written to a file
with a path specified under `urlfile` in the `results` directory.

```{r bam_urls}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        bampaths <- getColumn(sal, step = "bwa_alignment", "outfiles", column = "samtools_sort_bam")
        symLink2bam(
            sysargs = bampaths, htmldir = c("~/.html/", "somedir/"),
            urlbase = "http://cluster.hpcc.ucr.edu/~tgirke/",
            urlfile = "./results/IGVurl.txt"
        )
    },
    step_name = "bam_urls",
    dependency = "bwa_alignment",
    run_step = "optional"
)

```

### Variant calling

The following performs variant calling with `GATK` and `BCFtools` on a single 
machine by `runWF` function for each sample sequentially. If a cluster compute
is available, running in parallel mode on a compute cluster can be performed by 
`runWF`, making available the resources and choose `run_session = "compute"`. 

Not all users have a cluster system, so here to demonstrate an example of variant calling 
workflow, only single-machine commands are shown. For cluster jobs, please refer 
to our main [vignette](http://www.bioconductor.org/packages/release/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html). 

In addition, the user would choose only one variant caller here rather than 
running several ones. However, the workflow manager allows keeping multiple 
options available for running the analysis.

#### Variant calling with `GATK`

The following steps are based on `GATK 4.1.1.0` [Best Practice](https://software.broadinstitute.org/gatk/best-practices/). 
There are 10 individual steps where the user can choose where to jump in and where to skip. 
All scripts are located at `param/cwl/gatk`. `BQSR` (Base Quality Score Recalibration) 
and `VQSR` (Variant Quality Score Recalibration) are very specific 
to a limited species like human, so this workflow does not support these steps. 

#### Step 1: `fastq` to `ubam`

Convert `fastq` files to `bam` files to prepare for the following step. It is very 
important to specific your sequencing platform, default is `illumina`. User need 
to change `param/cwl/gatk/gatk_fastq2ubam.cwl` if the platform is different. Platform information 
is needed for the variant caller in later steps to correct calling parameters.

```{r fastq2ubam}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "fastq2ubam",
    targets = "targetsPE.txt",
    wf_file = "gatk/workflow_gatk_fastq2ubam.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(
        FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"
    ),
    dependency = c("faidx_index")
)

```

#### Step 2: Merge `bam` and `ubam`

This step merges a `bam` and `ubam` and creates a third `bam` file that contains 
alignment information and remaining information that was removed by the aligner like `BWA`. 
The removed information is essential for variant statistics calculation. Previous steps are 
recommended, but variant calling can still be performed without these steps.

```{r merge_bam}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "merge_bam",
    targets = c("bwa_alignment", "fastq2ubam"),
    wf_file = "gatk/workflow_gatk_mergebams.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(
        bwa_men_sam = "_bwasam_",
        ubam = "_ubam_",
        SampleName = "_SampleName_"
    ),
    rm_targets_col = c("preprocessReads_1", "preprocessReads_2"),
    dependency = c("bwa_alignment", "fastq2ubam")
)

```


#### Step 3: Sort `bam` files by genomic coordinates

Sort `bam` files by genomic coordinates.

```{r sort}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "sort",
    targets = "merge_bam",
    wf_file = "gatk/workflow_gatk_sort.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(merge_bam = "_mergebam_", SampleName = "_SampleName_"),
    rm_targets_col = c(
        "bwa_men_sam", "ubam", "SampleName_fastq2ubam",
        "Factor_fastq2ubam", "SampleLong_fastq2ubam",
        "Experiment_fastq2ubam", "Date_fastq2ubam"
    ),
    dependency = c("merge_bam")
)

```

#### Step 4: Mark duplicates

Mark PCR artifacts in sequencing. A `duplicate_metrics` file will also be produced 
by this step, but will not be used for the next step. This file is just for the user 
to check duplicates status summary.

```{r mark_dup}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "mark_dup",
    targets = "sort",
    wf_file = "gatk/workflow_gatk_markduplicates.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(sort_bam = "_sort_", SampleName = "_SampleName_"),
    rm_targets_col = c("merge_bam"),
    dependency = c("sort")
)

```

#### Step 5: Fixing tags

Takes the `bam` from the last step and calculates the NM, MD, and UQ tags. 
These tags are important for variant calling and filtering. 
This step is recommended but can be skipped.  

```{r fix_tag}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "fix_tag",
    targets = "mark_dup",
    wf_file = "gatk/workflow_gatk_fixtag.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(mark_bam = "_mark_", SampleName = "_SampleName_"),
    rm_targets_col = c("sort_bam"),
    dependency = c("mark_dup")
)

```

Up till this step, sample preprocess is done. All analysis ready `BAM` files and 
their index `.bai` files are created. Individual and cohort calling by 
`HaplotypeCaller` is performed from the next step.

#### Step 6: HaplotypeCaller `gvcf`

The `HaplotypeCaller` is running a **gvcf** mode in this step. G stands for 'genomic'. 
The file not only contains variant sites information but also non-variant sites information; 
thus, at the following step, the cohort caller can use this information to validate the true variants.

```{r hap_caller}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "hap_caller",
    targets = "fix_tag",
    wf_file = "gatk/workflow_gatk_haplotypecaller.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(fixtag_bam = "_fixed_", SampleName = "_SampleName_"),
    rm_targets_col = c("mark_bam"),
    dependency = c("fix_tag")
)

```

#### Step 7: Import all `gvcfs`

It is recommended to import all **gvcfs** to a 
[TileDB](https://github.com/Intel-HLS/GenomicsDB/wiki) database for fast cohort 
variant calling at the following step. Note: if you are working with non-diploid data, 
use `CombineGVCFs` function from `GATK` and change the `gvcf_db_folder` parameter 
in `param/cwl/gatk/gatk.yaml` to be your combined **gvcf** file path.

**Important**: Make sure all samples' `*.g.vcf.gz` files are in the results folder,
also the `tbi index` files also should be there.

```{r import}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "import",
    targets = NULL, dir = FALSE,
    wf_file = "gatk/workflow_gatk_genomicsDBImport.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = c("hap_caller")
)

```

#### Step 8: Cohort calling of `gvcf`

Assess variants by information from all `gvcfs`. A collective `vcf` called 
`samples.vcf.gz` is created by default naming.

```{r call_variants}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "call_variants",
    targets = NULL, dir = FALSE,
    wf_file = "gatk/workflow_gatk_genotypeGVCFs.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = c("import")
)

```

#### Step 9: Cohort hard filter variants

Variant Quality Score Recalibration (VQSR) is not included in this workflow. 
Variants are hard filtered together.
See this [Post](https://gatkforums.broadinstitute.org/gatk/discussion/2806/howto-apply-hard-filters-to-a-call-set) for parameters for hard filtering. Change these settings in `param/cwl/gak/gatk_variantFiltration.sh` if needed. VQSR requires a large quantity of 
samples to be training data before you can do filtering. Read this 
[post](https://gatk.broadinstitute.org/hc/en-us/articles/360035531612-Variant-Quality-Score-Recalibration-VQSR-) for more information.

```{r filter}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "filter",
    targets = NULL, dir = FALSE,
    wf_file = "gatk/workflow_gatk_variantFiltration.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    dependency = c("call_variants")
)

```

#### Step 10: Extract variant

After cohort calling, filtering, all variants for all samples are stored in one big file. 
Extract variants for each sample and save them separately (only variants that have 
passed the filters are stored).

```{r create_vcf}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "create_vcf",
    targets = "hap_caller",
    wf_file = "gatk/workflow_gatk_select_variant.cwl",
    input_file = "gatk/gatk.yaml",
    dir_path = "param/cwl",
    inputvars = c(SampleName = "_SampleName_"),
    dependency = c("hap_caller", "filter")
)

```


### Variant calling with `BCFtools`

Alternative option with `BCFtool`: 

The following runs the variant calling with `BCFtools`. This tool takes `BWA` 
aligned `BAM` files, sort, mark duplicates by `samtools` and finally call variants 
by `BCFtools`. 
**For legacy reasons we keep this option.**

```{r create_vcf_BCFtool}
#| eval: false
#| spr: true
appendStep(sal) <- SYSargsList(
    step_name = "create_vcf_BCFtool",
    targets = "bwa_alignment", dir = TRUE,
    wf_file = "workflow-bcftools/workflow_bcftools.cwl",
    input_file = "workflow-bcftools/bcftools.yml",
    dir_path = "param/cwl",
    inputvars = c(bwa_men_sam = "_bwasam_", SampleName = "_SampleName_"),
    rm_targets_col = c("preprocessReads_1", "preprocessReads_2"),
    dependency = "bwa_alignment",
    run_step = "optional"
)

```

Variant calling ends here. Downstream analysis starts from the next section.

### Inspect VCF file 

Scripts of downstream analysis are stored in `param/cwl/varseq_downstream`.

**optional**: This step is not included in the default workflow. After successfully 
execute the entire workflow, users may load individual vcf files to R for 
other analysis like below.

VCF files can be imported into R with the `readVcf` function. 
Both `VCF` and `VRanges` objects provide convenient data structure for
working with variant data (_e.g._ SNP quality filtering). 

This step is not included in the default workflow steps, but can be useful to 
inspect individual sample's raw variants.

```{r inspect_vcf}
#| eval: false
library(VariantAnnotation)
vcf_raw <- getColumn(sal, "create_vcf")
vcf <- readVcf(vcf_raw[1], "A. thaliana")
vcf
vr <- as(vcf, "VRanges")
vr

```

### Filter variants

The function `filterVars` filters VCF files based on user definable
quality parameters. It sequentially imports each VCF file into R, applies the
filtering on an internally generated `VRanges` object and then writes
the results to a new subsetted VCF file. The filter parameters are passed on to
the corresponding argument as a character string. The function applies this
filter to the internally generated `VRanges` object using the standard
subsetting syntax for two dimensional objects such as: `vr[filter, ]`.

#### Filter variants called by `GATK` 

The below example filters for variants that are supported by `>=x`
reads and >=80% of them support the called variants. In addition, all
variants need to pass `>=x` of the soft filters recorded in the VCF
files generated by GATK. Since the toy data used for this workflow is
very small, the chosen settings are unreasonabley relaxed. A more
reasonable filter setting is given in the line below (here commented
out).

There is already some cohort filtering in GATK step 10. Some additional hard 
filtering is provided here. This step is included here, but in a real analysis,
you may skip this step.

For real samples, use following filters:
`filter <- "totalDepth(vr) >= 20 & (altDepth(vr) / totalDepth(vr) >= 0.8)"`

```{r filter_vcf}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        vcf_raw <- getColumn(sal, "create_vcf")
        library(VariantAnnotation)
        filter <- "totalDepth(vr) >= 2 & (altDepth(vr) / totalDepth(vr) >= 0.8)"
        vcf_filter <- suppressWarnings(filterVars(vcf_raw, filter, organism = "A. thaliana", out_dir = "results/vcf_filter"))
        # dump the filtered path variable to running enviornment so 
        # other sysArg steps can get its values
        updateColumn(sal, 'create_vcf', "outfiles") <- data.frame(vcf_filter=vcf_filter)
    },
    step_name = "filter_vcf",
    dependency = "create_vcf"
)

```


#### Filter variants called by `BCFtools`  

The following shows how to filter the VCF files generated by `BCFtools` using
similar parameter settings as in the previous filtering of the GATK
results.

```{r filter_vcf_BCFtools}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        vcf_raw <- getColumn(sal, step = "create_vcf_BCFtool", 
                             position = "outfiles", column = "bcftools_call")
        library(VariantAnnotation)
        filter <- "rowSums(vr) >= 2 & (rowSums(vr[,3:4])/rowSums(vr[,1:4]) >= 0.8)"
        vcf_filter_bcf <- suppressWarnings(filterVars(vcf_raw, filter, organism = "A. thaliana", out_dir = "results/vcf_filter_BCFtools", varcaller = "bcftools"))
        
        updateColumn(sal, 'create_vcf', "outfiles") <- data.frame(vcf_filter_bcf=vcf_filter_bcf)
    },
    step_name = "filter_vcf_BCFtools",
    dependency = "create_vcf_BCFtool",
    run_step = "optional"
)

```

Check filtering outcome for one sample

This mini step can be used to compare `vcfs` files before and after filtering.
This can be used once the workflow has been run, and make sure "filter_vcf" is 
done, since it is an optional step.

```{r check_filter}
#| eval: false
copyEnvir(sal, "vcf_raw", globalenv())
copyEnvir(sal, "vcf_filter", globalenv())
length(as(readVcf(vcf_raw[1], genome = "Ath"), "VRanges")[, 1])
length(as(readVcf(vcf_filter[1], genome = "Ath"), "VRanges")[, 1])

```

### Annotate filtered variants

The function `variantReport` generates a variant report using
utilities provided by the `VariantAnnotation` package. The report for
each sample is written to a tabular file containing genomic context annotations
(_e.g._ coding or non-coding SNPs, amino acid changes, IDs of affected
genes, etc.) along with confidence statistics for each variant. The CWL
file `param/cwl/varseq_downstream/annotate.cwl` defines the paths to the input 
and output files which are stored in a `SYSargs2` instance. 

#### Basics of annotating variants

This step can be run after running the default workflow, not included in the default. 

Variants overlapping with common annotation features can be identified with `locateVariants`.

```{r annotate_basics}
#| eval: false
library("GenomicFeatures")
# comment the next line if optional step "filter_vcf" is included
vcf_filter <- getColumn(sal, "create_vcf")
# uncomment the next line if optional step "filter_vcf" is included
# copyEnvir(sal, "vcf_filter", globalenv())
txdb <- loadDb("./data/tair10.sqlite")
vcf <- readVcf(vcf_filter[1], "A. thaliana")
locateVariants(vcf, txdb, CodingVariants())

```

Synonymous/non-synonymous variants of coding sequences are computed by the 
`predictCoding` function for variants overlapping with coding regions.

```{r annotate_basics_non}
#| eval: false
fa <- FaFile("data/tair10.fasta")
predictCoding(vcf, txdb, seqSource = fa)

```

#### Annotate filtered variants `GATK` or `BCFtools`

**required**

```{r annotate_vcf}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        # get the filtered vcf path from R running environment
        copyEnvir(sal, "vcf_filter", globalenv())  
        library("GenomicFeatures")
        txdb <- loadDb("./data/tair10.sqlite")
        fa <- FaFile("data/tair10.fasta")
        vcf_anno <- suppressMessages(suppressWarnings(variantReport(vcf_filter, txdb = txdb, fa = fa, organism = "A. thaliana", out_dir = "results/vcf_anno")))
    },
    step_name = "annotate_vcf",
    dependency = "filter_vcf"
)

```

View annotation result for single sample

```{r view_annotation}
#| eval: false
copyEnvir(sal, "vcf_anno", globalenv())
read.delim(vcf_anno[1])[38:40, ]

```

### Combine annotation results among samples

To simplify comparisons among samples, the `combineVarReports`
function combines all variant annotation reports referenced in a
`SYSargs2` instance (here `args`). At the same time the function
allows to consider only certain feature types of interest. For instance, the
below setting `filtercol=c(Consequence="nonsynonymous")` will include
only nonsysynonymous variances listed in the `Consequence` column of
the annotation reports. To omit filtering, one can use the setting
`filtercol="All"`.

#### Combine results 

**required**

```{r combine_var}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        combineDF <- combineVarReports(vcf_anno, filtercol = c(Consequence = "nonsynonymous"))
        write.table(combineDF, "./results/combineDF_nonsyn.tsv", quote = FALSE, row.names = FALSE, sep = "\t")
    },
    step_name = "combine_var",
    dependency = "annotate_vcf"
)

```

### Summary statistics of variants

The `varSummary` function counts the number of variants for each feature type
included in the annotation reports.

### Summary of variants

**required**

```{r summary_var}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        write.table(varSummary(vcf_anno), "./results/variantStats.tsv", quote = FALSE, col.names = NA, sep = "\t")
    },
    step_name = "summary_var",
    dependency = "combine_var"
)

```

### Venn diagram of variants

**Optional** but included in the default

The venn diagram utilities defined by the `systemPipeR` package can be used to
identify common and unique variants reported for different samples
and/or variant callers. The below generates a 3-way venn diagram
comparing 3 samples for each of the two variant callers.

```{r venn_diagram}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        ## make a list of first three samples
        varlist <- sapply(names(vcf_anno[1:3]), function(x) as.character(read.delim(vcf_anno[x])$VARID))
        vennset <- overLapper(varlist, type = "vennsets")
        png("./results/vennplot_var.png")
        vennPlot(list(vennset), mymain = "Venn Plot of First 3 Samples", mysub = "", colmode = 2, ccol = c("red", "blue"))
        dev.off()
    },
    step_name = "venn_diagram",
    dependency = "annotate_vcf"
)

```

![Figure 3: Venn Diagram for 3 samples from GATK and BCFtools.](results/vennplot_var.png)


### Plot variants programmatically 

**Optional** but included in default

The following plots a selected variant with `ggbio`.

In this example, the input `BAM` file is from the `GATK` step 5, analysis ready bam. 
You can use other aligned `BAMs` as well, but make sure they are indexed. The `VCF` 
file is taken from `Inspect VCF file ` section or you can load your own vcf.

```{r plot_variant}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        # get the filtered vcf path from R running environment
        copyEnvir(sal, "vcf_filter", globalenv())  
        library(ggbio)
        library(VariantAnnotation)
        mychr <- "ChrM"
        mystart <- 19000
        myend <- 21000
        bams <- getColumn(sal, "fix_tag")
        vcf <- suppressWarnings(readVcf(vcf_filter["M6B"], "A. thaliana"))
        ga <- readGAlignments(bams["M6B"], use.names = TRUE, param = ScanBamParam(which = GRanges(mychr, IRanges(mystart, myend))))
        p1 <- autoplot(ga, geom = "rect")
        p2 <- autoplot(ga, geom = "line", stat = "coverage")
        p3 <- autoplot(vcf[seqnames(vcf) == mychr], type = "fixed") +
            xlim(mystart, myend) +
            theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank())
        p4 <- autoplot(loadDb("./data/tair10.sqlite"), which = GRanges(mychr, IRanges(mystart, myend)), names.expr = "gene_id")
        p1_4 <- tracks(Reads = p1, Coverage = p2, Variant = p3, Transcripts = p4, heights = c(0.3, 0.2, 0.1, 0.35)) + ylab("")
        ggbio::ggsave(p1_4, filename = "./results/plot_variant.png", units = "in")
    },
    step_name = "plot_variant",
    dependency = "filter_vcf"
)

```

![Figure 4: Plot variants programmatically.](results/plot_variant.png){fig-align="center" width=50%}

### Workflow Information

```{r sessionInfo}
#| eval: false
#| spr: true
appendStep(sal) <- LineWise(
    code = {
        sessionInfo()
        }, 
    step_name = "sessionInfo", 
    dependency = "plot_variant")

```

## Running workflow

### Interactive job submissions in a single machine

For running the workflow, `runWF` function will execute 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)

```

### Parallelization on clusters

Alternatively, the computation can be greatly accelerated by processing many files 
in parallel using several compute nodes of a cluster, where a scheduling/queuing
system is used for load balancing. 

The `resources` list object provides the number of independent parallel cluster 
processes defined under the `Njobs` element in the list. The following example 
will run 18 processes in parallel using each 4 CPU cores. 
If the resources available on a cluster allow running all 18 processes at the 
same time, then the shown sample submission will utilize in a total of 72 CPU cores.

Note, `runWF` can be used with most queueing systems as it is based on utilities 
from the `batchtools` package, which supports the use of template files (_`*.tmpl`_)
for defining the run parameters of different schedulers. To run the following 
code, one needs to have both a `conffile` (see _`.batchtools.conf.R`_ samples [here](https://mllg.github.io/batchtools/)) 
and a `template` file (see _`*.tmpl`_ samples [here](https://github.com/mllg/batchtools/tree/master/inst/templates)) 
for the queueing available on a system. The following example uses the sample 
`conffile` and `template` files for the Slurm scheduler provided by this package. 

The resources can be appended when the step is generated, or it is possible to 
add these resources later, as the following example using the `addResources` 
function:

```{r runWF_cluster}
#| eval: false
# wall time in mins, memory in MB
resources <- list(conffile=".batchtools.conf.R",
                  template="batchtools.slurm.tmpl", 
                  Njobs=18, 
                  walltime=120,
                  ntasks=1,
                  ncpus=4, 
                  memory=1024,
                  partition = "short"
                  )
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)

```

### Visualize workflow

_`systemPipeR`_ workflows instances can be visualized with the `plotWF` function.

```{r plotWF}
#| eval: false
plotWF(sal, rstudio = TRUE)

```

### Checking workflow status

To check the summary of the workflow, we can use:

```{r statusWF}
#| eval: false
sal
statusWF(sal)

```

### Accessing logs report

_`systemPipeR`_ compiles all the workflow execution logs in one central location,
making it easier to check any standard output (`stdout`) or standard error
(`stderr`) for any command-line tools used on the workflow or the R code stdout.

```{r logsWF}
#| eval: false
sal <- renderLogs(sal)

```

### 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("trimmomatic/0.39", "samtools/1.14", "gatk/4.2.0.0", "bcftools/1.15", 
          "bwa/0.7.17"), sep = "\n")
}

```

### Report 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}
sessionInfo()

```

## Funding

This project was supported by funds from the National Institutes of
Health (NIH) and the National Science Foundation (NSF).

## References


