---
title: "systemPipeR: Workflow and Visualization Toolkit" 
author: "Author: Daniela Cassol, Le Zhang, and Thomas Girke"
date: "Last update: `r format(Sys.time(), '%d %B, %Y')`" 
output:
  BiocStyle::html_document:
    toc_float: true
    code_folding: show
package: systemPipeR
vignette: |
  %\VignetteEncoding{UTF-8}
  %\VignetteIndexEntry{systemPipeR: Workflow and Visualization Toolkit}
  %\VignetteEngine{knitr::rmarkdown}
fontsize: 14pt
bibliography: bibtex.bib
editor_options: 
  markdown: 
    wrap: 80
  chunk_output_type: console
---

<!--
- Compile from command-line
Rscript -e "rmarkdown::render('systemPipeR.Rmd', c('BiocStyle::html_document'), clean=F); knitr::knit('systemPipeR.Rmd', tangle=FALSE)"
-->

```{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 = 100, 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 = 100), tidy = FALSE)
```

```{r setting, eval=TRUE, echo=FALSE}
if (file.exists("spr_project")) unlink("spr_project", recursive = TRUE)
is_win <- Sys.info()[['sysname']] != "Windows"
```

```{r load_library, eval=TRUE, include=FALSE}
suppressPackageStartupMessages({
    library(systemPipeR)
})
```

# Introduction

[_`systemPipeR`_](http://www.bioconductor.org/packages/devel/bioc/html/systemPipeR.html)
provides flexible utilities for designing, building and running automated end-to-end
analysis workflows for a wide range of research applications, including
next-generation sequencing (NGS) experiments, such as RNA-Seq, ChIP-Seq,
VAR-Seq and Ribo-Seq [@H_Backman2016-bt]. Important features include a uniform
workflow interface across different data analysis applications, automated
report generation, and support for running both R and command-line software,
such as NGS aligners or peak/variant callers, on local computers or compute
clusters (Figure \@ref(fig:utilities)). The latter supports interactive job submissions and batch submissions to queuing systems of clusters. 

_`systemPipeR`_ has been designed to improve the reproducibility of large-scale data
analysis projects while substantially reducing the time it takes to analyze 
complex omics data sets. It provides a uniform workflow interface and management
system that allows the users to run selected workflow steps, as well as customize and design entirely new workflows. Additionally, the package take advantage of central community S4 
classes of the Bioconductor ecosystem, and enhances them with command-line software support.

The main motivation and advantages of using _`systemPipeR`_ for complex data analysis tasks are:

1. Design of complex workflows involving multiple R/Bioconductor packages
2. Common workflow interface for different applications
3. User-friendly access to widely used Bioconductor utilities 
4. Support of command-line software from within R
5. Reduced complexity of using compute clusters from R 
6. Accelerated runtime of workflows via parallelization on computer systems with multiple CPU cores and/or multiple nodes
6. Improved reproducibility by automating the generation of analysis reports 

```{r utilities, eval=TRUE, warning= FALSE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Relevant features in `systemPipeR`. Workflow design concepts are illustrated under (A). Examples of `systemPipeR's` visualization functionalities are given under (B).", warning=FALSE}
knitr::include_graphics("images/utilities.png")
```

A central concept for designing workflows within the _`systemPipeR`_
environment is the use of workflow management containers. Workflow management
containers facilitate the design and execution of complex data
analysis steps. For its command-line interface _`systemPipeR`_ adopts the
widely used [Common Workflow Language](https://www.commonwl.org/) (CWL)
[@Amstutz2016-ka]. The interface to CWL is established by _`systemPipeR's`_
workflow control class called _`SYSargsList`_ (see Figure \@ref(fig:general)). This 
design offers many advantages such as: (i) options to run workflows either
entirely from within R, from various command-line wrappers (e.g., *cwl-runner*)
or from other languages (*, e.g.,* Bash or Python). Apart from providing
support for both command-line and R/Bioconductor software, the package provides
utilities for containerization, parallel evaluations on computer clusters and
automated generation of interactive analysis reports.

```{r general, warning= FALSE, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Overview of `systemPipeR` workflows management instances. A) A typical analysis workflow requires multiple software tools (red), as well the description of the input (green) and output files, including analysis reports (purple). B) `systemPipeR` provides multiple utilities to design and build a workflow, allowing multi-instance, integration of R code and command-line software, a simple and efficient annotation system, that allows automatic control of the input and output data, and multiple resources to manage the entire workflow. c) Options are provided to execute single or multiple workflow steps, while enabling scalability, checkpoints, and generation of technical and scientific reports.", warning=FALSE}

knitr::include_graphics("images/general.png")
```

An important feature of _`systemPipeR's`_ CWL interface is that it provides two
options to run command-line tools and workflows based on CWL. First, one can
run CWL in its native way via an R-based wrapper utility for *cwl-runner* or
*cwl-tools* (CWL-based approach). Second, one can run workflows using CWL's
command-line and workflow instructions from within R (R-based approach). In the
latter case the same CWL workflow definition files (*e.g.* `*.cwl` and `*.yml`)
are used but rendered and executed entirely with R functions defined by
_`systemPipeR`_, and thus use CWL mainly as a command-line and workflow
definition format rather than execution software to run workflows. In this regard
_`systemPipeR`_ also provides several convenience functions that are useful for
designing and debugging workflows, such as a command-line rendering function to
retrieve the exact command-line strings for each data set and processing step
prior to running a command-line.

## Workflow Management with _`SYSargsList`_ 

The `SYSargsList` S4 class is a list-like container that stores the paths to
all input and output files along with the corresponding parameters used in each
analysis step (see Figure \@ref(fig:sysargslistImage)). `SYSargsList` instances are constructed from an
optional targets files, and two CWL parameter files including `*.cwl` and
`*.yml` (for details, see below). When running preconfigured NGS workflows, the
only input the user needs to provide is the initial targets file containing the
paths to the input files (e.g., FASTQ) and experiment design information, such
as sample labels and biological replicates. Subsequent targets instances are
created automatically, based on the connectivity establish between each
workflow step. _`SYSargsList`_ containers store all information required for
one or multiple steps. This establishes central control for running, monitoring
and debugging complex workflows from start to finish. 

```{r sysargslistImage, warning= FALSE, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Workflow steps with input/output file operations are controlled by the _`SYSargsList`_ container. Each of its components (_`SYSargs2`_) are constructed from an optional *targets* and two *param* files. Alternatively, _`LineWise`_ instances containing pure R code can be used.", warning=FALSE}

knitr::include_graphics("images/SYSargsList.png")
```

## Command-line software support

_`systemPipeR`_ adopted the widely used community standard Common Workflow 
Language (CWL) [@Amstutz2016-ka] for describing command-line tools and workflows
in a declarative, generic, and reproducible manner. CWL specifications are 
text-based files and can be structured using YAML (https://yaml.org/) syntax. 
Therefore, the description files can be easily accessed and are readable. 
The significant advantage of adopting CWL as a standard description of 
command-line tools within \textit{systemPipeR} is the flexibility of workflow 
reusability for different computing environments and workflow frameworks in the 
community, improving reproducibility, portability, and shareability between
collaborators and community.

Following the CWL Command Line Tool Description Specification 
(https://www.commonwl.org/v1.2/CommandLineTool.html), the basic elements of the 
CWL tool description are defined in two files. Figure \@ref(fig:sprandCWL)A-B illustrate 
the “hello world” example. The main file contains all the information necessary 
to build the command-line that will be executed, specifying the input, expected
output files, and arguments for the command-line (Figure \@ref(fig:sprandCWL)A). 
The second file is optional yet provides flexibility to assign values to 
parameters required to make the input or output objects when building the command-line. 

```{r sprandCWL, warning=FALSE, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Example of 'Hello World' message using CWL syntax and demonstrating the connectivity with `systemPipeR`. (A) This file describes the command-line tool, here using 'echo' command. (B) This file describes all the parameters variables connected with the tool specification file. Here the reference value of the input parameter can be specific or can be filled dynamically, adding a variable that connects with the targets files from `systemPipeR`. (C) `SYSargsList` function provides the 'inputvars' arguments to build the connectivity between the CWL description of parameters and the targets files. The argument requires a named vector where each vector element is required to be defined in the CWL description of parameters file (B), and the names of the elements are needed to match the column names defined in the targets file (D).", warning=FALSE}

knitr::include_graphics("images/SPR_CWL_hello.png")
```

# Getting Started

## Installation

[_`systemPipeR`_](http://www.bioconductor.org/packages/devel/bioc/html/systemPipeR.html) 
environment can be installed from the R console using the [_`BiocManager::install`_](https://cran.r-project.org/web/packages/BiocManager/index.html) 
command. The associated data package [_`systemPipeRdata`_](http://www.bioconductor.org/packages/devel/data/experiment/html/systemPipeRdata.html) 
can be installed the same way. The latter is a helper package for generating _`systemPipeR`_ 
workflow environments with a single command containing all parameter files and 
sample data required to quickly test and run workflows. 

```{r install, eval=FALSE}
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager")
BiocManager::install("systemPipeR")
BiocManager::install("systemPipeRdata")
```

Please note that if you desire to use a third-party command-line tool, the particular
tool and dependencies need to be installed and exported in your PATH. 
See [details](#third-party-software-tools). 

## Loading package and documentation

```{r documentation, eval=FALSE}
library("systemPipeR") # Loads the package
library(help="systemPipeR") # Lists package info
vignette("systemPipeR") # Opens vignette
```

## How to get help for systemPipeR

All questions about the package or any particular function should be posted to 
the Bioconductor support site [https://support.bioconductor.org](https://support.bioconductor.org). 

Please add the "_`systemPipeR`_" tag to your question, and the package authors will 
automatically receive an alert. 

We appreciate receiving reports of bugs in the functions or documentation and
suggestions for improvement. For that, please consider opening an issue at
[GitHub](https://github.com/tgirke/systemPipeR/issues/new). 

# Project structure

_`systemPipeR`_ expects a project directory structure that consists of a directory
where users may store all the raw data, the results directory that will be reserved 
for all the outfiles files or new output folders, and the parameters directory. 

This structure allows reproducibility and collaboration across the data science 
team since internally relative paths are used. Users could transfer this project 
to a different location and still be able to run the entire workflow. Also, it 
increases efficiency and data management once the raw data is kept in a separate 
folder and avoids duplication. 

### Load sample data and workflow templates

The mini sample FASTQ files used by this overview vignette as well as the
associated workflow reporting vignettes can be loaded via the
_`systemPipeRdata`_ package as shown below. The chosen data set
[`SRP010938`](http://www.ncbi.nlm.nih.gov/sra/?term=SRP010938) obtains 18
paired-end (PE) read sets from _Arabidposis thaliana_ [@Howard2013-fq]. To
minimize processing time during testing, each FASTQ file has been subsetted to
90,000-100,000 randomly sampled PE reads that map to the first 100,000
nucleotides of each chromosome of the _A. thalina_ genome. The corresponding
reference genome sequence (FASTA) and its GFF annotation files (provided in the
same download) have been truncated accordingly. This way the entire test sample
data set requires less than 200MB disk storage space. A PE read set has been
chosen for this test data set for flexibility, because it can be used for
testing both types of analysis routines requiring either SE (single-end) reads
or PE reads. 

The following generates a fully populated _`systemPipeR`_ workflow environment 
(here for RNA-Seq) in the current working directory of an R session. At this time 
the package includes workflow templates for RNA-Seq, ChIP-Seq, VAR-Seq, and Ribo-Seq. 
Templates for additional NGS applications will be provided in the future.

## Directory Structure {#directory-structure}

[_`systemPipeRdata`_](http://bioconductor.org/packages/devel/data/experiment/html/systemPipeRdata.html), helper package, provides pre-configured workflows, reporting 
templates, and sample data loaded as demonstrated below. With a single command, 
the package allows creating the workflow environment containing the structure 
described here (see Figure \@ref(fig:dir)).

```{r generate_workenvir, eval=FALSE}
genWorkenvir(workflow="rnaseq")
setwd("rnaseq")
```

Directory names are indicated in <span style="color:grey">***green***</span>.
Users can change this structure as needed, but need to adjust the code in their 
workflows accordingly. 

* <span style="color:green">_**workflow/**_</span> (*e.g.* *myproject/*) 
    + This is the root directory of the R session running the workflow.
    + Run script ( *\*.Rmd*) and sample annotation (*targets.txt*) files are located here.
    + Note, this directory can have any name (*e.g.* <span style="color:green">_**myproject**_</span>). Changing its name does not require any modifications in the run script(s).
  + **Important subdirectories**: 
    + <span style="color:green">_**param/**_</span> 
        + <span style="color:green">_**param/cwl/**_</span>: This subdirectory stores all the parameter and configuration files. To organize workflows, each can have its own subdirectory, where all `*.cwl` and `*input.yml` files need to be in the same subdirectory. 
    + <span style="color:green">_**data/**_ </span>
        + Raw data (*e.g.* FASTQ files)
        + FASTA file of reference (*e.g.* reference genome)
        + Annotation files
        + Metadata
        + etc.
    + <span style="color:green">_**results/**_</span>
        + Analysis results are usually written to this directory, including: alignment, variant and peak files (BAM, VCF, BED); tabular result files; and image/plot files.
        + Note, the user has the option to organize results files for a given sample and analysis step in a separate subdirectory.

```{r dir, eval=TRUE, echo=FALSE, warning= FALSE, out.width="100%", fig.align = "center", fig.cap= "*systemPipeR's* preconfigured directory structure.", warning=FALSE}
knitr::include_graphics("images/spr_project.png")  
```

The following parameter files are included in each workflow template:

1. *`targets.txt`*: initial one provided by user; downstream *`targets_*.txt`* files are generated automatically
2. *`*.param/cwl`*: defines parameter for input/output file operations, *e.g.*:
    + *`hisat2/hisat2-mapping-se.cwl`* 
    + *`hisat2/hisat2-mapping-se.yml`*
3. Configuration files for computer cluster environments (skip on single machines):
    + *`.batchtools.conf.R`*: defines the type of scheduler for *`batchtools`* pointing to template file of cluster, and located in user's home directory
    + *`batchtools.*.tmpl`*: specifies parameters of scheduler used by a system, *e.g.* Torque, SGE, Slurm, etc.

## Structure of initial _`targets`_ file

The _`targets`_ file defines all input files (_e.g._ FASTQ, BAM, BCF) and sample 
comparisons of an analysis workflow. It can, also, store any number of descriptive 
information for each sample used in the workflow. 

The following shows the format of a sample _`targets`_ file included in the
package. It also can be viewed and downloaded 
from _`systemPipeR`'s_ GitHub repository [here](https://github.com/tgirke/systemPipeR/blob/master/inst/extdata/targets.txt). 
Please note that here it is represented a tabular file, however _`systemPipeR`_ can 
import the inputs information from a `YAML` files, as well as 
`SummarizedExperiment` object. For more details on how to create custom `targets`, 
please find here. 

Users should note here, the usage of targets files is optional when using
_`systemPipeR's`_ new workflow management interface. They can be replaced by a standard YAML
input file used by CWL. Since for organizing experimental variables targets
files are extremely useful and user-friendly. Thus, we encourage users to keep using 
them. 

### Structure of _`targets`_ file for single-end (SE) samples

In a target file with a single type of input files, here FASTQ files of 
single-end (SE) reads, the first column describe the path and the second column
represents a unique `id` name for each sample. The third column called `Factor` 
represents the biological replicates. All subsequent columns are additional 
information, and any number of extra columns can be added as needed.

```{r targetsSE, eval=TRUE}
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") 
showDF(read.delim(targetspath, comment.char = "#"))
```

To work with custom data, users need to generate a _`targets`_ file containing 
the paths to their own FASTQ files and then provide under _`targetspath`_ the
path to the corresponding _`targets`_ file. 

### Structure of _`targets`_ file for paired-end (PE) samples

For paired-end (PE) samples, the structure of the targets file is similar, where
users need to provide two FASTQ path columns: *`FileName1`* and *`FileName2`* 
with the paths to the PE FASTQ files. 

```{r targetsPE, eval=TRUE}
targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))
```

### Sample comparisons

Sample comparisons are defined in the header lines of the _`targets`_ file 
starting with '``# <CMP>``'. 

```{r comment_lines, echo=TRUE}
readLines(targetspath)[1:4]
```

The function _`readComp`_ imports the comparison information and stores it in a 
_`list`_. Alternatively, _`readComp`_ can obtain the comparison information from 
the corresponding _`SYSargsList`_ step (see below). Note, these header lines are 
optional. They are mainly useful for controlling comparative analyses according 
to certain biological expectations, such as identifying differentially expressed
genes in RNA-Seq experiments based on simple pair-wise comparisons.
 
```{r targetscomp, eval=TRUE}
readComp(file = targetspath, format = "vector", delim = "-")
```

## Downstream targets files description

After the step which required the initial targets file information, the downstream 
targets files are created automatically (see Figure \@ref(fig:targetsFig)).
Each step that uses the previous step outfiles as an input, the new targets input 
will be managed internally by the workflow instances, establishing connectivity 
among the steps in the workflow. 
_`systemPipeR`_ provides features to automatically and systematically build this 
connection, providing security that all the samples will be managed efficiently 
and reproducibly.

```{r targetsFig, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "_`systemPipeR`_ automatically creates the downstream `targets` files based on the previous steps outfiles. A) Usually, users provide the initial `targets` files, and this step will generate some outfiles, as demonstrated on B. Then, those files are used to build the new `targets` files as inputs in the next step. _`systemPipeR`_ (C) manages this connectivity among the steps automatically for the users.", warning=FALSE}
knitr::include_graphics("images/targets_con.png")  
```

## Structure of the new parameters files

The parameters and configuration required for running command-line software are 
provided by the widely used community standard [Common Workflow Language](https://www.commonwl.org/) (CWL) 
[@Amstutz2016-ka], which describes parameters analysis workflows in a generic 
and reproducible manner. For R-based workflow steps, param files are not required. 
For a complete overview of the CWL syntax, please see the [section](#cwl) below. 
Also, we have a dedicated section explain how to _`systemPipeR`_ establish the 
connection between the CWL parameters files and the targets files. Please see [here](#cwl_targets). 

```{r cleaning1, eval=TRUE, include=FALSE}
if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE)
## NOTE: Removing previous project create in the quick starts section
```

# Quick Start

This section will demonstrate how to build a basic workflow. The main features 
will be briefly illustrated here. The following section will discuss all the 
alternatives to design and build the workflow and the support features available 
in _`systemPipeR`_.

- Load sample data and directory structure

[*systemPipeRdata*](http://bioconductor.org/packages/release/data/experiment/html/systemPipeRdata.html) package is a helper package to generate a fully populated [*systemPipeR*](http://bioconductor.org/packages/release/bioc/html/systemPipeR.html)
workflow environment in the current working directory with a single command. 
All the instruction for generating the workflow are provide in the *systemPipeRdata* vignette [here](http://www.bioconductor.org/packages/devel/data/experiment/vignettes/systemPipeRdata/inst/doc/systemPipeRdata.html#1_Introduction). 

```{r genNew_wf, eval=TRUE}
systemPipeRdata::genWorkenvir(workflow = "new", mydirname = "spr_project")
```

```{r, eval=FALSE, warning=FALSE}
setwd("spr_project")
```

Typically, the user wants to record here the sources and versions of the
reference genome sequence along with the corresponding annotations. In
the provided sample data set all data inputs are stored in a `data`
subdirectory and all results will be written to a separate `results` directory,
while the `Rmarkdown` workflow file and the `targets` file are expected to be 
located in the parent directory. 

```{r setting_dir, include=FALSE, warning=FALSE}
knitr::opts_knit$set(root.dir = 'spr_project')
```

- Create a project

_`systemPipeR`_ workflows can be designed and built from start to finish with a 
single command, importing from an R Markdown file or stepwise in interactive 
mode from the R console. 

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

```{r SPRproject_ex, eval=TRUE}
sal <- SPRproject() 
```

- Build workflow from RMarkdown file

The workflow can be created by adding all the analysis steps interact or importing 
the steps from an R Markdown file.
After initializing the project with `SPRproject` function, `importWF` function
will scan and import all the R chunks from the R Markdown file and build all
the workflow instances. Then, each R chuck in the file will be converted in a workflow step. 

```{r importwf_example, eval=TRUE}
sal <- importWF(sal, 
                file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"),
                verbose = FALSE)
```

This workflow contains five steps. First, it loads _`systemPipeR`_  library, then 
exports `iris` data set, separating each species' data in a specific file. 
In the third step, each one of those files will be compressed using `gzip` 
command line, and later these files will be decompressed. Even though each step 
is a straightforward process, it is a simple demonstration of the connectivity 
between command-line and R-based analysis steps. Finally, the last step will
calculate a mean for sepal length and width and petal length and width, respectively,
and plot these statistics.

- 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 run_echo, eval=is_win}
sal <- runWF(sal)
```

- Visualize workflow

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

```{r plot_echo, eval=TRUE}
plotWF(sal, width = "80%", rstudio = TRUE)
```

- Checking workflow status

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

```{r status_echo, eval=TRUE}
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 logs_echo, eval=FALSE}
sal <- renderLogs(sal)
```

```{r cleaning2, eval=TRUE, include=FALSE}
if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE)
## NOTE: Removing previous project create in the quick starts section
```

# Project initialization

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

```{r SPRproject1, eval=TRUE}
sal <- SPRproject(projPath = getwd(), overwrite = TRUE) 
```

Internally, `SPRproject` function will create a hidden folder called `.SPRproject`,
by default, to store all the log files.
A `YAML` file, here called `SYSargsList.yml`, has been created, which initially
contains the basic location of the project structure; however, every time the 
workflow object `sal` is updated in R, the new information will also be store in this 
flat-file database for easy recovery.
If you desire different names for the logs folder and the `YAML` file, these can 
be modified as follows:

```{r SPRproject_logs, eval=FALSE}
sal <- SPRproject(logs.dir= ".SPRproject", sys.file=".SPRproject/SYSargsList.yml") 
```

Also, this function will check and/or create the basic folder structure if missing, 
which means `data`, `param`, and `results` folder, as described [here](#directory-structure). 
If the user wants to use a different names for these directories, can be specified 
as follows:

```{r SPRproject_dir, eval=FALSE}
sal <- SPRproject(data = "data", param = "param", results = "results") 
```

It is possible to separate all the R objects created within the workflow analysis 
from the current environment. `SPRproject` function provides the option to create 
a new environment, and in this way, it is not overwriting any object you may want
to have at your current section. 

```{r SPRproject_env, eval=FALSE}
sal <- SPRproject(envir = new.env()) 
```

In this stage, the object `sal` is a empty container, except for the project information. The project information can be accessed by the `projectInfo` method:

```{r projectInfo, eval=TRUE}
sal
projectInfo(sal)
```

Also, the `length` function will return how many steps this workflow contains, 
and in this case, it is empty, as follow:

```{r length, eval=TRUE}
length(sal)
```

# Workflow Design 

_`systemPipeR`_ workflows can be designed and built from start to finish with a single command, importing from an R Markdown file or stepwise in interactive mode from the R console. 
In the [next section](#appendstep), we will demonstrate how to build the workflow in an interactive mode, and in the [following section](#importWF), we will show how to build from a file. 

New workflows are constructed, or existing ones modified, by connecting each step
via `appendStep` method. Each `SYSargsList` instance contains instructions needed 
for processing a set of input files with a specific command-line and the paths to 
the corresponding outfiles generated.

The constructor function `Linewise` is used to build the R code-based step. 
For more details about this S4 class container, see [here](#linewise). 

## Build workflow interactive {#appendstep}

This tutorial shows a straightforward example for describing and explaining all main features available within systemPipeR to design, build, manage, run, and visualize the workflow. In summary, we are exporting a dataset to multiple files, compressing and decompressing each one of the files, importing to R, and finally performing a statistical analysis. 

In the previous section, we initialize the project by building the `sal` object.
Until this moment, the container has no steps:

```{r sal_check, eval=TRUE}
sal
```

Next, we need to populate the object created with the first step in the
workflow.

### Adding the first step 

The first step is R code based, and we are splitting the `iris` dataset by `Species`
and for each `Species` will be saved on file. Please note that this code will
not be executed now; it is just store in the container for further execution. 

This constructor function requires the `step_name` and the R-based code under 
the `code` argument. 
The R code should be enclosed by braces (`{}`) and separated by a new line. 

```{r, firstStep_R, eval=TRUE}
appendStep(sal) <- LineWise(code = {
                              mapply(function(x, y) write.csv(x, y),
                                     split(iris, factor(iris$Species)),
                                     file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv"))
                                     ) 
                            },
                            step_name = "export_iris")
```

For a brief overview of the workflow, we can check the object as follows:

```{r show, eval=TRUE}
sal
```

Also, for printing and double-check the R code in the step, we can use the 
`codeLine` method:

```{r codeLine, eval=TRUE}
codeLine(sal)
```

### Adding more steps

Next, an example of how to compress the exported files using 
[`gzip`](https://www.gnu.org/software/gzip/) command-line. 

The constructor function creates an `SYSargsList` S4 class object using data from
three input files:

    - CWL command-line specification file (`wf_file` argument);
    - Input variables (`input_file` argument);
    - Targets file (`targets` argument).

In CWL, files with the extension `.cwl` define the parameters of a chosen
command-line step or workflow, while files with the extension `.yml` define the
input variables of command-line steps. 

The `targets` file is optional for workflow steps lacking `input` files. The connection 
between `input` variables and the `targets` file is defined under the `inputvars` 
argument. It is required a `named vector`, where each element name needs to match
with column names in the `targets` file, and the value must match the names of 
the `input` variables defined in the `*.yml` files (see Figure \@ref(fig:sprandCWL)). 

A detailed description of the dynamic between `input` variables and `targets` 
files can be found [here](#cwl_targets). 
In addition, the CWL syntax overview can be found [here](#cwl). 

Besides all the data form `targets`, `wf_file`, `input_file` and `dir_path` arguments,
`SYSargsList` constructor function options include: 

  - `step_name`: a unique *name* for the step. This is not mandatory; however, 
    it is highly recommended. If no name is provided, a default `step_x`, where
    `x` reflects the step index, will be added. 
  - `dir`: this option allows creating an exclusive subdirectory for the step 
    in the workflow. All the outfiles and log files for this particular step will 
    be generated in the respective folders.
  - `dependency`: after the first step, all the additional steps appended to 
    the workflow require the information of the dependency tree. 

The `appendStep<-` method is used to append a new step in the workflow.

```{r gzip_secondStep, eval=TRUE}
targetspath <- system.file("extdata/cwl/gunzip", "targets_gunzip.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "gzip", 
                      targets = targetspath, dir = TRUE,
                      wf_file = "gunzip/workflow_gzip.cwl", input_file = "gunzip/gzip.yml",
                      dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                      inputvars = c(FileName = "_FILE_PATH_", SampleName = "_SampleName_"), 
                      dependency = "export_iris")
```

Note: This will not work if the `gzip` is not available on your system 
(installed and exported to PATH) and may only work on Windows systems using PowerShell. 

For a overview of the workflow, we can check the object as follows:

```{r}
sal
```

Note that we have two steps, and it is expected three files from the second step.
Also, the workflow status is *Pending*, which means the workflow object is 
rendered in R; however, we did not execute the workflow yet. 
In addition to this summary, it can be observed this step has three command lines. 

For more details about the command-line rendered for each target file, it can be 
checked as follows: 

```{r}
cmdlist(sal, step = "gzip")
```

#### Using the `outfiles` for the next step

For building this step, all the previous procedures are being used to append the 
next step. However, here, we can observe power features that build the 
connectivity between steps in the workflow.

In this example, we would like to use the outfiles from *gzip* Step, as
input from the next step, which is the *gunzip*. In this case, let's look at the 
outfiles from the first step:

```{r}
outfiles(sal)
```

The column we want to use is "gzip_file". For the argument `targets` in the 
`SYSargsList` function, it should provide the name of the correspondent step in
the Workflow and which `outfiles` you would like to be incorporated in the next 
step. 
The argument `inputvars` allows the connectivity between `outfiles` and the 
new `targets` file. Here, the name of the previous `outfiles` should be provided 
it. Please note that all `outfiles` column names must be unique.

It is possible to keep all the original columns from the `targets` files or remove
some columns for a clean `targets` file.
The argument `rm_targets_col` provides this flexibility, where it is possible to
specify the names of the columns that should be removed. If no names are passing
here, the new columns will be appended. 

```{r gunzip, eval=TRUE}
appendStep(sal) <- SYSargsList(step_name = "gunzip", 
                      targets = "gzip", dir = TRUE,
                      wf_file = "gunzip/workflow_gunzip.cwl", input_file = "gunzip/gunzip.yml",
                      dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                      inputvars = c(gzip_file = "_FILE_PATH_", SampleName = "_SampleName_"), 
                      rm_targets_col = "FileName", 
                      dependency = "gzip")
```

We can check the targets automatically create for this step, 
based on the previous `outfiles`:

```{r targetsWF_3, eval=TRUE}
targetsWF(sal[3])
```

We can also check all the expected `outfiles` for this particular step, as follows:

```{r outfiles_2, eval=TRUE}
outfiles(sal[3])
```

Now, we can observe that the third step has been added and contains one substep.

```{r}
sal
```

In addition, we can access all the command lines for each one of the substeps. 

```{r, eval=TRUE}
cmdlist(sal["gzip"], targets = 1)
```

#### Getting data from a workflow instance 

The final step in this simple workflow is an R code step. For that, we are using
the `LineWise` constructor function as demonstrated above. 

One interesting feature showed here is the `getColumn` method that allows 
extracting the information for a workflow instance. Those files can be used in
an R code, as demonstrated below. 

```{r getColumn, eval=TRUE}
getColumn(sal, step = "gunzip", 'outfiles')
```

```{r, iris_stats, eval=TRUE}
appendStep(sal) <- LineWise(code = {
                    df <- lapply(getColumn(sal, step = "gunzip", 'outfiles'), function(x) read.delim(x, sep = ",")[-1])
                    df <- do.call(rbind, df)
                    stats <- data.frame(cbind(mean = apply(df[,1:4], 2, mean), sd = apply(df[,1:4], 2, sd)))
                    stats$species <- rownames(stats)
                    
                    plot <- ggplot2::ggplot(stats, ggplot2::aes(x = species, y = mean, fill = species)) + 
                      ggplot2::geom_bar(stat = "identity", color = "black", position = ggplot2::position_dodge()) +
                      ggplot2::geom_errorbar(ggplot2::aes(ymin = mean-sd, ymax = mean+sd), width = .2, position = ggplot2::position_dodge(.9)) 
                    },
                    step_name = "iris_stats", 
                    dependency = "gzip")
```

## Build workflow from a {R Markdown} {#importWF}

The precisely same workflow can be created by importing the steps from an 
R Markdown file.
As demonstrated above, it is required to initialize the project with `SPRproject` function. 

`importWF` function will scan and import all the R chunk from the R Markdown file 
and build all the workflow instances. Then, each R chuck in the file will be 
converted in a workflow step. 

```{r importWF_rmd, eval=TRUE}
sal_rmd <- SPRproject(logs.dir = ".SPRproject_rmd") 

sal_rmd <- importWF(sal_rmd, 
                file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"))
```

Let's explore the workflow to check the steps:

```{r importWF_details}
stepsWF(sal_rmd)
dependency(sal_rmd)
codeLine(sal_rmd)
targetsWF(sal_rmd)
```

### Rules to create the R Markdown to import as workflow

To include a particular code chunk from the R Markdown file in the workflow
analysis, please use the following code chunk options:

    - `spr=TRUE'`: for code chunks with step workflow. 

For example: 

> *```{r step_1, eval=TRUE, spr=TRUE}*

> *```{r step_2, eval=FALSE, spr=TRUE}*

`importWF` function can ignore `eval` option in code chunk, and in this case, 
both of the examples steps above would be incorporated in the workflow. 

For `spr = TRUE`, the last object assigned must to be the `SYSargsList`, for example:

```{r fromFile_example_rules_cmd, eval=FALSE}
targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "Example", 
                      targets = targetspath, 
                      wf_file = "example/example.cwl", input_file = "example/example.yml", 
                      dir_path = system.file("extdata/cwl", package = "systemPipeR"), 
                      inputvars = c(Message = "_STRING_", SampleName = "_SAMPLE_"))
```

OR 

```{r fromFile_example_rules_r, eval=FALSE}
appendStep(sal) <- LineWise(code = {
                              library(systemPipeR)
                            },
                            step_name = "load_lib")
```

Also, note that all the required files or objects to generate one particular step 
must be defined in an imported R code. 
The motivation for this is that when R Markdown files are imported, the `spr = TRUE` 
flag will be evaluated, append, and stored in the workflow control class as 
the `SYSargsList` object. 

The workflow object name used in the R Markdown (e.g. `appendStep(sal)`) needs to 
be the same when the `importWF` function is used. It is important to keep consistency. 
If different object names are used, when running the workflow, you can see a error, 
like *Error: <objectName> object not found.*.

# Running the workflow

For running the workflow, `runWF` function will execute all the command lines 
store in the workflow container.

```{r runWF, eval=is_win}
sal <- runWF(sal)
```

This essential function allows the user to choose one or multiple steps to be 
executed using the `steps` argument. However, it is necessary to follow the 
workflow dependency graph. If a selected step depends on a previous step(s) that
was not executed, the execution will fail. 

```{r runWF_error, eval=FALSE}
sal <- runWF(sal, steps = c(1,3))
```

Also, it allows forcing the execution of the steps, even if the status of the 
step is `'Success'` and all the expected `outfiles` exists. 
Another feature of the `runWF` function is ignoring all the warnings 
and errors and running the workflow by the arguments `warning.stop` and 
`error.stop`, respectively.

```{r runWF_force, eval=FALSE}
sal <- runWF(sal, force = TRUE, warning.stop = FALSE, error.stop = TRUE)
```

When the project was initialized by `SPRproject` function, it was created an 
environment for all objects created during the workflow execution. This 
environment can be accessed as follows:

```{r runWF_env, eval=FALSE}
viewEnvir(sal)
```

The workflow execution allows to save this environment for future recovery:

```{r runWF_saveenv, eval=FALSE}
sal <- runWF(sal, saveEnv = TRUE)
```

## Workflow status

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

```{r show_statusWF, eval=TRUE}
sal
```

To access more details about the workflow instances, we can use the `statusWF` method:

```{r statusWF, eval=TRUE}
statusWF(sal)
```

## Parallelization on clusters

This section of the tutorial provides an introduction to the usage of the 
_`systemPipeR`_ features on a cluster.

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}
resources <- list(conffile=".batchtools.conf.R",
                  template="batchtools.slurm.tmpl", 
                  Njobs=18, 
                  walltime=120, ## minutes
                  ntasks=1,
                  ncpus=4, 
                  memory=1024, ## Mb
                  partition = "short"
                  )
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)
```

Note: The example is submitting the jog to `short` partition. If you desire to 
use a different partition, please adjust accordingly.

# Visualize workflow

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

This function will make a plot of selected workflow instance and the following 
information is displayed on the plot:

    - Workflow structure (dependency graphs between different steps); 
    - Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`; 
    - Sample status and statistics; 
    - Workflow timing: running duration time. 

If no argument is provided, the basic plot will automatically detect width, 
height, layout, plot method, branches, _etc_. 

```{r, eval=TRUE}
plotWF(sal, show_legend = TRUE, width = "80%", rstudio = TRUE)
```

For more details about the `plotWF` function, please see [here](#plotWF).

# Technical 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`.
Also, the workflow plot is appended at the beginning of the report, making it 
easier to click on the respective step.

```{r, eval=FALSE}
sal <- renderLogs(sal)
```

# Exported the workflow 

_`systemPipeR`_ workflow management system allows to translate and export the 
workflow build interactively to R Markdown format or an executable bash script.
This feature advances the reusability of the workflow, as well as the flexibility
for workflow execution.

## R Markdown file

`sal2rmd` function takes an `SYSargsList` workflow container and translates it to 
SPR workflow template R markdown format. This file can be imported with the 
`importWF` function, as demonstrated above.

```{r, eval=FALSE}
sal2rmd(sal)
```

## Bash script

`sal2bash` function takes an `SYSargsList` workflow container and translates 
it to an executable bash script, so one can run the workflow without loading 
`SPR` or using an R console.


```{r, eval=FALSE}
sal2bash(sal)
```

It will be generated on the project root an executable bash script, called by
default the `spr_wf.sh`. Also, a directory `./spr_wf` will be created and store 
all the R scripts based on the workflow steps. Please note that this function will 
"collapse" adjacent R steps into one file as much as possible.

# Project Resume and Restart

If you desire to resume or restart a project that has been initialized in the past, 
`SPRproject` function allows this operation.

With the `resume` option, it is possible to load the `SYSargsList` object in R and 
resume the analysis. Please, make sure to provide the `logs.dir` location, and the 
corresponded `YAML` file name, if the default names were not used when the project was created.

```{r SPR_resume, eval=FALSE}
sal <- SPRproject(resume = TRUE, logs.dir = ".SPRproject", 
                  sys.file = ".SPRproject/SYSargsList.yml") 
```

If you choose to save the environment in the last analysis, you can recover all 
the files created in that particular section. `SPRproject` function allows this 
with `load.envir` argument. Please note that the environment was saved only with
you run the workflow in the last section (`runWF()`).

```{r resume_load, eval=FALSE}
sal <- SPRproject(resume = TRUE, load.envir = TRUE) 
```

After loading the workflow at your current section, you can check the objects 
created in the old environment and decide if it is necessary to copy them to the
current environment.

```{r envir, eval=FALSE}
viewEnvir(sal)
copyEnvir(sal, list="plot", new.env = globalenv())
```

The `resume` option will keep all previous logs in the folder; however, if you desire to 
clean the execution (delete all the log files) history and restart the workflow, 
the `restart=TRUE` option can be used.

```{r restart_load, eval=FALSE}
sal <- SPRproject(restart = TRUE, load.envir = FALSE) 
```

The last and more drastic option from `SYSproject` function is to `overwrite` the
logs and the `SYSargsList` object. This option will delete the hidden folder and the 
information on the `SYSargsList.yml` file. This will not delete any parameter
file nor any results it was created in previous runs. Please use with caution. 

```{r SPR_overwrite, eval=FALSE}
sal <- SPRproject(overwrite = TRUE) 
```

# Exploring workflow instances {#sysargslist}

_`systemPipeR`_ provide several accessor methods and useful functions to explore 
`SYSargsList` workflow object.

## Accessor Methods

Several accessor methods are available that are named after the slot names of
the `SYSargsList` workflow object.

```{r}
names(sal)
```

- Check the length of the workflow:

```{r}
length(sal)
```

- Check the steps of the workflow:

```{r}
stepsWF(sal)
```

- Checking the command-line for each target sample:

`cmdlist()` method printing the system commands for running command-line
software as specified by a given `*.cwl` file combined with the paths to the
input samples (*e.g.* FASTQ files) provided by a `targets` file. The example below
shows the `cmdlist()` output for running `gzip` and `gunzip` on the first sample. 
Evaluating the output of `cmdlist()` can be very helpful for designing
and debugging `*.cwl` files of new command-line software or changing the
parameter settings of existing ones.

```{r}
cmdlist(sal, step = c(2,3), targets = 1)
```

- Check the workflow status:

```{r}
statusWF(sal)
```

- Check the workflow targets files:

```{r}
targetsWF(sal[2])
```

- Checking the expected outfiles files:

The `outfiles` components of `SYSargsList` define the expected outfiles files 
for each step in the workflow, some of which are the input for the next workflow step.

```{r}
outfiles(sal[2])
```

- Check the workflow dependencies:

```{r}
dependency(sal)
```

- Check the sample comparisons:

Sample comparisons are defined in the header lines of the `targets` file 
starting with '``# <CMP>``'. This information can be accessed as follows:

```{r, eval=FALSE}
targetsheader(sal, step = "Quality")
```

- Get the workflow steps names: 

```{r}
stepName(sal)
```

- Get the Sample Id for on particular step:

```{r}
SampleName(sal, step = "gzip")
SampleName(sal, step = "iris_stats")
```

- Get the `outfiles` or `targets` column files: 

```{r}
getColumn(sal, "outfiles", step = "gzip", column = "gzip_file")
getColumn(sal, "targetsWF", step = "gzip", column = "FileName")
```

- Get the R code for a `LineWise` step:

```{r}
codeLine(sal, step = "export_iris")
```

- View all the objects in the running environment:

```{r}
viewEnvir(sal)
```

- Copy one or multiple objects from the running environment to a new environment:

```{r}
copyEnvir(sal, list = c("plot"), new.env = globalenv(), silent = FALSE)
```

- Accessing the `*.yml` data

```{r}
yamlinput(sal, step = "gzip")
```

## Subsetting the workflow details

- The `SYSargsList` class and its subsetting operator `[`:

```{r}
sal[1]
sal[1:3]
sal[c(1,3)]
```

- The `SYSargsList` class and its subsetting by steps and input samples: 

```{r}
sal_sub <- subset(sal, subset_steps = c( 2,3), input_targets = ("SE"), keep_steps = TRUE)
stepsWF(sal_sub)
targetsWF(sal_sub)
outfiles(sal_sub)
```

- The `SYSargsList` class and its operator `+`

```{r, eval=FALSE}
sal[1] + sal[2] + sal[3]
```

## Replacement Methods

- Update a `input` parameter in the workflow

```{r, eval=TRUE}
sal_c <- sal
## check values
yamlinput(sal_c, step = "gzip")
## check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
## Replace
yamlinput(sal_c, step = "gzip", paramName = "ext") <- "txt.gz"

## check NEW values
yamlinput(sal_c, step = "gzip")
## Check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
```

- Append and Replacement methods for R Code Steps

```{r, sal_lw_rep, eval=TRUE}
appendCodeLine(sal_c, step = "export_iris", after = 1) <- "log_cal_100 <- log(100)"
codeLine(sal_c, step = "export_iris")

replaceCodeLine(sal_c, step="export_iris", line = 2) <- LineWise(code={
                    log_cal_100 <- log(50)
                    })
codeLine(sal_c, step = 1)
```

For more details about the `LineWise` class, please see [below](#linewise).

- Rename a Step

```{r}
renameStep(sal_c, step = 1) <- "newStep"
renameStep(sal_c, c(1, 2)) <- c("newStep2", "newIndex")
sal_c
names(outfiles(sal_c))
names(targetsWF(sal_c))
dependency(sal_c)
```

- Replace a Step

```{r, eval=FALSE}
sal_test <- sal[c(1,2)]
replaceStep(sal_test, step = 1, step_name = "gunzip" ) <- sal[3]
sal_test
```

Note: Please use this method with attention, because it can disrupt all 
the dependency graphs. 

- Removing a Step

```{r}
sal_test <- sal[-2]
sal_test
```

# CWL syntax {#cwl}

This section will introduce how CWL describes command-line tools and the 
specification and terminology of each file. 
For complete documentation, please check the CommandLineTools documentation [here](https://www.commonwl.org/v1.2/CommandLineTool.html) 
and [here](https://www.commonwl.org/v1.2/Workflow.html) for Workflows and the user guide [here](https://www.commonwl.org/user_guide/). 

CWL command-line specifications are written in [YAML](http://yaml.org/) format.

In CWL, files with the extension `.cwl` define the parameters of a chosen 
command-line step or workflow, while files with the extension `.yml` define 
the input variables of command-line steps. 

## CWL `CommandLineTool`

`CommandLineTool` by CWL definition is a standalone process, with no interaction 
if other programs, execute a program, and produce output. 

Let's explore the `*.cwl` file: 

```{r}
dir_path <- system.file("extdata/cwl", package = "systemPipeR")
cwl <- yaml::read_yaml(file.path(dir_path, "example/example.cwl"))
```

- The `cwlVersion` component shows the CWL specification version used by the document. 
- The `class` component shows this document describes a `CommandLineTool.`
Note that CWL has another `class`, called `Workflow` which represents a union of one 
or more command-line tools together. 

```{r}
cwl[1:2]
```

- `baseCommand` component provides the name of the software that we desire to execute.

```{r}
cwl[3]
```

- The `inputs` section provides the input information to run the tool. Important 
components of this section are: 
    - `id`: each input has an id describing the input name;
    - `type`: describe the type of input value (string, int, long, float, double, 
    File, Directory or Any);
    - `inputBinding`: It is optional. This component indicates if the input 
    parameter should appear on the command-line. If this component is missing 
    when describing an input parameter, it will not appear in the command-line 
    but can be used to build the command-line.

```{r}
cwl[4]
```

- The `outputs` section should provide a list of the expected outputs after running the command-line tools. Important 
components of this section are: 
    - `id`: each input has an id describing the output name;
    - `type`: describe the type of output value (string, int, long, float, double, 
    File, Directory, Any or `stdout`);
    - `outputBinding`: This component defines how to set the outputs values. The `glob` component will define the name of the output value. 

```{r}
cwl[5]
```

- `stdout`: component to specify a `filename` to capture standard output.
Note here we are using a syntax that takes advantage of the inputs section, 
using results_path parameter and also the `SampleName` to construct the output `filename.` 

```{r}
cwl[6]
```

## CWL `Workflow`

`Workflow` class in CWL is defined by multiple process steps, where can have 
interdependencies between the steps, and the output for one step can be used as 
input in the further steps.

```{r}
cwl.wf <- yaml::read_yaml(file.path(dir_path, "example/workflow_example.cwl"))
```

- The `cwlVersion` component shows the CWL specification version used by the document. 
- The `class` component shows this document describes a `Workflow`.

```{r}
cwl.wf[1:2]
```

- The `inputs` section describes the inputs of the workflow.

```{r}
cwl.wf[3]
```

- The `outputs` section describes the outputs of the workflow.

```{r}
cwl.wf[4]
```

- The `steps` section describes the steps of the workflow. In this simple example, 
we demonstrate one step.

```{r}
cwl.wf[5]
```

## CWL Input Parameter 

Next, let's explore the *.yml* file, which provide the input parameter values for all
the components we describe above. 

For this simple example, we have three parameters defined:

```{r}
yaml::read_yaml(file.path(dir_path, "example/example_single.yml"))
```

Note that if we define an input component in the *.cwl* file, this value needs 
to be also defined here in the *.yml* file. 

# How to connect CWL description files within _`systemPipeR`_ {#cwl_targets}

This section will demonstrate how to connect CWL parameters files to create 
workflows. In addition, we will show how the workflow can be easily scalable 
with _`systemPipeR`_.

`SYSargsList` container stores all the information and instructions needed for processing 
a set of input files with a single or many command-line steps within a workflow 
(i.e. several components of the software or several independent software tools). 
The `SYSargsList` object is created and fully populated with the `SYSargsList` construct
function. 
Full documentation of `SYSargsList` management instances can be found [here](#sysargslist)
and [here](#appendstep).

The following imports a `.cwl` file (here `example.cwl`) for running the `echo Hello World!` 
example.

```{r fromFile, eval=TRUE}
HW <- SYSargsList(wf_file = "example/workflow_example.cwl", 
                  input_file = "example/example_single.yml", 
                  dir_path = system.file("extdata/cwl", package = "systemPipeR"))
HW
cmdlist(HW)
```

However, we are limited to run just one command-line or one sample in this example. 
To scale the command-line over many samples, a simple solution offered by `systemPipeR` 
is to provide a `variable` for each of the parameters that we want to run with multiple samples. 

Let's explore the example:

```{r}
yml <- yaml::read_yaml(file.path(dir_path, "example/example.yml"))
yml
```

For the `message` and `SampleName` parameter, we are passing a variable connecting 
with a third file called `targets.` 

Now, let's explore the `targets` file structure:

```{r}
targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR")
read.delim(targetspath, comment.char = "#")
```

The `targets` file defines all input files or values and sample ids of an analysis workflow. 
For this example, we have defined a string message for the `echo` command-line tool,
in the first column that will be evaluated, and the second column is the 
`SampleName` id for each one of the messages.
Any number of additional columns can be added as needed.

Users should note here, the usage of `targets` files is optional when using 
`systemPipeR's` new CWL interface. Since for organizing experimental variables targets 
files are extremely useful and user-friendly. Thus, we encourage users to keep using them.

### How to connect the parameter files and `targets` file information?

The constructor function creates an `SYSargsList` S4 class object connecting three input files:

  - CWL command-line specification file (`wf_file` argument);
  - Input variables (`input_file` argument);
  - Targets file (`targets` argument).
    
As demonstrated above, the latter is optional for workflow steps lacking input files. 
The connection between input variables (here defined by `input_file` argument) 
and the `targets` file are defined under the `inputvars` argument. 
A named vector is required, where each element name needs to match with column 
names in the `targets` file, and the value must match the names of the *.yml* 
variables. This is used to replace the CWL variable and construct all the command-line
for that particular step. 

The variable pattern `_XXXX_` is used to distinguish CWL variables that target 
columns will replace. This pattern is recommended for consistency and easy identification
but not enforced.

The following imports a `.cwl` file (same example demonstrated above) for running
the `echo Hello World` example. However, now we are connecting the variable defined 
on the `.yml` file with the `targets` file inputs.

```{r fromFile_example, eval=TRUE}
HW_mul <- SYSargsList(step_name = "echo", 
                      targets=targetspath, 
                      wf_file="example/workflow_example.cwl", input_file="example/example.yml", 
                      dir_path = dir_path, 
                      inputvars = c(Message = "_STRING_", SampleName = "_SAMPLE_"))
HW_mul
cmdlist(HW_mul)
```

# Creating the CWL param files from the command-line

Users need to define the command-line in a pseudo-bash script format:

```{r cmd, eval=TRUE}
command <- "
    hisat2 \
    -S <F, out: ./results/M1A.sam> \
    -x <F: ./data/tair10.fasta> \
     -k <int: 1> \
    -min-intronlen <int: 30> \
    -max-intronlen <int: 3000> \
    -threads <int: 4> \
    -U <F: ./data/SRR446027_1.fastq.gz>
"
```

## Define prefix and defaults

- First line is the base command. Each line is an argument with its default value.

- For argument lines (starting from the second line), any word before the first 
  space with leading `-` or `--` in each will be treated as a prefix, like `-S` or 
  `--min`. Any line without this first word will be treated as no prefix. 
  
- All defaults are placed inside `<...>`.

- First argument is the input argument type. `F` for "File", "int", "string" are unchanged.

- Optional: use the keyword `out` followed the type with a `,` comma separation to 
  indicate if this argument is also an CWL output.
  
- Then, use `:` to separate keywords and default values, any non-space value after the `:`
  will be treated as the default value. 
  
- If any argument has no default value, just a flag, like `--verbose`, there is no need to add any `<...>`

## `createParam` Function

`createParam` function requires the `string` as defined above as an input. 

First of all, the function will print the three components of the `cwl` file:
    - `BaseCommand`: Specifies the program to execute. 
    - `Inputs`: Defines the input parameters of the process.
    - `Outputs`: Defines the parameters representing the output of the process.
    
The four component is the original command-line.

If in interactive mode, the function will verify that everything is correct and 
will ask you to proceed. Here, the user can answer "no" and provide more 
information at the string level. Another question is to save the param created here.

If running the workflow in non-interactive mode, the `createParam` function will 
consider "yes" and returning the container.

```{r}
cmd <- createParam(command, writeParamFiles = FALSE)
```

If the user chooses not to save the `param` files on the above operation, 
it can use the `writeParamFiles` function.

```{r saving, eval=FALSE}
writeParamFiles(cmd, overwrite = TRUE)
```

## How to access and edit param files

### Print a component

```{r}
printParam(cmd, position = "baseCommand") ## Print a baseCommand section
printParam(cmd, position = "outputs")
printParam(cmd, position = "inputs", index = 1:2) ## Print by index
printParam(cmd, position = "inputs", index = -1:-2) ## Negative indexing printing to exclude certain indices in a position
```

### Subsetting the command-line

```{r}
cmd2 <- subsetParam(cmd, position = "inputs", index = 1:2, trim = TRUE)
cmdlist(cmd2)

cmd2 <- subsetParam(cmd, position = "inputs", index = c("S", "x"), trim = TRUE)
cmdlist(cmd2)
```

### Replacing a existing argument in the command-line

```{r}
cmd3 <- replaceParam(cmd, "base", index = 1, replace = list(baseCommand = "bwa"))
cmdlist(cmd3)
```

```{r}
new_inputs <- new_inputs <- list(
    "new_input1" = list(type = "File", preF="-b", yml ="myfile"),
    "new_input2" = "-L <int: 4>"
)
cmd4 <- replaceParam(cmd, "inputs", index = 1:2, replace = new_inputs)
cmdlist(cmd4)
```

### Adding new arguments

```{r}
newIn <- new_inputs <- list(
    "new_input1" = list(type = "File", preF="-b1", yml ="myfile1"),
    "new_input2" = list(type = "File", preF="-b2", yml ="myfile2"),
    "new_input3" = "-b3 <F: myfile3>"
)
cmd5 <- appendParam(cmd, "inputs", index = 1:2, append = new_inputs)
cmdlist(cmd5)

cmd6 <- appendParam(cmd, "inputs", index = 1:2, after=0, append = new_inputs)
cmdlist(cmd6)
```

### Editing `output` param

```{r}
new_outs <- list(
    "sam_out" = "<F: $(inputs.results_path)/test.sam>"
) 
cmd7 <- replaceParam(cmd, "outputs", index = 1, replace = new_outs)
output(cmd7) 
```

### Internal Check

```{r sysargs2, eval=TRUE}
cmd <- "
    hisat2 \
    -S <F, out: _SampleName_.sam> \
    -x <F: ./data/tair10.fasta> \
    -k <int: 1> \
    -min-intronlen <int: 30> \
    -max-intronlen <int: 3000> \
    -threads <int: 4> \
    -U <F: _FASTQ_PATH1_>
"
WF <- createParam(cmd, overwrite = TRUE, writeParamFiles = TRUE, confirm = TRUE) 
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
WF_test <- loadWorkflow(targets = targetspath, wf_file="hisat2.cwl",
                   input_file="hisat2.yml", dir_path = "param/cwl/hisat2/")
WF_test <- renderWF(WF_test, inputvars = c(FileName = "_FASTQ_PATH1_"))
WF_test
cmdlist(WF_test)[1:2]
```

# Visualize Workflow - Full details {#plotWF}

## Color and text

On the plot, different colors and numbers indicate different status. This information 
can be found also in the plot legends.

**Shapes:**

- circular steps: pure R code steps
- rounded squares steps: `sysargs` steps, steps that will invoke command-line calls
- blue colored steps and arrows: main branch (see [main branch](#main-branch) section below)

**Step colors**

- black: pending steps
- <span class="text-success">Green</span>: successful steps
- <span class="text-danger">Green</span>: failed steps

**Number and colors**

<br>There are 4 numbers in the second row of each step, separated by `/`

- <span class="text-success">First No.</span>: number of passed samples
- <span class="text-warning">Second No.</span>: number of warning samples
- <span class="text-danger">Third No.</span>: number of erros samples
- <span style="color:blue;">Forth No.</span>: number of total samples

**Duration**

<br>This is shown after the sample information, as how long it took to run this step. 
Units are a few seconds (**s**), some minutes (**m**), or some hours (**h**). 

## on hover

When the mouse is hovering on each step, detailed information will be displayed.

## logging

The workflow steps will also become clickable if `in_log = TRUE`. This will create links 
for each step that navigate to corresponding log section in the SPR 
[workflow log file](change to page that introduce the log file). Normally this option 
is handled by SPR log file generating function to create this plot on top of the log file,
so when a certain step is click, it will navigate to the detailed section down the page. 
Here is only an example to demo how the plot can be clickable (will not navigate you to 
anywhere). Visit [this page](link) to see a real example. 

```{r}
plotWF(sal, in_log = TRUE)
```

## Plot Method

The default plotting method is `svg`. It means the plot is generated by `svg` embedding. 
Sometimes certain browsers may not display `svg` correctly. In this case, the other option 
is to use `png` to embed the plot. However, you will **lose hovering, clicking and some** 
**responsiveness** (plot auto resizing ability) of the plot. 

```{r}
plotWF(sal, plot_method = "png")
```

## Rstudio

By default, even if you are working inside Rstudio the plot is **not displayed in Rstudio viewer**. 
This is because the workflow steps will be too small inside Rstudio viewer too 
see the details. We recommend to view it in a larger space, so by default it will 
open up your web browser to display it. You can enforce `rstudio = TRUE` to see it 
in Rstudio Viewer. 

```{r}
plotWF(sal, rstudio = TRUE)
```

## Responsiveness

This is a term often used in web development. It means will the plot resize itself 
if the user resize the document window? By default, `plotWF` will be responsive, 
meaning it will fit current window container size and adjust the size once the window 
size has changed. To always display the full sized plot, use `responsive = FALSE`,
useful for embedding the plot in a full-screen mode.

<div style="height: 500px">
```{r}
plotWF(sal, responsive = FALSE)
```
</div>

For the plot above, you need to scroll to see the plot.

## Layout

There a few different layout you can choose. There is no best layout. It all depends 
on the workflow structure you have. The default is `compact` but we recommend you 
to try different layouts to find the best fitting one. 

- `compact`: try to plot steps as close as possible.
- `vertical`: main branch will be placed vertically and side branches will be placed
  on the same horizontal level and sub steps of side branches will be placed
  vertically.
- `horizontal`: main branch is placed horizontally and side branches and sub
  steps will be placed vertically.
- `execution`: a linear plot to show the workflow execution order of all steps.

**vertical**

```{r}
plotWF(sal, layout = "vertical", height = "600px")
```
The plot is very long, use `height` to make it smaller.

**horizontal**

```{r}
plotWF(sal, layout = "horizontal")
```

**execution**

```{r}
plotWF(sal, layout = "execution", height = "600px", responsive = FALSE)
```

The plot is very long but if we use `height` to limit to a smaller size, details are 
hard to see. Then it will be good to use `height` and `responsive = FALSE` together. 

## Main branch 

From the plots above, you can see that there are many steps which do not connect to any 
other steps in downstream. These dead-ends are called ending steps. If we connect the first step,
steps in between and these ending step, this will become a branch. Imagine the workflow is 
a top-bottom tree structure and the root is the first step. Therefore, there are 
many possible ways to connect the workflow. For the convenience of plotting, we 
introduce a concept of _"main branch"_, meaning one of the possible connecting 
strategies that will be placed at the center of the plot. Other steps that are not 
in this major branch will surround this major space. This  _"main branch"_ will not 
affect how a workflow is run, but an algorithm to **compute the best visualization**. 

This main branch will not impact the `compact` layout so much but will have a huge 
effect on `horizontal` and `vertical` layouts. 

The algorithm in `plotWF` will automatically choose a best branch for 
you by default. In simple words, it favors: (a). branches that connect first and last step;
(b). as long as possible. 

You can also choose a branch you want by `branch_method = "choose"`. It will first 
list all possible branches, and then give you a prompt to ask for your favorite branch. 
Here, for rendering the Rmarkdown, we cannot have a prompt, so we use a second argument 
in combination, `branch_no = x` to directly choose a branch and skip the prompt. Also, 
we use the `verbose = TRUE` to imitate the branch listing in console. In a real case, 
you only need `branch_method = "choose"`.

To have the main branch marked, `mark_main_branch = TRUE` must be added (default `FALSE`). 
Watch closely how the plot change by choosing different branches. Here we use `vertical`
layout to demo. Remember, the main branch is marked in blue.

```{r collapse=TRUE}
plotWF(sal,  mark_main_branch = TRUE, layout = "vertical", branch_method = "choose", branch_no = 1, verbose = FALSE)
```

## Legends 

The legend can also be removed by `show_legend = FALSE`

```{r}
plotWF(sal, show_legend = FALSE, height = "500px")
```

## Output formats 

There are current three output formats: `"html"` and `"dot"`, `"dot_print"`. If first 
two were chosen, you also need provide a path `out_path` to save the file. 

- html: a single html file contains the plot.
- dot: a DOT script file with the code to reproduce the plot in a [graphiz](https://graphviz.org/)
  DOT engine. 
- dot_print: directly cat the dot script to console. 

### HTML
HTML format is very useful if you want to view the plot later or share it to other 
people. This format is also helpful when you are working on a remote computer 
cluster. To view the workflow plot, a browser device (viewer) must be available, 
but often time this is not the case for computer clusters. When you plot a workflow 
and see the message `"Couldn't get a file descriptor referring to the console"`,
it means your computer (cluster) does not have a browser device. Saving to 
HTML format is the best option. 

```{r}
plotWF(sal, out_format = "html", out_path = "example_out.html")
file.exists("example_out.html")
```

### DOT
Saving workflow plot to `dot` format allows one to reproduce the plot with 
the Graphviz language. 
```{r}
plotWF(sal, out_format = "dot", out_path = "example_out.dot")
cat(readLines("example_out.dot")[1:5], sep = "\n")
```

### DOT print
Instead of saving the Graphviz plotting code to a file, this option directly 
prints out the code on console. If you have a Graphviz plotting device in hand,
simply copy and paste the code to that engine to reproduce the plot. For example, 
use our [Workflow Plot Editor](https://systempipe.org/sp/spr/viz_editor).

```{r eval=FALSE}
plotWF(sal, out_format = "dot_print") 
```

### Save to a static image file

Some users may want to save the plot to a static image, like `.png` format. We will 
need do some extra work to save the file. The reason we cannot directly save it to 
a png file is the plot is generated in real-time by a browser javascript engine. It 
requires one type of javascript engine, like Chrome, MS Edge, Viewer in Rstudio,
to render the plot before we can see it.

#### Interactive 

- With the `plot_ctr = TRUE` (default) option, a plot control panel is displayed 
  on the top-left corner. One can choose from different formats like _png_, _jpg_, 
  _svg_ to download them from the webpage. To enable these buttons, internet connection is required.
  The underlying Javascript libraries are download on the fly. Please make sure 
  internet is connected. There are known conflicts of underlying web format
  creating libraries with R markdown web libraries, so some of these buttons may not work
  inside R markdown as you are seeing in this vignette right now. However, they 
  should work properly if the workflow plot is saved to a stand-alone HTML file.
- If you are working in Rstudio, you can also use the `export` button in the viewer to save 
an image file. 

  
#### Non-interactive

If you cannot have an interactive session, like submitting a job to a cluster, 
but still want the png, we recommend to use the {[webshot2](https://github.com/wch/webshot)}
package to screenshot the plot. It runs headless Chrome in the back-end (which has a javascript engine).

Install the package 

```{r eval=FALSE}
remotes::install_github("rstudio/webshot2")
```

Save to html first

```{r eval=FALSE}
plotWF(sal, out_format = "html", out_path = "example_out.html")
file.exists("example_out.html")
```

Use `webshot2` to save the image

```{r eval=FALSE}
webshot2::webshot("example_out.html", "example_out.png")
```


## plotWF in Shiny
`plotWF` is able to render plots in Shiny apps. Please refer to the help document of
`renderPlotwf`, `plotwfOutput`, and check out systemPipeShiny [demos](https://systempipe.org/sps/#demo).

# Inner Classes

`SYSargsList` steps are can be defined with two inner classes, `SYSargs2` and
`LineWise`. Next, more details on both classes.

## `SYSargs2` Class {#sysargs2}

*`SYSargs2`* workflow control class, an S4 class, is a list-like container where
each instance stores all the input/output paths and parameter components
required for a particular data analysis step. *`SYSargs2`* instances are
generated by two constructor functions, *loadWF* and *renderWF*, using as data
input *targets* or *yaml* files as well as two *cwl* parameter files (for
details see below).

In CWL, files with the extension *`.cwl`* define the parameters of a chosen
command-line step or workflow, while files with the extension *`.yml`* define
the input variables of command-line steps. Note, input variables provided by a
*targets* file can be passed on to a *`SYSargs2`* instance via the *inputvars*
argument of the *renderWF* function.

The following imports a *`.cwl`* file (here *`hisat2-mapping-se.cwl`*) for
running the short read aligner HISAT2 [@Kim2015-ve]. For more details about the
file structure and how to design or customize our own software tools, please
check `systemPipeR and CWL` pipeline.

```{r sysargs2_cwl_structure, echo = FALSE, eval=FALSE}
hisat2.cwl <- system.file("extdata", "cwl/hisat2/hisat2-mapping-se.cwl", package = "systemPipeR")
yaml::read_yaml(hisat2.cwl)
```

```{r sysargs2_yaml_structure, echo = FALSE, eval=FALSE}
hisat2.yml <- system.file("extdata", "cwl/hisat2/hisat2-mapping-se.yml", package = "systemPipeR")
yaml::read_yaml(hisat2.yml)
```

The *loadWF* and *renderWF* functions render the proper command-line strings for
each sample and software tool.

```{r SYSargs2_structure, eval=TRUE}
library(systemPipeR)
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
dir_path <- system.file("extdata/cwl", package = "systemPipeR")
WF <- loadWF(targets = targetspath, wf_file = "hisat2/hisat2-mapping-se.cwl",
                   input_file = "hisat2/hisat2-mapping-se.yml",
                   dir_path = dir_path)

WF <- renderWF(WF, inputvars = c(FileName = "_FASTQ_PATH1_", 
                                 SampleName = "_SampleName_"))
```

Several accessor methods are available that are named after the slot names of
the *`SYSargs2`* object.

```{r names_WF, eval=TRUE}
names(WF)
```

Of particular interest is the *`cmdlist()`* method. It constructs the system
commands for running command-line software as specified by a given *`.cwl`* file
combined with the paths to the input samples (*e.g.* FASTQ files) provided by a
*`targets`* file. The example below shows the *`cmdlist()`* output for running
HISAT2 on the first SE read sample. Evaluating the output of *`cmdlist()`* can
be very helpful for designing and debugging *`.cwl`* files of new command-line
software or changing the parameter settings of existing ones.

```{r cmdlist, eval=TRUE}
cmdlist(WF)[1]
```

The output components of *`SYSargs2`* define the expected output files for each
step in the workflow; some of which are the input for the next workflow step,
here next *`SYSargs2`* instance.

```{r output_WF, eval=TRUE}
output(WF)[1]
```

The targets components of `SYSargs2` object can be accessed by the targets
method. Here, for single-end (SE) samples, the structure of the targets file is
defined by:

- `FileName`: specify the FASTQ files path;
- `SampleName`: Unique IDs for each sample;
- `Factor`: ID for each treatment or condition.

```{r, targets_WF, eval=TRUE}
targets(WF)[1]
as(WF, "DataFrame")
```

Please note, to work with custom data, users need to generate a *`targets`* file
containing the paths to their own FASTQ files and then provide under
*`targetspath`* the path to the corresponding *`targets`* file.

In addition, if the [Environment Modules](http://modules.sourceforge.net/) is
available, it is possible to define which module should be loaded, as shown
here:

```{r, module_WF, eval=TRUE}
modules(WF)
```

Additional information can be accessed, as the parameters files location and the
`inputvars` provided to generate the object.

```{r, other_WF, eval=FALSE}
files(WF)
inputvars(WF)
```

## LineWise Class {#linewise}

`LineWise` was designed to store all the R code chunk when an RMarkdown file is
imported as a workflow.

```{r lw, eval=TRUE}
rmd <- system.file("extdata", "spr_simple_lw.Rmd", package = "systemPipeR")
sal_lw <- SPRproject(overwrite = TRUE)
sal_lw <- importWF(sal_lw, rmd, verbose = FALSE)
codeLine(sal_lw)
```

- Coerce methods available:

```{r, lw_coerce, eval=TRUE}
lw <- stepsWF(sal_lw)[[2]]
## Coerce
ll <- as(lw, "list")
class(ll)
lw <- as(ll, "LineWise")
lw
```

- Access details

```{r, lw_access, eval=TRUE}
length(lw)
names(lw)
codeLine(lw)
codeChunkStart(lw)
rmdPath(lw)
```

- Subsetting

```{r, lw_sub, eval=TRUE}
l <- lw[2]
codeLine(l)
l_sub <- lw[-2]
codeLine(l_sub)
```

- Replacement methods

```{r, lw_rep, eval=TRUE}
replaceCodeLine(lw, line = 2) <- "5+5"
codeLine(lw)
appendCodeLine(lw, after = 0) <- "6+7"
codeLine(lw)
```

- Replacement methods for `SYSargsList`

```{r, sal_rep_append, eval=FALSE}
replaceCodeLine(sal_lw, step = 2, line = 2) <- LineWise(code={
                                                             "5+5"
                                                                })
codeLine(sal_lw, step = 2)

appendCodeLine(sal_lw, step = 2) <- "66+55"
codeLine(sal_lw, step = 2)

appendCodeLine(sal_lw, step = 1, after = 1) <- "66+55"
codeLine(sal_lw, step = 1)
```

## Workflow design structure using *`SYSargs`*: Previous version

Instances of this S4 object class are constructed by the *`systemArgs`* function
from two simple tabular files: a *`targets`* file and a *`param`* file. The
latter is optional for workflow steps lacking command-line software. Typically,
a *`SYSargs`* instance stores all sample-level inputs as well as the paths to
the corresponding outputs generated by command-line- or R-based software
generating sample-level output files, such as read preprocessors
(trimmed/filtered FASTQ files), aligners (SAM/BAM files), variant callers
(VCF/BCF files) or peak callers (BED/WIG files). Each sample level input/output
operation uses its own *`SYSargs`* instance. The outpaths of *`SYSargs`* usually
define the sample inputs for the next *`SYSargs`* instance. This connectivity is
established by writing the outpaths with the *`writeTargetsout`* function to a
new *`targets`* file that serves as input to the next *`systemArgs`* call.
Typically, the user has to provide only the initial *`targets`* file. All
downstream *`targets`* files are generated automatically. By chaining several
*`SYSargs`* steps together one can construct complex workflows involving many
sample-level input/output file operations with any combination of command-line
or R-based software.

<center><img src="images/SystemPipeR_Workflow.png"></center>

```{r, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Workflow design structure of *`systemPipeR`* using previous version of *`SYSargs`*", warning=FALSE}
# knitr::include_graphics(system.file("extdata/images", "SystemPipeR_Workflow.png", package = "systemPipeR"))
```

# Third-party software tools {#third-party-software-tools}

Current, *systemPipeR* provides the _`param`_ file templates for third-party 
software tools. Please check the listed software tools. 

```{r table_tools, echo=FALSE, message=FALSE}
library(magrittr)
SPR_software <- system.file("extdata", "SPR_software.csv", package = "systemPipeR")
software <- read.delim(SPR_software, sep = ",", comment.char = "#")
colors <- colorRampPalette((c("darkseagreen", "indianred1")))(length(unique(software$Category)))
id <- as.numeric(c((unique(software$Category))))
software %>%
  dplyr::mutate(Step = kableExtra::cell_spec(Step, color = "white", bold = TRUE,
    background = factor(Category, id, colors)
  )) %>%
   dplyr::select(Tool, Description, Step) %>%
  dplyr::arrange(Tool) %>% 
  kableExtra::kable(escape = FALSE, align = "c", col.names = c("Tool Name",	"Description", "Step")) %>%
  kableExtra::kable_styling(c("striped", "hover", "condensed"), full_width = TRUE) %>%
  kableExtra::scroll_box(width = "80%", height = "500px")
```

Remember, if you desire to run any of these tools, make sure to have the 
respective software installed on your system and configure in the `PATH`. 
You can check as follows: 

```{r test_tool_path, eval=FALSE}
tryCMD(command="gzip") 
```

```{r cleaning3, eval=TRUE, include=FALSE}
if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE)
## NOTE: Removing previous project create in the quick starts section
```

# Workflow commom steps overview

## Project initialization

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

```{r SPRproject2, eval=FALSE}
sal <- SPRproject() 
```

### Required packages and resources

The `systemPipeR` package needs to be loaded [@H_Backman2016-bt].

```{r load_SPR, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            library(systemPipeR)
                            }, 
                            step_name = "load_SPR")
```

## Read Preprocessing

### 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 adapter 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 adapter trimming with
the `trimLRPatterns` function from the `Biostrings` package.

Here, we are appending this step to 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 = system.file("extdata/cwl", package = "systemPipeR"),
    inputvars = c(
        FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"
    ),
    dependency = c("load_SPR"))
```

After the preprocessing 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 `appendStep` function is
automatically handling this connectivity between steps. Please check the next 
step for more details.

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. 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 custom 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 custom 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
```

### Preprocessing with TrimGalore!

[TrimGalore!](http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) is 
a wrapper tool to consistently apply quality and adapter trimming to fastq files, 
with some extra functionality for removing Reduced Representation Bisulfite-Seq 
(RRBS) libraries. 

```{r trimGalore, eval=FALSE, spr=TRUE}
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "trimGalore", 
                               targets = targetspath, dir = TRUE,
                               wf_file = "trim_galore/trim_galore-se.cwl", 
                               input_file = "trim_galore/trim_galore-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"), 
                               dependency = "load_SPR", 
                               run_step = "optional")
```


### Preprocessing with Trimmomatic

[Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) software [@Bolger2014-yr] 
performs a variety of useful trimming tasks for Illumina paired-end and single
ended data. Here, an example of how to perform this task using parameters template
files for trimming FASTQ files.

```{r trimmomatic, eval=FALSE, spr=TRUE}
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "trimmomatic", 
                               targets = targetspath, dir = TRUE,
                               wf_file = "trimmomatic/trimmomatic-se.cwl", 
                               input_file = "trimmomatic/trimmomatic-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"), 
                               dependency = "load_SPR", 
                               run_step = "optional")
```

## 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 are
written to a PDF file named `fastqReport.pdf`.

```{r fastq_report, eval=FALSE, message=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(code = {
                fastq <- getColumn(sal, step = "preprocessing", "targetsWF", column = 1)
                fqlist <- seeFastq(fastq = fastq, batchsize = 10000, klength = 8)
                pdf("./results/fastqReport.pdf", height = 18, width = 4*length(fqlist))
                seeFastqPlot(fqlist)
                dev.off()
                }, step_name = "fastq_report", 
                dependency = "preprocessing")
``` 

<center><img src="results/fastqReport.png"></center>
<div align="center">**Figure 1:** FASTQ quality report </div></br>

## NGS Alignment software

After quality control, the sequence reads can be aligned to a reference genome or 
transcriptome database. The following sessions present some NGS sequence alignment 
software. Select the most accurate aligner and determining the optimal parameter 
for your custom data set project.

For all the following examples, it is necessary to install the respective software 
and export the `PATH` accordingly. 

### Alignment with `HISAT2` 

The following steps will demonstrate how to use the short read aligner `Hisat2`
[@Kim2015-ve] in both interactive job submissions and batch submissions to
queuing systems of clusters using the _`systemPipeR's`_ new CWL command-line interface.

- Build `Hisat2` index.

```{r hisat_index, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "hisat_index", 
                               targets = NULL, dir = FALSE,
                               wf_file = "hisat2/hisat2-index.cwl", 
                               input_file = "hisat2/hisat2-index.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = NULL, 
                               dependency = "preprocessing")
```

The parameter settings of the aligner are defined in the `workflow_hisat2-se.cwl` 
and `workflow_hisat2-se.yml` files. The following shows how to construct the 
corresponding *SYSargsList* object, and append to *sal* workflow.

- Alignment with `HISAT2` and `SAMtools`

It possible to build an workflow with `HISAT2` and `SAMtools`.

```{r hisat_mapping_samtools, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "hisat_mapping", 
                               targets = "preprocessing", dir = TRUE,
                               wf_file = "workflow-hisat2/workflow_hisat2-se.cwl", 
                               input_file = "workflow-hisat2/workflow_hisat2-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), 
                               dependency = c("hisat_index"), 
                               run_session = "remote")
```

### Alignment with `Tophat2`

The NGS reads of this project can also be aligned against the reference genome 
sequence using `Bowtie2/TopHat2` [@Kim2013-vg; @Langmead2012-bs]. 

- Build _`Bowtie2`_ index.

```{r bowtie_index, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "bowtie_index", 
                               targets = NULL, dir = FALSE,
                               wf_file = "bowtie2/bowtie2-index.cwl", 
                               input_file = "bowtie2/bowtie2-index.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = NULL, 
                               dependency = "preprocessing", 
                               run_step = "optional")
```


The parameter settings of the aligner are defined in the `workflow_tophat2-mapping.cwl` 
and `tophat2-mapping-pe.yml` files. The following shows how to construct the 
corresponding *SYSargsList* object, using the outfiles from the `preprocessing` step.

```{r tophat2_mapping, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "tophat2_mapping", 
                               targets = "preprocessing", dir = TRUE,
                               wf_file = "tophat2/workflow_tophat2-mapping-se.cwl", 
                               input_file = "tophat2/tophat2-mapping-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), 
                               dependency = c("bowtie_index"), 
                               run_session = "remote", 
                               run_step = "optional")
```


### Alignment with _`Bowtie2`_ (_e.g._ for miRNA profiling)

The following example runs _`Bowtie2`_ as a single process without submitting it to a cluster.

```{r bowtie2_mapping, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "bowtie2_mapping", 
                               targets = "preprocessing", dir = TRUE,
                               wf_file = "bowtie2/workflow_bowtie2-mapping-se.cwl", 
                               input_file = "bowtie2/bowtie2-mapping-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), 
                               dependency = c("bowtie_index"), 
                               run_session = "remote", 
                               run_step = "optional")
```

### Alignment with _`BWA-MEM`_ (_e.g._ for VAR-Seq)

The following example runs BWA-MEM as a single process without submitting it to a cluster.

- Build the index:

```{r bwa_index, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "bwa_index", 
                               targets = NULL, dir = FALSE,
                               wf_file = "bwa/bwa-index.cwl", 
                               input_file = "bwa/bwa-index.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = NULL, 
                               dependency = "preprocessing", 
                               run_step = "optional")
```

- Prepare the alignment step:

```{r bwa_mapping, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "bwa_mapping", 
                               targets = "preprocessing", dir = TRUE,
                               wf_file = "bwa/bwa-se.cwl", 
                               input_file = "bwa/bwa-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), 
                               dependency = c("bwa_index"), 
                               run_session = "remote", 
                               run_step = "optional")
```

### Alignment with _`Rsubread`_ (_e.g._ for RNA-Seq)

The following example shows how one can use within the \Rpackage{systemPipeR} environment the R-based aligner \Rpackage{Rsubread}, allowing running from R or command-line. 

- Build the index:

```{r rsubread_index, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "rsubread_index", 
                               targets = NULL, dir = FALSE,
                               wf_file = "rsubread/rsubread-index.cwl", 
                               input_file = "rsubread/rsubread-index.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = NULL, 
                               dependency = "preprocessing", 
                               run_step = "optional")
```

- Prepare the alignment step:

```{r rsubread_mapping, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "rsubread", 
                               targets = "preprocessing", dir = TRUE,
                               wf_file = "rsubread/rsubread-mapping-se.cwl", 
                               input_file = "rsubread/rsubread-mapping-se.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), 
                               dependency = c("rsubread_index"), 
                               run_session = "remote", 
                               run_step = "optional")
```


### Alignment with _`gsnap`_ (_e.g._ for VAR-Seq and RNA-Seq)

Another R-based short read aligner is _`gsnap`_ from the _`gmapR`_ package [@Wu2010-iq]. 
The code sample below introduces how to run this aligner on multiple nodes of a compute cluster.

- Build the index:

```{r gsnap_index, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "gsnap_index", 
                               targets = NULL, dir = FALSE,
                               wf_file = "gsnap/gsnap-index.cwl", 
                               input_file = "gsnap/gsnap-index.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                               inputvars = NULL, 
                               dependency = "preprocessing", 
                               run_step = "optional")
```

- Prepare the alignment step:

```{r gsnap_mapping, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(step_name = "gsnap", 
                               targets = "targetsPE.txt", dir = TRUE,
                               wf_file = "gsnap/gsnap-mapping-pe.cwl", 
                               input_file = "gsnap/gsnap-mapping-pe.yml", 
                               dir_path = system.file("extdata/cwl", package = "systemPipeR"),
                                inputvars = c(FileName1 = "_FASTQ_PATH1_", 
                                              FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"), 
                               dependency = c("gsnap_index"), 
                               run_session = "remote", 
                               run_step = "optional")
```


## Create symbolic links for viewing BAM files in IGV

The genome browser IGV supports reading of indexed/sorted BAM files via web URLs. This way it can be avoided to create unnecessary copies of these large files. To enable this approach, an HTML directory with Http access needs to be available in the user account (_e.g._ _`home/publichtml`_) of a system. If this is not the case then the BAM files need to be moved or copied to the system where IGV runs. In the following, _`htmldir`_ defines the path to the HTML directory with http access where the symbolic links to the BAM files will be stored. The corresponding URLs will be written to a text file specified under the `_urlfile`_ argument. 

```{r bam_IGV, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
    code = {
        bampaths <- getColumn(sal, step = "hisat2_mapping", "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_IGV",
    dependency = "hisat2_mapping",
    run_step = "optional"
)
```

## Read counting for mRNA profiling experiments

Create _`txdb`_ (needs to be done only once).

```{r create_txdb, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(code = {
                            library(GenomicFeatures)
                            txdb <- makeTxDbFromGFF(file="data/tair10.gff", format="gff", 
                                                    dataSource="TAIR", organism="Arabidopsis thaliana")
                            saveDb(txdb, file="./data/tair10.sqlite")
                            }, 
                            step_name = "create_txdb", 
                            dependency = "hisat_mapping")
```


The following performs read counting with _`summarizeOverlaps`_ in parallel mode with multiple cores. 

```{r read_counting, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            library(BiocParallel)
                            txdb <- loadDb("./data/tair10.sqlite")
                            eByg <- exonsBy(txdb, by="gene")
                            outpaths <- getColumn(sal, step = "hisat_mapping", 'outfiles', column = 2)
                            bfl <- BamFileList(outpaths, yieldSize=50000, index=character())
                            multicoreParam <- MulticoreParam(workers=4); register(multicoreParam); registered()
                            counteByg <- bplapply(bfl, function(x) summarizeOverlaps(eByg, x, mode="Union",
                                                                                     ignore.strand=TRUE,
                                                                                     inter.feature=TRUE,
                                                                                     singleEnd=TRUE))
                            # Note: for strand-specific RNA-Seq set 'ignore.strand=FALSE' and for PE data set 'singleEnd=FALSE'
                            countDFeByg <- sapply(seq(along=counteByg), 
                                                  function(x) assays(counteByg[[x]])$counts)
                            rownames(countDFeByg) <- names(rowRanges(counteByg[[1]]))
                            colnames(countDFeByg) <- names(bfl)
                            rpkmDFeByg <- apply(countDFeByg, 2, function(x) returnRPKM(counts=x, ranges=eByg))
                            write.table(countDFeByg, "results/countDFeByg.xls", 
                                        col.names=NA, quote=FALSE, sep="\t")
                            write.table(rpkmDFeByg, "results/rpkmDFeByg.xls", 
                                        col.names=NA, quote=FALSE, sep="\t")
                            }, 
                            step_name = "read_counting", 
                            dependency = "create_txdb")
```

Please note, in addition to read counts this step generates RPKM normalized expression values. For most statistical differential expression or abundance analysis methods, such as _`edgeR`_ or _`DESeq2`_, the raw count values should be used as input. The usage of RPKM values should be restricted to specialty applications required by some users, _e.g._ manually comparing the expression levels of different genes or features. 

#### Read and alignment count stats

Generate a table of read and alignment counts for all samples. 

```{r align_stats, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            read_statsDF <- alignStats(args)
                            write.table(read_statsDF, "results/alignStats.xls", 
                                        row.names = FALSE, quote = FALSE, sep = "\t")
                            }, 
                            step_name = "align_stats", 
                            dependency = "hisat_mapping")
```


The following shows the first four lines of the sample alignment stats file 
provided by the _`systemPipeR`_ package. For simplicity the number of PE reads 
is multiplied here by 2 to approximate proper alignment frequencies where each 
read in a pair is counted. 

```{r align_stats2, eval=TRUE}
read.table(system.file("extdata", "alignStats.xls", package = "systemPipeR"), header = TRUE)[1:4,]
```

## Read counting for miRNA profiling experiments

Download `miRNA` genes from `miRBase`.

```{r read_counting_mirna, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            system("wget https://www.mirbase.org/ftp/CURRENT/genomes/ath.gff3 -P ./data/")
                            gff <- rtracklayer::import.gff("./data/ath.gff3")
                            gff <- split(gff, elementMetadata(gff)$ID)
                            bams <- getColumn(sal, step = "bowtie2_mapping", 'outfiles', column = 2)
                            bfl <- BamFileList(bams, yieldSize=50000, index=character())
                            countDFmiR <- summarizeOverlaps(gff, bfl, mode="Union",
                                                            ignore.strand = FALSE, inter.feature = FALSE) 
                            countDFmiR <- assays(countDFmiR)$counts
                            # Note: inter.feature=FALSE important since pre and mature miRNA ranges overlap
                            rpkmDFmiR <- apply(countDFmiR, 2, function(x) returnRPKM(counts = x, ranges = gff))
                            write.table(assays(countDFmiR)$counts, "results/countDFmiR.xls", 
                                        col.names=NA, quote=FALSE, sep="\t")
                            write.table(rpkmDFmiR, "results/rpkmDFmiR.xls", col.names=NA, quote=FALSE, sep="\t")
                            }, 
                            step_name = "read_counting_mirna", 
                            dependency = "bowtie2_mapping")
```


## Correlation analysis of samples

The following computes the sample-wise Spearman correlation coefficients from the _`rlog`_ (regularized-logarithm) transformed expression values generated with the _`DESeq2`_ package. After transformation to a distance matrix, hierarchical clustering is performed with the _`hclust`_ function and the result is plotted as a dendrogram ([sample\_tree.pdf](sample_tree.png)). 

```{r sample_tree_rlog, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            library(DESeq2, warn.conflicts=FALSE, quietly=TRUE)
                            library(ape, warn.conflicts=FALSE)
                            countDFpath <- system.file("extdata", "countDFeByg.xls", package="systemPipeR")
                            countDF <- as.matrix(read.table(countDFpath))
                            colData <- data.frame(row.names = targetsWF(sal)[[2]]$SampleName,  
                                                  condition=targetsWF(sal)[[2]]$Factor)
                            dds <- DESeqDataSetFromMatrix(countData = countDF, colData = colData, 
                                                          design = ~ condition)
                            d <- cor(assay(rlog(dds)), method = "spearman")
                            hc <- hclust(dist(1-d))
                            plot.phylo(as.phylo(hc), type = "p", edge.col = 4, edge.width = 3,
                                       show.node.label = TRUE, no.margin = TRUE)
                            }, 
                            step_name = "sample_tree_rlog", 
                            dependency = "read_counting")
```

<center><img src="results/sample_tree.png"></center>

<div align="center">**Figure 2:** Correlation dendrogram of samples for _`rlog`_ values. </div></br>

## DEG analysis with _`edgeR`_

The following _`run_edgeR`_ function is a convenience wrapper for
identifying differentially expressed genes (DEGs) in batch mode with
_`edgeR`_'s GML method [@Robinson2010-uk] for any number of
pairwise sample comparisons specified under the _`cmp`_ argument. Users
are strongly encouraged to consult the 
[_`edgeR`_](\href{http://www.bioconductor.org/packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf) vignette 
for more detailed information on this topic and how to properly run _`edgeR`_ 
on data sets with more complex experimental designs. 

```{r edger, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
                            targets <- read.delim(targetspath, comment = "#")
                            cmp <- readComp(file = targetspath, format = "matrix", delim = "-")
                            countDFeBygpath <- system.file("extdata", "countDFeByg.xls", package = "systemPipeR")
                            countDFeByg <- read.delim(countDFeBygpath, row.names = 1)
                            edgeDF <- run_edgeR(countDF = countDFeByg, targets = targets, cmp = cmp[[1]],
                                                independent = FALSE, mdsplot = "")
                            DEG_list <- filterDEGs(degDF = edgeDF, filter = c(Fold = 2, FDR = 10))
                            }, 
                            step_name = "edger", 
                            dependency = "read_counting")
```

Filter and plot DEG results for up and down-regulated genes. Because of the small size of the toy data set used by this vignette, the _FDR_ value has been set to a relatively high threshold (here 10%). More commonly used _FDR_ cutoffs are 1% or 5%. The definition of '_up_' and '_down_' is given in the corresponding help file. To open it, type _`?filterDEGs`_ in the R console. 

<center><img src="results/DEGcounts.png"></center>

<div align="center">**Figure 3:** Up and down regulated DEGs identified by _`edgeR`_. </div></br>

## DEG analysis with _`DESeq2`_ 

The following _`run_DESeq2`_ function is a convenience wrapper for
identifying DEGs in batch mode with _`DESeq2`_ [@Love2014-sh] for any number of
pairwise sample comparisons specified under the _`cmp`_ argument. Users
are strongly encouraged to consult the 
[_`DESeq2`_](http://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.pdf) vignette
for more detailed information on this topic and how to properly run _`DESeq2`_ 
on data sets with more complex experimental designs. 

```{r deseq2, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            degseqDF <- run_DESeq2(countDF=countDFeByg, targets=targets, cmp=cmp[[1]],
                                                   independent=FALSE)
                            DEG_list2 <- filterDEGs(degDF=degseqDF, filter=c(Fold=2, FDR=10))
                            }, 
                            step_name = "deseq2", 
                            dependency = "read_counting")
```

## Venn Diagrams

The function _`overLapper`_ can compute Venn intersects for large numbers of sample sets (up to 20 or more) and _`vennPlot`_ can plot 2-5 way Venn diagrams. A useful feature is the possibility to combine the counts from several Venn comparisons with the same number of sample sets in a single Venn diagram (here for 4 up and down DEG sets).

```{r vennplot, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            vennsetup <- overLapper(DEG_list$Up[6:9], type="vennsets")
                            vennsetdown <- overLapper(DEG_list$Down[6:9], type="vennsets")
                            vennPlot(list(vennsetup, vennsetdown), mymain="", mysub="", 
                                     colmode=2, ccol=c("blue", "red"))
                            }, 
                            step_name = "vennplot", 
                            dependency = "edger")
```

<center><img src="results/vennplot.png"></center>
<div align="center">**Figure 4:** Venn Diagram for 4 Up and Down DEG Sets. </div></br>

## GO term enrichment analysis of DEGs

### Obtain gene-to-GO mappings

The following shows how to obtain gene-to-GO mappings from _`biomaRt`_ (here for _A. thaliana_) and how to organize them for the downstream GO term enrichment analysis. Alternatively, the gene-to-GO mappings can be obtained for many organisms from Bioconductor's  _`*.db`_ genome annotation packages or GO annotation files provided by various genome databases. For each annotation, this relatively slow preprocessing step needs to be performed only once. Subsequently, the preprocessed data can be loaded with the _`load`_ function as shown in the next subsection. 

```{r get_go_biomart, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            library("biomaRt")
                            listMarts() # To choose BioMart database
                            listMarts(host="plants.ensembl.org")
                            m <- useMart("plants_mart", host="https://plants.ensembl.org")
                            listDatasets(m)
                            m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org")
                            listAttributes(m) # Choose data types you want to download
                            go <- getBM(attributes=c("go_id", "tair_locus", "namespace_1003"), mart=m)
                            go <- go[go[,3]!="",]; go[,3] <- as.character(go[,3])
                            go[go[,3]=="molecular_function", 3] <- "F"
                            go[go[,3]=="biological_process", 3] <- "P"
                            go[go[,3]=="cellular_component", 3] <- "C"
                            go[1:4,]
                            dir.create("./data/GO")
                            write.table(go, "data/GO/GOannotationsBiomart_mod.txt", 
                                        quote=FALSE, row.names=FALSE, col.names=FALSE, sep="\t")
                            catdb <- makeCATdb(myfile="data/GO/GOannotationsBiomart_mod.txt",
                                               lib=NULL, org="", colno=c(1,2,3), idconv=NULL)
                            save(catdb, file="data/GO/catdb.RData")
                            }, 
                            step_name = "get_go_biomart", 
                            dependency = "edger")
```

### Batch GO term enrichment analysis

Apply the enrichment analysis to the DEG sets obtained in the above differential expression analysis. Note, in the following example the _FDR_ filter is set here to an unreasonably high value, simply because of the small size of the toy data set used in this vignette. Batch enrichment analysis of many gene sets is performed with the _`GOCluster_Report`_ function. When _`method="all"`_, it returns all GO terms passing the p-value cutoff specified under the _`cutoff`_ arguments. When _`method="slim"`_, it returns only the GO terms specified under the _`myslimv`_ argument. The given example shows how one can obtain such a GO slim vector from BioMart for a specific organism. 

```{r go_enrichment, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            load("data/GO/catdb.RData")
                            DEG_list <- filterDEGs(degDF=edgeDF, filter=c(Fold=2, FDR=50), plot=FALSE)
                            up_down <- DEG_list$UporDown; names(up_down) <- paste(names(up_down), "_up_down", sep="")
                            up <- DEG_list$Up; names(up) <- paste(names(up), "_up", sep="")
                            down <- DEG_list$Down; names(down) <- paste(names(down), "_down", sep="")
                            DEGlist <- c(up_down, up, down)
                            DEGlist <- DEGlist[sapply(DEGlist, length) > 0]
                            BatchResult <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="all",
                                                            id_type="gene", CLSZ=2, cutoff=0.9,
                                                            gocats=c("MF", "BP", "CC"), recordSpecGO=NULL)
                            library("biomaRt")
                            m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org")
                            goslimvec <- as.character(getBM(attributes=c("goslim_goa_accession"), mart=m)[,1])
                            BatchResultslim <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="slim",
                                                                id_type="gene", myslimv=goslimvec, CLSZ=10,
                                                                cutoff=0.01, gocats=c("MF", "BP", "CC"),
                                                                recordSpecGO=NULL)
                            gos <- BatchResultslim[grep("M6-V6_up_down", BatchResultslim$CLID), ]
                            gos <- BatchResultslim
                            pdf("GOslimbarplotMF.pdf", height=8, width=10); goBarplot(gos, gocat="MF"); dev.off()
                            goBarplot(gos, gocat="BP")
                            goBarplot(gos, gocat="CC")
                            }, 
                            step_name = "go_enrichment", 
                            dependency = "get_go_biomart")
```

### Plot batch GO term results

The _`data.frame`_ generated by _`GOCluster_Report`_ can be plotted with the _`goBarplot`_ function. Because of the variable size of the sample sets, it may not always be desirable to show the results from different DEG sets in the same bar plot. Plotting single sample sets is achieved by subsetting the input data frame as shown in the first line of the following example. 

<center><img src="results/GOslimbarplotMF.png"></center>
<div align="center">**Figure 5:** GO Slim Barplot for MF Ontology.</div></br>


## Clustering and heat maps

The following example performs hierarchical clustering on the _`rlog`_ transformed expression matrix subsetted by the DEGs identified in the 
above differential expression analysis. It uses a Pearson correlation-based distance measure and complete linkage for cluster join.

```{r hierarchical_clustering, message=FALSE, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise({
                            library(pheatmap)
                            geneids <- unique(as.character(unlist(DEG_list[[1]])))
                            y <- assay(rlog(dds))[geneids, ]
                            pdf("heatmap1.pdf")
                            pheatmap(y, scale="row", clustering_distance_rows="correlation",
                                     clustering_distance_cols="correlation")
                            dev.off()
                            }, 
                            step_name = "hierarchical_clustering", 
                            dependency = c("sample_tree_rlog", "edgeR"))
```


<center><img src="results/heatmap1.png"></center>
<div align="center">**Figure 7:** Heat map with hierarchical clustering dendrograms of DEGs.</div></br>

## Visualize workflow

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

This function will make a plot of selected workflow instance and the following information is displayed on the plot:

- Workflow structure (dependency graphs between different steps); 
- Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`; 
- Sample status and statistics; 
- Workflow timing: running duration time. 

If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, etc.

```{r}
plotWF(sal, show_legend = TRUE, width = "80%")
```

## 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_steps, 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_steps, eval=FALSE}
resources <- list(conffile=".batchtools.conf.R",
                  template="batchtools.slurm.tmpl", 
                  Njobs=18, 
                  walltime=120, ## minutes
                  ntasks=1,
                  ncpus=4, 
                  memory=1024, ## Mb
                  partition = "short"
                  )
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)
```

### Checking workflow status

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

```{r statusWF_steps, 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_steps, eval=FALSE}
sal <- renderLogs(sal)
```

```{r quiting, eval=TRUE, echo=FALSE}
knitr::opts_knit$set(root.dir = '../')
unlink("rnaseq", recursive = TRUE)
```

# Version information

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

# Funding

This project is funded by NSF award [ABI-1661152](https://www.nsf.gov/awardsearch/showAward?AWD_ID=1661152). 

# References
