1 Introduction

systemPipeR is a multipurpose data analysis workflow environment that unifies R with command-line tools (H Backman and Girke 2016). It enables scientists to analyze many types of large- or small-scale data on local or distributed computer systems with a high level of reproducibility, scalability and portability (Figure 1). At its core is a command-line interface (CLI) that adopts the Common Workflow Language (CWL, Crusoe et al. (2021)). This design allows users to choose for each analysis step the optimal R or command-line software. It supports both end-to-end and partial execution of workflows with built-in restart functionalities. Efficient management of complex analysis tasks is accomplished by a flexible workflow control container class (SYSargsList). Handling of large numbers of input samples and experimental designs is facilitated by consistent sample annotation mechanisms. As a multi-purpose workflow toolkit, systemPipeR enables users to run existing workflows, customize them or design entirely new ones while taking advantage of widely adopted data structures within the Bioconductor ecosystem. Another important core functionality is the generation of reproducible scientific analysis and technical reports. For result interpretation, systemPipeR offers a wide range of plotting functionality, while an associated Shiny App offers many useful functionalities for interactive result exploration.

Relevant features in `systemPipeR`. Workflow design concepts are illustrated under (A). Examples of `systemPipeR's` visualization functionalities are given under (B).

Figure 1: Relevant features in systemPipeR
Workflow design concepts are illustrated under (A). Examples of systemPipeR's visualization functionalities are given under (B).

systemPipeR's CWL interface 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 step prior to running a command-line. Auto-generation of CWL parameter files is also supported, where users can simply provide the command-line strings new software to a function and the corresponding *.cwl and *.yml are generated for them.

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.

Figure 2: 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.

1.1 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 3). 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.

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.

Figure 3: 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.

1.2 Command-line software support

systemPipeR adopted the widely used community standard Common Workflow Language (CWL) (Amstutz et al. 2016) 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 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 4A-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 4A). 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.

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).

Figure 4: 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).

2 Getting Started

2.1 Installation

systemPipeR environment can be installed from the R console using the BiocManager::install command. The associated data package systemPipeRdata 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.

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.

2.2 Loading package and documentation

library("systemPipeR") # Loads the package
library(help="systemPipeR") # Lists package info
vignette("systemPipeR") # Opens vignette

2.3 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.

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.

3 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.

3.0.1 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 obtains 18 paired-end (PE) read sets from Arabidposis thaliana (Howard et al. 2013). 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.

3.1 Directory Structure

systemPipeRdata, 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 5).

genWorkenvir(workflow="rnaseq")
setwd("rnaseq")

Directory names are indicated in green. Users can change this structure as needed, but need to adjust the code in their workflows accordingly.

  • workflow/ (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. myproject). Changing its name does not require any modifications in the run script(s).
    • Important subdirectories:
      • param/
        • param/cwl/: 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.
      • data/
        • Raw data (e.g. FASTQ files)
        • FASTA file of reference (e.g. reference genome)
        • Annotation files
        • Metadata
        • etc.
      • results/
        • 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.
*systemPipeR's* preconfigured directory structure.

Figure 5: systemPipeR’s preconfigured directory structure

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.

3.2 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. 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.

3.2.1 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.

targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") 
showDF(read.delim(targetspath, comment.char = "#"))
## Loading required namespace: DT

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.

3.2.2 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.

targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))

3.2.3 Sample comparisons

Sample comparisons are defined in the header lines of the targets file starting with ‘# <CMP>’.

readLines(targetspath)[1:4]
## [1] "# Project ID: Arabidopsis - Pseudomonas alternative splicing study (SRA: SRP010938; PMID: 24098335)"                                                                              
## [2] "# The following line(s) allow to specify the contrasts needed for comparative analyses, such as DEG identification. All possible comparisons can be specified with 'CMPset: ALL'."
## [3] "# <CMP> CMPset1: M1-A1, M1-V1, A1-V1, M6-A6, M6-V6, A6-V6, M12-A12, M12-V12, A12-V12"                                                                                             
## [4] "# <CMP> CMPset2: ALL"

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.

readComp(file = targetspath, format = "vector", delim = "-")
## $CMPset1
## [1] "M1-A1"   "M1-V1"   "A1-V1"   "M6-A6"   "M6-V6"   "A6-V6"   "M12-A12" "M12-V12" "A12-V12"
## 
## $CMPset2
##  [1] "M1-A1"   "M1-V1"   "M1-M6"   "M1-A6"   "M1-V6"   "M1-M12"  "M1-A12"  "M1-V12"  "A1-V1"  
## [10] "A1-M6"   "A1-A6"   "A1-V6"   "A1-M12"  "A1-A12"  "A1-V12"  "V1-M6"   "V1-A6"   "V1-V6"  
## [19] "V1-M12"  "V1-A12"  "V1-V12"  "M6-A6"   "M6-V6"   "M6-M12"  "M6-A12"  "M6-V12"  "A6-V6"  
## [28] "A6-M12"  "A6-A12"  "A6-V12"  "V6-M12"  "V6-A12"  "V6-V12"  "M12-A12" "M12-V12" "A12-V12"

3.3 Downstream targets files description

After the step which required the initial targets file information, the downstream targets files are created automatically (see Figure 6). 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.

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

Figure 6: 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.

3.4 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 (CWL) (Amstutz et al. 2016), 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 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.

4 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 package is a helper package to generate a fully populated systemPipeR 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.

systemPipeRdata::genWorkenvir(workflow = "new", mydirname = "spr_project")
## [1] "Generated spr_project directory. Next run in new directory, the R code from *.Rmd template interactively. Alternatively, workflows can be exectued with a single command as instructed in the vignette."
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.

  • 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:

sal <- SPRproject() 
## Creating directory '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject'
## Creating file '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject/SYSargsList.yml'
  • 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.

sal <- importWF(sal, 
                file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"),
                verbose = FALSE)
## Now check if required tools are installed 
## Check if they are in path:
## Checking path for gzip
## PASS
## Checking path for gunzip
## PASS
##   step_name   tool in_path
## 1      gzip   gzip    TRUE
## 2    gunzip gunzip    TRUE
## All required tools in PATH, skip module check. If you want to check modules use `listCmdModules`

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.

sal <- runWF(sal)
## Running Step:  load_library 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
## Step Status:  Success 
## Running Step:  export_iris 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
## Step Status:  Success 
## Running Step:  gzip 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==============================                                                            |  33%
  |                                                                                                
  |============================================================                              |  67%
  |                                                                                                
  |==========================================================================================| 100%
## ---- Summary ---- 
##    Targets Total_Files Existing_Files Missing_Files    gzip
## SE      SE           1              1             0 Success
## VE      VE           1              1             0 Success
## VI      VI           1              1             0 Success
## 
## Step Status:  Success 
## Running Step:  gunzip 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==============================                                                            |  33%
  |                                                                                                
  |============================================================                              |  67%
  |                                                                                                
  |==========================================================================================| 100%
## ---- Summary ---- 
##    Targets Total_Files Existing_Files Missing_Files  gunzip
## SE      SE           1              1             0 Success
## VE      VE           1              1             0 Success
## VI      VI           1              1             0 Success
## 
## Step Status:  Success 
## Running Step:  stats 
## Running Session: Management

## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
## Step Status:  Success 
## Done with workflow running, now consider rendering logs & reports
## To render logs, run:    sal <- renderLogs(sal)
## From command-line:      Rscript -e "sal = systemPipeR::SPRproject(resume = TRUE); sal = systemPipeR::renderLogs(sal)"
## To render reports, run: sal <- renderReport(sal)
## From command-line:      Rscript -e "sal= s ystemPipeR::SPRproject(resume = TRUE); sal = systemPipeR::renderReport(sal)"
## This message is displayed once per R session
  • Visualize workflow

systemPipeR workflows instances can be visualized with the plotWF function.

plotWF(sal, width = "80%", rstudio = TRUE)
  • Checking workflow status

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

statusWF(sal)
## $load_library
## DataFrame with 1 row and 2 columns
##           Step      Status
##    <character> <character>
## 1 load_library     Success
## 
## $export_iris
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1 export_iris     Success
## 
## $gzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files     gzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $gunzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files   gunzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $stats
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1       stats     Success
  • 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.

sal <- renderLogs(sal)

5 Project initialization

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

sal <- SPRproject(projPath = getwd(), overwrite = TRUE) 
## Creating directory '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject'
## Creating file '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject/SYSargsList.yml'

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:

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. If the user wants to use a different names for these directories, can be specified as follows:

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.

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:

sal
## Instance of 'SYSargsList': 
##  No workflow steps added
projectInfo(sal)
## $project
## [1] "/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project"
## 
## $data
## [1] "data"
## 
## $param
## [1] "param"
## 
## $results
## [1] "results"
## 
## $logsDir
## [1] ".SPRproject"
## 
## $sysargslist
## [1] ".SPRproject/SYSargsList.yml"

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

length(sal)
## [1] 0

6 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, we will demonstrate how to build the workflow in an interactive mode, and in the following section, 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.

6.1 Build workflow interactive

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:

sal
## Instance of 'SYSargsList': 
##  No workflow steps added

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

6.1.1 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.

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:

sal
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Pending
## 

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

codeLine(sal)
## export_iris
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))

6.1.2 Adding more steps

Next, an example of how to compress the exported files using 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 4).

A detailed description of the dynamic between input variables and targets files can be found here. In addition, the CWL syntax overview can be found here.

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.

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:

sal
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Pending
##        2. gzip --> Status: Pending 
##            Total Files: 3 | Existing: 0 | Missing: 3 
##          2.1. gzip
##              cmdlist: 3 | Pending: 3
## 

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:

cmdlist(sal, step = "gzip")
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c  results/setosa.csv > results/SE.csv.gz"
## 
## 
## $gzip$VE
## $gzip$VE$gzip
## [1] "gzip -c  results/versicolor.csv > results/VE.csv.gz"
## 
## 
## $gzip$VI
## $gzip$VI$gzip
## [1] "gzip -c  results/virginica.csv > results/VI.csv.gz"

6.1.2.1 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:

outfiles(sal)
## $export_iris
## DataFrame with 0 rows and 0 columns
## 
## $gzip
## DataFrame with 3 rows and 1 column
##            gzip_file
##          <character>
## SE results/SE.csv.gz
## VE results/VE.csv.gz
## VI results/VI.csv.gz

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.

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:

targetsWF(sal[3])
## $gunzip
## DataFrame with 3 rows and 2 columns
##            gzip_file  SampleName
##          <character> <character>
## SE results/SE.csv.gz          SE
## VE results/VE.csv.gz          VE
## VI results/VI.csv.gz          VI

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

outfiles(sal[3])
## $gunzip
## DataFrame with 3 rows and 1 column
##       gunzip_file
##       <character>
## SE results/SE.csv
## VE results/VE.csv
## VI results/VI.csv

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

sal
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Pending
##        2. gzip --> Status: Pending 
##            Total Files: 3 | Existing: 0 | Missing: 3 
##          2.1. gzip
##              cmdlist: 3 | Pending: 3
##        3. gunzip --> Status: Pending 
##            Total Files: 3 | Existing: 0 | Missing: 3 
##          3.1. gunzip
##              cmdlist: 3 | Pending: 3
## 

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

cmdlist(sal["gzip"], targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c  results/setosa.csv > results/SE.csv.gz"

6.1.2.2 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.

getColumn(sal, step = "gunzip", 'outfiles')
##               SE               VE               VI 
## "results/SE.csv" "results/VE.csv" "results/VI.csv"
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")

6.2 Build workflow from a {R Markdown}

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.

sal_rmd <- SPRproject(logs.dir = ".SPRproject_rmd") 
## Creating directory '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject_rmd'
## Creating file '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject_rmd/SYSargsList.yml'
sal_rmd <- importWF(sal_rmd, 
                file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"))
## Reading Rmd file
## 
##  ---- Actions ----
## Checking chunk eval values
## Checking chunk SPR option
## Ignore non-SPR chunks: 17
## Parse chunk code
## Checking preprocess code for each step
## No preprocessing code for SPR steps found
## Now importing step 'load_library' 
## Now importing step 'export_iris' 
## Now importing step 'gzip' 
## Now importing step 'gunzip' 
## Now importing step 'stats' 
## Now back up current Rmd file as template for `renderReport`
## Template for renderReport is stored at 
##  /tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject_rmd/workflow_template.Rmd 
##  Edit this file manually is not recommended 
## Now check if required tools are installed 
## Check if they are in path:
## Checking path for gzip
## PASS
## Checking path for gunzip
## PASS
##   step_name   tool in_path
## 1      gzip   gzip    TRUE
## 2    gunzip gunzip    TRUE
## All required tools in PATH, skip module check. If you want to check modules use `listCmdModules`Import  done

Let’s explore the workflow to check the steps:

stepsWF(sal_rmd)
## $load_library
## Instance of 'LineWise'
##     Code Chunk length: 1
## 
## $export_iris
## Instance of 'LineWise'
##     Code Chunk length: 1
## 
## $gzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 3 (SE...VI), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 3, output: 3
##       cmdlist: 3
##    Sub Steps:
##       1. gzip (rendered: TRUE)
## 
## 
## 
## $gunzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 3 (SE...VI), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 3, output: 3
##       cmdlist: 3
##    Sub Steps:
##       1. gunzip (rendered: TRUE)
## 
## 
## 
## $stats
## Instance of 'LineWise'
##     Code Chunk length: 5
dependency(sal_rmd)
## $load_library
## [1] NA
## 
## $export_iris
## [1] "load_library"
## 
## $gzip
## [1] "export_iris"
## 
## $gunzip
## [1] "gzip"
## 
## $stats
## [1] "gunzip"
codeLine(sal_rmd)
## gzip AND gunzip step have been dropped because it is not a LineWise object.
## load_library
##     library(systemPipeR)
## export_iris
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
## stats
##     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 = 0.2, position = ggplot2::position_dodge(0.9))
targetsWF(sal_rmd)
## $load_library
## DataFrame with 0 rows and 0 columns
## 
## $export_iris
## DataFrame with 0 rows and 0 columns
## 
## $gzip
## DataFrame with 3 rows and 2 columns
##                  FileName  SampleName
##               <character> <character>
## SE     results/setosa.csv          SE
## VE results/versicolor.csv          VE
## VI  results/virginica.csv          VI
## 
## $gunzip
## DataFrame with 3 rows and 2 columns
##            gzip_file  SampleName
##          <character> <character>
## SE results/SE.csv.gz          SE
## VE results/VE.csv.gz          VE
## VI results/VI.csv.gz          VI
## 
## $stats
## DataFrame with 0 rows and 0 columns

6.2.1 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:

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

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: object not found..

7 Running the workflow

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

sal <- runWF(sal)
## Running Step:  export_iris 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
## Step Status:  Success 
## Running Step:  gzip 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==============================                                                            |  33%
  |                                                                                                
  |============================================================                              |  67%
  |                                                                                                
  |==========================================================================================| 100%
## ---- Summary ---- 
##    Targets Total_Files Existing_Files Missing_Files    gzip
## SE      SE           1              1             0 Success
## VE      VE           1              1             0 Success
## VI      VI           1              1             0 Success
## 
## Step Status:  Success 
## Running Step:  gunzip 
## Running Session: Management 
## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==============================                                                            |  33%
  |                                                                                                
  |============================================================                              |  67%
  |                                                                                                
  |==========================================================================================| 100%
## ---- Summary ---- 
##    Targets Total_Files Existing_Files Missing_Files  gunzip
## SE      SE           1              1             0 Success
## VE      VE           1              1             0 Success
## VI      VI           1              1             0 Success
## 
## Step Status:  Success 
## Running Step:  iris_stats 
## Running Session: Management

## 
  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
## Step Status:  Success

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.

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.

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:

viewEnvir(sal)

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

sal <- runWF(sal, saveEnv = TRUE)

7.1 Workflow status

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

sal
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Success
##        2. gzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          2.1. gzip
##              cmdlist: 3 | Success: 3
##        3. gunzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          3.1. gunzip
##              cmdlist: 3 | Success: 3
##        4. iris_stats --> Status: Success
## 

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

statusWF(sal)
## $export_iris
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1 export_iris     Success
## 
## $gzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files     gzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $gunzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files   gunzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $iris_stats
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1  iris_stats     Success

7.2 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) and a template file (see *.tmpl samples here) 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:

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.

8 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.

plotWF(sal, show_legend = TRUE, width = "80%", rstudio = TRUE)

For more details about the plotWF function, please see here.

9 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.

sal <- renderLogs(sal)

10 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.

10.1 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.

sal2rmd(sal)

10.2 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.

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.

11 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.

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()).

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.

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.

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.

sal <- SPRproject(overwrite = TRUE) 

12 Exploring workflow instances

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

12.1 Accessor Methods

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

names(sal)
## [1] "stepsWF"            "statusWF"           "targetsWF"          "outfiles"          
## [5] "SE"                 "dependency"         "targets_connection" "projectInfo"       
## [9] "runInfo"
  • Check the length of the workflow:
length(sal)
## [1] 4
  • Check the steps of the workflow:
stepsWF(sal)
## $export_iris
## Instance of 'LineWise'
##     Code Chunk length: 1
## 
## $gzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 3 (SE...VI), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 3, output: 3
##       cmdlist: 3
##    Sub Steps:
##       1. gzip (rendered: TRUE)
## 
## 
## 
## $gunzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 3 (SE...VI), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 3, output: 3
##       cmdlist: 3
##    Sub Steps:
##       1. gunzip (rendered: TRUE)
## 
## 
## 
## $iris_stats
## Instance of 'LineWise'
##     Code Chunk length: 5
  • 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.

cmdlist(sal, step = c(2,3), targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c  results/setosa.csv > results/SE.csv.gz"
## 
## 
## 
## $gunzip
## $gunzip$SE
## $gunzip$SE$gunzip
## [1] "gunzip -c  ./results/gzip/SE.csv.gz > results/SE.csv"
  • Check the workflow status:
statusWF(sal)
## $export_iris
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1 export_iris     Success
## 
## $gzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files     gzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $gunzip
## DataFrame with 3 rows and 5 columns
##        Targets Total_Files Existing_Files Missing_Files   gunzip
##    <character>   <numeric>      <numeric>     <numeric> <matrix>
## SE          SE           1              1             0  Success
## VE          VE           1              1             0  Success
## VI          VI           1              1             0  Success
## 
## $iris_stats
## DataFrame with 1 row and 2 columns
##          Step      Status
##   <character> <character>
## 1  iris_stats     Success
  • Check the workflow targets files:
targetsWF(sal[2])
## $gzip
## DataFrame with 3 rows and 2 columns
##                  FileName  SampleName
##               <character> <character>
## SE     results/setosa.csv          SE
## VE results/versicolor.csv          VE
## VI  results/virginica.csv          VI
  • 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.

outfiles(sal[2])
## $gzip
## DataFrame with 3 rows and 1 column
##                gzip_file
##              <character>
## 1 ./results/gzip/SE.cs..
## 2 ./results/gzip/VE.cs..
## 3 ./results/gzip/VI.cs..
  • Check the workflow dependencies:
dependency(sal)
## $export_iris
## [1] NA
## 
## $gzip
## [1] "export_iris"
## 
## $gunzip
## [1] "gzip"
## 
## $iris_stats
## [1] "gzip"
  • 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:

targetsheader(sal, step = "Quality")
  • Get the workflow steps names:
stepName(sal)
## [1] "export_iris" "gzip"        "gunzip"      "iris_stats"
  • Get the Sample Id for on particular step:
SampleName(sal, step = "gzip")
## [1] "SE" "VE" "VI"
SampleName(sal, step = "iris_stats")
## NULL
  • Get the outfiles or targets column files:
getColumn(sal, "outfiles", step = "gzip", column = "gzip_file")
##                         SE                         VE                         VI 
## "./results/gzip/SE.csv.gz" "./results/gzip/VE.csv.gz" "./results/gzip/VI.csv.gz"
getColumn(sal, "targetsWF", step = "gzip", column = "FileName")
##                       SE                       VE                       VI 
##     "results/setosa.csv" "results/versicolor.csv"  "results/virginica.csv"
  • Get the R code for a LineWise step:
codeLine(sal, step = "export_iris")
## export_iris
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
  • View all the objects in the running environment:
viewEnvir(sal)
## <environment: 0x56339dbfdc10>
## [1] "df"    "plot"  "stats"
  • Copy one or multiple objects from the running environment to a new environment:
copyEnvir(sal, list = c("plot"), new.env = globalenv(), silent = FALSE)
## <environment: 0x56339dbfdc10>
## Copying to 'new.env': 
## plot
  • Accessing the *.yml data
yamlinput(sal, step = "gzip")
## $file
## $file$class
## [1] "File"
## 
## $file$path
## [1] "_FILE_PATH_"
## 
## 
## $SampleName
## [1] "_SampleName_"
## 
## $ext
## [1] "csv.gz"
## 
## $results_path
## $results_path$class
## [1] "Directory"
## 
## $results_path$path
## [1] "./results"

12.2 Subsetting the workflow details

  • The SYSargsList class and its subsetting operator [:
sal[1]
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Success
## 
sal[1:3]
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Success
##        2. gzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          2.1. gzip
##              cmdlist: 3 | Success: 3
##        3. gunzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          3.1. gunzip
##              cmdlist: 3 | Success: 3
## 
sal[c(1,3)]
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Success
##        2. gunzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          2.1. gunzip
##              cmdlist: 3 | Success: 3
## 
  • The SYSargsList class and its subsetting by steps and input samples:
sal_sub <- subset(sal, subset_steps = c( 2,3), input_targets = ("SE"), keep_steps = TRUE)
stepsWF(sal_sub)
## $export_iris
## Instance of 'LineWise'
##     Code Chunk length: 1
## 
## $gzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 1 (SE...SE), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 1, output: 1
##       cmdlist: 1
##    Sub Steps:
##       1. gzip (rendered: TRUE)
## 
## 
## 
## $gunzip
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 1 (SE...SE), targetsheader: 1 (lines)
##       modules: 0
##       wf: 1, clt: 1, yamlinput: 4 (inputs)
##       input: 1, output: 1
##       cmdlist: 1
##    Sub Steps:
##       1. gunzip (rendered: TRUE)
## 
## 
## 
## $iris_stats
## Instance of 'LineWise'
##     Code Chunk length: 5
targetsWF(sal_sub)
## $export_iris
## DataFrame with 0 rows and 0 columns
## 
## $gzip
## DataFrame with 1 row and 2 columns
##              FileName  SampleName
##           <character> <character>
## SE results/setosa.csv          SE
## 
## $gunzip
## DataFrame with 1 row and 2 columns
##                 gzip_file  SampleName
##               <character> <character>
## SE ./results/gzip/SE.cs..          SE
## 
## $iris_stats
## DataFrame with 0 rows and 0 columns
outfiles(sal_sub)
## $export_iris
## DataFrame with 0 rows and 0 columns
## 
## $gzip
## DataFrame with 1 row and 1 column
##                gzip_file
##              <character>
## 1 ./results/gzip/SE.cs..
## 
## $gunzip
## DataFrame with 1 row and 1 column
##              gunzip_file
##              <character>
## 1 ./results/gunzip/SE...
## 
## $iris_stats
## DataFrame with 0 rows and 0 columns
  • The SYSargsList class and its operator +
sal[1] + sal[2] + sal[3]

12.3 Replacement Methods

  • Update a input parameter in the workflow
sal_c <- sal
## check values
yamlinput(sal_c, step = "gzip")
## $file
## $file$class
## [1] "File"
## 
## $file$path
## [1] "_FILE_PATH_"
## 
## 
## $SampleName
## [1] "_SampleName_"
## 
## $ext
## [1] "csv.gz"
## 
## $results_path
## $results_path$class
## [1] "Directory"
## 
## $results_path$path
## [1] "./results"
## check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c  results/setosa.csv > results/SE.csv.gz"
## Replace
yamlinput(sal_c, step = "gzip", paramName = "ext") <- "txt.gz"

## check NEW values
yamlinput(sal_c, step = "gzip")
## $file
## $file$class
## [1] "File"
## 
## $file$path
## [1] "_FILE_PATH_"
## 
## 
## $SampleName
## [1] "_SampleName_"
## 
## $ext
## [1] "txt.gz"
## 
## $results_path
## $results_path$class
## [1] "Directory"
## 
## $results_path$path
## [1] "./results"
## Check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c  results/setosa.csv > results/SE.txt.gz"
  • Append and Replacement methods for R Code Steps
appendCodeLine(sal_c, step = "export_iris", after = 1) <- "log_cal_100 <- log(100)"
codeLine(sal_c, step = "export_iris")
## export_iris
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
##     log_cal_100 <- log(100)
replaceCodeLine(sal_c, step="export_iris", line = 2) <- LineWise(code={
                    log_cal_100 <- log(50)
                    })
codeLine(sal_c, step = 1)
## export_iris
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
##     log_cal_100 <- log(50)

For more details about the LineWise class, please see below.

  • Rename a Step
renameStep(sal_c, step = 1) <- "newStep"
renameStep(sal_c, c(1, 2)) <- c("newStep2", "newIndex")
sal_c
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. newStep2 --> Status: Success
##        2. newIndex --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          2.1. gzip
##              cmdlist: 3 | Success: 3
##        3. gunzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          3.1. gunzip
##              cmdlist: 3 | Success: 3
##        4. iris_stats --> Status: Success
## 
names(outfiles(sal_c))
## [1] "newStep2"   "newIndex"   "gunzip"     "iris_stats"
names(targetsWF(sal_c))
## [1] "newStep2"   "newIndex"   "gunzip"     "iris_stats"
dependency(sal_c)
## $newStep2
## [1] NA
## 
## $newIndex
## [1] "newStep2"
## 
## $gunzip
## [1] "newIndex"
## 
## $iris_stats
## [1] "newIndex"
  • Replace a Step
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
sal_test <- sal[-2]
sal_test
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. export_iris --> Status: Success
##        2. gunzip --> Status: Success 
##            Total Files: 3 | Existing: 3 | Missing: 0 
##          2.1. gunzip
##              cmdlist: 3 | Success: 3
##        3. iris_stats --> Status: Success
## 

13 CWL syntax

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 and here for Workflows and the user guide here.

CWL command-line specifications are written in YAML 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.

13.1 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:

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.
cwl[1:2]
## $cwlVersion
## [1] "v1.0"
## 
## $class
## [1] "CommandLineTool"
  • baseCommand component provides the name of the software that we desire to execute.
cwl[3]
## $baseCommand
## [1] "echo"
  • 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.
cwl[4]
## $inputs
## $inputs$message
## $inputs$message$type
## [1] "string"
## 
## $inputs$message$inputBinding
## $inputs$message$inputBinding$position
## [1] 1
## 
## 
## 
## $inputs$SampleName
## $inputs$SampleName$type
## [1] "string"
## 
## 
## $inputs$results_path
## $inputs$results_path$type
## [1] "Directory"
  • 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.
cwl[5]
## $outputs
## $outputs$string
## $outputs$string$type
## [1] "stdout"
  • 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.
cwl[6]
## $stdout
## [1] "$(inputs.results_path.basename)/$(inputs.SampleName).txt"

13.2 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.

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.
cwl.wf[1:2]
## $class
## [1] "Workflow"
## 
## $cwlVersion
## [1] "v1.0"
  • The inputs section describes the inputs of the workflow.
cwl.wf[3]
## $inputs
## $inputs$message
## [1] "string"
## 
## $inputs$SampleName
## [1] "string"
## 
## $inputs$results_path
## [1] "Directory"
  • The outputs section describes the outputs of the workflow.
cwl.wf[4]
## $outputs
## $outputs$string
## $outputs$string$outputSource
## [1] "echo/string"
## 
## $outputs$string$type
## [1] "stdout"
  • The steps section describes the steps of the workflow. In this simple example, we demonstrate one step.
cwl.wf[5]
## $steps
## $steps$echo
## $steps$echo$`in`
## $steps$echo$`in`$message
## [1] "message"
## 
## $steps$echo$`in`$SampleName
## [1] "SampleName"
## 
## $steps$echo$`in`$results_path
## [1] "results_path"
## 
## 
## $steps$echo$out
## [1] "[string]"
## 
## $steps$echo$run
## [1] "example/example.cwl"

13.3 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:

yaml::read_yaml(file.path(dir_path, "example/example_single.yml"))
## $message
## [1] "Hello World!"
## 
## $SampleName
## [1] "M1"
## 
## $results_path
## $results_path$class
## [1] "Directory"
## 
## $results_path$path
## [1] "./results"

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

14 How to connect CWL description files within systemPipeR

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 and here.

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

HW <- SYSargsList(wf_file = "example/workflow_example.cwl", 
                  input_file = "example/example_single.yml", 
                  dir_path = system.file("extdata/cwl", package = "systemPipeR"))
HW
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. Step_x --> Status: Pending 
##            Total Files: 1 | Existing: 0 | Missing: 1 
##          1.1. echo
##              cmdlist: 1 | Pending: 1
## 
cmdlist(HW)
## $Step_x
## $Step_x$defaultid
## $Step_x$defaultid$echo
## [1] "echo Hello World! > results/M1.txt"

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:

yml <- yaml::read_yaml(file.path(dir_path, "example/example.yml"))
yml
## $message
## [1] "_STRING_"
## 
## $SampleName
## [1] "_SAMPLE_"
## 
## $results_path
## $results_path$class
## [1] "Directory"
## 
## $results_path$path
## [1] "./results"

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:

targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR")
read.delim(targetspath, comment.char = "#")
##               Message SampleName
## 1        Hello World!         M1
## 2          Hello USA!         M2
## 3 Hello Bioconductor!         M3

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.

14.0.1 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.

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
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. echo --> Status: Pending 
##            Total Files: 3 | Existing: 0 | Missing: 3 
##          1.1. echo
##              cmdlist: 3 | Pending: 3
## 
cmdlist(HW_mul)
## $echo
## $echo$M1
## $echo$M1$echo
## [1] "echo Hello World! > results/M1.txt"
## 
## 
## $echo$M2
## $echo$M2$echo
## [1] "echo Hello USA! > results/M2.txt"
## 
## 
## $echo$M3
## $echo$M3$echo
## [1] "echo Hello Bioconductor! > results/M3.txt"

15 Creating the CWL param files from the command-line

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

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>
"

15.1 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 <...>

15.2 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.

cmd <- createParam(command, writeParamFiles = FALSE)
## *****BaseCommand*****
## hisat2 
## *****Inputs*****
## S:
##     type: File
##     preF: -S
##     yml: ./results/M1A.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## k:
##     type: int
##     preF: -k
##     yml: 1
## min-intronlen:
##     type: int
##     preF: -min-intronlen
##     yml: 30
## max-intronlen:
##     type: int
##     preF: -max-intronlen
##     yml: 3000
## threads:
##     type: int
##     preF: -threads
##     yml: 4
## U:
##     type: File
##     preF: -U
##     yml: ./data/SRR446027_1.fastq.gz
## *****Outputs*****
## output1:
##     type: File
##     value: ./results/M1A.sam
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz

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

writeParamFiles(cmd, overwrite = TRUE)

15.3 How to access and edit param files

15.3.2 Subsetting the command-line

cmd2 <- subsetParam(cmd, position = "inputs", index = 1:2, trim = TRUE)
## *****Inputs*****
## S:
##     type: File
##     preF: -S
##     yml: ./results/M1A.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta
cmdlist(cmd2)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta"
cmd2 <- subsetParam(cmd, position = "inputs", index = c("S", "x"), trim = TRUE)
## *****Inputs*****
## S:
##     type: File
##     preF: -S
##     yml: ./results/M1A.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta
cmdlist(cmd2)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta"

15.3.3 Replacing a existing argument in the command-line

cmd3 <- replaceParam(cmd, "base", index = 1, replace = list(baseCommand = "bwa"))
## Replacing baseCommand
## *****BaseCommand*****
## bwa 
## *****Parsed raw command line*****
## bwa -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd3)
## $defaultid
## $defaultid$hisat2
## [1] "bwa -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
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)
## Replacing inputs
## *****Inputs*****
## new_input1:
##     type: File
##     preF: -b
##     yml: myfile
## new_input2:
##     type: int
##     preF: -L
##     yml: 4
## k:
##     type: int
##     preF: -k
##     yml: 1
## min-intronlen:
##     type: int
##     preF: -min-intronlen
##     yml: 30
## max-intronlen:
##     type: int
##     preF: -max-intronlen
##     yml: 3000
## threads:
##     type: int
##     preF: -threads
##     yml: 4
## U:
##     type: File
##     preF: -U
##     yml: ./data/SRR446027_1.fastq.gz
## *****Parsed raw command line*****
## hisat2 -b myfile -L 4 -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd4)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -b myfile -L 4 -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"

15.3.4 Adding new arguments

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)
## Replacing inputs
## *****Inputs*****
## S:
##     type: File
##     preF: -S
##     yml: ./results/M1A.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## k:
##     type: int
##     preF: -k
##     yml: 1
## min-intronlen:
##     type: int
##     preF: -min-intronlen
##     yml: 30
## max-intronlen:
##     type: int
##     preF: -max-intronlen
##     yml: 3000
## threads:
##     type: int
##     preF: -threads
##     yml: 4
## U:
##     type: File
##     preF: -U
##     yml: ./data/SRR446027_1.fastq.gz
## new_input1:
##     type: File
##     preF: -b1
##     yml: myfile1
## new_input2:
##     type: File
##     preF: -b2
##     yml: myfile2
## new_input3:
##     type: File
##     preF: -b3
##     yml: myfile3
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz -b1 myfile1 -b2 myfile2 -b3 myfile3
cmdlist(cmd5)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz -b1 myfile1 -b2 myfile2 -b3 myfile3"
cmd6 <- appendParam(cmd, "inputs", index = 1:2, after=0, append = new_inputs)
## Replacing inputs
## *****Inputs*****
## new_input1:
##     type: File
##     preF: -b1
##     yml: myfile1
## new_input2:
##     type: File
##     preF: -b2
##     yml: myfile2
## new_input3:
##     type: File
##     preF: -b3
##     yml: myfile3
## S:
##     type: File
##     preF: -S
##     yml: ./results/M1A.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## k:
##     type: int
##     preF: -k
##     yml: 1
## min-intronlen:
##     type: int
##     preF: -min-intronlen
##     yml: 30
## max-intronlen:
##     type: int
##     preF: -max-intronlen
##     yml: 3000
## threads:
##     type: int
##     preF: -threads
##     yml: 4
## U:
##     type: File
##     preF: -U
##     yml: ./data/SRR446027_1.fastq.gz
## *****Parsed raw command line*****
## hisat2 -b1 myfile1 -b2 myfile2 -b3 myfile3 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd6)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -b1 myfile1 -b2 myfile2 -b3 myfile3 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"

15.3.5 Editing output param

new_outs <- list(
    "sam_out" = "<F: $(inputs.results_path)/test.sam>"
) 
cmd7 <- replaceParam(cmd, "outputs", index = 1, replace = new_outs)
## Replacing outputs
## *****Outputs*****
## sam_out:
##     type: File
##     value: $(inputs.results_path)/test.sam
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
output(cmd7) 
## $defaultid
## $defaultid$hisat2
## [1] "./results/test.sam"

15.3.6 Internal Check

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) 
## *****BaseCommand*****
## hisat2 
## *****Inputs*****
## S:
##     type: File
##     preF: -S
##     yml: _SampleName_.sam
## x:
##     type: File
##     preF: -x
##     yml: ./data/tair10.fasta
## k:
##     type: int
##     preF: -k
##     yml: 1
## min-intronlen:
##     type: int
##     preF: -min-intronlen
##     yml: 30
## max-intronlen:
##     type: int
##     preF: -max-intronlen
##     yml: 3000
## threads:
##     type: int
##     preF: -threads
##     yml: 4
## U:
##     type: File
##     preF: -U
##     yml: _FASTQ_PATH1_
## *****Outputs*****
## output1:
##     type: File
##     value: _SampleName_.sam
## *****Parsed raw command line*****
## hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U _FASTQ_PATH1_ 
##   Written content of 'commandLine' to file: 
##  param/cwl/hisat2/hisat2.cwl 
##   Written content of 'commandLine' to file: 
##  param/cwl/hisat2/hisat2.yml
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
## Instance of 'SYSargs2':
##    Slot names/accessors: 
##       targets: 18 (M1A...V12B), targetsheader: 4 (lines)
##       modules: 1
##       wf: 0, clt: 1, yamlinput: 9 (inputs)
##       input: 18, output: 18
##       cmdlist: 18
##    Sub Steps:
##       1. hisat2 (rendered: TRUE)
cmdlist(WF_test)[1:2]
## $M1A
## $M1A$hisat2
## [1] "hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
## 
## 
## $M1B
## $M1B$hisat2
## [1] "hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446028_1.fastq.gz"

16 Inner Classes

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

16.1 SYSargs2 Class

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 (Kim, Langmead, and Salzberg 2015). For more details about the file structure and how to design or customize our own software tools, please check systemPipeR and CWL pipeline.

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

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.

names(WF)
##  [1] "targets"           "targetsheader"     "modules"           "wf"               
##  [5] "clt"               "yamlinput"         "cmdlist"           "input"            
##  [9] "output"            "files"             "inputvars"         "cmdToCwl"         
## [13] "status"            "internal_outfiles"

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.

cmdlist(WF)[1]
## $M1A
## $M1A$`hisat2-mapping-se`
## [1] "hisat2 -S ./results/M1A.sam  -x ./data/tair10.fasta  -k 1  --min-intronlen 30  --max-intronlen 3000  -U ./data/SRR446027_1.fastq.gz --threads 4"

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.

output(WF)[1]
## $M1A
## $M1A$`hisat2-mapping-se`
## [1] "./results/M1A.sam"

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.
targets(WF)[1]
## $M1A
## $M1A$FileName
## [1] "./data/SRR446027_1.fastq.gz"
## 
## $M1A$SampleName
## [1] "M1A"
## 
## $M1A$Factor
## [1] "M1"
## 
## $M1A$SampleLong
## [1] "Mock.1h.A"
## 
## $M1A$Experiment
## [1] 1
## 
## $M1A$Date
## [1] "23-Mar-2012"
as(WF, "DataFrame")
## DataFrame with 18 rows and 6 columns
##                   FileName  SampleName      Factor  SampleLong  Experiment        Date
##                <character> <character> <character> <character> <character> <character>
## 1   ./data/SRR446027_1.f..         M1A          M1   Mock.1h.A           1 23-Mar-2012
## 2   ./data/SRR446028_1.f..         M1B          M1   Mock.1h.B           1 23-Mar-2012
## 3   ./data/SRR446029_1.f..         A1A          A1    Avr.1h.A           1 23-Mar-2012
## 4   ./data/SRR446030_1.f..         A1B          A1    Avr.1h.B           1 23-Mar-2012
## 5   ./data/SRR446031_1.f..         V1A          V1    Vir.1h.A           1 23-Mar-2012
## ...                    ...         ...         ...         ...         ...         ...
## 14  ./data/SRR446040_1.f..        M12B         M12  Mock.12h.B           1 23-Mar-2012
## 15  ./data/SRR446041_1.f..        A12A         A12   Avr.12h.A           1 23-Mar-2012
## 16  ./data/SRR446042_1.f..        A12B         A12   Avr.12h.B           1 23-Mar-2012
## 17  ./data/SRR446043_1.f..        V12A         V12   Vir.12h.A           1 23-Mar-2012
## 18  ./data/SRR446044_1.f..        V12B         V12   Vir.12h.B           1 23-Mar-2012

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 is available, it is possible to define which module should be loaded, as shown here:

modules(WF)
##        module1 
## "hisat2/2.1.0"

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

files(WF)
inputvars(WF)

16.2 LineWise Class

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

rmd <- system.file("extdata", "spr_simple_lw.Rmd", package = "systemPipeR")
sal_lw <- SPRproject(overwrite = TRUE)
## Recreating directory '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject'
## Creating file '/tmp/RtmpUy9Tuw/Rbuild3aef6f23436e71/systemPipeR/vignettes/spr_project/.SPRproject/SYSargsList.yml'
sal_lw <- importWF(sal_lw, rmd, verbose = FALSE)
## Now check if required tools are installed 
## There is no commandline (SYSargs) step in this workflow, skip.
codeLine(sal_lw)
## firstStep
##     mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
## secondStep
##     setosa <- read.delim("results/setosa.csv", sep = ",")
##     versicolor <- read.delim("results/versicolor.csv", sep = ",")
##     virginica <- read.delim("results/virginica.csv", sep = ",")
  • Coerce methods available:
lw <- stepsWF(sal_lw)[[2]]
## Coerce
ll <- as(lw, "list")
class(ll)
## [1] "list"
lw <- as(ll, "LineWise")
lw
## Instance of 'LineWise'
##     Code Chunk length: 3
  • Access details
length(lw)
## [1] 3
names(lw)
## [1] "codeLine"       "codeChunkStart" "stepName"       "dependency"     "status"        
## [6] "files"          "runInfo"
codeLine(lw)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## versicolor <- read.delim("results/versicolor.csv", sep = ",")
## virginica <- read.delim("results/virginica.csv", sep = ",")
codeChunkStart(lw)
## integer(0)
rmdPath(lw)
## character(0)
  • Subsetting
l <- lw[2]
codeLine(l)
## versicolor <- read.delim("results/versicolor.csv", sep = ",")
l_sub <- lw[-2]
codeLine(l_sub)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## virginica <- read.delim("results/virginica.csv", sep = ",")
  • Replacement methods
replaceCodeLine(lw, line = 2) <- "5+5"
codeLine(lw)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## 5 + 5
## virginica <- read.delim("results/virginica.csv", sep = ",")
appendCodeLine(lw, after = 0) <- "6+7"
codeLine(lw)
## 6 + 7
## setosa <- read.delim("results/setosa.csv", sep = ",")
## 5 + 5
## virginica <- read.delim("results/virginica.csv", sep = ",")
  • Replacement methods for SYSargsList
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)

16.3 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.

17 Third-party software tools

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

Tool Name Description Step
bwa BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome.  Alignment
Bowtie2 Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. Alignment
FASTX-Toolkit FASTX-Toolkit is a collection of command line tools for Short-Reads FASTA/FASTQ files preprocessing. Read Preprocessing
TransRate Transrate is software for de-novo transcriptome assembly quality analysis. Quality
Gsnap GSNAP is a genomic short-read nucleotide alignment program. Alignment
Samtools Samtools is a suite of programs for interacting with high-throughput sequencing data. Post-processing
Trimmomatic Trimmomatic is a flexible read trimming tool for Illumina NGS data. Read Preprocessing
Rsubread Rsubread is a Bioconductor software package that provides high-performance alignment and read counting functions for RNA-seq reads. Alignment
Picard Picard is a set of command line tools for manipulating high-throughput sequencing (HTS) data and formats such as SAM/BAM/CRAM and VCF. Manipulating HTS data
Busco BUSCO assesses genome assembly and annotation completeness with Benchmarking Universal Single-Copy Orthologs. Quality
Hisat2 HISAT2 is a fast and sensitive alignment program for mapping NGS reads (both DNA and RNA) to reference genomes. Alignment
Tophat2 TopHat is a fast splice junction mapper for RNA-Seq reads. Alignment
GATK Variant Discovery in High-Throughput Sequencing Data. Variant Discovery
Trim_galore Trim Galore is a wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files. Read Preprocessing
TransDecoder TransDecoder identifies candidate coding regions within transcript sequences. Find Coding Regions
Trinotate Trinotate is a comprehensive annotation suite designed for automatic functional annotation of transcriptomes. Transcriptome Functional Annotation
STAR STAR is an ultrafast universal RNA-seq aligner. Alignment
Trinity Trinity assembles transcript sequences from Illumina RNA-Seq data. denovo Transcriptome Assembly
MACS2 MACS2 identifies transcription factor binding sites in ChIP-seq data. Peak calling
Kallisto kallisto is a program for quantifying abundances of transcripts from RNA-Seq data. Read counting
BCFtools BCFtools is a program for variant calling and manipulating files in the Variant Call Format (VCF) and its binary counterpart BCF. Variant Discovery
Bismark Bismark is a program to map bisulfite treated sequencing reads to a genome of interest and perform methylation calls in a single step. Bisulfite mapping
Fastqc FastQC is a quality control tool for high throughput sequence data. Quality
Blast BLAST finds regions of similarity between biological sequences. Blast

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.

There are a few ways to check if the required tools/modules are installed. The easiest way is automatically performed for users by calling the importWF function. At the end of the import, all required tools and modules are automatically listed and checked for users.

There are a few other methods that one could use to perform the tool validation, please read details on our website, the Before running section.

18 Workflow commom steps overview

18.1 Project initialization

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

sal <- SPRproject() 

18.1.1 Required packages and resources

The systemPipeR package needs to be loaded (H Backman and Girke 2016).

appendStep(sal) <- LineWise({
                            library(systemPipeR)
                            }, 
                            step_name = "load_SPR")

18.2 Read Preprocessing

18.2.1 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.

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:

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:

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

18.2.2 Preprocessing with TrimGalore!

TrimGalore! 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.

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")

18.2.3 Preprocessing with Trimmomatic

Trimmomatic software (Bolger, Lohse, and Usadel 2014) 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.

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")

18.3 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.

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")
Figure 1: FASTQ quality report


18.4 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.

18.4.1 Alignment with HISAT2

The following steps will demonstrate how to use the short read aligner Hisat2 (Kim, Langmead, and Salzberg 2015) 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.
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.

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")

18.4.2 Alignment with Tophat2

The NGS reads of this project can also be aligned against the reference genome sequence using Bowtie2/TopHat2 (Kim et al. 2013; Langmead and Salzberg 2012).

  • Build Bowtie2 index.
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.

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")

18.4.3 Alignment with Bowtie2 (e.g. for miRNA profiling)

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

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")

18.4.4 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:
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:
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")

18.4.5 Alignment with Rsubread (e.g. for RNA-Seq)

The following example shows how one can use within the environment the R-based aligner , allowing running from R or command-line.

  • Build the index:
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:
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")

18.4.6 Alignment with gsnap (e.g. for VAR-Seq and RNA-Seq)

Another R-based short read aligner is gsnap from the gmapR package (Wu and Nacu 2010). The code sample below introduces how to run this aligner on multiple nodes of a compute cluster.

  • Build the index:
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:
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")

18.6 Read counting for mRNA profiling experiments

Create txdb (needs to be done only once).

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.

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.

18.6.0.1 Read and alignment count stats

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

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.

read.table(system.file("extdata", "alignStats.xls", package = "systemPipeR"), header = TRUE)[1:4,]
##   FileName Nreads2x Nalign Perc_Aligned Nalign_Primary Perc_Aligned_Primary
## 1      M1A   192918 177961     92.24697         177961             92.24697
## 2      M1B   197484 159378     80.70426         159378             80.70426
## 3      A1A   189870 176055     92.72397         176055             92.72397
## 4      A1B   188854 147768     78.24457         147768             78.24457

18.7 Read counting for miRNA profiling experiments

Download miRNA genes from miRBase.

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")

18.8 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).

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")
Figure 2: Correlation dendrogram of samples for rlog values.


18.9 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 (Robinson, McCarthy, and Smyth 2010) for any number of pairwise sample comparisons specified under the cmp argument. Users are strongly encouraged to consult the edgeR vignette for more detailed information on this topic and how to properly run edgeR on data sets with more complex experimental designs.

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.

Figure 3: Up and down regulated DEGs identified by edgeR.


18.10 DEG analysis with DESeq2

The following run_DESeq2 function is a convenience wrapper for identifying DEGs in batch mode with DESeq2 (Love, Huber, and Anders 2014) for any number of pairwise sample comparisons specified under the cmp argument. Users are strongly encouraged to consult the DESeq2 vignette for more detailed information on this topic and how to properly run DESeq2 on data sets with more complex experimental designs.

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")

18.11 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).

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")
Figure 4: Venn Diagram for 4 Up and Down DEG Sets.


18.12 GO term enrichment analysis of DEGs

18.12.1 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.

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")

18.12.2 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.

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")

18.12.3 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.

Figure 5: GO Slim Barplot for MF Ontology.


18.13 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.

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"))
Figure 7: Heat map with hierarchical clustering dendrograms of DEGs.


18.14 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.

plotWF(sal, show_legend = TRUE, width = "80%")

To check more details of plotWF, visit our website.

18.15 Running workflow

18.15.1 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.

sal <- runWF(sal)

18.15.2 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) and a template file (see *.tmpl samples here) 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:

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)

18.15.3 Checking workflow status

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

sal
statusWF(sal)

18.15.4 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.

sal <- renderLogs(sal)

19 Version information

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB             
##  [4] LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] magrittr_2.0.3              systemPipeR_2.8.0           ShortRead_1.60.0           
##  [4] GenomicAlignments_1.38.0    SummarizedExperiment_1.32.0 Biobase_2.62.0             
##  [7] MatrixGenerics_1.14.0       matrixStats_1.0.0           BiocParallel_1.36.0        
## [10] Rsamtools_2.18.0            Biostrings_2.70.0           XVector_0.42.0             
## [13] GenomicRanges_1.54.0        GenomeInfoDb_1.38.0         IRanges_2.36.0             
## [16] S4Vectors_0.40.0            BiocGenerics_0.48.0         BiocStyle_2.30.0           
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0        viridisLite_0.4.2       dplyr_1.1.3             farver_2.1.1           
##  [5] bitops_1.0-7            fastmap_1.1.1           RCurl_1.98-1.12         digest_0.6.33          
##  [9] lifecycle_1.0.3         ellipsis_0.3.2          compiler_4.3.1          rlang_1.1.1            
## [13] sass_0.4.7              tools_4.3.1             utf8_1.2.4              yaml_2.3.7             
## [17] systemPipeRdata_2.5.1   knitr_1.44              S4Arrays_1.2.0          labeling_0.4.3         
## [21] htmlwidgets_1.6.2       interp_1.1-4            DelayedArray_0.28.0     xml2_1.3.5             
## [25] RColorBrewer_1.1-3      abind_1.4-5             withr_2.5.1             hwriter_1.3.2.1        
## [29] grid_4.3.1              fansi_1.0.5             latticeExtra_0.6-30     colorspace_2.1-0       
## [33] ggplot2_3.4.4           scales_1.2.1            cli_3.6.1               rmarkdown_2.25         
## [37] crayon_1.5.2            generics_0.1.3          remotes_2.4.2.1         rstudioapi_0.15.0      
## [41] httr_1.4.7              cachem_1.0.8            stringr_1.5.0           zlibbioc_1.48.0        
## [45] rvest_1.0.3             parallel_4.3.1          BiocManager_1.30.22     vctrs_0.6.4            
## [49] webshot_0.5.5           Matrix_1.6-1.1          jsonlite_1.8.7          bookdown_0.36          
## [53] systemfonts_1.0.5       jpeg_0.1-10             magick_2.8.1            crosstalk_1.2.0        
## [57] jquerylib_0.1.4         glue_1.6.2              codetools_0.2-19        DT_0.30                
## [61] stringi_1.7.12          gtable_0.3.4            deldir_1.0-9            munsell_0.5.0          
## [65] tibble_3.2.1            pillar_1.9.0            htmltools_0.5.6.1       GenomeInfoDbData_1.2.11
## [69] R6_2.5.1                evaluate_0.22           kableExtra_1.3.4        lattice_0.22-5         
## [73] png_0.1-8               bslib_0.5.1             Rcpp_1.0.11             svglite_2.1.2          
## [77] SparseArray_1.2.0       xfun_0.40               pkgconfig_2.0.3

20 Funding

This project is funded by NSF award ABI-1661152.

References

Amstutz, Peter, Michael R Crusoe, Nebojša Tijanić, Brad Chapman, John Chilton, Michael Heuer, Andrey Kartashov, et al. 2016. “Common Workflow Language, V1.0,” July. https://doi.org/10.6084/m9.figshare.3115156.v2.

Bolger, Anthony M, Marc Lohse, and Bjoern Usadel. 2014. “Trimmomatic: A Flexible Trimmer for Illumina Sequence Data.” Bioinformatics 30 (15): 2114–20.

Crusoe, Michael R, Sanne Abeln, Alexandru Iosup, Peter Amstutz, John Chilton, Nebojša Tijanić, Hervé Ménager, Stian Soiland-Reyes, Bogdan Gavrilovic, and Carole Goble. 2021. “Methods Included: Standardizing Computational Reuse and Portability with the Common Workflow Language,” May. http://arxiv.org/abs/2105.07028.

H Backman, Tyler W, and Thomas Girke. 2016. “systemPipeR: NGS workflow and report generation environment.” BMC Bioinformatics 17 (1): 388. https://doi.org/10.1186/s12859-016-1241-0.

Howard, Brian E, Qiwen Hu, Ahmet Can Babaoglu, Manan Chandra, Monica Borghi, Xiaoping Tan, Luyan He, et al. 2013. “High-Throughput RNA Sequencing of Pseudomonas-Infected Arabidopsis Reveals Hidden Transcriptome Complexity and Novel Splice Variants.” PLoS One 8 (10): e74183. https://doi.org/10.1371/journal.pone.0074183.

Kim, Daehwan, Ben Langmead, and Steven L Salzberg. 2015. “HISAT: A Fast Spliced Aligner with Low Memory Requirements.” Nat. Methods 12 (4): 357–60.

Kim, Daehwan, Geo Pertea, Cole Trapnell, Harold Pimentel, Ryan Kelley, and Steven L Salzberg. 2013. “TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions.” Genome Biol. 14 (4): R36. https://doi.org/10.1186/gb-2013-14-4-r36.

Langmead, Ben, and Steven L Salzberg. 2012. “Fast Gapped-Read Alignment with Bowtie 2.” Nat. Methods 9 (4): 357–59. https://doi.org/10.1038/nmeth.1923.

Love, Michael, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biol. 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

Robinson, M D, D J McCarthy, and G K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.

Wu, T D, and S Nacu. 2010. “Fast and SNP-tolerant Detection of Complex Variants and Splicing in Short Reads.” Bioinformatics 26 (7): 873–81. https://doi.org/10.1093/bioinformatics/btq057.