--- title: "systemPipeR: Workflow Environment for Data Analysis and Report Generation" author: "Author: Le Zhang, Daniela Cassol, and Thomas Girke" date: "Last update: `r format(Sys.time(), '%d %B, %Y')`" output: BiocStyle::html_document: toc_float: true code_folding: show package: systemPipeR vignette: | %\VignetteEncoding{UTF-8} %\VignetteIndexEntry{Overview} %\VignetteEngine{knitr::rmarkdown} fontsize: 14pt bibliography: bibtex.bib editor_options: markdown: wrap: 80 chunk_output_type: console --- ```{css, echo=FALSE} pre code { white-space: pre !important; overflow-x: scroll !important; word-break: keep-all !important; word-wrap: initial !important; } ``` ```{r style, echo = FALSE, results = 'asis'} BiocStyle::markdown() options(width = 100, max.print = 1000) knitr::opts_chunk$set( eval = as.logical(Sys.getenv("KNITR_EVAL", "TRUE")), cache = as.logical(Sys.getenv("KNITR_CACHE", "TRUE")), tidy.opts = list(width.cutoff = 100), tidy = FALSE) ``` ```{r setting, eval=TRUE, echo=FALSE} if (file.exists("spr_project")) unlink("spr_project", recursive = TRUE) is_win <- Sys.info()[['sysname']] != "Windows" ``` ```{r load_library, eval=TRUE, include=FALSE} suppressPackageStartupMessages({ library(systemPipeR) }) ``` # Introduction [_`systemPipeR`_](http://www.bioconductor.org/packages/devel/bioc/html/systemPipeR.html) is a versatile workflow environment for data analysis that integrates R with command-line (CL) software [@H_Backman2016-bt]. This platform allows scientists to analyze diverse data types on personal or distributed computer systems. It ensures a high level of reproducibility, scalability, and portability (Figure \@ref(fig:utilities)). Central to `systemPipeR` is a CL interface (CLI) that adopts the Common Workflow Language [CWL, @Crusoe2021-iq]. Using this CLI, users can select the optimal R or CL software for each analysis step. The platform supports end-to-end and partial execution of workflows, with built-in restart capabilities. A workflow control container class manages analysis tasks of varying complexity. Standardized processing routines for metadata facilitate the handling of large numbers of input samples and complex experimental designs. As a multipurpose workflow management toolkit, `systemPipeR` enables users to run existing workflows, customize them, or create entirely new ones while leveraging widely adopted data structures within the Bioconductor ecosystem. Another key aspect of `systemPipeR` is its ability to generate reproducible scientific analysis and technical reports. For result interpretation, it offers a range of graphics functionalities. Additionally, an associated Shiny App provides various interactive features for result exploration, and enhancing the user experience. ```{r utilities, eval=TRUE, warning= FALSE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Important functionalities of systemPipeR. (A) Illustration of workflow design concepts, and (B) examples of visualization functionalities for NGS data.", warning=FALSE} knitr::include_graphics("images/utilities.png") ``` ## Workflow management container A central component of `systemPipeR` is `SYSargsList` or short `SAL`, a container for workflow management. This S4 class stores all relevant information for running and monitoring workflows. It captures the connectivity between workflow steps, the paths to their input and output data, and pertinent parameter values used in each step (see Figure \@ref(fig:sysargslistImage)). `SAL` instances can be constructed from a specific metadata table, referred to as targets file, R code and/or CWL parameter files (details provided below). For preconfigured NGS workflows, users need to provide only a targets file and the initial input data specified in that file (such as FASTQ files). The targets file provides information about the input data and additional metadata describing the experimental design, such as sample labels, replicate information, and other relevant details. As the workflow progresses, subsequent input and output data are generated by the individual analysis steps, that can be retrieved as descendent targets instances. ```{r sysargslistImage, warning= FALSE, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Workflow management class. Workflows are defined and managed by the `SYSargsList` (`SAL`) control class. Components of `SAL` include `SYSargs2` and/or `LineWise` for defining CL- and R-based workflow steps, respectively. The former are constructed from a `targets` and two CWL *param* files, and the latter comprises mainly R code.", warning=FALSE} knitr::include_graphics("images/SYSargsList.png") ``` ## CL interface (CLI) _`systemPipeR`_ adopts the Common Workflow Language (CWL) [@Amstutz2016-ka], which is a widely used community standard for describing CL tools and workflows in a declarative, generic, and reproducible manner. CWL specifications are text-based YAML (https://yaml.org/) files that are straightforward to create and to modify. Integrating CWL in `systemPipeR` enhances the sharability, standardization, extensibility and portability of data analysis workflows. Following the [CWL Specifications](https://www.commonwl.org/v1.2/CommandLineTool.html), the basic description for executing a command-line software are defined by two files: a cwl step definition file and a yml configuration file. Figure \@ref(fig:sprandCWL) illustrates the utilitity of the two files using “Hello World” as an example. The cwl file (A) defines the CL tool (C) along with its parameters, and the yml file (B) assigns values to the corresponding parameters. For convenience, parameter values can be provided by the targets file (D, see above), and automatically passed on to the corresponding parameters in the yml file. The usage of a targets file greatly simplifies the operation of the system for users, because a tabular metadata file is intuitive to maintain, and it eliminates the need of modifying the more complex cwl and yml files directly. ```{r sprandCWL, warning=FALSE, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Parameter files. Illustration how the different fields in cwl, yml and targets files are connected to assemble command-line calls, here for 'Hello World' example.", warning=FALSE} knitr::include_graphics("images/SPR_CWL_hello.png") ``` ## Additional important functionalities The package also provides several convenience functions that are useful for designing and testing workflows, such as a command-line rendering function that assembles from the parameter files (cwl, yml and targets) the exact command-line strings for each step prior to running a command-line tool. Auto-generation of CWL parameter files is also supported. Here, users can simply provide the command-line strings for a command-line software of interest to a rendering function that generates the corresponding `*.cwl` and `*.yml` files for them. # Getting started ## Installation The [_`systemPipeR`_](http://www.bioconductor.org/packages/devel/bioc/html/systemPipeR.html) package can be installed from the R console using the [_`BiocManager::install`_](https://cran.r-project.org/web/packages/BiocManager/index.html) command. The associated data package [_`systemPipeRdata`_](http://www.bioconductor.org/packages/devel/data/experiment/html/systemPipeRdata.html) can be installed the same way. The latter is a helper package for generating _`systemPipeR`_ workflow test instances with a single command. These instances contain all parameter files and sample data required to quickly test and run workflows. ```{r install, eval=FALSE} if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install("systemPipeR") BiocManager::install("systemPipeRdata") ``` For a workflow to run successfully, all CL tools used by a workflow need to be installed and executable on the user's system, where the analysis will be performed (details provided [below](#third-party-software-tools)). ## Five minute quick-start The following demonstrates how to initialize, run and monitor workflows, and subsequently create analysis reports. __1. Create workflow environment.__ The chosen example uses the `genWorenvir` function from the `systemPipeRdata` package to create an RNA-Seq workflow environment that is fully populated with test data, including FASTQ files, reference genome and annotation data. After this, the user's R session needs to be directed into the resulting `rnaseq` directory (here with `setwd`). ```{r eval=FALSE} systemPipeRdata::genWorkenvir(workflow = "rnaseq") setwd("rnaseq") ``` __2. Initialize project and import workflow from `Rmd` template.__ New workflow instances are created with the `SPRproject` function. When calling this function, a project directory with the default name `.SPRproject` is created within the workflow directory. Progress information and log files of a workflow run will be stored in this directory. After this, workflow steps can be loaded into `sal` one-by-one, or all at once with the `importWF` function. The latter reads all steps of a from a workflow Rmd file (here `systemPipeRNAseq.Rmd`) defining the analysis steps. ```{r eval=FALSE} library(systemPipeR) # Initialize workflow project sal <- SPRproject() ## Creating directory '/home/myuser/systemPipeR/rnaseq/.SPRproject' ## Creating file '/home/myuser/systemPipeR/rnaseq/.SPRproject/SYSargsList.yml' sal <- importWF(sal, file_path = "systemPipeRNAseq.Rmd") # import into sal the WF steps defined by chosen Rmd file ``` ![](https://systempipe.org/sp/spr/sp_run/listCmdTools.png) The `importWF` function also checks the availability of required CL software. All dependency CL software needs to be installed and exported to a user's `PATH`. In the given example, the CL tools `trimmomatic`, `hisat2-build`, `hisat2`, and `samtools` are not available. This is expected here since this vignette was rendered on Bioconductor's package build system where custom command-line software is not necessarily installed. __3. Status summary.__ An overview of the workflow steps and their status information can be printed by typing `sal`. For space reasons, the following shows only the first three steps of the RNA-Seq workflow. At this stage all workflow steps are listed as pending since none of them have been executed yet. ```{r eval=FALSE} sal ## Instance of 'SYSargsList': ## WF Steps: ## 1. load_SPR --> Status: Pending ## 2. preprocessing --> Status: Pending ## Total Files: 36 | Existing: 0 | Missing: 36 ## 2.1. preprocessReads-pe ## cmdlist: 18 | Pending: 18 ## 3. trimming --> Status: Pending ## Total Files: 72 | Existing: 0 | Missing: 72 ## ... ``` __4. Run workflow.__ Next, one can execute the entire workflow from start to finish. The `steps` argument of `runWF` can be used to run only selected steps. For details, consult the help file with `?runWF`. During the run detailed status information will be provided for each workflow step. ```{r eval=FALSE} sal <- runWF(sal) ``` After completing all or only some steps, the status of workflow steps can always be checked with the summary print function. If a workflow step was completed, its status will change from `Pending` to `Success` or `Failed`. ```{r eval=FALSE} sal ``` ![](https://systempipe.org/sp/spr/sp_run/runwf.png) __5. Workflow topology graph.__ Workflows can be displayed as topology graphs using the `plotWF` function (not evaluated here). The run status information for each step and various other details are embedded in these graphs. Examples of the workflow plot are available in the [visualize workflow section](#visualize-workflow) below. ```{r eval=FALSE} plotWF(sal) ``` __6. Report generation.__ The `renderReport` and `renderLogs` function can be used for generating scientific and technical reports, respectively. ```{r eval=FALSE} # Scietific report sal <- renderReport(sal) # Technical (log) report sal <- renderLogs(sal) ``` # Project structure The root directory of `systemPipeR` projects contains by default three user facing sub-directories: `data`, `results` and `param`. A fourth sub-directory is a log directory with the default name `.SPRproject` that is created when initializing a workflow run with the `SPRproject` function (see above). Users can change the recommended directory structure, but will need to adjust in some cases the code in their workflows. Just adding directories is possible without requiring changes to the workflows. The following directory tree summarizes the expected content in each default directory (names given in ***green***). * _**workflow/**_ + This is the root directory of a workflow. It can have any name and includes the following files: + Workflow *Rmd* and metadata targets file(s) + Optionally, configuration files for computer clusters, such as `.batchtools.conf.R` and `tmpl` files for `batchtools` and `BiocParallel`. + Additional files can be added as needed. + Default sub-directories: + _**param/**_ + CWL parameter files are organized by CL tools (under _**cwl/**_), each with its own sub-directory that contains the corresponding `cwl` and `yml` files. Previous versions of parameter files are stored in a separate sub-directory. + _**data/**_ + Raw input and/or assay data (*e.g.* FASTQ files) + Reference data, including genome sequences, annotation files, databases, etc. + Any number of sub-directories can be added to organize the data under this directory. + Other input data + _**results/**_ + Analysis results are written to this directory. Examples include tables, plots, or NGS results such as alignment (BAM), variant (VCF), peak (BED) files. + Any number of sub-directories can be created to organize the analysis results under this directory. + _**.SPRproject/**_ + Hidden log directory (name starts with a dot) created by `SPRproject` function at the beginning of a workflow run. + Run status information and log files of a workflow run are stored here. ## 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 two _`targets`_ file examples included in the package. They can also be viewed and downloaded from _`systemPipeR`'s_ GitHub repository [here](https://github.com/tgirke/systemPipeR/blob/master/inst/extdata/targets.txt). As an alternative to using targets files, `YAML` files can be used instead. Since organizing experimental variables in tabular files is straightforward, the following sections of this vignette focus on the usage of targets files. This also integrates well with the widely used `SummarizedExperiment` object class. ### Structure of _`targets`_ file for single-end (SE) samples In a `targets` file with a single type of input files, here FASTQ files of single-end (SE) reads, the first three columns are mandatory including their column names, while it is four mandatory columns for FASTQ files of PE reads. All subsequent columns are optional and any number of additional columns can be added as needed. The columns in `targets` files are expected to be tab separated (TSV format). The `SampleName` column contains usually short labels for referencing samples (here FASTQ files) across many workflow steps (_e.g._ plots and column titles). Importantly, the labels used in the `SampleName` column need to be unique, while technical or biological replicates are indicated by the same values under the `Factor` column. For readability and transparency, it is useful to use here a short, consistent and informative syntax for naming samples and replicates. This is important since the values used under the `SampleName` and `Factor` columns are intended to be used as labels for naming columns or plotting features in downstream analysis steps. ```{r targetsSE, eval=TRUE} targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") showDF(read.delim(targetspath, comment.char = "#")) ``` To work with custom data, users need to generate a `targets` file containing the paths to their own FASTQ files and then provide under `targetspath` the path to the corresponding `targets` file. ### Structure of _`targets`_ file for paired-end (PE) samples For paired-end (PE) samples, the structure of the targets file is similar. The main difference is that `targets` files for PE data have two FASTQ path columns (here `FileName1` and `FileName2`) each containing the paths to the corresponding PE FASTQ files. ```{r targetsPE, eval=TRUE} targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR") showDF(read.delim(targetspath, comment.char = "#")) ``` ### Sample comparisons Sample comparisons of comparative experiments, such as differentially expressed genes (DEGs), can be specified in the header lines starting with the tag '``# ``' of a `targets` file. Their usage is optional, but useful for controlling comparative analyses according to certain biological expectations, such as identifying DEGs in RNA-Seq experiments based on simple pair-wise comparisons. ```{r comment_lines, echo=TRUE} readLines(targetspath)[1:4] ``` The function `readComp` imports the comparison information and stores it in a `list`. Alternatively, `readComp` can obtain the comparison information from a `SYSargsList` instance containing the `targets` file information (see below). ```{r targetscomp, eval=TRUE} readComp(file = targetspath, format = "vector", delim = "-") ``` ## Downstream targets files description After the step which required the initial targets file information, the downstream targets files are created automatically (see Figure \@ref(fig:targetsFig)). Each step that uses the previous step outfiles as an input, the new targets input will be managed internally by the workflow instances, establishing connectivity among the steps in the workflow. _`systemPipeR`_ provides features to automatically and systematically build this connection, providing security that all the samples will be managed efficiently and reproducibly. ```{r targetsFig, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "_`systemPipeR`_ automatically creates the downstream `targets` files based on the previous steps outfiles. A) Usually, users provide the initial `targets` files, and this step will generate some outfiles, as demonstrated on B. Then, those files are used to build the new `targets` files as inputs in the next step. _`systemPipeR`_ (C) manages this connectivity among the steps automatically for the users.", warning=FALSE} knitr::include_graphics("images/targets_con.png") ``` ## Structure of the new parameters files The parameters and configuration required for running command-line software are provided by the widely used community standard [Common Workflow Language](https://www.commonwl.org/) (CWL) [@Amstutz2016-ka], which describes parameters analysis workflows in a generic and reproducible manner. For R-based workflow steps, param files are not required. For a complete overview of the CWL syntax, please see the [section](#cwl) below. Also, we have a dedicated section explain how to _`systemPipeR`_ establish the connection between the CWL parameters files and the targets files. Please see [here](#cwl_targets). ```{r cleaning1, eval=TRUE, include=FALSE} if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE) ## NOTE: Removing previous project create in the quick starts section ``` # Project initialization To create a workflow within _`systemPipeR`_, we can start by defining an empty container and checking the directory structure: ```{r SPRproject1, eval=TRUE} sal <- SPRproject(projPath = getwd(), overwrite = TRUE) ``` Internally, `SPRproject` function will create a hidden folder called `.SPRproject`, by default, to store all the log files. A `YAML` file, here called `SYSargsList.yml`, has been created, which initially contains the basic location of the project structure; however, every time the workflow object `sal` is updated in R, the new information will also be store in this flat-file database for easy recovery. If you desire different names for the logs folder and the `YAML` file, these can be modified as follows: ```{r SPRproject_logs, eval=FALSE} sal <- SPRproject(logs.dir= ".SPRproject", sys.file=".SPRproject/SYSargsList.yml") ``` Also, this function will check and/or create the basic folder structure if missing, which means `data`, `param`, and `results` folder, as described [here](#directory-structure). If the user wants to use a different names for these directories, can be specified as follows: ```{r SPRproject_dir, eval=FALSE} sal <- SPRproject(data = "data", param = "param", results = "results") ``` It is possible to separate all the R objects created within the workflow analysis from the current environment. `SPRproject` function provides the option to create a new environment, and in this way, it is not overwriting any object you may want to have at your current section. ```{r SPRproject_env, eval=FALSE} sal <- SPRproject(envir = new.env()) ``` In this stage, the object `sal` is a empty container, except for the project information. The project information can be accessed by the `projectInfo` method: ```{r projectInfo, eval=TRUE} sal projectInfo(sal) ``` Also, the `length` function will return how many steps this workflow contains, and in this case, it is empty, as follow: ```{r length, eval=TRUE} length(sal) ``` # Workflow Design _`systemPipeR`_ workflows can be designed and built from start to finish with a single command, importing from an R Markdown file or stepwise in interactive mode from the R console. In the [next section](#appendstep), we will demonstrate how to build the workflow in an interactive mode, and in the [following section](#importWF), we will show how to build from a file. New workflows are constructed, or existing ones modified, by connecting each step via `appendStep` method. Each `SYSargsList` instance contains instructions needed for processing a set of input files with a specific command-line and the paths to the corresponding outfiles generated. The constructor function `Linewise` is used to build the R code-based step. For more details about this S4 class container, see [here](#linewise). ## Build workflow interactive {#appendstep} This tutorial shows a straightforward example for describing and explaining all main features available within systemPipeR to design, build, manage, run, and visualize the workflow. In summary, we are exporting a dataset to multiple files, compressing and decompressing each one of the files, importing to R, and finally performing a statistical analysis. In the previous section, we initialize the project by building the `sal` object. Until this moment, the container has no steps: ```{r sal_check, eval=TRUE} sal ``` Next, we need to populate the object created with the first step in the workflow. ### Adding the first step The first step is R code based, and we are splitting the `iris` dataset by `Species` and for each `Species` will be saved on file. Please note that this code will not be executed now; it is just store in the container for further execution. This constructor function requires the `step_name` and the R-based code under the `code` argument. The R code should be enclosed by braces (`{}`) and separated by a new line. ```{r, firstStep_R, eval=TRUE} appendStep(sal) <- LineWise(code = { mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")) ) }, step_name = "export_iris") ``` For a brief overview of the workflow, we can check the object as follows: ```{r show, eval=TRUE} sal ``` Also, for printing and double-check the R code in the step, we can use the `codeLine` method: ```{r codeLine, eval=TRUE} codeLine(sal) ``` ### Adding more steps Next, an example of how to compress the exported files using [`gzip`](https://www.gnu.org/software/gzip/) command-line. The constructor function creates an `SYSargsList` S4 class object using data from three input files: - CWL command-line specification file (`wf_file` argument); - Input variables (`input_file` argument); - Targets file (`targets` argument). In CWL, files with the extension `.cwl` define the parameters of a chosen command-line step or workflow, while files with the extension `.yml` define the input variables of command-line steps. The `targets` file is optional for workflow steps lacking `input` files. The connection between `input` variables and the `targets` file is defined under the `inputvars` argument. It is required a `named vector`, where each element name needs to match with column names in the `targets` file, and the value must match the names of the `input` variables defined in the `*.yml` files (see Figure \@ref(fig:sprandCWL)). A detailed description of the dynamic between `input` variables and `targets` files can be found [here](#cwl_targets). In addition, the CWL syntax overview can be found [here](#cwl). Besides all the data form `targets`, `wf_file`, `input_file` and `dir_path` arguments, `SYSargsList` constructor function options include: - `step_name`: a unique *name* for the step. This is not mandatory; however, it is highly recommended. If no name is provided, a default `step_x`, where `x` reflects the step index, will be added. - `dir`: this option allows creating an exclusive subdirectory for the step in the workflow. All the outfiles and log files for this particular step will be generated in the respective folders. - `dependency`: after the first step, all the additional steps appended to the workflow require the information of the dependency tree. The `appendStep<-` method is used to append a new step in the workflow. ```{r gzip_secondStep, eval=TRUE} targetspath <- system.file("extdata/cwl/gunzip", "targets_gunzip.txt", package = "systemPipeR") appendStep(sal) <- SYSargsList(step_name = "gzip", targets = targetspath, dir = TRUE, wf_file = "gunzip/workflow_gzip.cwl", input_file = "gunzip/gzip.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName = "_FILE_PATH_", SampleName = "_SampleName_"), dependency = "export_iris") ``` Note: This will not work if the `gzip` is not available on your system (installed and exported to PATH) and may only work on Windows systems using PowerShell. For a overview of the workflow, we can check the object as follows: ```{r} sal ``` Note that we have two steps, and it is expected three files from the second step. Also, the workflow status is *Pending*, which means the workflow object is rendered in R; however, we did not execute the workflow yet. In addition to this summary, it can be observed this step has three command lines. For more details about the command-line rendered for each target file, it can be checked as follows: ```{r} cmdlist(sal, step = "gzip") ``` #### Using the `outfiles` for the next step For building this step, all the previous procedures are being used to append the next step. However, here, we can observe power features that build the connectivity between steps in the workflow. In this example, we would like to use the outfiles from *gzip* Step, as input from the next step, which is the *gunzip*. In this case, let's look at the outfiles from the first step: ```{r} outfiles(sal) ``` The column we want to use is "gzip_file". For the argument `targets` in the `SYSargsList` function, it should provide the name of the correspondent step in the Workflow and which `outfiles` you would like to be incorporated in the next step. The argument `inputvars` allows the connectivity between `outfiles` and the new `targets` file. Here, the name of the previous `outfiles` should be provided it. Please note that all `outfiles` column names must be unique. It is possible to keep all the original columns from the `targets` files or remove some columns for a clean `targets` file. The argument `rm_targets_col` provides this flexibility, where it is possible to specify the names of the columns that should be removed. If no names are passing here, the new columns will be appended. ```{r gunzip, eval=TRUE} appendStep(sal) <- SYSargsList(step_name = "gunzip", targets = "gzip", dir = TRUE, wf_file = "gunzip/workflow_gunzip.cwl", input_file = "gunzip/gunzip.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(gzip_file = "_FILE_PATH_", SampleName = "_SampleName_"), rm_targets_col = "FileName", dependency = "gzip") ``` We can check the targets automatically create for this step, based on the previous `outfiles`: ```{r targetsWF_3, eval=TRUE} targetsWF(sal[3]) ``` We can also check all the expected `outfiles` for this particular step, as follows: ```{r outfiles_2, eval=TRUE} outfiles(sal[3]) ``` Now, we can observe that the third step has been added and contains one substep. ```{r} sal ``` In addition, we can access all the command lines for each one of the substeps. ```{r, eval=TRUE} cmdlist(sal["gzip"], targets = 1) ``` #### Getting data from a workflow instance The final step in this simple workflow is an R code step. For that, we are using the `LineWise` constructor function as demonstrated above. One interesting feature showed here is the `getColumn` method that allows extracting the information for a workflow instance. Those files can be used in an R code, as demonstrated below. ```{r getColumn, eval=TRUE} getColumn(sal, step = "gunzip", 'outfiles') ``` ```{r, iris_stats, eval=TRUE} appendStep(sal) <- LineWise(code = { df <- lapply(getColumn(sal, step = "gunzip", 'outfiles'), function(x) read.delim(x, sep = ",")[-1]) df <- do.call(rbind, df) stats <- data.frame(cbind(mean = apply(df[,1:4], 2, mean), sd = apply(df[,1:4], 2, sd))) stats$species <- rownames(stats) plot <- ggplot2::ggplot(stats, ggplot2::aes(x = species, y = mean, fill = species)) + ggplot2::geom_bar(stat = "identity", color = "black", position = ggplot2::position_dodge()) + ggplot2::geom_errorbar(ggplot2::aes(ymin = mean-sd, ymax = mean+sd), width = .2, position = ggplot2::position_dodge(.9)) }, step_name = "iris_stats", dependency = "gzip") ``` ## Build workflow from a {R Markdown} {#importWF} The precisely same workflow can be created by importing the steps from an R Markdown file. As demonstrated above, it is required to initialize the project with `SPRproject` function. `importWF` function will scan and import all the R chunk from the R Markdown file and build all the workflow instances. Then, each R chuck in the file will be converted in a workflow step. ```{r importWF_rmd, eval=TRUE} sal_rmd <- SPRproject(logs.dir = ".SPRproject_rmd") sal_rmd <- importWF(sal_rmd, file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR")) ``` Let's explore the workflow to check the steps: ```{r importWF_details} stepsWF(sal_rmd) dependency(sal_rmd) codeLine(sal_rmd) targetsWF(sal_rmd) ``` ### Rules to create the R Markdown to import as workflow To include a particular code chunk from the R Markdown file in the workflow analysis, please use the following code chunk options: - `spr=TRUE'`: for code chunks with step workflow. For example: > *```{r step_1, eval=TRUE, spr=TRUE}* > *```{r step_2, eval=FALSE, spr=TRUE}* `importWF` function can ignore `eval` option in code chunk, and in this case, both of the examples steps above would be incorporated in the workflow. For `spr = TRUE`, the last object assigned must to be the `SYSargsList`, for example: ```{r fromFile_example_rules_cmd, eval=FALSE} targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR") appendStep(sal) <- SYSargsList(step_name = "Example", targets = targetspath, wf_file = "example/example.cwl", input_file = "example/example.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(Message = "_STRING_", SampleName = "_SAMPLE_")) ``` OR ```{r fromFile_example_rules_r, eval=FALSE} appendStep(sal) <- LineWise(code = { library(systemPipeR) }, step_name = "load_lib") ``` Also, note that all the required files or objects to generate one particular step must be defined in an imported R code. The motivation for this is that when R Markdown files are imported, the `spr = TRUE` flag will be evaluated, append, and stored in the workflow control class as the `SYSargsList` object. The workflow object name used in the R Markdown (e.g. `appendStep(sal)`) needs to be the same when the `importWF` function is used. It is important to keep consistency. If different object names are used, when running the workflow, you can see a error, like *Error: object not found.*. # Running the workflow For running the workflow, `runWF` function will execute all the command lines store in the workflow container. ```{r runWF, eval=is_win} sal <- runWF(sal) ``` This essential function allows the user to choose one or multiple steps to be executed using the `steps` argument. However, it is necessary to follow the workflow dependency graph. If a selected step depends on a previous step(s) that was not executed, the execution will fail. ```{r runWF_error, eval=FALSE} sal <- runWF(sal, steps = c(1,3)) ``` Also, it allows forcing the execution of the steps, even if the status of the step is `'Success'` and all the expected `outfiles` exists. Another feature of the `runWF` function is ignoring all the warnings and errors and running the workflow by the arguments `warning.stop` and `error.stop`, respectively. ```{r runWF_force, eval=FALSE} sal <- runWF(sal, force = TRUE, warning.stop = FALSE, error.stop = TRUE) ``` When the project was initialized by `SPRproject` function, it was created an environment for all objects created during the workflow execution. This environment can be accessed as follows: ```{r runWF_env, eval=FALSE} viewEnvir(sal) ``` The workflow execution allows to save this environment for future recovery: ```{r runWF_saveenv, eval=FALSE} sal <- runWF(sal, saveEnv = TRUE) ``` ## Workflow status To check the summary of the workflow, we can use: ```{r show_statusWF, eval=TRUE} sal ``` To access more details about the workflow instances, we can use the `statusWF` method: ```{r statusWF, eval=TRUE} statusWF(sal) ``` ## Parallelization on clusters This section of the tutorial provides an introduction to the usage of the _`systemPipeR`_ features on a cluster. Alternatively, the computation can be greatly accelerated by processing many files in parallel using several compute nodes of a cluster, where a scheduling/queuing system is used for load balancing. The `resources` list object provides the number of independent parallel cluster processes defined under the `Njobs` element in the list. The following example will run 18 processes in parallel using each 4 CPU cores. If the resources available on a cluster allow running all 18 processes at the same time, then the shown sample submission will utilize in a total of 72 CPU cores. Note, `runWF` can be used with most queueing systems as it is based on utilities from the `batchtools` package, which supports the use of template files (_`*.tmpl`_) for defining the run parameters of different schedulers. To run the following code, one needs to have both a `conffile` (see _`.batchtools.conf.R`_ samples [here](https://mllg.github.io/batchtools/)) and a `template` file (see _`*.tmpl`_ samples [here](https://github.com/mllg/batchtools/tree/master/inst/templates)) for the queueing available on a system. The following example uses the sample `conffile` and `template` files for the Slurm scheduler provided by this package. The resources can be appended when the step is generated, or it is possible to add these resources later, as the following example using the `addResources` function: ```{r runWF_cluster, eval=FALSE} resources <- list(conffile=".batchtools.conf.R", template="batchtools.slurm.tmpl", Njobs=18, walltime=120, ## minutes ntasks=1, ncpus=4, memory=1024, ## Mb partition = "short" ) sal <- addResources(sal, c("hisat2_mapping"), resources = resources) sal <- runWF(sal) ``` Note: The example is submitting the jog to `short` partition. If you desire to use a different partition, please adjust accordingly. # Visualize workflow _`systemPipeR`_ workflows instances can be visualized with the `plotWF` function. This function will make a plot of selected workflow instance and the following information is displayed on the plot: - Workflow structure (dependency graphs between different steps); - Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`; - Sample status and statistics; - Workflow timing: running duration time. If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, _etc_. ```{r, eval=TRUE} plotWF(sal, show_legend = TRUE, width = "80%", rstudio = TRUE) ``` For more details about the `plotWF` function, please see [here](#plotWF). # Technical report _`systemPipeR`_ compiles all the workflow execution logs in one central location, making it easier to check any standard output (`stdout`) or standard error (`stderr`) for any command-line tools used on the workflow or the R code `stdout`. Also, the workflow plot is appended at the beginning of the report, making it easier to click on the respective step. ```{r, eval=FALSE} sal <- renderLogs(sal) ``` # Exported the workflow _`systemPipeR`_ workflow management system allows to translate and export the workflow build interactively to R Markdown format or an executable bash script. This feature advances the reusability of the workflow, as well as the flexibility for workflow execution. ## R Markdown file `sal2rmd` function takes an `SYSargsList` workflow container and translates it to SPR workflow template R markdown format. This file can be imported with the `importWF` function, as demonstrated above. ```{r, eval=FALSE} sal2rmd(sal) ``` ## Bash script `sal2bash` function takes an `SYSargsList` workflow container and translates it to an executable bash script, so one can run the workflow without loading `SPR` or using an R console. ```{r, eval=FALSE} sal2bash(sal) ``` It will be generated on the project root an executable bash script, called by default the `spr_wf.sh`. Also, a directory `./spr_wf` will be created and store all the R scripts based on the workflow steps. Please note that this function will "collapse" adjacent R steps into one file as much as possible. # Project Resume and Restart If you desire to resume or restart a project that has been initialized in the past, `SPRproject` function allows this operation. With the `resume` option, it is possible to load the `SYSargsList` object in R and resume the analysis. Please, make sure to provide the `logs.dir` location, and the corresponded `YAML` file name, if the default names were not used when the project was created. ```{r SPR_resume, eval=FALSE} sal <- SPRproject(resume = TRUE, logs.dir = ".SPRproject", sys.file = ".SPRproject/SYSargsList.yml") ``` If you choose to save the environment in the last analysis, you can recover all the files created in that particular section. `SPRproject` function allows this with `load.envir` argument. Please note that the environment was saved only with you run the workflow in the last section (`runWF()`). ```{r resume_load, eval=FALSE} sal <- SPRproject(resume = TRUE, load.envir = TRUE) ``` After loading the workflow at your current section, you can check the objects created in the old environment and decide if it is necessary to copy them to the current environment. ```{r envir, eval=FALSE} viewEnvir(sal) copyEnvir(sal, list="plot", new.env = globalenv()) ``` The `resume` option will keep all previous logs in the folder; however, if you desire to clean the execution (delete all the log files) history and restart the workflow, the `restart=TRUE` option can be used. ```{r restart_load, eval=FALSE} sal <- SPRproject(restart = TRUE, load.envir = FALSE) ``` The last and more drastic option from `SYSproject` function is to `overwrite` the logs and the `SYSargsList` object. This option will delete the hidden folder and the information on the `SYSargsList.yml` file. This will not delete any parameter file nor any results it was created in previous runs. Please use with caution. ```{r SPR_overwrite, eval=FALSE} sal <- SPRproject(overwrite = TRUE) ``` # Exploring workflow instances {#sysargslist} _`systemPipeR`_ provide several accessor methods and useful functions to explore `SYSargsList` workflow object. ## Accessor Methods Several accessor methods are available that are named after the slot names of the `SYSargsList` workflow object. ```{r} names(sal) ``` - Check the length of the workflow: ```{r} length(sal) ``` - Check the steps of the workflow: ```{r} stepsWF(sal) ``` - Checking the command-line for each target sample: `cmdlist()` method printing the system commands for running command-line software as specified by a given `*.cwl` file combined with the paths to the input samples (*e.g.* FASTQ files) provided by a `targets` file. The example below shows the `cmdlist()` output for running `gzip` and `gunzip` on the first sample. Evaluating the output of `cmdlist()` can be very helpful for designing and debugging `*.cwl` files of new command-line software or changing the parameter settings of existing ones. ```{r} cmdlist(sal, step = c(2,3), targets = 1) ``` - Check the workflow status: ```{r} statusWF(sal) ``` - Check the workflow targets files: ```{r} targetsWF(sal[2]) ``` - Checking the expected outfiles files: The `outfiles` components of `SYSargsList` define the expected outfiles files for each step in the workflow, some of which are the input for the next workflow step. ```{r} outfiles(sal[2]) ``` - Check the workflow dependencies: ```{r} dependency(sal) ``` - Check the sample comparisons: Sample comparisons are defined in the header lines of the `targets` file starting with '``# ``'. This information can be accessed as follows: ```{r, eval=FALSE} targetsheader(sal, step = "Quality") ``` - Get the workflow steps names: ```{r} stepName(sal) ``` - Get the Sample Id for on particular step: ```{r} SampleName(sal, step = "gzip") SampleName(sal, step = "iris_stats") ``` - Get the `outfiles` or `targets` column files: ```{r} getColumn(sal, "outfiles", step = "gzip", column = "gzip_file") getColumn(sal, "targetsWF", step = "gzip", column = "FileName") ``` - Get the R code for a `LineWise` step: ```{r} codeLine(sal, step = "export_iris") ``` - View all the objects in the running environment: ```{r} viewEnvir(sal) ``` - Copy one or multiple objects from the running environment to a new environment: ```{r} copyEnvir(sal, list = c("plot"), new.env = globalenv(), silent = FALSE) ``` - Accessing the `*.yml` data ```{r} yamlinput(sal, step = "gzip") ``` ## Subsetting the workflow details - The `SYSargsList` class and its subsetting operator `[`: ```{r} sal[1] sal[1:3] sal[c(1,3)] ``` - The `SYSargsList` class and its subsetting by steps and input samples: ```{r} sal_sub <- subset(sal, subset_steps = c( 2,3), input_targets = ("SE"), keep_steps = TRUE) stepsWF(sal_sub) targetsWF(sal_sub) outfiles(sal_sub) ``` - The `SYSargsList` class and its operator `+` ```{r, eval=FALSE} sal[1] + sal[2] + sal[3] ``` ## Replacement Methods - Update a `input` parameter in the workflow ```{r, eval=TRUE} sal_c <- sal ## check values yamlinput(sal_c, step = "gzip") ## check on command-line cmdlist(sal_c, step = "gzip", targets = 1) ## Replace yamlinput(sal_c, step = "gzip", paramName = "ext") <- "txt.gz" ## check NEW values yamlinput(sal_c, step = "gzip") ## Check on command-line cmdlist(sal_c, step = "gzip", targets = 1) ``` - Append and Replacement methods for R Code Steps ```{r, sal_lw_rep, eval=TRUE} appendCodeLine(sal_c, step = "export_iris", after = 1) <- "log_cal_100 <- log(100)" codeLine(sal_c, step = "export_iris") replaceCodeLine(sal_c, step="export_iris", line = 2) <- LineWise(code={ log_cal_100 <- log(50) }) codeLine(sal_c, step = 1) ``` For more details about the `LineWise` class, please see [below](#linewise). - Rename a Step ```{r} renameStep(sal_c, step = 1) <- "newStep" renameStep(sal_c, c(1, 2)) <- c("newStep2", "newIndex") sal_c names(outfiles(sal_c)) names(targetsWF(sal_c)) dependency(sal_c) ``` - Replace a Step ```{r, eval=FALSE} sal_test <- sal[c(1,2)] replaceStep(sal_test, step = 1, step_name = "gunzip" ) <- sal[3] sal_test ``` Note: Please use this method with attention, because it can disrupt all the dependency graphs. - Removing a Step ```{r} sal_test <- sal[-2] sal_test ``` # CWL syntax {#cwl} This section will introduce how CWL describes command-line tools and the specification and terminology of each file. For complete documentation, please check the CommandLineTools documentation [here](https://www.commonwl.org/v1.2/CommandLineTool.html) and [here](https://www.commonwl.org/v1.2/Workflow.html) for Workflows and the user guide [here](https://www.commonwl.org/user_guide/). CWL command-line specifications are written in [YAML](http://yaml.org/) format. In CWL, files with the extension `.cwl` define the parameters of a chosen command-line step or workflow, while files with the extension `.yml` define the input variables of command-line steps. ## CWL `CommandLineTool` `CommandLineTool` by CWL definition is a standalone process, with no interaction if other programs, execute a program, and produce output. Let's explore the `*.cwl` file: ```{r} dir_path <- system.file("extdata/cwl", package = "systemPipeR") cwl <- yaml::read_yaml(file.path(dir_path, "example/example.cwl")) ``` - The `cwlVersion` component shows the CWL specification version used by the document. - The `class` component shows this document describes a `CommandLineTool.` Note that CWL has another `class`, called `Workflow` which represents a union of one or more command-line tools together. ```{r} cwl[1:2] ``` - `baseCommand` component provides the name of the software that we desire to execute. ```{r} cwl[3] ``` - The `inputs` section provides the input information to run the tool. Important components of this section are: - `id`: each input has an id describing the input name; - `type`: describe the type of input value (string, int, long, float, double, File, Directory or Any); - `inputBinding`: It is optional. This component indicates if the input parameter should appear on the command-line. If this component is missing when describing an input parameter, it will not appear in the command-line but can be used to build the command-line. ```{r} cwl[4] ``` - The `outputs` section should provide a list of the expected outputs after running the command-line tools. Important components of this section are: - `id`: each input has an id describing the output name; - `type`: describe the type of output value (string, int, long, float, double, File, Directory, Any or `stdout`); - `outputBinding`: This component defines how to set the outputs values. The `glob` component will define the name of the output value. ```{r} cwl[5] ``` - `stdout`: component to specify a `filename` to capture standard output. Note here we are using a syntax that takes advantage of the inputs section, using results_path parameter and also the `SampleName` to construct the output `filename.` ```{r} cwl[6] ``` ## CWL `Workflow` `Workflow` class in CWL is defined by multiple process steps, where can have interdependencies between the steps, and the output for one step can be used as input in the further steps. ```{r} cwl.wf <- yaml::read_yaml(file.path(dir_path, "example/workflow_example.cwl")) ``` - The `cwlVersion` component shows the CWL specification version used by the document. - The `class` component shows this document describes a `Workflow`. ```{r} cwl.wf[1:2] ``` - The `inputs` section describes the inputs of the workflow. ```{r} cwl.wf[3] ``` - The `outputs` section describes the outputs of the workflow. ```{r} cwl.wf[4] ``` - The `steps` section describes the steps of the workflow. In this simple example, we demonstrate one step. ```{r} cwl.wf[5] ``` ## CWL Input Parameter Next, let's explore the *.yml* file, which provide the input parameter values for all the components we describe above. For this simple example, we have three parameters defined: ```{r} yaml::read_yaml(file.path(dir_path, "example/example_single.yml")) ``` Note that if we define an input component in the *.cwl* file, this value needs to be also defined here in the *.yml* file. # How to connect CWL description files within _`systemPipeR`_ {#cwl_targets} This section will demonstrate how to connect CWL parameters files to create workflows. In addition, we will show how the workflow can be easily scalable with _`systemPipeR`_. `SYSargsList` container stores all the information and instructions needed for processing a set of input files with a single or many command-line steps within a workflow (i.e. several components of the software or several independent software tools). The `SYSargsList` object is created and fully populated with the `SYSargsList` construct function. Full documentation of `SYSargsList` management instances can be found [here](#sysargslist) and [here](#appendstep). The following imports a `.cwl` file (here `example.cwl`) for running the `echo Hello World!` example. ```{r fromFile, eval=TRUE} HW <- SYSargsList(wf_file = "example/workflow_example.cwl", input_file = "example/example_single.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR")) HW cmdlist(HW) ``` However, we are limited to run just one command-line or one sample in this example. To scale the command-line over many samples, a simple solution offered by `systemPipeR` is to provide a `variable` for each of the parameters that we want to run with multiple samples. Let's explore the example: ```{r} yml <- yaml::read_yaml(file.path(dir_path, "example/example.yml")) yml ``` For the `message` and `SampleName` parameter, we are passing a variable connecting with a third file called `targets.` Now, let's explore the `targets` file structure: ```{r} targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR") read.delim(targetspath, comment.char = "#") ``` The `targets` file defines all input files or values and sample ids of an analysis workflow. For this example, we have defined a string message for the `echo` command-line tool, in the first column that will be evaluated, and the second column is the `SampleName` id for each one of the messages. Any number of additional columns can be added as needed. Users should note here, the usage of `targets` files is optional when using `systemPipeR's` new CWL interface. Since for organizing experimental variables targets files are extremely useful and user-friendly. Thus, we encourage users to keep using them. ### How to connect the parameter files and `targets` file information? The constructor function creates an `SYSargsList` S4 class object connecting three input files: - CWL command-line specification file (`wf_file` argument); - Input variables (`input_file` argument); - Targets file (`targets` argument). As demonstrated above, the latter is optional for workflow steps lacking input files. The connection between input variables (here defined by `input_file` argument) and the `targets` file are defined under the `inputvars` argument. A named vector is required, where each element name needs to match with column names in the `targets` file, and the value must match the names of the *.yml* variables. This is used to replace the CWL variable and construct all the command-line for that particular step. The variable pattern `_XXXX_` is used to distinguish CWL variables that target columns will replace. This pattern is recommended for consistency and easy identification but not enforced. The following imports a `.cwl` file (same example demonstrated above) for running the `echo Hello World` example. However, now we are connecting the variable defined on the `.yml` file with the `targets` file inputs. ```{r fromFile_example, eval=TRUE} HW_mul <- SYSargsList(step_name = "echo", targets=targetspath, wf_file="example/workflow_example.cwl", input_file="example/example.yml", dir_path = dir_path, inputvars = c(Message = "_STRING_", SampleName = "_SAMPLE_")) HW_mul cmdlist(HW_mul) ``` # Creating the CWL param files from the command-line Users need to define the command-line in a pseudo-bash script format: ```{r cmd, eval=TRUE} command <- " hisat2 \ -S \ -x \ -k \ -min-intronlen \ -max-intronlen \ -threads \ -U " ``` ## Define prefix and defaults - First line is the base command. Each line is an argument with its default value. - For argument lines (starting from the second line), any word before the first space with leading `-` or `--` in each will be treated as a prefix, like `-S` or `--min`. Any line without this first word will be treated as no prefix. - All defaults are placed inside `<...>`. - First argument is the input argument type. `F` for "File", "int", "string" are unchanged. - Optional: use the keyword `out` followed the type with a `,` comma separation to indicate if this argument is also an CWL output. - Then, use `:` to separate keywords and default values, any non-space value after the `:` will be treated as the default value. - If any argument has no default value, just a flag, like `--verbose`, there is no need to add any `<...>` ## `createParam` Function `createParam` function requires the `string` as defined above as an input. First of all, the function will print the three components of the `cwl` file: - `BaseCommand`: Specifies the program to execute. - `Inputs`: Defines the input parameters of the process. - `Outputs`: Defines the parameters representing the output of the process. The four component is the original command-line. If in interactive mode, the function will verify that everything is correct and will ask you to proceed. Here, the user can answer "no" and provide more information at the string level. Another question is to save the param created here. If running the workflow in non-interactive mode, the `createParam` function will consider "yes" and returning the container. ```{r} cmd <- createParam(command, writeParamFiles = FALSE) ``` If the user chooses not to save the `param` files on the above operation, it can use the `writeParamFiles` function. ```{r saving, eval=FALSE} writeParamFiles(cmd, overwrite = TRUE) ``` ## How to access and edit param files ### Print a component ```{r} printParam(cmd, position = "baseCommand") ## Print a baseCommand section printParam(cmd, position = "outputs") printParam(cmd, position = "inputs", index = 1:2) ## Print by index printParam(cmd, position = "inputs", index = -1:-2) ## Negative indexing printing to exclude certain indices in a position ``` ### Subsetting the command-line ```{r} cmd2 <- subsetParam(cmd, position = "inputs", index = 1:2, trim = TRUE) cmdlist(cmd2) cmd2 <- subsetParam(cmd, position = "inputs", index = c("S", "x"), trim = TRUE) cmdlist(cmd2) ``` ### Replacing a existing argument in the command-line ```{r} cmd3 <- replaceParam(cmd, "base", index = 1, replace = list(baseCommand = "bwa")) cmdlist(cmd3) ``` ```{r} new_inputs <- new_inputs <- list( "new_input1" = list(type = "File", preF="-b", yml ="myfile"), "new_input2" = "-L " ) cmd4 <- replaceParam(cmd, "inputs", index = 1:2, replace = new_inputs) cmdlist(cmd4) ``` ### Adding new arguments ```{r} newIn <- new_inputs <- list( "new_input1" = list(type = "File", preF="-b1", yml ="myfile1"), "new_input2" = list(type = "File", preF="-b2", yml ="myfile2"), "new_input3" = "-b3 " ) cmd5 <- appendParam(cmd, "inputs", index = 1:2, append = new_inputs) cmdlist(cmd5) cmd6 <- appendParam(cmd, "inputs", index = 1:2, after=0, append = new_inputs) cmdlist(cmd6) ``` ### Editing `output` param ```{r} new_outs <- list( "sam_out" = "" ) cmd7 <- replaceParam(cmd, "outputs", index = 1, replace = new_outs) output(cmd7) ``` ### Internal Check ```{r sysargs2, eval=TRUE} cmd <- " hisat2 \ -S \ -x \ -k \ -min-intronlen \ -max-intronlen \ -threads \ -U " WF <- createParam(cmd, overwrite = TRUE, writeParamFiles = TRUE, confirm = TRUE) targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") WF_test <- loadWorkflow(targets = targetspath, wf_file="hisat2.cwl", input_file="hisat2.yml", dir_path = "param/cwl/hisat2/") WF_test <- renderWF(WF_test, inputvars = c(FileName = "_FASTQ_PATH1_")) WF_test cmdlist(WF_test)[1:2] ``` # Inner Classes `SYSargsList` steps are can be defined with two inner classes, `SYSargs2` and `LineWise`. Next, more details on both classes. ## `SYSargs2` Class {#sysargs2} *`SYSargs2`* workflow control class, an S4 class, is a list-like container where each instance stores all the input/output paths and parameter components required for a particular data analysis step. *`SYSargs2`* instances are generated by two constructor functions, *loadWF* and *renderWF*, using as data input *targets* or *yaml* files as well as two *cwl* parameter files (for details see below). In CWL, files with the extension *`.cwl`* define the parameters of a chosen command-line step or workflow, while files with the extension *`.yml`* define the input variables of command-line steps. Note, input variables provided by a *targets* file can be passed on to a *`SYSargs2`* instance via the *inputvars* argument of the *renderWF* function. The following imports a *`.cwl`* file (here *`hisat2-mapping-se.cwl`*) for running the short read aligner HISAT2 [@Kim2015-ve]. For more details about the file structure and how to design or customize our own software tools, please check `systemPipeR and CWL` pipeline. ```{r sysargs2_cwl_structure, echo = FALSE, eval=FALSE} hisat2.cwl <- system.file("extdata", "cwl/hisat2/hisat2-mapping-se.cwl", package = "systemPipeR") yaml::read_yaml(hisat2.cwl) ``` ```{r sysargs2_yaml_structure, echo = FALSE, eval=FALSE} hisat2.yml <- system.file("extdata", "cwl/hisat2/hisat2-mapping-se.yml", package = "systemPipeR") yaml::read_yaml(hisat2.yml) ``` The *loadWF* and *renderWF* functions render the proper command-line strings for each sample and software tool. ```{r SYSargs2_structure, eval=TRUE} library(systemPipeR) targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") dir_path <- system.file("extdata/cwl", package = "systemPipeR") WF <- loadWF(targets = targetspath, wf_file = "hisat2/hisat2-mapping-se.cwl", input_file = "hisat2/hisat2-mapping-se.yml", dir_path = dir_path) WF <- renderWF(WF, inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_")) ``` Several accessor methods are available that are named after the slot names of the *`SYSargs2`* object. ```{r names_WF, eval=TRUE} names(WF) ``` Of particular interest is the *`cmdlist()`* method. It constructs the system commands for running command-line software as specified by a given *`.cwl`* file combined with the paths to the input samples (*e.g.* FASTQ files) provided by a *`targets`* file. The example below shows the *`cmdlist()`* output for running HISAT2 on the first SE read sample. Evaluating the output of *`cmdlist()`* can be very helpful for designing and debugging *`.cwl`* files of new command-line software or changing the parameter settings of existing ones. ```{r cmdlist, eval=TRUE} cmdlist(WF)[1] ``` The output components of *`SYSargs2`* define the expected output files for each step in the workflow; some of which are the input for the next workflow step, here next *`SYSargs2`* instance. ```{r output_WF, eval=TRUE} output(WF)[1] ``` The targets components of `SYSargs2` object can be accessed by the targets method. Here, for single-end (SE) samples, the structure of the targets file is defined by: - `FileName`: specify the FASTQ files path; - `SampleName`: Unique IDs for each sample; - `Factor`: ID for each treatment or condition. ```{r, targets_WF, eval=TRUE} targets(WF)[1] as(WF, "DataFrame") ``` Please note, to work with custom data, users need to generate a *`targets`* file containing the paths to their own FASTQ files and then provide under *`targetspath`* the path to the corresponding *`targets`* file. In addition, if the [Environment Modules](http://modules.sourceforge.net/) is available, it is possible to define which module should be loaded, as shown here: ```{r, module_WF, eval=TRUE} modules(WF) ``` Additional information can be accessed, as the parameters files location and the `inputvars` provided to generate the object. ```{r, other_WF, eval=FALSE} files(WF) inputvars(WF) ``` ## LineWise Class {#linewise} `LineWise` was designed to store all the R code chunk when an RMarkdown file is imported as a workflow. ```{r lw, eval=TRUE} rmd <- system.file("extdata", "spr_simple_lw.Rmd", package = "systemPipeR") sal_lw <- SPRproject(overwrite = TRUE) sal_lw <- importWF(sal_lw, rmd, verbose = FALSE) codeLine(sal_lw) ``` - Coerce methods available: ```{r, lw_coerce, eval=TRUE} lw <- stepsWF(sal_lw)[[2]] ## Coerce ll <- as(lw, "list") class(ll) lw <- as(ll, "LineWise") lw ``` - Access details ```{r, lw_access, eval=TRUE} length(lw) names(lw) codeLine(lw) codeChunkStart(lw) rmdPath(lw) ``` - Subsetting ```{r, lw_sub, eval=TRUE} l <- lw[2] codeLine(l) l_sub <- lw[-2] codeLine(l_sub) ``` - Replacement methods ```{r, lw_rep, eval=TRUE} replaceCodeLine(lw, line = 2) <- "5+5" codeLine(lw) appendCodeLine(lw, after = 0) <- "6+7" codeLine(lw) ``` - Replacement methods for `SYSargsList` ```{r, sal_rep_append, eval=FALSE} replaceCodeLine(sal_lw, step = 2, line = 2) <- LineWise(code={ "5+5" }) codeLine(sal_lw, step = 2) appendCodeLine(sal_lw, step = 2) <- "66+55" codeLine(sal_lw, step = 2) appendCodeLine(sal_lw, step = 1, after = 1) <- "66+55" codeLine(sal_lw, step = 1) ``` ## Workflow design structure using *`SYSargs`*: Previous version Instances of this S4 object class are constructed by the *`systemArgs`* function from two simple tabular files: a *`targets`* file and a *`param`* file. The latter is optional for workflow steps lacking command-line software. Typically, a *`SYSargs`* instance stores all sample-level inputs as well as the paths to the corresponding outputs generated by command-line- or R-based software generating sample-level output files, such as read preprocessors (trimmed/filtered FASTQ files), aligners (SAM/BAM files), variant callers (VCF/BCF files) or peak callers (BED/WIG files). Each sample level input/output operation uses its own *`SYSargs`* instance. The outpaths of *`SYSargs`* usually define the sample inputs for the next *`SYSargs`* instance. This connectivity is established by writing the outpaths with the *`writeTargetsout`* function to a new *`targets`* file that serves as input to the next *`systemArgs`* call. Typically, the user has to provide only the initial *`targets`* file. All downstream *`targets`* files are generated automatically. By chaining several *`SYSargs`* steps together one can construct complex workflows involving many sample-level input/output file operations with any combination of command-line or R-based software.
```{r, eval=TRUE, echo=FALSE, out.width="100%", fig.align = "center", fig.cap= "Workflow design structure of *`systemPipeR`* using previous version of *`SYSargs`*", warning=FALSE} # knitr::include_graphics(system.file("extdata/images", "SystemPipeR_Workflow.png", package = "systemPipeR")) ``` # Third-party software tools {#third-party-software-tools} Current, *systemPipeR* provides the _`param`_ file templates for third-party software tools. Please check the listed software tools. ```{r table_tools, echo=FALSE, message=FALSE} library(magrittr) SPR_software <- system.file("extdata", "SPR_software.csv", package = "systemPipeR") software <- read.delim(SPR_software, sep = ",", comment.char = "#") colors <- colorRampPalette((c("darkseagreen", "indianred1")))(length(unique(software$Category))) id <- as.numeric(c((unique(software$Category)))) software %>% dplyr::mutate(Step = kableExtra::cell_spec(Step, color = "white", bold = TRUE, background = factor(Category, id, colors) )) %>% dplyr::select(Tool, Description, Step) %>% dplyr::arrange(Tool) %>% kableExtra::kable(escape = FALSE, align = "c", col.names = c("Tool Name", "Description", "Step")) %>% kableExtra::kable_styling(c("striped", "hover", "condensed"), full_width = TRUE) %>% kableExtra::scroll_box(width = "80%", height = "500px") ``` Remember, if you desire to run any of these tools, make sure to have the respective software installed on your system and configure in the `PATH`. 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](https://systempipe.org/sp/spr/sp_run/step_run/#before-running). ```{r cleaning3, eval=TRUE, include=FALSE} if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE) ## NOTE: Removing previous project create in the quick starts section ``` # Workflow commom steps overview ## Project initialization To create a Workflow within _`systemPipeR`_, we can start by defining an empty container and checking the directory structure: ```{r SPRproject2, eval=FALSE} sal <- SPRproject() ``` ### Required packages and resources The `systemPipeR` package needs to be loaded [@H_Backman2016-bt]. ```{r load_SPR, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ library(systemPipeR) }, step_name = "load_SPR") ``` ## Read Preprocessing ### Preprocessing with `preprocessReads` function The function `preprocessReads` allows to apply predefined or custom read preprocessing functions to all FASTQ files referenced in a `SYSargsList` container, such as quality filtering or adapter trimming routines. Internally, `preprocessReads` uses the `FastqStreamer` function from the `ShortRead` package to stream through large FASTQ files in a memory-efficient manner. The following example performs adapter trimming with the `trimLRPatterns` function from the `Biostrings` package. Here, we are appending this step to the `SYSargsList` object created previously. All the parameters are defined on the `preprocessReads/preprocessReads-pe.yml` file. ```{r preprocessing, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList( step_name = "preprocessing", targets = "targetsPE.txt", dir = TRUE, wf_file = "preprocessReads/preprocessReads-pe.cwl", input_file = "preprocessReads/preprocessReads-pe.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c( FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_" ), dependency = c("load_SPR")) ``` After the preprocessing step, the `outfiles` files can be used to generate the new targets files containing the paths to the trimmed FASTQ files. The new targets information can be used for the next workflow step instance, _e.g._ running the NGS alignments with the trimmed FASTQ files. The `appendStep` function is automatically handling this connectivity between steps. Please check the next step for more details. The following example shows how one can design a custom read _'preprocessReads'_ function using utilities provided by the `ShortRead` package, and then run it in batch mode with the _'preprocessReads'_ function. Here, it is possible to replace the function used on the `preprocessing` step and modify the `sal` object. Because it is a custom function, it is necessary to save the part in the R object, and internally the `preprocessReads.doc.R` is loading the custom function. If the R object is saved with a different name (here `"param/customFCT.RData"`), please replace that accordingly in the `preprocessReads.doc.R`. Please, note that this step is not added to the workflow, here just for demonstration. First, we defined the custom function in the workflow: ```{r custom_preprocessing_function, eval=FALSE} appendStep(sal) <- LineWise( code = { filterFct <- function(fq, cutoff = 20, Nexceptions = 0) { qcount <- rowSums(as(quality(fq), "matrix") <= cutoff, na.rm = TRUE) # Retains reads where Phred scores are >= cutoff with N exceptions fq[qcount <= Nexceptions] } save(list = ls(), file = "param/customFCT.RData") }, step_name = "custom_preprocessing_function", dependency = "preprocessing" ) ``` After, we can edit the input parameter: ```{r editing_preprocessing, message=FALSE, eval=FALSE} yamlinput(sal, "preprocessing")$Fct yamlinput(sal, "preprocessing", "Fct") <- "'filterFct(fq, cutoff=20, Nexceptions=0)'" yamlinput(sal, "preprocessing")$Fct ## check the new function cmdlist(sal, "preprocessing", targets = 1) ## check if the command line was updated with success ``` ### Preprocessing with TrimGalore! [TrimGalore!](http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) is a wrapper tool to consistently apply quality and adapter trimming to fastq files, with some extra functionality for removing Reduced Representation Bisulfite-Seq (RRBS) libraries. ```{r trimGalore, eval=FALSE, spr=TRUE} targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") appendStep(sal) <- SYSargsList(step_name = "trimGalore", targets = targetspath, dir = TRUE, wf_file = "trim_galore/trim_galore-se.cwl", input_file = "trim_galore/trim_galore-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"), dependency = "load_SPR", run_step = "optional") ``` ### Preprocessing with Trimmomatic [Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) software [@Bolger2014-yr] performs a variety of useful trimming tasks for Illumina paired-end and single ended data. Here, an example of how to perform this task using parameters template files for trimming FASTQ files. ```{r trimmomatic, eval=FALSE, spr=TRUE} targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") appendStep(sal) <- SYSargsList(step_name = "trimmomatic", targets = targetspath, dir = TRUE, wf_file = "trimmomatic/trimmomatic-se.cwl", input_file = "trimmomatic/trimmomatic-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"), dependency = "load_SPR", run_step = "optional") ``` ## FASTQ quality report The following `seeFastq` and `seeFastqPlot` functions generate and plot a series of useful quality statistics for a set of FASTQ files, including per cycle quality box plots, base proportions, base-level quality trends, relative k-mer diversity, length, and occurrence distribution of reads, number of reads above quality cutoffs and mean quality distribution. The results are written to a PDF file named `fastqReport.pdf`. ```{r fastq_report, eval=FALSE, message=FALSE, spr=TRUE} appendStep(sal) <- LineWise(code = { fastq <- getColumn(sal, step = "preprocessing", "targetsWF", column = 1) fqlist <- seeFastq(fastq = fastq, batchsize = 10000, klength = 8) pdf("./results/fastqReport.pdf", height = 18, width = 4*length(fqlist)) seeFastqPlot(fqlist) dev.off() }, step_name = "fastq_report", dependency = "preprocessing") ```
**Figure 1:** FASTQ quality report

## NGS Alignment software After quality control, the sequence reads can be aligned to a reference genome or transcriptome database. The following sessions present some NGS sequence alignment software. Select the most accurate aligner and determining the optimal parameter for your custom data set project. For all the following examples, it is necessary to install the respective software and export the `PATH` accordingly. ### Alignment with `HISAT2` The following steps will demonstrate how to use the short read aligner `Hisat2` [@Kim2015-ve] in both interactive job submissions and batch submissions to queuing systems of clusters using the _`systemPipeR's`_ new CWL command-line interface. - Build `Hisat2` index. ```{r hisat_index, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "hisat_index", targets = NULL, dir = FALSE, wf_file = "hisat2/hisat2-index.cwl", input_file = "hisat2/hisat2-index.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = NULL, dependency = "preprocessing") ``` The parameter settings of the aligner are defined in the `workflow_hisat2-se.cwl` and `workflow_hisat2-se.yml` files. The following shows how to construct the corresponding *SYSargsList* object, and append to *sal* workflow. - Alignment with `HISAT2` and `SAMtools` It possible to build an workflow with `HISAT2` and `SAMtools`. ```{r hisat_mapping_samtools, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "hisat_mapping", targets = "preprocessing", dir = TRUE, wf_file = "workflow-hisat2/workflow_hisat2-se.cwl", input_file = "workflow-hisat2/workflow_hisat2-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("hisat_index"), run_session = "remote") ``` ### Alignment with `Tophat2` The NGS reads of this project can also be aligned against the reference genome sequence using `Bowtie2/TopHat2` [@Kim2013-vg; @Langmead2012-bs]. - Build _`Bowtie2`_ index. ```{r bowtie_index, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "bowtie_index", targets = NULL, dir = FALSE, wf_file = "bowtie2/bowtie2-index.cwl", input_file = "bowtie2/bowtie2-index.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = NULL, dependency = "preprocessing", run_step = "optional") ``` The parameter settings of the aligner are defined in the `workflow_tophat2-mapping.cwl` and `tophat2-mapping-pe.yml` files. The following shows how to construct the corresponding *SYSargsList* object, using the outfiles from the `preprocessing` step. ```{r tophat2_mapping, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "tophat2_mapping", targets = "preprocessing", dir = TRUE, wf_file = "tophat2/workflow_tophat2-mapping-se.cwl", input_file = "tophat2/tophat2-mapping-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("bowtie_index"), run_session = "remote", run_step = "optional") ``` ### Alignment with _`Bowtie2`_ (_e.g._ for miRNA profiling) The following example runs _`Bowtie2`_ as a single process without submitting it to a cluster. ```{r bowtie2_mapping, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "bowtie2_mapping", targets = "preprocessing", dir = TRUE, wf_file = "bowtie2/workflow_bowtie2-mapping-se.cwl", input_file = "bowtie2/bowtie2-mapping-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("bowtie_index"), run_session = "remote", run_step = "optional") ``` ### Alignment with _`BWA-MEM`_ (_e.g._ for VAR-Seq) The following example runs BWA-MEM as a single process without submitting it to a cluster. - Build the index: ```{r bwa_index, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "bwa_index", targets = NULL, dir = FALSE, wf_file = "bwa/bwa-index.cwl", input_file = "bwa/bwa-index.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = NULL, dependency = "preprocessing", run_step = "optional") ``` - Prepare the alignment step: ```{r bwa_mapping, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "bwa_mapping", targets = "preprocessing", dir = TRUE, wf_file = "bwa/bwa-se.cwl", input_file = "bwa/bwa-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("bwa_index"), run_session = "remote", run_step = "optional") ``` ### Alignment with _`Rsubread`_ (_e.g._ for RNA-Seq) The following example shows how one can use within the \Rpackage{systemPipeR} environment the R-based aligner \Rpackage{Rsubread}, allowing running from R or command-line. - Build the index: ```{r rsubread_index, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "rsubread_index", targets = NULL, dir = FALSE, wf_file = "rsubread/rsubread-index.cwl", input_file = "rsubread/rsubread-index.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = NULL, dependency = "preprocessing", run_step = "optional") ``` - Prepare the alignment step: ```{r rsubread_mapping, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "rsubread", targets = "preprocessing", dir = TRUE, wf_file = "rsubread/rsubread-mapping-se.cwl", input_file = "rsubread/rsubread-mapping-se.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("rsubread_index"), run_session = "remote", run_step = "optional") ``` ### Alignment with _`gsnap`_ (_e.g._ for VAR-Seq and RNA-Seq) Another R-based short read aligner is _`gsnap`_ from the _`gmapR`_ package [@Wu2010-iq]. The code sample below introduces how to run this aligner on multiple nodes of a compute cluster. - Build the index: ```{r gsnap_index, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "gsnap_index", targets = NULL, dir = FALSE, wf_file = "gsnap/gsnap-index.cwl", input_file = "gsnap/gsnap-index.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = NULL, dependency = "preprocessing", run_step = "optional") ``` - Prepare the alignment step: ```{r gsnap_mapping, eval=FALSE, spr=TRUE} appendStep(sal) <- SYSargsList(step_name = "gsnap", targets = "targetsPE.txt", dir = TRUE, wf_file = "gsnap/gsnap-mapping-pe.cwl", input_file = "gsnap/gsnap-mapping-pe.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"), dependency = c("gsnap_index"), run_session = "remote", run_step = "optional") ``` ## Create symbolic links for viewing BAM files in IGV The genome browser IGV supports reading of indexed/sorted BAM files via web URLs. This way it can be avoided to create unnecessary copies of these large files. To enable this approach, an HTML directory with Http access needs to be available in the user account (_e.g._ _`home/publichtml`_) of a system. If this is not the case then the BAM files need to be moved or copied to the system where IGV runs. In the following, _`htmldir`_ defines the path to the HTML directory with http access where the symbolic links to the BAM files will be stored. The corresponding URLs will be written to a text file specified under the `_urlfile`_ argument. ```{r bam_IGV, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise( code = { bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam") symLink2bam( sysargs = bampaths, htmldir = c("~/.html/", "somedir/"), urlbase = "http://cluster.hpcc.ucr.edu/~tgirke/", urlfile = "./results/IGVurl.txt") }, step_name = "bam_IGV", dependency = "hisat2_mapping", run_step = "optional" ) ``` ## Read counting for mRNA profiling experiments Create _`txdb`_ (needs to be done only once). ```{r create_txdb, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise(code = { library(txdbmaker) txdb <- makeTxDbFromGFF(file="data/tair10.gff", format="gff", dataSource="TAIR", organism="Arabidopsis thaliana") saveDb(txdb, file="./data/tair10.sqlite") }, step_name = "create_txdb", dependency = "hisat_mapping") ``` The following performs read counting with _`summarizeOverlaps`_ in parallel mode with multiple cores. ```{r read_counting, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ library(BiocParallel) txdb <- loadDb("./data/tair10.sqlite") eByg <- exonsBy(txdb, by="gene") outpaths <- getColumn(sal, step = "hisat_mapping", 'outfiles', column = 2) bfl <- BamFileList(outpaths, yieldSize=50000, index=character()) multicoreParam <- MulticoreParam(workers=4); register(multicoreParam); registered() counteByg <- bplapply(bfl, function(x) summarizeOverlaps(eByg, x, mode="Union", ignore.strand=TRUE, inter.feature=TRUE, singleEnd=TRUE)) # Note: for strand-specific RNA-Seq set 'ignore.strand=FALSE' and for PE data set 'singleEnd=FALSE' countDFeByg <- sapply(seq(along=counteByg), function(x) assays(counteByg[[x]])$counts) rownames(countDFeByg) <- names(rowRanges(counteByg[[1]])) colnames(countDFeByg) <- names(bfl) rpkmDFeByg <- apply(countDFeByg, 2, function(x) returnRPKM(counts=x, ranges=eByg)) write.table(countDFeByg, "results/countDFeByg.xls", col.names=NA, quote=FALSE, sep="\t") write.table(rpkmDFeByg, "results/rpkmDFeByg.xls", col.names=NA, quote=FALSE, sep="\t") }, step_name = "read_counting", dependency = "create_txdb") ``` Please note, in addition to read counts this step generates RPKM normalized expression values. For most statistical differential expression or abundance analysis methods, such as _`edgeR`_ or _`DESeq2`_, the raw count values should be used as input. The usage of RPKM values should be restricted to specialty applications required by some users, _e.g._ manually comparing the expression levels of different genes or features. #### Read and alignment count stats Generate a table of read and alignment counts for all samples. ```{r align_stats, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ read_statsDF <- alignStats(args) write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE, sep = "\t") }, step_name = "align_stats", dependency = "hisat_mapping") ``` The following shows the first four lines of the sample alignment stats file provided by the _`systemPipeR`_ package. For simplicity the number of PE reads is multiplied here by 2 to approximate proper alignment frequencies where each read in a pair is counted. ```{r align_stats2, eval=TRUE} read.table(system.file("extdata", "alignStats.xls", package = "systemPipeR"), header = TRUE)[1:4,] ``` ## Read counting for miRNA profiling experiments Download `miRNA` genes from `miRBase`. ```{r read_counting_mirna, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ system("wget https://www.mirbase.org/ftp/CURRENT/genomes/ath.gff3 -P ./data/") gff <- rtracklayer::import.gff("./data/ath.gff3") gff <- split(gff, elementMetadata(gff)$ID) bams <- getColumn(sal, step = "bowtie2_mapping", 'outfiles', column = 2) bfl <- BamFileList(bams, yieldSize=50000, index=character()) countDFmiR <- summarizeOverlaps(gff, bfl, mode="Union", ignore.strand = FALSE, inter.feature = FALSE) countDFmiR <- assays(countDFmiR)$counts # Note: inter.feature=FALSE important since pre and mature miRNA ranges overlap rpkmDFmiR <- apply(countDFmiR, 2, function(x) returnRPKM(counts = x, ranges = gff)) write.table(assays(countDFmiR)$counts, "results/countDFmiR.xls", col.names=NA, quote=FALSE, sep="\t") write.table(rpkmDFmiR, "results/rpkmDFmiR.xls", col.names=NA, quote=FALSE, sep="\t") }, step_name = "read_counting_mirna", dependency = "bowtie2_mapping") ``` ## Correlation analysis of samples The following computes the sample-wise Spearman correlation coefficients from the _`rlog`_ (regularized-logarithm) transformed expression values generated with the _`DESeq2`_ package. After transformation to a distance matrix, hierarchical clustering is performed with the _`hclust`_ function and the result is plotted as a dendrogram ([sample\_tree.pdf](sample_tree.png)). ```{r sample_tree_rlog, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ library(DESeq2, warn.conflicts=FALSE, quietly=TRUE) library(ape, warn.conflicts=FALSE) countDFpath <- system.file("extdata", "countDFeByg.xls", package="systemPipeR") countDF <- as.matrix(read.table(countDFpath)) colData <- data.frame(row.names = targetsWF(sal)[[2]]$SampleName, condition=targetsWF(sal)[[2]]$Factor) dds <- DESeqDataSetFromMatrix(countData = countDF, colData = colData, design = ~ condition) d <- cor(assay(rlog(dds)), method = "spearman") hc <- hclust(dist(1-d)) plot.phylo(as.phylo(hc), type = "p", edge.col = 4, edge.width = 3, show.node.label = TRUE, no.margin = TRUE) }, step_name = "sample_tree_rlog", dependency = "read_counting") ```
**Figure 2:** Correlation dendrogram of samples for _`rlog`_ values.

## DEG analysis with _`edgeR`_ The following _`run_edgeR`_ function is a convenience wrapper for identifying differentially expressed genes (DEGs) in batch mode with _`edgeR`_'s GML method [@Robinson2010-uk] for any number of pairwise sample comparisons specified under the _`cmp`_ argument. Users are strongly encouraged to consult the [_`edgeR`_](\href{http://www.bioconductor.org/packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf) vignette for more detailed information on this topic and how to properly run _`edgeR`_ on data sets with more complex experimental designs. ```{r edger, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") targets <- read.delim(targetspath, comment = "#") cmp <- readComp(file = targetspath, format = "matrix", delim = "-") countDFeBygpath <- system.file("extdata", "countDFeByg.xls", package = "systemPipeR") countDFeByg <- read.delim(countDFeBygpath, row.names = 1) edgeDF <- run_edgeR(countDF = countDFeByg, targets = targets, cmp = cmp[[1]], independent = FALSE, mdsplot = "") DEG_list <- filterDEGs(degDF = edgeDF, filter = c(Fold = 2, FDR = 10)) }, step_name = "edger", dependency = "read_counting") ``` Filter and plot DEG results for up and down-regulated genes. Because of the small size of the toy data set used by this vignette, the _FDR_ value has been set to a relatively high threshold (here 10%). More commonly used _FDR_ cutoffs are 1% or 5%. The definition of '_up_' and '_down_' is given in the corresponding help file. To open it, type _`?filterDEGs`_ in the R console.
**Figure 3:** Up and down regulated DEGs identified by _`edgeR`_.

## DEG analysis with _`DESeq2`_ The following _`run_DESeq2`_ function is a convenience wrapper for identifying DEGs in batch mode with _`DESeq2`_ [@Love2014-sh] for any number of pairwise sample comparisons specified under the _`cmp`_ argument. Users are strongly encouraged to consult the [_`DESeq2`_](http://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.pdf) vignette for more detailed information on this topic and how to properly run _`DESeq2`_ on data sets with more complex experimental designs. ```{r deseq2, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ degseqDF <- run_DESeq2(countDF=countDFeByg, targets=targets, cmp=cmp[[1]], independent=FALSE) DEG_list2 <- filterDEGs(degDF=degseqDF, filter=c(Fold=2, FDR=10)) }, step_name = "deseq2", dependency = "read_counting") ``` ## Venn Diagrams The function _`overLapper`_ can compute Venn intersects for large numbers of sample sets (up to 20 or more) and _`vennPlot`_ can plot 2-5 way Venn diagrams. A useful feature is the possibility to combine the counts from several Venn comparisons with the same number of sample sets in a single Venn diagram (here for 4 up and down DEG sets). ```{r vennplot, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ vennsetup <- overLapper(DEG_list$Up[6:9], type="vennsets") vennsetdown <- overLapper(DEG_list$Down[6:9], type="vennsets") vennPlot(list(vennsetup, vennsetdown), mymain="", mysub="", colmode=2, ccol=c("blue", "red")) }, step_name = "vennplot", dependency = "edger") ```
**Figure 4:** Venn Diagram for 4 Up and Down DEG Sets.

## GO term enrichment analysis of DEGs ### Obtain gene-to-GO mappings The following shows how to obtain gene-to-GO mappings from _`biomaRt`_ (here for _A. thaliana_) and how to organize them for the downstream GO term enrichment analysis. Alternatively, the gene-to-GO mappings can be obtained for many organisms from Bioconductor's _`*.db`_ genome annotation packages or GO annotation files provided by various genome databases. For each annotation, this relatively slow preprocessing step needs to be performed only once. Subsequently, the preprocessed data can be loaded with the _`load`_ function as shown in the next subsection. ```{r get_go_biomart, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ library("biomaRt") listMarts() # To choose BioMart database listMarts(host="plants.ensembl.org") m <- useMart("plants_mart", host="https://plants.ensembl.org") listDatasets(m) m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org") listAttributes(m) # Choose data types you want to download go <- getBM(attributes=c("go_id", "tair_locus", "namespace_1003"), mart=m) go <- go[go[,3]!="",]; go[,3] <- as.character(go[,3]) go[go[,3]=="molecular_function", 3] <- "F" go[go[,3]=="biological_process", 3] <- "P" go[go[,3]=="cellular_component", 3] <- "C" go[1:4,] dir.create("./data/GO") write.table(go, "data/GO/GOannotationsBiomart_mod.txt", quote=FALSE, row.names=FALSE, col.names=FALSE, sep="\t") catdb <- makeCATdb(myfile="data/GO/GOannotationsBiomart_mod.txt", lib=NULL, org="", colno=c(1,2,3), idconv=NULL) save(catdb, file="data/GO/catdb.RData") }, step_name = "get_go_biomart", dependency = "edger") ``` ### Batch GO term enrichment analysis Apply the enrichment analysis to the DEG sets obtained in the above differential expression analysis. Note, in the following example the _FDR_ filter is set here to an unreasonably high value, simply because of the small size of the toy data set used in this vignette. Batch enrichment analysis of many gene sets is performed with the _`GOCluster_Report`_ function. When _`method="all"`_, it returns all GO terms passing the p-value cutoff specified under the _`cutoff`_ arguments. When _`method="slim"`_, it returns only the GO terms specified under the _`myslimv`_ argument. The given example shows how one can obtain such a GO slim vector from BioMart for a specific organism. ```{r go_enrichment, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ load("data/GO/catdb.RData") DEG_list <- filterDEGs(degDF=edgeDF, filter=c(Fold=2, FDR=50), plot=FALSE) up_down <- DEG_list$UporDown; names(up_down) <- paste(names(up_down), "_up_down", sep="") up <- DEG_list$Up; names(up) <- paste(names(up), "_up", sep="") down <- DEG_list$Down; names(down) <- paste(names(down), "_down", sep="") DEGlist <- c(up_down, up, down) DEGlist <- DEGlist[sapply(DEGlist, length) > 0] BatchResult <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="all", id_type="gene", CLSZ=2, cutoff=0.9, gocats=c("MF", "BP", "CC"), recordSpecGO=NULL) library("biomaRt") m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org") goslimvec <- as.character(getBM(attributes=c("goslim_goa_accession"), mart=m)[,1]) BatchResultslim <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="slim", id_type="gene", myslimv=goslimvec, CLSZ=10, cutoff=0.01, gocats=c("MF", "BP", "CC"), recordSpecGO=NULL) gos <- BatchResultslim[grep("M6-V6_up_down", BatchResultslim$CLID), ] gos <- BatchResultslim pdf("GOslimbarplotMF.pdf", height=8, width=10); goBarplot(gos, gocat="MF"); dev.off() goBarplot(gos, gocat="BP") goBarplot(gos, gocat="CC") }, step_name = "go_enrichment", dependency = "get_go_biomart") ``` ### Plot batch GO term results The _`data.frame`_ generated by _`GOCluster_Report`_ can be plotted with the _`goBarplot`_ function. Because of the variable size of the sample sets, it may not always be desirable to show the results from different DEG sets in the same bar plot. Plotting single sample sets is achieved by subsetting the input data frame as shown in the first line of the following example.
**Figure 5:** GO Slim Barplot for MF Ontology.

## Clustering and heat maps The following example performs hierarchical clustering on the _`rlog`_ transformed expression matrix subsetted by the DEGs identified in the above differential expression analysis. It uses a Pearson correlation-based distance measure and complete linkage for cluster join. ```{r hierarchical_clustering, message=FALSE, eval=FALSE, spr=TRUE} appendStep(sal) <- LineWise({ library(pheatmap) geneids <- unique(as.character(unlist(DEG_list[[1]]))) y <- assay(rlog(dds))[geneids, ] pdf("heatmap1.pdf") pheatmap(y, scale="row", clustering_distance_rows="correlation", clustering_distance_cols="correlation") dev.off() }, step_name = "hierarchical_clustering", dependency = c("sample_tree_rlog", "edgeR")) ```
**Figure 7:** Heat map with hierarchical clustering dendrograms of DEGs.

## Visualize workflow systemPipeR workflows instances can be visualized with the `plotWF` function. This function will make a plot of selected workflow instance and the following information is displayed on the plot: - Workflow structure (dependency graphs between different steps); - Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`; - Sample status and statistics; - Workflow timing: running duration time. If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, etc. ```{r} plotWF(sal, show_legend = TRUE, width = "80%") ``` To check more details of `plotWF`, visit [our website](https://systempipe.org/sp/spr/sp_run/step_vis/). ## Running workflow ### Interactive job submissions in a single machine For running the workflow, `runWF` function will execute all the steps store in the workflow container. The execution will be on a single machine without submitting to a queuing system of a computer cluster. ```{r runWF_steps, eval=FALSE} sal <- runWF(sal) ``` ### Parallelization on clusters Alternatively, the computation can be greatly accelerated by processing many files in parallel using several compute nodes of a cluster, where a scheduling/queuing system is used for load balancing. The `resources` list object provides the number of independent parallel cluster processes defined under the `Njobs` element in the list. The following example will run 18 processes in parallel using each 4 CPU cores. If the resources available on a cluster allow running all 18 processes at the same time, then the shown sample submission will utilize in a total of 72 CPU cores. Note, `runWF` can be used with most queueing systems as it is based on utilities from the `batchtools` package, which supports the use of template files (_`*.tmpl`_) for defining the run parameters of different schedulers. To run the following code, one needs to have both a `conffile` (see _`.batchtools.conf.R`_ samples [here](https://mllg.github.io/batchtools/)) and a `template` file (see _`*.tmpl`_ samples [here](https://github.com/mllg/batchtools/tree/master/inst/templates)) for the queueing available on a system. The following example uses the sample `conffile` and `template` files for the Slurm scheduler provided by this package. The resources can be appended when the step is generated, or it is possible to add these resources later, as the following example using the `addResources` function: ```{r runWF_cluster_steps, eval=FALSE} resources <- list(conffile=".batchtools.conf.R", template="batchtools.slurm.tmpl", Njobs=18, walltime=120, ## minutes ntasks=1, ncpus=4, memory=1024, ## Mb partition = "short" ) sal <- addResources(sal, c("hisat2_mapping"), resources = resources) sal <- runWF(sal) ``` ### Checking workflow status To check the summary of the workflow, we can use: ```{r statusWF_steps, eval=FALSE} sal statusWF(sal) ``` ### Accessing logs report _`systemPipeR`_ compiles all the workflow execution logs in one central location, making it easier to check any standard output (`stdout`) or standard error (`stderr`) for any command-line tools used on the workflow or the R code stdout. ```{r logsWF_steps, eval=FALSE} sal <- renderLogs(sal) ``` ```{r quiting, eval=TRUE, echo=FALSE} knitr::opts_knit$set(root.dir = '../') unlink("rnaseq", recursive = TRUE) ``` # Version information ```{r sessionInfo} sessionInfo() ``` # Funding This project is funded by NSF award [ABI-1661152](https://www.nsf.gov/awardsearch/showAward?AWD_ID=1661152). # References