```{css, echo=FALSE} pre code { white-space: pre !important; overflow-x: scroll !important; word-break: keep-all !important; word-wrap: initial !important; }
```r 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)
if (file.exists("spr_project")) unlink("spr_project", recursive = TRUE) is_win <- Sys.info()[['sysname']] != "Windows"
suppressPackageStartupMessages({ library(systemPipeR) })
systemPipeR
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.
knitr::include_graphics("images/utilities.png")
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 each analysis step in 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)). Typically, SAL
instances are constructed
from an intial metadata targets table, R code and CWL parameter files for each
R- and CL-based analysis step in workflows (details provided below).
For preconfigured workflows, users only need to provide their input data (such as FASTQ
files) and the corresponding metadata in a targets file. The latter describes the
experimental design, defines sample labels, replicate information, and other
relevant information.
knitr::include_graphics("images/SYSargsList.png")
systemPipeR
adopts the Common Workflow Language
(CWL), which is a widely used community
standard for describing CL tools and workflows in a declarative, generic, and
reproducible manner [@Amstutz2016-ka]. CWL specifications are human-readable
YAML 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, the basic description for executing a CL
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
parameters of CL tool or workflow (C), and the yml file (B) assigns the input
variables to the corresponding parameters. For convenience, in systemPipeR
parameter
values can be provided by a 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. The structure of
targets
files is explained in the corresponding section
below. A detailed overview of the CWL syntax is provided in
the CWL syntax section below, and the details for connecting the input
information in targets
with CWL parameters are described
here.
knitr::include_graphics("images/SPR_CWL_hello.png")
systemPipeRdata
, a companion package to systemPipeR
, offers a collection of
workflow templates that are ready to use. With a single command, users can
easily load these templates onto their systems. Once loaded, users have the
flexibility to utilize the templates as they are or modify them as needed. More
in-depth information can be found in the main vignette of systemPipeRdata,
which can be accessed
here.
The package also provides several convenience
functions that are useful for designing and testing workflows, such as a
CL rendering function that assembles from the parameter files (cwl, yml and
targets) the exact CL strings for each step prior to running a CL tool.
Auto-generation of CWL parameter files is also supported. Here, users can simply
provide the CL strings for a CL software of interest to a rendering function that generates
the corresponding *.cwl
and *.yml
files for them. And auto-conversion of workflows to
executable Bash scripts.
The systemPipeR
package can be installed from the R console using the BiocManager::install
command. The associated systemPipeRdata
package can be installed the same way. The latter is a data 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.
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 a user's system, where the analysis will be performed (details provided below).
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 a small test data set, 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
). A list of available workflow templates
is available in the vignette of the systemPipeRdata
package here.
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 from a workflow Rmd file (here systemPipeRNAseq.Rmd
)
defining the analysis steps.
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 ## The following print statements, issued during the import, are shortened for brevity ## Import messages for first 3 of 20 steps total ## Parse chunk code ## Now importing step 'load_SPR' ## Now importing step 'preprocessing' ## Now importing step 'trimming' ## Now importing step '...' ## ... ## Now check if required CL tools are installed ## Messages for 4 of 7 CL tools total ## step_name tool in_path ## 1 trimming trimmomatic TRUE ## 2 hisat2_index hisat2-build TRUE ## 3 hisat2_mapping hisat2 TRUE ## 4 hisat2_mapping samtools TRUE ## ...
The importWF
function also checks the availability of the R packages and CL
software tools used by a workflow. 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 listed. If the in_path
column shows FALSE
for
any of them, then the missing CL software needs to be installed and made available in a user's
PATH
prior to running the workflow. Note, the shown availability table of CL tools can
also be returned with listCmdTools(sal, check_path=TRUE)
, and the availability of individual CL
tools can be checked with tryCL
, e.g. for hisat2
use: tryCL(command = "hisat2")
.
3. Status summary. An overview of the workflow steps and their status
information can be returned by typing sal
. For space reasons, the following
shows only the first 3 of a total of 20 steps of the RNA-Seq workflow. At this
stage all workflow steps are listed as pending since none of them have been executed yet.
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. - 20. not shown here for brevity
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.
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
.
sal
knitr::include_graphics("images/runwf.png")
5. Workflow topology graph. Workflows can be displayed as topology graphs
using the plotWF
function. The run status information for each step and various
other details are embedded in these graphs. Additional details are provided in the visualize workflow
section below.
plotWF(sal)
knitr::include_graphics("images/plotWF.png")
6. Report generation. The renderReport
and renderLogs
function can be used
for generating scientific and technical reports, respectively. Alternatively, scientific
reports can be generated with the render
function of the rmarkdown
package.
# Scietific report sal <- renderReport(sal) rmarkdown::render("systemPipeRNAseq.Rmd", clean = TRUE, output_format = "BiocStyle::html_document") # Technical (log) report sal <- renderLogs(sal)
The root directory of systemPipeR
workflows contains by default three user
facing sub-directories: data
, results
and param
. A fourth sub-directory is
a hidden 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 to the default structure is possible without requiring changes
to the workflows. The following directory tree summarizes the expected content in
each default directory (names given in green).
.batchtools.conf.R
and tmpl
files for batchtools
and BiocParallel
. cwl
and yml
files. Previous versions of parameter files are stored in a separate sub-directory. SPRproject
function at the beginning of a workflow run. It is a hidden directory because its name starts with a dot. targets
file {#targets-files}A targets
file defines the input files (e.g. FASTQ, BAM, BCF) and
sample comparisons used in a data analysis workflow. It can also store any number of
additional descriptive information for each sample. How the input
information is passed on from a targets
file to the CWL parameter files is
introduced above, and additional details are given below. 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.
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. Their usage also
integrates well with the widely used SummarizedExperiment
object class.
Descendant targets files can be extracted for each step with input/output operations where the output of the previous step(s) serves as input to the current step, and the output of the current step becomes the input of the next step. This connectivity among input/output operations is automatically tracked throughout workflows. This way it is straightforward to start workflows at different processing stages. For instance, one can intialize an RNA-Seq workflow at the stage of raw sequence files (FASTQ), alignment files (BAM) or a precomputed read count table.
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 provided under the SampleName
and Factor
columns are intended to
be used as labels for naming the columns or plotting features in downstream
analysis steps.
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.
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.
targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR") showDF(read.delim(targetspath, comment.char = "#"))
If needed, sample comparisons of comparative experiments, such as differentially expressed genes (DEGs), can be
specified in the header lines of a targets
file that start with a # <CMP>
tag.
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.
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).
readComp(file = targetspath, format = "vector", delim = "-")
if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE) ## NOTE: Removes previous project created in the quick-start section
A systemPipeR
workflow instance is initialized with the SPRproject
function. This function
call creates an empty SAL
container instance and at the same time a linked project
log directory that acts as a flat-file database of a workflow. A YAML file is automatically
included in the project directory that specifies the basic location of the workflow project.
Every time the SAL
container is updated in R with a new workflow step or a modification
to an existing step, the changes are automatically recorded in the flat-file database. This
is important for tracking the run status of workflows and providing restart functionality for
workflows.
sal <- SPRproject()
If overwrite
is set to TRUE
, a new project log directory will be created and any existing
one deleted. This option should be used with caution. It is mainly useful when developing
and testing workflows, but should be avoided in production runs of workflows.
sal <- SPRproject(projPath = getwd(), overwrite = TRUE)
The function checks whether the expected workflow directories (see here) exist, and will create them if any of them is missing. If needed users can change the default names of these directories as shown.
sal <- SPRproject(data = "data", param = "param", results = "results")
Similarly, the default names of the log directory and YAML
file can be changed.
sal <- SPRproject(logs.dir= ".SPRproject", sys.file=".SPRproject/SYSargsList.yml")
It is also possible to use for all workflow steps a dedicated R environment that is separate from the current environment. This way R objects generated by workflow steps will not overwrite objects with the same names in the current environment.
sal <- SPRproject(envir = new.env())
At this stage, sal
is an empty SAL
(SYSargsList
) container that only contains
the basic information about the project's directory structure that can be accessed with
projectInfo
.
sal
projectInfo(sal)
The number of workflow steps stored in a SAL
object can be returned with the length
function. At this stage
it returns zero since no workflow steps have been loaded into sal
yet.
length(sal)
Workflows in systemPipeR
can be constructed stepwise in interactive
mode by evaluating the code of individual workflow steps in the R
console one-by-one. Alternatively, one can import all steps of a workflow with
a single import command at once, either from an R script or an R
Markdown workflow file. For the purpose of explaining the details about
constructing and connecting different types of workflow steps, this tutorial
section introduces first the interactive approach. After this the automated
import of entire workflows with many steps is explained where the individual
steps are defined the same way. In all cases, workflow steps are loaded to a
SAL
workflow container with the proper connectivity information using
systemPipeR's
appendStep
method where steps can be comprised of R code or
CL calls.
The following demonstrates how to design, load and run workflows using a simple data processing routine as an example. This mini workflow will export a test dataset to multiple files, compress/decompress the exported files, import them back into R, and then perform a simple statistical analysis and plot the results.
The sal
object of the new workflow project (directory named.SPRproject
) was
intialized in the previous section. At this point this sal
instance contains
no data analysis steps since none have been loaded so far.
sal
Next, workflow steps will be added to sal
.
The first step in the chosen example workflow comprises R code that will be
stored in a LineWise
object. It is constructed with the LineWise
function,
and then appended to sal
with the appendStep<-
method. The R code of an
analysis step is assigned to the code
argument of the LineWise
function. In this
assignment the R code has to be enclosed by braces ({...}
) and separted from
them by new lines. Additionally, the workflow step should be given a descriptive name
under the step_name
argument. Step names are required to be unique throughout
workflows. During the construction of workflow steps, the included R code will
not be executed. The execution of workflow steps is explained in a separate
section below.
In the given code example, the iris
dataset is split by the species
names under the Species
column, and then the resulting data.frames
are
exported to three tabular files.
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")
After adding the R code, sal
contains now one workflow step.
sal
To extract the code of an R step stored in a SAL
object, the codeLine
method can be used.
codeLine(sal)
CL steps are stored as SYSargs2
objects that are constructed with the
SYSargsList
function, and then appended to sal
with the appendStep<-
method. As outlined in the introduction (see here), CL steps
are defined by two CWL parameter files (yml
configuration and cwl
step
definition files) and an optional targets
file. How parameter values in the
targets
file are passed on to the corresponding entries in the yml
file, is
defined by a named vector
that is assigned to the inputvars
argument of the
SYSargsList
function. A parameter connection is established if a name assigned to
inputvars
has matching column and element names in the targets
and yml
files,
respectively (Fig \@ref(fig:sprandCWL)). More details about parameter passing and CWL
syntax are provied below (see here and here).
The most important other arguments of the SYSargsList
function are listed below. For more
information, users want to consult the function's help with ?SYSargsList
.
step_name
: a unique name for the step. If no name is provided, a default
step_x
name will be assigned, where x
is the step index.dir
: if TRUE
(default) all output files generated by a workflow step will be written to a
subdirectory with the same name as step_name
. This is useful for organizing result files. dependency
: assign here the name of the step the current step depends on. This is mandatory
for all steps in a workflow, except the first one. The dependency tree of a workflow is
based on the dependency connections among steps.In the specific example code given below, a CL step is added to the workflow where the
gzip
software is used to compress the
files that were generated in the previous step.
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")
After adding the above CL step, sal
contains now two steps.
sal
The individual CL calls, that will be executed by the gzip
step, can be rendered and viewed
with the cmdlist
function. Under the targets
argument one can subset the CL calls to
specific samples by assigning the corresponding names or index numbers.
cmdlist(sal, step = "gzip") # cmdlist(sal, step = "gzip", targets=c("SE"))
In many use cases the output files, generated by an upstream workflow step, serve as input
to a downstream step. To establish these input/output connections, the names (paths) of the
output files generated by each step needs to be accessible. This information
can be extracted from SAL
objects with the outfiles
accessor method as shown below.
# outfiles(sal) # output files of all steps in sal outfiles(sal)['gzip'] # output files of 'gzip' step # colnames(outfiles(sal)$gzip) # returns column name passed on to `inputvars`
Note, the names of this and other important accessor methods for 'SAL' objects
can be looked up conveniently with names(sal)
(for more details see here).
In the chosen workflow example, the output files (here compressed gz
files), that
were generated by the previous gzip
step, will be uncompressed in the current step with the
gunzip
software. The corresponding input files for the gunzip
step are listed under the
gzip_file
column above. For defining the gunzip
step, the values 'gzip' and 'gzip_file'
will be used under the targets
and inputvars
arguments of the SYSargsList
function,
respectively. The argument rm_targets_col
allows to drop columns in the targets
instance of the new step. The remaining parameters settings are similar to those in the
previous step.
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")
After adding the above new step, sal
contains now a third step.
sal
The targets
instance of the new step can be returned with the targetsWF
method
where the output files from the previous step are listed under the first column (input).
targetsWF(sal['gunzip'])
As before, the output files of the new step can be returned with outfiles
.
outfiles(sal['gunzip'])
Finally, the corresponding CL calls of the new step can be returned with the cmdlist
function (here for first entry).
cmdlist(sal["gunzip"], targets = 1)
The final step in this sample workflow is an R step that uses the files from a previous
step as input. In this case the getColumn
method is used to obtain the paths to the files
generated in a previous step, which is in the given example the 'gunzip' step..
getColumn(sal, step = "gunzip", 'outfiles')
In this R step, the tabular files generated in the previous gunzip
CL step
are imported into R and row appended to a single data.frame
. Next the
column-wise mean values are calculated for the first four columns.
Subsequently, the results are plotted as a bar diagram with error bars.
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$size <- rownames(stats) plot <- ggplot2::ggplot(stats, ggplot2::aes(x = size, y = mean, fill = size)) + 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")
This is the final step of this demonstration resulting in a sal
workflow container with
a total of four steps.
sal
The above process of loading workflow steps one-by-one into a SAL
workflow
container can be easily automated by storing the step definitions in an R or
Rmd script, and then importing them from there into an R session.
1. Loading workflows from an R script. For importing workflow steps from an
R script, the code of the workflow steps needs to be stored in an R script
from where it can be imported with R's source
command. Applied to
the above workflow example (see here), this means nothing else
than saving the code of the four workflow steps to an R script where each step is declared
with the standard CL or R step syntax: appendStep(sal) <- SYSargsList/LineWise(...)
.
At the beginning of the R script one has to load the systemPipeR
library, and
initialize a new workflow project and associated SAL
container with SPRproject()
.
After sourcing the R script from R, the fully populated SAL
container will be
loaded into that session, and the workflow is ready to be executed (see below).
2. Loading workflows from an R Markdown file. As an alternative to plain R
scripts, R Markdown (Rmd) scripts provide a more adaptable solution for
defining workflows. An Rmd file can be converted into various publication-ready
formats, such as HTML or PDF. These formats can incorporate not only the
analysis code but also the results the code generates, including tables and figures.
This approach enables the creation of reproducible analysis reports for
workflows. This reporting feature is crucial for reproducibility,
documentation, and visual interpretation of the analysis results. The following illustrates this
approach for the same four workflow steps used in the previous section here,
that is included in an Rmd file of the systemPipeR
package. Note, the path to this Rmd file
is retrieved with R's system.file
function.
Prior to importing the workflow from an Rmd file, it is required to initialize for it a new
workflow project with the SPRproject
function. Next, the importWF
function is used to scan
the Rmd file for code chunks that define workflow steps, and subsequently import them in to the
SAL
workflow container of the project.
sal_rmd <- SPRproject(logs.dir = ".SPRproject_rmd") sal_rmd <- importWF(sal_rmd, file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"))
After the import, the new sal_rmd
workflow container, that is fully populated with all four workflow
steps from before, can be inspected with several accessor functions (not
evaluated here). Additional details about accessor functions are provided here.
sal_rmd stepsWF(sal_rmd) dependency(sal_rmd) cmdlist(sal_rmd) codeLine(sal_rmd) targetsWF(sal_rmd) statusWF(sal_rmd)
In standard R Markdown (Rmd) files, code chunks are enclosed by new lines
starting with three backticks. The backtick line at the start of a code chunk
is followed by braces that can contain arguments controlling the code chunk's
behavior. To formally declare a workflow step in an R Markdown file's argument
line, systemPipeR
introduces a special argument named spr
. When
using importWF
to scan an R Markdown file, only code chunks with spr=TRUE
in
their argument line will be recognized as workflow steps and loaded into the
provided SAL
workflow container. This design allows for the inclusion of
standard code chunks not part of a workflow and renders them as usual. Here are
two examples of argument settings that will both result in the inclusion of the
corresponding code chunk as a workflow step since spr
is set to TRUE
in both
cases. Notably, in one case, the standard R Markdown argument eval
is assigned
FALSE
, preventing the rmarkdown::render
function from evaluating the
corresponding code chunk.
Examples: workflow code chunks are declared by spr
flag in their argument line:
In addition to including spr = TRUE
, the actual code of workflow steps has additional
requirements. First, the last assignment in a code chunk of a workflow step needs to be an
appendStep
of SAL
using SYSargsList
or LineWise
for CL or R code, respectively. This
requirement is met if there are no other assignments outside of appnedStep
. Second,
R workflow steps need to be largely self contained by generating and/or loading the dependencies
required to execute the code. Third, in most cases the name of a SAL
container should remain
the same throughout a workflow. This avoids errors such as: 'Error:
Example of last assignment in a CL step.
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_"))
Example of last assignment in an R step.
appendStep(sal) <- LineWise(code = { library(systemPipeR) }, step_name = "load_lib")
In systemPipeR
, the runWF
function serves as the primary tool for executing
workflows. It is responsible for running the code specified in the steps of a
populated SAL
workflow container. The following runWF
command will run the
test workflow from above from start to finish. This test workflow was first assembled step-by-step,
allowing for a thorough examination of its behavior. Subsequently, the same workflow
was imported from an Rmd file to demonstrate how to auto-load all steps of a workflow
at once into a SAL
container. Please refer to the provided link here
for more information about this process.
sal <- runWF(sal)
The runWF
function allows to choose one or multiple steps to be
executed via its steps
argument. When using partial workflow executions, it is important
to pay attention to the requirements of the dependency graph of a workflow. If a selected step
depends on one or more previous steps, that have not been executed and completed yet,
then the execution of the chosen step(s) will not be possible.
sal <- runWF(sal, steps = c(1,3))
Importantly, by default, already completed workflow steps with a status of 'Success
' (for
example, all output files exist) will not be repeated unnecessarily unless one explicitly sets
the force parameter to TRUE. Skipping such steps can save time, particularly
when optimizing workflows or adding new samples to previously completed runs.
Additionally, one may find it useful in certain situations to ignore warnings or
errors without terminating workflow runs. This behavior can be enabled by setting
warning.stop=TRUE
and/or error.stop=TRUE
.
sal <- runWF(sal, force = TRUE, warning.stop = FALSE, error.stop = TRUE)
When starting a new workflow project with the SPRproject
function, a new R environment
will be initialized that stores the objects generated by the workflow steps. The content
of this R environment can be inspected with the viewEnvir
function.
viewEnvir(sal)
The runWF
function saves the new R environment to an rds
file under .SPRproject
when saveEnv=TRUE
, which
is done by default. For additional details, users want to consult the corresponding help document
with ?runWF
.
sal <- runWF(sal, saveEnv = TRUE)
A status summary of the executed workflows can be returned by typing sal
.
sal
Several accessor functions can be used to retrieve additional information about workflows and their run status. The code box below lists these functions, omitting their output for brevity. Although some of these functions have been introduced above already, they are included here again for easy reference. Additional, details on these functions can be found here.
stepsWF(sal) dependency(sal) cmdlist(sal) codeLine(sal) targetsWF(sal) statusWF(sal) projectInfo(sal)
While SAL
objects are autosaved when working with workflows, it
can be sometimes safer to explicity save the object before closing R.
sal <- write_SYSargsList(sal)
Some computing systems, such as HPC clusters, allow users to load software via
an Environment Modules system into their
PATH
. If a module system is available, the function module
allows to
interact with it from within R. Specific actions are controlled by values
passed on to the action_type
argument of the module
function, such as
loading and unloading software with load
and unload
, respectively.
Additionally, dedicated functions are provided for certain actions. The
following code examples are not evaluated since they will only work on systems where
an Environment Modules software is installed. A full list of actions and
additional functions for Environment Modules can be accessed with ?module
.
module(action_type="load", module_name="hisat2") moduleload("hisat2") # alternative command moduleunload("hisat2") modulelist() # list software loaded into current session moduleAvail() # list all software available in module system
Note, the module load/unload actions can be defined in the R/Rmd workflow
scripts or in the CWL parameter files. The listCmdModules
function can be
used, to list the names and versions of all software tools that are loaded via
Environment Modules in each step of a SAL
workflow container. Independent of
the usage of an Environment Modules system, all CL software used by each step
in a workflow can be listed with listCmdTools
. The output of both fumction
calls is not shown below for the same reason as in the previous code chunk.
listCmdModules(sal) listCmdTools(sal)
The processing time of computationally expensive steps can be greatly accelerated by
processing many input files in parallel using several CPUs and/or computer nodes
of an HPC or cloud system, where a scheduling system is used for load balancing.
To simplify for users the configuration and execution of workflow steps in serial or parallel mode,
systemPipeR
uses for both the same runWF
function. Parallelization simply
requires appending of the parallelization parameters to the settings of the corresponding workflow
steps each requesting the computing resources specified by the user, such as
the number of CPU cores, RAM and run time. These resource settings are
stored in the corresponding workflow step of the SAL
workflow container.
After adding the parallelization parameters, runWF
will execute the chosen steps
in parallel mode as instructed.
The following example applies to an alignment step of an RNA-Seq workflow. The
above demonstration workflow is not used here since it is too simple to benefit
from parallel processing. In the chosen alignment example, the parallelization
parameters are added to the alignment step (here hisat2_mapping
) of SAL
via
a resources
list. The given parameter settings will run 18 processes (Njobs
) in
parallel using for each 4 CPU cores (ncpus
), thus utilizing a total of 72 CPU
cores. The runWF
function can be used with most queueing systems as it is based on
utilities defined by the batchtools
package, which supports the use of
template files (*.tmpl
) for defining the run parameters of different
schedulers. In the given example below, a conffile
(see
.batchtools.conf.R
samples here) and
a template
file (see *.tmpl
samples
here) need to be present
on the highest level of a user's workflow project. The following example uses the sample
conffile
and template
files for the Slurm scheduler that are both provided by this
package.
The resources
list can be added to analysis steps when a workflow is loaded into SAL
.
Alternatively, one can add the resource settings with the addResources
function
to any step of a pre-populated SAL
container afterwards. For workflow steps with the same resource
requirements, one can add them to several steps at once with a single call to addResources
by
specifying multiple step names under the step
argument.
resources <- list(conffile=".batchtools.conf.R", template="batchtools.slurm.tmpl", Njobs=18, walltime=120, ## in minutes ntasks=1, ncpus=4, memory=1024, ## in Mb partition = "short" ) sal <- addResources(sal, step=c("hisat2_mapping"), resources = resources) sal <- runWF(sal)
The above example will submit via runWF(sal)
the hisat2_mapping step
to a partition (queue) called short
on an HPC cluster. Users need to adjust this and
other parameters, that are defined in the resources
list, to their cluster environment .
Workflows instances can be visualized as topology graphs with the plotWF
function.
The resulting plot includes the following information.
If no layout parameters are provided, then plotWF
will automatically detect reasonable settings
for a user's system, including width, height, layout, plot method, branch styles and others.
plotWF(sal, show_legend = TRUE, width = "80%", rstudio = TRUE)
For more details about the plotWF
function, please visit its help with ?plotWF
.
systemPipeR
produces two report types: Scientific and Technical. The
Scientific Report resembles a scientific publication detailing data analysis,
results, interpretation information, and user-provided text. The Technical
Report provides logging information useful for assessing workflow steps and
troubleshooting problems.
After a workflow run, systemPipeR's
renderReport
or rmarkdown's
render
function can be used to generate Scientific Reports in HTML, PDF or other
formats. The former uses the final SAL
instance as input, and the latter the
underlying Rmd source file. The resulting reports mimic research papers by combining
user-generated text with analysis results, creating reproducible analysis
reports. This reporting infrastructure offers support for citations,
auto-generated bibliographies, code chunks with syntax highlighting, and inline
evaluation of variables to update text content. Tables and figures in a report
can be automatically updated when the document is rebuilt or workflows are
rerun, ensuring data components are always current. This automation increases
reproducibility and saves time creating Scientific Reports. Furthermore, the
workflow topology maps described earlier can be incorporated into Scientific
Reports, enabling integration between Scientific and Technical Reports.
sal <- renderReport(sal) rmarkdown::render("my.Rmd", clean = TRUE, output_format = "BiocStyle::html_document")
Note, my.Rmd
in the last code line needs to be replaced with the name (path) of
the source Rmd
file used for generating the SAL
workflow container.
The package collects technical information about workflow runs in a project’s
log directory (default name: .SPRproject
). After partial or full completion
of a workflow, the logging information of a run is used by the renderLog
function to generate a Technical Report in HTML or other formats. The report
includes software execution commands, warnings and errors messages of each
workflow step. Easy visual navigation of Technical Reports is provided by
including an interactive instance of the corresponding workflow topology graph.
The technical details in these reports help assess the success of each workflow
step and facilitate troubleshooting.
sal <- renderLogs(sal)
The SAL workflow containers of systemPipeR
provide versatile conversion and
export options to both Rmd and executable Bash scripts. This feature not only
enhances the portability and reusability of workflows across different systems
but also promotes transparency, enabling efficient testing and
troubleshooting.
A populated SAL
workflow container can be converted to an Rmd file using the
sal2rmd
function. If needed, this Rmd
file can be used to construct a SAL
workflow container with the importWF
function as introduced above. This
functionality is useful for building templates of workflow Rmds and sharing
them with other systems.
sal2rmd(sal)
The sal2bash
function converts and exports workflows stored in SAL containers
into executable Bash scripts. This enables users to run their workflows as Bash
scripts from the command line. The function takes a SAL container as input and
generates a spr_wf.sh
file in the project's root directory as output.
Additionally, it creates a spr_bash
directory that stores all R-based workflow
steps as separate R scripts. To minimize the number of R scripts needed, the
function combines adjacent R steps into a single file.
sal2bash(sal)
The ability to restart existing workflows projects is important for continuing analyses that could not be completed, or to make changes without repeating already completed steps. Two main options are provided to restart workflows. Another option is provided that resets workflows to the very beginning, which effectively deletes the previous environment.
1. The resume=TRUE
option will initialize the latest instance of a SAL
object stored in the logs.dir
including its log files. When this option is used, a workflow can be continued where it was left off,
for example after closing and restarting R from the same directory on the same system. If the project was created
with custom directory and/or file names, then those names need to be specified under the log.dir
and sys.file
arguments of the SPRproject
function, respectively, otherwise the default names will be used.
sal <- SPRproject(resume = TRUE)
If the R environment was saved, one can recover with load.envir=TRUE
all
objects that were created during the previous workflow run. The same is possible with
the restart
option. For more details, please consult the help for the runWF
function.
sal <- SPRproject(resume = TRUE, load.envir = TRUE)
After resuming the workflow with load.envir
enabled, one can inspect the objects
created in the old environment, and decide if it is necessary to copy them to the
current environment.
viewEnvir(sal) copyEnvir(sal, list="plot", new.env = globalenv())
2. The restart=TRUE
option will also use the latest instance of the SAL
object stored in
the logs.dir
, but the previous log files will be deleted.
sal <- SPRproject(restart = TRUE)
3. The overwrite=TRUE
option will reset the workflow project to the very beginning by deleting the
log.dir
directory (.SPRproject
) that stores the previous SAL
instance and all its log files. At the same time
a new and empty ‘SAL’ workflow container will be created. This option should be used with caution
since it will effectively delete the workflow environment. Output files written by the
workflow steps to the results
directory will not be deleted when this option is used.
sal <- SPRproject(overwrite = TRUE)
This section describes methods for accessing, subsetting and modifying SAL
workflow objects.
Workflows and their run status can be explored further using a range of
accessor functions for SAL
objects.
The number of steps in a workflow can be determined with the length
function.
length(sal)
The corresponding names of workflow steps can be returned with stepName
.
stepName(sal)
CL software used by each step in a workflow can be listed with listCmdTools
.
listCmdTools(sal)
Some computing systems (often HPC clusters) allow users to load CL software via
an Environment Modules system into their PATH.
If this is the case, then the exact verions of the software tools loaded via the
module system can be listed for SAL
and SYSargs2
objects with listCmdModules
and modules
, respectively. The example workflow used here
does not make use of Environment Modules, and thus there are no software tools
to list here.
listCmdModules(sal) modules(stepsWF(sal)$gzip)
For more information on how to work with Environment Modules in systemPipeR
, please
visit the help with ?module
, ?modules
and ?listCmdModules
.
Several accessor functions are named after the corresponding slot names in
SAL
objects. This makes it easy to look them up with names()
, and then
apply them to sal
as the only argument, such as runInfo(sal)
.
names(sal)
The individual workflow steps in a SAL
container are stored as SYSargs2
and LineWise
components. They can be returned with the stepsWF
function.
stepsWF(sal)
The accessor function of SYSargs2
and LineWise
objects can be returned similarly
(here for gzip
step).
names(stepsWF(sal)$gzip)
The statusWF
function returns a status summary for each step in a SAL
workflow instance.
statusWF(sal)
The targets
instances for each step in a workflow can be returned with targetsWF
. The
below applies it to the second step.
targetsWF(sal[2])
If a workflow contains sample comparisons, that have been specified in the header
lines of a targets file starting with a # <CMP> tag
, then they can be returned
with the targetsheader
functions. This does not apply to the current demo sal
instance, and thus the function returns NULL
. For more details, consult the targets
file section here.
targetsheader(sal, step = "Quality")
The outfiles
component of a SAL
object stores the paths to the expected outfiles files
for each step in a workflow. Some of them are the input for downstream workflow steps.
outfiles(sal[2])
The dependency
step(s) in a workflow can be obtained with the
dependency
function. This information is used to construct the toplogy
graph of a workflow (see here).
dependency(sal)
The sample names (IDs) stored in the corresponding column of a targets file
can be returned with the SampleName
function.
SampleName(sal, step = "gzip")
The getColumn
method can be used to obtain the paths to the files generated in a
specified step.
getColumn(sal, "outfiles", step = "gzip", column = "gzip_file") getColumn(sal, "targetsWF", step = "gzip", column = "FileName")
The yamlinput
function returns the parameters of a workflow step defined in the
corresponding yml file.
yamlinput(sal, step = "gzip")
The exact syntax for running CL software on each input data set in a workflow can
be returned with the cmdlist
function. The CL calls are assembled from the corresponding
yml
and cwl
, and an optional targets
file as described in the above CLI section
here. The example below shows the CL syntax for running gzip
and gunzip
on the first input sample. Evaluating the output of cmdlist
can
be very helpful for designing and debugging CWL parameter files to support new CL
software or changing their settings.
cmdlist(sal, step = c(2,3), targets = 1)
Similarly, the codeLine
function returns the R code of a LineWise
workflow step.
codeLine(sal, step = "export_iris")
The objects generated in a workflow's run environment can be accessed with viewEnvir
.
viewEnvir(sal)
If needed one or multiple objects can be copied from the run environment of a workflow to the current environment of an R session.
copyEnvir(sal, list = c("plot"), new.env = globalenv(), silent = FALSE)
The bracket operator can be used to subset workflow by steps. For instance, the current
instance of sal
has four steps, and sal[1:2]
will subset the workflow to the first two
steps.
length(sal) sal[1:2]
In addition to subsetting by steps, one can subset workflows by input samples. The following illustrates this for the first two input samples, but omits the function output for brevity.
sal_sub <- subset(sal, subset_steps = c(2,3), input_targets = c("SE", "VE"), keep_steps = TRUE) stepsWF(sal_sub) targetsWF(sal_sub) outfiles(sal_sub)
For appending workflow steps, one can use the +
operator.
sal[1] + sal[2] + sal[3]
Replacement methods are implemented to make adjustments to certain paramer settings and R code in workflow steps.
## create a copy of sal for testing sal_c <- sal ## view original value, here restricted to 'ext' slot yamlinput(sal_c, step = "gzip")$ext ## Replace value under 'ext' yamlinput(sal_c, step = "gzip", paramName = "ext") <- "txt.gz" ## view modified value, here restricted to 'ext' slot yamlinput(sal_c, step = "gzip")$ext ## Evaluate resulting CL call cmdlist(sal_c, step = "gzip", targets = 1)
Code lines can be added with appendCodeLine
to R steps (LineWise
) as shown in the
following example.
appendCodeLine(sal_c, step = "export_iris", after = 1) <- "log_cal_100 <- log(100)" codeLine(sal_c, step = "export_iris")
In addition, code lines can be replaced with the replaceCodeLine
function.
For additional details about the LineWise
class, please see the example
above and the detailed description of the LineWise
class
here.
replaceCodeLine(sal_c, step="export_iris", line = 2) <- LineWise(code={ log_cal_100 <- log(50) }) codeLine(sal_c, step = "export_iris")
Renaming of workflow steps is possible with the renameStep
function.
renameStep(sal_c, c(1, 2)) <- c("newStep2", "newIndex") sal_c names(outfiles(sal_c)) names(targetsWF(sal_c)) dependency(sal_c)
The replaceStep
function can be used to replace entire workflow steps. When
replacing workflow steps, it is important to maintain a valid dependency graph
among the affected steps.
sal_test <- sal[c(1,2)] replaceStep(sal_test, step = 1, step_name = "gunzip" ) <- sal[3] sal_test
If needed, workflow steps can be removed as follows.
sal_test <- sal[-2] sal_test
This section provides a concise overview of CWL
and its utilization within systemPipeR
. It covers fundamental CWL concepts, including
the CommandLineTool
and Workflow
classes for describing individual CL processes and
workflows. For further details, readers want to refer to CWL's comprehensive
CommandLineTool and
Workflow documentation, as well
as the examples provided in CWL's Beginner Tutorial
and User Guide. Additionally, familiarizing oneself
with CWL's YAML format
specifications can be beneficial.
As illustrated in the introduction (Fig 2), CWL files with the '.cwl
'
extension define the parameters of a specific CL step or workflow, while files
with the '.yml
' extension define their input values.
CommandLineTool
{#cwl-clt}A Command Line Tool (CommandLineTool
class) specifies a standalone process
that can be run by itself (without including interactions with other
programs), and has inputs and outputs.
The following inspects the basic structure of a '.cwl
' sample file for a CommandLineTool
that is provided by this package.
dir_path <- system.file("extdata/cwl", package = "systemPipeR") cwl <- yaml::read_yaml(file.path(dir_path, "example/example.cwl"))
Important components include:
1. cwlVersion
: version of CWL specification used by file.
2. class
: declares description of a CommandLineTool
.
cwl[1:2]
3. baseCommand
: name of CL tool.
cwl[3]
4. inputs
: defines input information to run the tool. This includes:
id
: each input has an id
including name.type
: type of input value; one of string
, int
, long
, float
, double
,
File
, Directory
or Any
.inputBinding
: indicates if the input parameter should appear in CL call. If
missing input will not appear in the CL call.cwl[4]
5.. outputs
: list of expected outputs after running the CL tool. Important components are:
id
: each input has an id
including name.type
: type of output value; one of string
, int
, long
, float
, double
,
File
, Directory
, Any
or stdout
);outputBinding
: defines how to set outputs values; glob
specifies output value's name. cwl[5]
6. stdout
: specifies filename
for standard output. Note, by default systemPipeR
constructs the output filename
from results_path
and SampleName
(see above).
cwl[6]
Workflow
{#cwl-wf}CWL's Workflow
class describes one or more workflow steps, declares
their interdependencies, and defines how CommandLineTools
are executed.
Its CWL file includes inputs, outputs, and steps.
The following illustrates the basic structure of a '.cwl
' sample file for a Workflow
that is provided by this package.
cwl.wf <- yaml::read_yaml(file.path(dir_path, "example/workflow_example.cwl"))
1. cwlVersion
: version of CWL specification used by file.
2. class
: declares description of a Workflow
that describes one or
more CommandLineTools
and their combined usage.
cwl.wf[1:2]
3. inputs
: defines the inputs of the workflow.
cwl.wf[3]
4. outputs
: defines the outputs of the workflow.
cwl.wf[4]
5. steps
: describes the steps of the workflow. The example below shows one step.
cwl.wf[5]
The .yml
file provides the input values for the parameters described above.
The following example includes input values for three parameters (message
,
SampleName
and results_path
).
yaml::read_yaml(file.path(dir_path, "example/example_single.yml"))
Note, the .yml
file needs to provide input values for each input parameter
specified in the corresponding .cwl
file (compare cwl[4]
above).
cwl
, yml
and targets
{#cwl_targets}This section illustrates how the parameters in CWL files (cwl
and yml
) are
interconnected to construct CL calls of steps, and subsequently assembled
to workflows.
A SAL
container (long name SYSargsList
) stores all information and instructions
needed for processing a set of inputs (incl. files) with a single or many CL steps within a workflow
The SAL
object is created and fully populated with the SYSargsList
constructor
function. More detailed documentation of SAL
workflow instances is available
here and here.
The following imports the .cwl
and .yml
files for running the echo Hello World!
example.
HW <- SYSargsList(wf_file = "example/workflow_example.cwl", input_file = "example/example_single.yml", dir_path = system.file("extdata/cwl", package = "systemPipeR")) HW cmdlist(HW)
The example provided is restricted to creating a CL call for a single input
(sample). To process multiple inputs, a straightforward approach is to assign
variables to the corresponding parameters instead of using fixed (hard-coded)
values. These variables can then be assigned the desired input values
iteratively, resulting in multiple CL calls, one for each input value. The
following illustrates this with an example, where the message
and SampleName
parameters are assigned variables that are labeled with tags of the form
_XXX_
. These variables will be assigned values stored in a targets
file.
yml <- yaml::read_yaml(file.path(dir_path, "example/example.yml")) yml
The content of the targets
file used for this example is shown below. The
general structure of targets
files is explained above.
targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR") read.delim(targetspath, comment.char = "#")
In the simple example given above the values stored under the Message
and
SampleName
columns of the targets file will be passed on to the corresponding
parameters with matching names in the yml
file, and from there to
the echo
command defined in the cwl
file (see here).
As mentioned previously, the usage of targets
files is optional in
systemPipeR
. Since targets
files provide an easy and efficient solution for
organizing experimental variables, their usage is highly encouraged and well
supported in systemPipeR
.
The SYSargsList
function constructs SAL
instances from the three parameter
files, that were introduced above (see Fig 3). The path to each file is assigned to its own
argument: wf_file
is assigned the path of a cwl
workflow file, input_file
the path of a yml
input file, and targets
the path of a targets
file. Additionally, a named
vector is provided under the inputvars
argument that defines which column
values in the targets
file are assigned to specific parameters in the yml
file. A parameter connection is established where a name in inputvars
has
matching column and parameter names in the targets
and yml
files,
respectively (Fig 3). A tagging syntax with the pattern _XXX_
is used to
indicate which parameters contain variables that will be assigned values from
the targets
file. The usage of this pattern is only recommended for
consistency and easy identification, but not enforced.
The SYSargslist
function call constructs the echo
commands (CL calls) based on the
parameters provided by the above described parameter file instances (cwl
, yml
and targets
)
as well as the variable mappings specified under the inputvars
argument.
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
The final CL calls (here echo
command) can be returned with the cmdlist
for
each string given under the Message
column of the targets
file. The values under
the SampleName
column are used to name the corresponding output files, each with a
txt
extension.
cmdlist(HW_mul)
To streamline the process of generating CWL parameter files (both cwl
and
yml
), users can simply provide the CL syntax for executing new software. This
action will automatically create the corresponding CWL parameter files, which
alleviates the need for manual creation of CWL files, reducing the
burden on users. This functionality is implemented in systemPipeR’s
createParam
function group.
To use this functionality, CL calls need to be provided in a pseudo-bash script format
and stored as a character vector
.
The following uses as example a CL call for the HISAT2 software.
hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
For the CL call above, the corresponding pseudo-bash syntax is given below.
Here, the CL string needs to be stored in a single slot of a character vector
,
here named command
. The formatting requirements for the CL string will be explained
next.
command <- " hisat2 \ -S <F, out: ./results/M1A.sam> \ -x <F: ./data/tair10.fasta> \ -k <int: 1> \ -min-intronlen <int: 30> \ -max-intronlen <int: 3000> \ -threads <int: 4> \ -U <F: ./data/SRR446027_1.fastq.gz> "
Format specifications for pseudo-bash syntax (Version 1)
\
at the end of a line.baseCommand
). It can include a subcommand, such as in git commit
where commit
is a subcommand. -S
and --min
. Values that should be assigned to
arguments are placed inside <...>
, also on the same line. Arguments and flags without
values lack this assignment. <TYPE:
, where TYPE
can be
F
for "File", "int", "string", etc. out
can be added after TYPE
separated by a comma.:
is used to separate keywords and default values. Any non-space value after the :
will be treated as a default value. Note, the above specifications are Version 1 (v1
) of the pseudo-bash syntax
used by the createParam
function below. There also is a Version 2 (v2
)
specification that supports additional features, but comes with more syntax
restrictions. Details on this are available in the help of the createParam
function.
createParam
FunctionThe createParam
function accepts as input a CL string that is formatted in the above
pseudo-bash syntax. As output it creates the corresponding CWL files (cwl
and yml
)
that will be written to the default directory: ./param/cwl/
. This path can be changed
under the file
argument. In addition, it constructs for the given CL string the corresponding
SYSargs2
object (here assigned to cmd
). The information printed as console output
contains the original CL string that is included for checking purposes. This CL string is
not included to the resulting CWL files.
cmd <- createParam(command, writeParamFiles = TRUE, overwrite=TRUE, confirm=TRUE)
Next, the cmdlist
can be used to check the correctness of the CL call defined
by the CWL parameter files generated by the createParam
command above.
cmdlist(cmd)
If the createParam
function is executed without creating the CWL parameter files right
away (argument setting writeParamFiles=FALSE
) then these files can be generated in a
separate step with writeParamFiles
.
writeParamFiles(cmd, overwrite = TRUE)
targets
fileThe following gives a more complete example where the CWL files are first created for a CL string,
and then loaded together with a targets
file into a SYSargs2
object. Next, the final CL calls
are assembled for each input sample with the renderWF
function. The final CL calls can then be
inspected with the cmdlist
function, where the below shows only the first 2 of a total of
18 CL calls for brevity.
command2 <- " hisat2 \ -S <F, out: _SampleName_.sam> \ -x <F: ./data/tair10.fasta> \ -k <int: 1> \ -min-intronlen <int: 30> \ -max-intronlen <int: 3000> \ -threads <int: 4> \ -U <F: _FASTQ_PATH1_> " WF <- createParam(command2, 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]
The following introduces several accessor and replacement functions that are useful for creating and revising CWL parameter files.
printParam(cmd, position = "baseCommand") ## Return baseCommand printParam(cmd, position = "outputs") ## Return outputs printParam(cmd, position = "inputs", index = 1:2) ## Return components by index printParam(cmd, position = "inputs", index = -1:-2) ## Negative index subsetting
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)
cmd3 <- replaceParam(cmd, "base", index = 1, replace = list(baseCommand = "bwa")) cmdlist(cmd3) ## Replacement changed baseCommand
new_inputs <- new_inputs <- list( "new_input1" = list(type = "File", preF="-b", yml ="myfile"), "new_input2" = "-L <int: 4>" ) cmd4 <- replaceParam(cmd, "inputs", index = 1:2, replace = new_inputs) cmdlist(cmd4)
newIn <- new_inputs <- list( "new_input1" = list(type = "File", preF="-b1", yml ="myfile1"), "new_input2" = list(type = "File", preF="-b2", yml ="myfile2"), "new_input3" = "-b3 <F: myfile3>" ) cmd5 <- appendParam(cmd, "inputs", index = 1:2, append = new_inputs) cmdlist(cmd5) cmd6 <- appendParam(cmd, "inputs", index = 1:2, after=0, append = new_inputs) cmdlist(cmd6)
output
parametersnew_outs <- list( "sam_out" = "<F: $(inputs.results_path)/test.sam>" ) cmd7 <- replaceParam(cmd, "outputs", index = 1, replace = new_outs) output(cmd7)
The workflow steps of SAL
(synonym SYSargsList
) containers are composed of SYSargs2
and/or LineWise
objects. These two classes are introduced here in more detail.
SYSargs2
class {#sysargs2}The SYSargs2
class stores workflow steps that run CL software. An instance of
SYSargs2
stores all the input/output paths and parameter components necessary
for executing a specific CL data analysis step. SYSargs2
instances are
created using two constructor functions: loadWF
and renderWF
. These
functions make use of a targets
(or yml
) and the two CWL parameter files
cwl
and yml
. The structure and content for the CWL files are described
above. The following creates a SYSargs2
instance using the cwl
and
yml
files for running the RNA-Seq read aligner HISAT2 [@Kim2015-ve]. Note,
when using the SYSargsList
method for constructing workflow steps
(see here), then the user will not need to run the loadWF
and renderWF
directly.
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_"))
In addition to SAL
objects (see here), the cmdlist
function accepts
SYSargs2
to constructs CL calls based on the parameter inputs imported from the
corresponding targets
, yml
and cwl
files.
cmdlist(WF)[1]
Several accessor methods are available that are named after the slot names of
SYSargs2
objects.
names(WF)
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,
e.g. a downstream SYSargs2
instance.
output(WF)[1]
The targets
method allows access to the targets
component of a SYSargs2
object. Refer to the information provided above for an
explanation of the targets
file structure.
targets(WF)[1] as(WF, "DataFrame")
If CL software is loaded via an Environment Modules system
into a user's PATH
, then this information can be accessed in SYSargs2
objects as shown
below. For more details on working with Environment Modules, see here.
modules(WF)
Additional accessible information includes the location of the parameters files,
inputvars
(see here) and more.
files(WF) inputvars(WF)
To define R code as workflow steps, the LineWise
class is used. The syntax
for declaring lines of R code as workflow steps in R or Rmd files is introduced
in the workflow design section. The following showcases
additional utilities for LineWise
objects.
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 a LineWise
object to a list
object and vice versa.
lw <- stepsWF(sal_lw)[[2]] ## Coerce ll <- as(lw, "list") class(ll) lw <- as(ll, "LineWise") lw
Accessing basic information of LineWise
objects.
length(lw) names(lw) codeLine(lw) codeChunkStart(lw) rmdPath(lw)
Subsetting LineWise
objects.
l <- lw[2] codeLine(l) l_sub <- lw[-2] codeLine(l_sub)
Replacement methods for changing R code in LineWise
objects.
replaceCodeLine(lw, line = 2) <- "5+5" codeLine(lw) appendCodeLine(lw, after = 0) <- "6+7" codeLine(lw)
For comparison, similar replacement methods are available for SAL
objects. They have been
covered above.
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)
Here is a partial list of CL software for which systemPipeR
includes CWL
parameter file templates. Notably, with the newly added auto-creation feature
for CWL files (see here), generating CWL parameter files for most CL
tools has become straightforward. Thus, maintaining and extending this list will
not be necessary anymore.
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, Step, Description) %>% dplyr::arrange(Tool) %>% kableExtra::kable("html", escape = FALSE, col.names = c("Tool Name", "Description", "Step")) %>% kableExtra::kable_material(c("striped", "hover", "condensed")) %>% kableExtra::scroll_box(width = "80%", height = "500px")
To run any of the tools mentioned, users must ensure that the necessary
software is installed on their system and added to the PATH
. There are
several methods to verify if the required tools/modules are installed. The
easiest method is automatically executed for users when they call the importWF
function, or just tryCL(<base_command>)
. In the print message of importWF
, all
necessary tools and modules are automatically listed and checked for users.
For additional tool validation methods, please refer to these instructions:
Five Minute Tutorial, Environment Modules, and
Managing Workflows.
if (file.exists(".SPRproject")) unlink(".SPRproject", recursive = TRUE) ## NOTE: Removing previous project create in the quick starts section
This section presents various data analysis functionalities that are valuable for many workflows. Some of these functionalities are R functions, while others are CWL interfaces to widely used CL software. A few of them are included for historical reasons.
To work with the following examples a new workflow project needs to be created.
The below includes the overwrite=TRUE
setting to remove any already
project directory.
sal <- SPRproject(projPath = getwd(), overwrite = TRUE)
The first step in the new workflow project is to load the systemPipeR
package.
appendStep(sal) <- LineWise({ library(systemPipeR) }, step_name = "load_SPR")
Importantly, in order to use the individual appendStep
operations below, one has
to pay attention to the step dependencies.
preprocessReads
functionThe function preprocessReads
allows to apply predefined or custom
read preprocessing functions to the 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.
In this step, the preprocessing parameters defined in the corresponding
*.pe.cwl
and *.pe.yml
files are added to a previously created SAL
object.
This preprocessing step is crucial for preparing the reads for further
analysis.
targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR") appendStep(sal) <- SYSargsList( step_name = "preprocessing", targets = targetspath, 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 preprocessReads
function. Here, it is possible to replace the function used on the
preprocessing
step and modify the corresponding sal
object. Because it is a
custom function, it is necessary to save this part in the R object, and
internally the preprocessReads.doc.R
script, that is stored in the param
directory
of the workflow templates, is loading the custom function. If the R
object is saved with a different name (here "param/customFCT.RData"
), one has
to adjust the corresponding path in the preprocessReads.doc.R
script.
First, the custom function is defined.
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 this the input parameters can be edited as shown here.
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
TrimGalore! is a wrapper tool for Cutadapt and FastQC to consistently apply quality and adapter trimming to fastq files.
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")
Trimmomatic software [@Bolger2014-yr] performs a variety of useful trimming tasks for Illumina paired-end and single ended reads. The following is an example of how to perform this task.
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")
The following seeFastq
and seeFastqPlot
functions generate and plot a series of useful
quality statistics for a set of FASTQ files, including per cycle quality box
plots, base proportions, base-level quality trends, relative k-mer
diversity, length, and occurrence distribution of reads, number of reads
above quality cutoffs and mean quality distribution. The results are
written to a PDF file named fastqReport.pdf
.
appendStep(sal) <- LineWise(code = { fastq <- getColumn(sal, step = "preprocessing", "targetsWF", column = 1) fqlist <- seeFastq(fastq = fastq, batchsize = 10000, klength = 8) pdf("./results/fastqReport.pdf", height = 18, width = 4*length(fqlist)) seeFastqPlot(fqlist) dev.off() }, step_name = "fastq_report", dependency = "preprocessing")
After quality control, the sequence reads can be aligned to a reference genome or transcriptome. The following gives examples for running several NGS read aligners.
HISAT2
The following steps demonstrate how to run the HISAT2
short read aligner
[@Kim2015-ve] from systemPipeR
.
To use an NGS aligner, one has to first index the reference genome. This is done
below with hisat2-build
.
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 append the alignment
step to the sal
workflow container. In this step several post-processing steps
with Samtools
are included to convert the SAM files, that were generated by HISAT2
,
to indexed and sorted BAM files. Those sub-steps are defined by the corresponding CWL workflow file
(see workflow_hisat2-se.cwl).
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(FileName1="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("hisat_index"), run_session = "compute")
Tophat2
The Bowtie2/Tophat2
suite is the predecessor of Hisat2
[@Kim2013-vg; @Langmead2012-bs].
How to run it via CWL is shown below.
First, the reference genome has to be indexed for Bowtie2
.
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")
Next, the alignment step is constructed with the parameter files workflow_tophat2-mapping.cwl
and tophat2-mapping-pe.yml
.
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")
Bowtie2
The following example runs Bowtie2
by itself (without Tophat2
or Hisat2
).
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")
BWA-MEM
The following example runs BWA-MEM, an aligner that is widely used for VAR-Seq experiments.
First, the reference genome has to be indexed for BWA-MEM
.
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")
Next, the reads can be aligned with BWA-MEM
.
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")
Rsubread
Rsubread
is an R package for processing short and long reads. It is well known for its
fast and accurate mapping performance of RNA-Seq reads.
First, the reference genome has to be indexed for Rsubread
.
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")
Next, the RNA-Seq reads can be aligned with Rsubread
.
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(FileName1="_FASTQ_PATH1_", SampleName="_SampleName_"), dependency = c("rsubread_index"), run_session = "compute", run_step = "optional")
gsnap
The gmapR
package provides an interface to work with the GSNAP
and GMAP
aligners from R [@Wu2010-iq].
First, the reference genome has to be indexed for GSNAP
.
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")
Next, the RNA-Seq reads are aligned with GSNAP
.
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")
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 https access needs to be
available in the user account (e.g. home/.html
) 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 https 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. To make the following code work, users need to change
the directory name (here somedir/
) and the username (here <username>
) to the
corresponding names on their system.
appendStep(sal) <- LineWise( code = { bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam") symLink2bam( sysargs = bampaths, htmldir = c("~/.html/", "somedir/"), urlbase = "https://cluster.hpcc.ucr.edu/<username>/", urlfile = "./results/IGVurl.txt") }, step_name = "bam_IGV", dependency = "hisat_mapping", run_step = "optional" )
Reads overlapping with annotation ranges of interest are counted for each
sample using the summarizeOverlaps
function [@Lawrence2013-kt].
First, the gene annotation ranges from a GFF file are stored in a TxDb
container for
efficient work with genomic features.
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")
Next, The read counting is preformed for exonic gene regions in a non-strand-specific manner while ignoring overlaps among different genes.
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.
The following provides an overview of the number of reads in each sample and how many of them aligned to the reference.
appendStep(sal) <- LineWise({ read_statsDF <- alignStats(args) write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE, sep = "\t") }, step_name = "align_stats", dependency = "hisat_mapping")
The following shows the first four lines of the sample alignment stats file
provided by the systemPipeR
package. For simplicity the number of PE reads
is multiplied here by 2 to approximate proper alignment frequencies where each
read in a pair is counted.
read.table(system.file("extdata", "alignStats.xls", package = "systemPipeR"), header = TRUE)[1:4,]
Example of downloading a GFF file for miRNA ranges from an organism of interest (here A. thaliana), and then use them for read counting, here RNA-Seq reads from the above steps.
appendStep(sal) <- LineWise({ system("wget https://www.mirbase.org/download/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")
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.
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")
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
vignette
for more detailed information on this topic and how to properly run edgeR
on data sets with more complex experimental designs.
appendStep(sal) <- LineWise({ targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR") targets <- read.delim(targetspath, comment = "#") cmp <- readComp(file = targetspath, format = "matrix", delim = "-") countDFeBygpath <- system.file("extdata", "countDFeByg.xls", package = "systemPipeR") countDFeByg <- read.delim(countDFeBygpath, row.names = 1) edgeDF <- run_edgeR(countDF = countDFeByg, targets = targets, cmp = cmp[[1]], independent = FALSE, mdsplot = "") DEG_list <- filterDEGs(degDF = edgeDF, filter = c(Fold = 2, FDR = 10)) }, step_name = "edger", dependency = "read_counting")
Filter and plot DEG results for up and down-regulated genes. Because of the
small size of the toy data set used by this vignette, the FDR cutoff 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.
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
vignette
for more detailed information on this topic and how to properly run DESeq2
on data sets with more complex experimental designs.
appendStep(sal) <- LineWise({ degseqDF <- run_DESeq2(countDF=countDFeByg, targets=targets, cmp=cmp[[1]], independent=FALSE) DEG_list2 <- filterDEGs(degDF=degseqDF, filter=c(Fold=2, FDR=10)) }, step_name = "deseq2", dependency = "read_counting")
The function overLapper
can compute Venn intersects for large numbers of
sample sets (up to 20 or more) and vennPlot
can plot 2-5 way Venn diagrams.
A useful feature is the possibility to combine the counts from several Venn
comparisons with the same number of sample sets in a single Venn diagram (here
for 4 up and down DEG sets).
appendStep(sal) <- LineWise({ vennsetup <- overLapper(DEG_list$Up[6:9], type="vennsets") vennsetdown <- overLapper(DEG_list$Down[6:9], type="vennsets") vennPlot(list(vennsetup, vennsetdown), mymain="", mysub="", colmode=2, ccol=c("blue", "red")) }, step_name = "vennplot", dependency = "edger")
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 step.
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")
Apply the enrichment analysis to the DEG sets obtained in the above
differential expression analysis. Note, in the following example the FDR
filter is set here to an unreasonably high value, simply because of the small
size of the toy data set used in this vignette. Batch enrichment analysis of
many gene sets is performed with the GOCluster_Report
function. When
method="all"
, it returns all GO terms passing the p-value cutoff specified
under the cutoff
arguments. When method="slim"
, it returns only the GO
terms specified under the myslimv
argument. The given example shows how one
can obtain such a GO slim vector from BioMart for a specific organism.
appendStep(sal) <- LineWise({ load("data/GO/catdb.RData") DEG_list <- filterDEGs(degDF=edgeDF, filter=c(Fold=2, FDR=50), plot=FALSE) up_down <- DEG_list$UporDown; names(up_down) <- paste(names(up_down), "_up_down", sep="") up <- DEG_list$Up; names(up) <- paste(names(up), "_up", sep="") down <- DEG_list$Down; names(down) <- paste(names(down), "_down", sep="") DEGlist <- c(up_down, up, down) DEGlist <- DEGlist[sapply(DEGlist, length) > 0] BatchResult <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="all", id_type="gene", CLSZ=2, cutoff=0.9, gocats=c("MF", "BP", "CC"), recordSpecGO=NULL) library("biomaRt") m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org") goslimvec <- as.character(getBM(attributes=c("goslim_goa_accession"), mart=m)[,1]) BatchResultslim <- GOCluster_Report(catdb=catdb, setlist=DEGlist, method="slim", id_type="gene", myslimv=goslimvec, CLSZ=10, cutoff=0.01, gocats=c("MF", "BP", "CC"), recordSpecGO=NULL) gos <- BatchResultslim[grep("M6-V6_up_down", BatchResultslim$CLID), ] gos <- BatchResultslim pdf("GOslimbarplotMF.pdf", height=8, width=10); goBarplot(gos, gocat="MF"); dev.off() goBarplot(gos, gocat="BP") goBarplot(gos, gocat="CC") }, step_name = "go_enrichment", dependency = "get_go_biomart")
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.
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 joining.
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"))
sessionInfo()
This project is funded by awards from the National Science Foundation (ABI-1661152], and the National Institute on Aging of the National Institutes of Health (U19AG023122).
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.