systemPipeR 2.11.3
systemPipeR
is a versatile workflow environment for data analysis that integrates R with command-line (CL) software (H Backman and Girke 2016). 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 1). Central to systemPipeR
is a CL interface (CLI) that adopts the Common Workflow Language (CWL, Crusoe et al. 2021). 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.
Figure 1: Important functionalities of systemPipeR
(A) Illustration of workflow design concepts, and (B) examples of visualization functionalities for NGS data.
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 2). 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.
Figure 2: Workflow management class
Workflows are defined and managed by the SYSargsList
(SAL
) control class. Components of SAL
include SYSargs2
and/or LineWise
for defining CL- and R-based workflow steps, respectively. The former are constructed from a targets
and two CWL param files, and the latter comprises mainly R code.
systemPipeR
adopts the Common Workflow Language (CWL) (Amstutz et al. 2016), which is a
widely used community standard for describing CL tools and workflows
in a declarative, generic, and reproducible manner. CWL specifications are
text-based YAML (https://yaml.org/) files that are straightforward to create and
to modify. Integrating CWL in systemPipeR
enhances the sharability, standardization,
extensibility and portability of data analysis workflows.
Following the CWL Specifications, the basic
description for executing a command-line software are defined by two files: a
cwl step definition file and a yml configuration file. Figure
3 illustrates the utilitity of the two files using “Hello World”
as an example. The cwl file (A) defines the CL tool (C) along with
its parameters, and the yml file (B) assigns values to the corresponding parameters.
For convenience, parameter values can be provided by 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.
Figure 3: Parameter files
Illustration how the different fields in cwl, yml and targets files are connected to assemble command-line calls, here for ‘Hello World’ example.
The package also provides several convenience
functions that are useful for designing and testing workflows, such as a
command-line rendering function that assembles from the parameter files (cwl, yml and
targets) the exact command-line strings for each step prior to running a command-line tool.
Auto-generation of CWL parameter files is also supported. Here, users can simply
provide the command-line strings for a command-line software of interest to a rendering function that generates
the corresponding *.cwl
and *.yml
files for them.
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
).
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
columns 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
Figure 4: Status check of workflow
The run status flags of each workflow step are given in its summary view.
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)
Figure 5: Toplogy graph of RNA-Seq workflow
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
fileA 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 = "#"))
## Loading required namespace: DT
To work with custom data, users need to generate a targets
file containing
the paths to their own FASTQ files and then provide under targetspath
the
path to the corresponding targets
file.
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]
## [1] "# Project ID: Arabidopsis - Pseudomonas alternative splicing study (SRA: SRP010938; PMID: 24098335)"
## [2] "# The following line(s) allow to specify the contrasts needed for comparative analyses, such as DEG identification. All possible comparisons can be specified with 'CMPset: ALL'."
## [3] "# <CMP> CMPset1: M1-A1, M1-V1, A1-V1, M6-A6, M6-V6, A6-V6, M12-A12, M12-V12, A12-V12"
## [4] "# <CMP> CMPset2: ALL"
The function readComp
imports the comparison information and stores it in a
list
. Alternatively, readComp
can obtain the comparison information from
a SYSargsList
instance containing the targets
file information (see below).
readComp(file = targetspath, format = "vector", delim = "-")
## $CMPset1
## [1] "M1-A1" "M1-V1" "A1-V1" "M6-A6" "M6-V6" "A6-V6" "M12-A12" "M12-V12" "A12-V12"
##
## $CMPset2
## [1] "M1-A1" "M1-V1" "M1-M6" "M1-A6" "M1-V6" "M1-M12" "M1-A12" "M1-V12" "A1-V1"
## [10] "A1-M6" "A1-A6" "A1-V6" "A1-M12" "A1-A12" "A1-V12" "V1-M6" "V1-A6" "V1-V6"
## [19] "V1-M12" "V1-A12" "V1-V12" "M6-A6" "M6-V6" "M6-M12" "M6-A12" "M6-V12" "A6-V6"
## [28] "A6-M12" "A6-A12" "A6-V12" "V6-M12" "V6-A12" "V6-V12" "M12-A12" "M12-V12" "A12-V12"
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)
## Creating directory: /tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/data
## Creating directory: /tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/param
## Creating directory '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject'
## Creating file '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject/SYSargsList.yml'
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
## Instance of 'SYSargsList':
## No workflow steps added
projectInfo(sal)
## $project
## [1] "/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes"
##
## $data
## [1] "data"
##
## $param
## [1] "param"
##
## $results
## [1] "results"
##
## $logsDir
## [1] ".SPRproject"
##
## $sysargslist
## [1] ".SPRproject/SYSargsList.yml"
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)
## [1] 0
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 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
## Instance of 'SYSargsList':
## No workflow steps added
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
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Pending
##
To extract the code of an R step stored in a SAL
object, the codeLine
method can be used.
codeLine(sal)
## export_iris
## mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
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 3). 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
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Pending
## 2. gzip --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 2.1. gzip
## cmdlist: 3 | Pending: 3
##
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")
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c results/setosa.csv > results/SE.csv.gz"
##
##
## $gzip$VE
## $gzip$VE$gzip
## [1] "gzip -c results/versicolor.csv > results/VE.csv.gz"
##
##
## $gzip$VI
## $gzip$VI$gzip
## [1] "gzip -c results/virginica.csv > results/VI.csv.gz"
# 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
## $gzip
## DataFrame with 3 rows and 1 column
## gzip_file
## <character>
## SE results/SE.csv.gz
## VE results/VE.csv.gz
## VI results/VI.csv.gz
# 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
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Pending
## 2. gzip --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 2.1. gzip
## cmdlist: 3 | Pending: 3
## 3. gunzip --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 3.1. gunzip
## cmdlist: 3 | Pending: 3
##
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'])
## $gunzip
## DataFrame with 3 rows and 2 columns
## gzip_file SampleName
## <character> <character>
## SE results/SE.csv.gz SE
## VE results/VE.csv.gz VE
## VI results/VI.csv.gz VI
As before, the output files of the new step can be returned with outfiles
.
outfiles(sal['gunzip'])
## $gunzip
## DataFrame with 3 rows and 1 column
## gunzip_file
## <character>
## SE results/SE.csv
## VE results/VE.csv
## VI results/VI.csv
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)
## $gunzip
## $gunzip$SE
## $gunzip$SE$gunzip
## [1] "gunzip -c results/SE.csv.gz > results/SE.csv"
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')
## SE VE VI
## "results/SE.csv" "results/VE.csv" "results/VI.csv"
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 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
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Pending
## 2. gzip --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 2.1. gzip
## cmdlist: 3 | Pending: 3
## 3. gunzip --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 3.1. gunzip
## cmdlist: 3 | Pending: 3
## 4. iris_stats --> Status: Pending
##
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")
## Creating directory '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject_rmd'
## Creating file '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject_rmd/SYSargsList.yml'
sal_rmd <- importWF(sal_rmd,
file_path = system.file("extdata", "spr_simple_wf.Rmd", package = "systemPipeR"))
## Reading Rmd file
##
## ---- Actions ----
## Checking chunk eval values
## Checking chunk SPR option
## Ignore non-SPR chunks: 17
## Parse chunk code
## Checking preprocess code for each step
## No preprocessing code for SPR steps found
## Now importing step 'load_library'
## Now importing step 'export_iris'
## Now importing step 'gzip'
## Now importing step 'gunzip'
## Now importing step 'stats'
## Now back up current Rmd file as template for `renderReport`
## Template for renderReport is stored at
## /tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject_rmd/workflow_template.Rmd
## Edit this file manually is not recommended
## Now check if required tools are installed
## Check if they are in path:
## Checking path for gzip
## PASS
## Checking path for gunzip
## PASS
## step_name tool in_path
## 1 gzip gzip TRUE
## 2 gunzip gunzip TRUE
## All required tools in PATH, skip module check. If you want to check modules use `listCmdModules`Import done
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).
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)
## Running Step: export_iris
## Running Session: Management
##
|
| | 0%
|
|==========================================================================================| 100%
## Step Status: Success
## Running Step: gzip
## Running Session: Management
##
|
| | 0%
|
|============================== | 33%
|
|============================================================ | 67%
|
|==========================================================================================| 100%
## ---- Summary ----
## Targets Total_Files Existing_Files Missing_Files gzip
## SE SE 1 1 0 Success
## VE VE 1 1 0 Success
## VI VI 1 1 0 Success
##
## Step Status: Success
## Running Step: gunzip
## Running Session: Management
##
|
| | 0%
|
|============================== | 33%
|
|============================================================ | 67%
|
|==========================================================================================| 100%
## ---- Summary ----
## Targets Total_Files Existing_Files Missing_Files gunzip
## SE SE 1 1 0 Success
## VE VE 1 1 0 Success
## VI VI 1 1 0 Success
##
## Step Status: Success
## Running Step: iris_stats
## Running Session: Management
##
|
| | 0%
|
|==========================================================================================| 100%
## Step Status: Success
## Done with workflow running, now consider rendering logs & reports
## To render logs, run: sal <- renderLogs(sal)
## From command-line: Rscript -e "sal = systemPipeR::SPRproject(resume = TRUE); sal = systemPipeR::renderLogs(sal)"
## To render reports, run: sal <- renderReport(sal)
## From command-line: Rscript -e "sal= s ystemPipeR::SPRproject(resume = TRUE); sal = systemPipeR::renderReport(sal)"
## This message is displayed once per R session
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
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Success
## 2. gzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 2.1. gzip
## cmdlist: 3 | Success: 3
## 3. gunzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 3.1. gunzip
## cmdlist: 3 | Success: 3
## 4. iris_stats --> Status: Success
##
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.
stepsWF(sal)
dependency(sal)
cmdlist(sal)
codeLine(sal)
targetsWF(sal)
statusWF(sal)
projectInfo(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 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.
- Workflow topology graph rendered based on dependencies among steps
- Workflow step status, e.g. Success, Error, Pending, Warnings
- Sample status and statistics
- Run time of individual steps
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, and interpretation information. 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 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
offer convenient conversion and
export options to Rmd and executable Bash scripts. This feature not only
enhances the portability and reusability of workflows across different systems
but also promotes transparency, facilitating 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)
If you desire to resume or restart a project that has been initialized in the past,
SPRproject
function allows this operation.
With the resume
option, it is possible to load the SYSargsList
object in R and
resume the analysis. Please, make sure to provide the logs.dir
location, and the
corresponded YAML
file name, if the default names were not used when the project was created.
sal <- SPRproject(resume = TRUE, logs.dir = ".SPRproject",
sys.file = ".SPRproject/SYSargsList.yml")
If you choose to save the environment in the last analysis, you can recover all
the files created in that particular section. SPRproject
function allows this
with load.envir
argument. Please note that the environment was saved only with
you run the workflow in the last section (runWF()
).
sal <- SPRproject(resume = TRUE, load.envir = TRUE)
After loading the workflow at your current section, you can check the objects created in the old environment and decide if it is necessary to copy them to the current environment.
viewEnvir(sal)
copyEnvir(sal, list="plot", new.env = globalenv())
The resume
option will keep all previous logs in the folder; however, if you desire to
clean the execution (delete all the log files) history and restart the workflow,
the restart=TRUE
option can be used.
sal <- SPRproject(restart = TRUE, load.envir = FALSE)
The last and more drastic option from SYSproject
function is to overwrite
the
logs and the SYSargsList
object. This option will delete the hidden folder and the
information on the SYSargsList.yml
file. This will not delete any parameter
file nor any results it was created in previous runs. Please use with caution.
sal <- SPRproject(overwrite = TRUE)
systemPipeR
provide several accessor methods and useful functions to explore
SYSargsList
workflow object.
Several accessor methods are available that are named after the slot names of
the SYSargsList
workflow object.
names(sal)
## [1] "stepsWF" "statusWF" "targetsWF" "outfiles"
## [5] "SE" "dependency" "targets_connection" "projectInfo"
## [9] "runInfo"
length(sal)
## [1] 4
stepsWF(sal)
## $export_iris
## Instance of 'LineWise'
## Code Chunk length: 1
##
## $gzip
## Instance of 'SYSargs2':
## Slot names/accessors:
## targets: 3 (SE...VI), targetsheader: 1 (lines)
## modules: 0
## wf: 1, clt: 1, yamlinput: 4 (inputs)
## input: 3, output: 3
## cmdlist: 3
## Sub Steps:
## 1. gzip (rendered: TRUE)
##
##
##
## $gunzip
## Instance of 'SYSargs2':
## Slot names/accessors:
## targets: 3 (SE...VI), targetsheader: 1 (lines)
## modules: 0
## wf: 1, clt: 1, yamlinput: 4 (inputs)
## input: 3, output: 3
## cmdlist: 3
## Sub Steps:
## 1. gunzip (rendered: TRUE)
##
##
##
## $iris_stats
## Instance of 'LineWise'
## Code Chunk length: 5
cmdlist()
method printing the system commands for running command-line
software as specified by a given *.cwl
file combined with the paths to the
input samples (e.g. FASTQ files) provided by a targets
file. The example below
shows the cmdlist()
output for running gzip
and gunzip
on the first sample.
Evaluating the output of cmdlist()
can be very helpful for designing
and debugging *.cwl
files of new command-line software or changing the
parameter settings of existing ones.
cmdlist(sal, step = c(2,3), targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c results/setosa.csv > results/SE.csv.gz"
##
##
##
## $gunzip
## $gunzip$SE
## $gunzip$SE$gunzip
## [1] "gunzip -c ./results/gzip/SE.csv.gz > results/SE.csv"
statusWF(sal)
## $export_iris
## DataFrame with 1 row and 2 columns
## Step Status
## <character> <character>
## 1 export_iris Success
##
## $gzip
## DataFrame with 3 rows and 5 columns
## Targets Total_Files Existing_Files Missing_Files gzip
## <character> <numeric> <numeric> <numeric> <matrix>
## SE SE 1 1 0 Success
## VE VE 1 1 0 Success
## VI VI 1 1 0 Success
##
## $gunzip
## DataFrame with 3 rows and 5 columns
## Targets Total_Files Existing_Files Missing_Files gunzip
## <character> <numeric> <numeric> <numeric> <matrix>
## SE SE 1 1 0 Success
## VE VE 1 1 0 Success
## VI VI 1 1 0 Success
##
## $iris_stats
## DataFrame with 1 row and 2 columns
## Step Status
## <character> <character>
## 1 iris_stats Success
targetsWF(sal[2])
## $gzip
## DataFrame with 3 rows and 2 columns
## FileName SampleName
## <character> <character>
## SE results/setosa.csv SE
## VE results/versicolor.csv VE
## VI results/virginica.csv VI
The outfiles
components of SYSargsList
define the expected outfiles files
for each step in the workflow, some of which are the input for the next workflow step.
outfiles(sal[2])
## $gzip
## DataFrame with 3 rows and 1 column
## gzip_file
## <character>
## 1 ./results/gzip/SE.cs..
## 2 ./results/gzip/VE.cs..
## 3 ./results/gzip/VI.cs..
dependency(sal)
## $export_iris
## [1] NA
##
## $gzip
## [1] "export_iris"
##
## $gunzip
## [1] "gzip"
##
## $iris_stats
## [1] "gzip"
Sample comparisons are defined in the header lines of the targets
file
starting with ‘# <CMP>
’. This information can be accessed as follows:
targetsheader(sal, step = "Quality")
stepName(sal)
## [1] "export_iris" "gzip" "gunzip" "iris_stats"
SampleName(sal, step = "gzip")
## [1] "SE" "VE" "VI"
SampleName(sal, step = "iris_stats")
## NULL
outfiles
or targets
column files:getColumn(sal, "outfiles", step = "gzip", column = "gzip_file")
## SE VE VI
## "./results/gzip/SE.csv.gz" "./results/gzip/VE.csv.gz" "./results/gzip/VI.csv.gz"
getColumn(sal, "targetsWF", step = "gzip", column = "FileName")
## SE VE VI
## "results/setosa.csv" "results/versicolor.csv" "results/virginica.csv"
LineWise
step:codeLine(sal, step = "export_iris")
## export_iris
## mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
viewEnvir(sal)
## <environment: 0x55abef9232d0>
## [1] "df" "plot" "stats"
copyEnvir(sal, list = c("plot"), new.env = globalenv(), silent = FALSE)
## <environment: 0x55abef9232d0>
## Copying to 'new.env':
## plot
*.yml
datayamlinput(sal, step = "gzip")
## $file
## $file$class
## [1] "File"
##
## $file$path
## [1] "_FILE_PATH_"
##
##
## $SampleName
## [1] "_SampleName_"
##
## $ext
## [1] "csv.gz"
##
## $results_path
## $results_path$class
## [1] "Directory"
##
## $results_path$path
## [1] "./results"
SYSargsList
class and its subsetting operator [
:sal[1]
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Success
##
sal[1:3]
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Success
## 2. gzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 2.1. gzip
## cmdlist: 3 | Success: 3
## 3. gunzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 3.1. gunzip
## cmdlist: 3 | Success: 3
##
sal[c(1,3)]
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Success
## 2. gunzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 2.1. gunzip
## cmdlist: 3 | Success: 3
##
SYSargsList
class and its subsetting by steps and input samples:sal_sub <- subset(sal, subset_steps = c( 2,3), input_targets = ("SE"), keep_steps = TRUE)
stepsWF(sal_sub)
## $export_iris
## Instance of 'LineWise'
## Code Chunk length: 1
##
## $gzip
## Instance of 'SYSargs2':
## Slot names/accessors:
## targets: 1 (SE...SE), targetsheader: 1 (lines)
## modules: 0
## wf: 1, clt: 1, yamlinput: 4 (inputs)
## input: 1, output: 1
## cmdlist: 1
## Sub Steps:
## 1. gzip (rendered: TRUE)
##
##
##
## $gunzip
## Instance of 'SYSargs2':
## Slot names/accessors:
## targets: 1 (SE...SE), targetsheader: 1 (lines)
## modules: 0
## wf: 1, clt: 1, yamlinput: 4 (inputs)
## input: 1, output: 1
## cmdlist: 1
## Sub Steps:
## 1. gunzip (rendered: TRUE)
##
##
##
## $iris_stats
## Instance of 'LineWise'
## Code Chunk length: 5
targetsWF(sal_sub)
## $export_iris
## DataFrame with 0 rows and 0 columns
##
## $gzip
## DataFrame with 1 row and 2 columns
## FileName SampleName
## <character> <character>
## SE results/setosa.csv SE
##
## $gunzip
## DataFrame with 1 row and 2 columns
## gzip_file SampleName
## <character> <character>
## SE ./results/gzip/SE.cs.. SE
##
## $iris_stats
## DataFrame with 0 rows and 0 columns
outfiles(sal_sub)
## $export_iris
## DataFrame with 0 rows and 0 columns
##
## $gzip
## DataFrame with 1 row and 1 column
## gzip_file
## <character>
## 1 ./results/gzip/SE.cs..
##
## $gunzip
## DataFrame with 1 row and 1 column
## gunzip_file
## <character>
## 1 ./results/gunzip/SE...
##
## $iris_stats
## DataFrame with 0 rows and 0 columns
SYSargsList
class and its operator +
sal[1] + sal[2] + sal[3]
input
parameter in the workflowsal_c <- sal
## check values
yamlinput(sal_c, step = "gzip")
## $file
## $file$class
## [1] "File"
##
## $file$path
## [1] "_FILE_PATH_"
##
##
## $SampleName
## [1] "_SampleName_"
##
## $ext
## [1] "csv.gz"
##
## $results_path
## $results_path$class
## [1] "Directory"
##
## $results_path$path
## [1] "./results"
## check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c results/setosa.csv > results/SE.csv.gz"
## Replace
yamlinput(sal_c, step = "gzip", paramName = "ext") <- "txt.gz"
## check NEW values
yamlinput(sal_c, step = "gzip")
## $file
## $file$class
## [1] "File"
##
## $file$path
## [1] "_FILE_PATH_"
##
##
## $SampleName
## [1] "_SampleName_"
##
## $ext
## [1] "txt.gz"
##
## $results_path
## $results_path$class
## [1] "Directory"
##
## $results_path$path
## [1] "./results"
## Check on command-line
cmdlist(sal_c, step = "gzip", targets = 1)
## $gzip
## $gzip$SE
## $gzip$SE$gzip
## [1] "gzip -c results/setosa.csv > results/SE.txt.gz"
appendCodeLine(sal_c, step = "export_iris", after = 1) <- "log_cal_100 <- log(100)"
codeLine(sal_c, step = "export_iris")
## export_iris
## mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
## log_cal_100 <- log(100)
replaceCodeLine(sal_c, step="export_iris", line = 2) <- LineWise(code={
log_cal_100 <- log(50)
})
codeLine(sal_c, step = 1)
## export_iris
## mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
## log_cal_100 <- log(50)
For more details about the LineWise
class, please see below.
renameStep(sal_c, step = 1) <- "newStep"
renameStep(sal_c, c(1, 2)) <- c("newStep2", "newIndex")
sal_c
## Instance of 'SYSargsList':
## WF Steps:
## 1. newStep2 --> Status: Success
## 2. newIndex --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 2.1. gzip
## cmdlist: 3 | Success: 3
## 3. gunzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 3.1. gunzip
## cmdlist: 3 | Success: 3
## 4. iris_stats --> Status: Success
##
names(outfiles(sal_c))
## [1] "newStep2" "newIndex" "gunzip" "iris_stats"
names(targetsWF(sal_c))
## [1] "newStep2" "newIndex" "gunzip" "iris_stats"
dependency(sal_c)
## $newStep2
## [1] NA
##
## $newIndex
## [1] "newStep2"
##
## $gunzip
## [1] "newIndex"
##
## $iris_stats
## [1] "newIndex"
sal_test <- sal[c(1,2)]
replaceStep(sal_test, step = 1, step_name = "gunzip" ) <- sal[3]
sal_test
Note: Please use this method with attention, because it can disrupt all the dependency graphs.
sal_test <- sal[-2]
sal_test
## Instance of 'SYSargsList':
## WF Steps:
## 1. export_iris --> Status: Success
## 2. gunzip --> Status: Success
## Total Files: 3 | Existing: 3 | Missing: 0
## 2.1. gunzip
## cmdlist: 3 | Success: 3
## 3. iris_stats --> Status: Success
##
This section will introduce how CWL describes command-line tools and the specification and terminology of each file. For complete documentation, please check the CommandLineTools documentation here and here for Workflows and the user guide here.
CWL command-line specifications are written in YAML format.
In CWL, files with the extension .cwl
define the parameters of a chosen
command-line step or workflow, while files with the extension .yml
define
the input variables of command-line steps.
CommandLineTool
CommandLineTool
by CWL definition is a standalone process, with no interaction
if other programs, execute a program, and produce output.
Let’s explore the *.cwl
file:
dir_path <- system.file("extdata/cwl", package = "systemPipeR")
cwl <- yaml::read_yaml(file.path(dir_path, "example/example.cwl"))
cwlVersion
component shows the CWL specification version used by the document.class
component shows this document describes a CommandLineTool.
Note that CWL has another class
, called Workflow
which represents a union of one
or more command-line tools together.cwl[1:2]
## $cwlVersion
## [1] "v1.0"
##
## $class
## [1] "CommandLineTool"
baseCommand
component provides the name of the software that we desire to execute.cwl[3]
## $baseCommand
## [1] "echo"
inputs
section provides the input information to run the tool. Important
components of this section are:
id
: each input has an id describing the input name;type
: describe the type of input value (string, int, long, float, double,
File, Directory or Any);inputBinding
: It is optional. This component indicates if the input
parameter should appear on the command-line. If this component is missing
when describing an input parameter, it will not appear in the command-line
but can be used to build the command-line.cwl[4]
## $inputs
## $inputs$message
## $inputs$message$type
## [1] "string"
##
## $inputs$message$inputBinding
## $inputs$message$inputBinding$position
## [1] 1
##
##
##
## $inputs$SampleName
## $inputs$SampleName$type
## [1] "string"
##
##
## $inputs$results_path
## $inputs$results_path$type
## [1] "Directory"
outputs
section should provide a list of the expected outputs after running the command-line tools. Important
components of this section are:
id
: each input has an id describing the output name;type
: describe the type of output value (string, int, long, float, double,
File, Directory, Any or stdout
);outputBinding
: This component defines how to set the outputs values. The glob
component will define the name of the output value.cwl[5]
## $outputs
## $outputs$string
## $outputs$string$type
## [1] "stdout"
stdout
: component to specify a filename
to capture standard output.
Note here we are using a syntax that takes advantage of the inputs section,
using results_path parameter and also the SampleName
to construct the output filename.
cwl[6]
## $stdout
## [1] "$(inputs.results_path.basename)/$(inputs.SampleName).txt"
Workflow
Workflow
class in CWL is defined by multiple process steps, where can have
interdependencies between the steps, and the output for one step can be used as
input in the further steps.
cwl.wf <- yaml::read_yaml(file.path(dir_path, "example/workflow_example.cwl"))
cwlVersion
component shows the CWL specification version used by the document.class
component shows this document describes a Workflow
.cwl.wf[1:2]
## $class
## [1] "Workflow"
##
## $cwlVersion
## [1] "v1.0"
inputs
section describes the inputs of the workflow.cwl.wf[3]
## $inputs
## $inputs$message
## [1] "string"
##
## $inputs$SampleName
## [1] "string"
##
## $inputs$results_path
## [1] "Directory"
outputs
section describes the outputs of the workflow.cwl.wf[4]
## $outputs
## $outputs$string
## $outputs$string$outputSource
## [1] "echo/string"
##
## $outputs$string$type
## [1] "stdout"
steps
section describes the steps of the workflow. In this simple example,
we demonstrate one step.cwl.wf[5]
## $steps
## $steps$echo
## $steps$echo$`in`
## $steps$echo$`in`$message
## [1] "message"
##
## $steps$echo$`in`$SampleName
## [1] "SampleName"
##
## $steps$echo$`in`$results_path
## [1] "results_path"
##
##
## $steps$echo$out
## [1] "[string]"
##
## $steps$echo$run
## [1] "example/example.cwl"
Next, let’s explore the .yml file, which provide the input parameter values for all the components we describe above.
For this simple example, we have three parameters defined:
yaml::read_yaml(file.path(dir_path, "example/example_single.yml"))
## $message
## [1] "Hello World!"
##
## $SampleName
## [1] "M1"
##
## $results_path
## $results_path$class
## [1] "Directory"
##
## $results_path$path
## [1] "./results"
Note that if we define an input component in the .cwl file, this value needs to be also defined here in the .yml file.
systemPipeR
This section will demonstrate how to connect CWL parameters files to create
workflows. In addition, we will show how the workflow can be easily scalable
with systemPipeR
.
SYSargsList
container stores all the information and instructions needed for processing
a set of input files with a single or many command-line steps within a workflow
(i.e. several components of the software or several independent software tools).
The SYSargsList
object is created and fully populated with the SYSargsList
construct
function.
Full documentation of SYSargsList
management instances can be found here
and here.
The following imports a .cwl
file (here example.cwl
) for running the echo Hello World!
example.
HW <- SYSargsList(wf_file = "example/workflow_example.cwl",
input_file = "example/example_single.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"))
HW
## Instance of 'SYSargsList':
## WF Steps:
## 1. Step_x --> Status: Pending
## Total Files: 1 | Existing: 0 | Missing: 1
## 1.1. echo
## cmdlist: 1 | Pending: 1
##
cmdlist(HW)
## $Step_x
## $Step_x$defaultid
## $Step_x$defaultid$echo
## [1] "echo Hello World! > results/M1.txt"
However, we are limited to run just one command-line or one sample in this example.
To scale the command-line over many samples, a simple solution offered by systemPipeR
is to provide a variable
for each of the parameters that we want to run with multiple samples.
Let’s explore the example:
yml <- yaml::read_yaml(file.path(dir_path, "example/example.yml"))
yml
## $message
## [1] "_STRING_"
##
## $SampleName
## [1] "_SAMPLE_"
##
## $results_path
## $results_path$class
## [1] "Directory"
##
## $results_path$path
## [1] "./results"
For the message
and SampleName
parameter, we are passing a variable connecting
with a third file called targets.
Now, let’s explore the targets
file structure:
targetspath <- system.file("extdata/cwl/example/targets_example.txt", package = "systemPipeR")
read.delim(targetspath, comment.char = "#")
## Message SampleName
## 1 Hello World! M1
## 2 Hello USA! M2
## 3 Hello Bioconductor! M3
The targets
file defines all input files or values and sample ids of an analysis workflow.
For this example, we have defined a string message for the echo
command-line tool,
in the first column that will be evaluated, and the second column is the
SampleName
id for each one of the messages.
Any number of additional columns can be added as needed.
Users should note here, the usage of targets
files is optional when using
systemPipeR's
new CWL interface. Since for organizing experimental variables targets
files are extremely useful and user-friendly. Thus, we encourage users to keep using them.
targets
file information?The constructor function creates an SYSargsList
S4 class object connecting three input files:
wf_file
argument);input_file
argument);targets
argument).As demonstrated above, the latter is optional for workflow steps lacking input files.
The connection between input variables (here defined by input_file
argument)
and the targets
file are defined under the inputvars
argument.
A named vector is required, where each element name needs to match with column
names in the targets
file, and the value must match the names of the .yml
variables. This is used to replace the CWL variable and construct all the command-line
for that particular step.
The variable pattern _XXXX_
is used to distinguish CWL variables that target
columns will replace. This pattern is recommended for consistency and easy identification
but not enforced.
The following imports a .cwl
file (same example demonstrated above) for running
the echo Hello World
example. However, now we are connecting the variable defined
on the .yml
file with the targets
file inputs.
HW_mul <- SYSargsList(step_name = "echo",
targets=targetspath,
wf_file="example/workflow_example.cwl", input_file="example/example.yml",
dir_path = dir_path,
inputvars = c(Message = "_STRING_", SampleName = "_SAMPLE_"))
HW_mul
## Instance of 'SYSargsList':
## WF Steps:
## 1. echo --> Status: Pending
## Total Files: 3 | Existing: 0 | Missing: 3
## 1.1. echo
## cmdlist: 3 | Pending: 3
##
cmdlist(HW_mul)
## $echo
## $echo$M1
## $echo$M1$echo
## [1] "echo Hello World! > results/M1.txt"
##
##
## $echo$M2
## $echo$M2$echo
## [1] "echo Hello USA! > results/M2.txt"
##
##
## $echo$M3
## $echo$M3$echo
## [1] "echo Hello Bioconductor! > results/M3.txt"
Users need to define the command-line in a pseudo-bash script format:
command <- "
hisat2 \
-S <F, out: ./results/M1A.sam> \
-x <F: ./data/tair10.fasta> \
-k <int: 1> \
-min-intronlen <int: 30> \
-max-intronlen <int: 3000> \
-threads <int: 4> \
-U <F: ./data/SRR446027_1.fastq.gz>
"
First line is the base command. Each line is an argument with its default value.
For argument lines (starting from the second line), any word before the first
space with leading -
or --
in each will be treated as a prefix, like -S
or
--min
. Any line without this first word will be treated as no prefix.
All defaults are placed inside <...>
.
First argument is the input argument type. F
for “File”, “int”, “string” are unchanged.
Optional: use the keyword out
followed the type with a ,
comma separation to
indicate if this argument is also an CWL output.
Then, use :
to separate keywords and default values, any non-space value after the :
will be treated as the default value.
If any argument has no default value, just a flag, like --verbose
, there is no need to add any <...>
createParam
FunctioncreateParam
function requires the string
as defined above as an input.
First of all, the function will print the three components of the cwl
file:
- BaseCommand
: Specifies the program to execute.
- Inputs
: Defines the input parameters of the process.
- Outputs
: Defines the parameters representing the output of the process.
The four component is the original command-line.
If in interactive mode, the function will verify that everything is correct and will ask you to proceed. Here, the user can answer “no” and provide more information at the string level. Another question is to save the param created here.
If running the workflow in non-interactive mode, the createParam
function will
consider “yes” and returning the container.
cmd <- createParam(command, writeParamFiles = FALSE)
## *****BaseCommand*****
## hisat2
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: ./data/SRR446027_1.fastq.gz
## *****Outputs*****
## output1:
## type: File
## value: ./results/M1A.sam
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
If the user chooses not to save the param
files on the above operation,
it can use the writeParamFiles
function.
writeParamFiles(cmd, overwrite = TRUE)
printParam(cmd, position = "baseCommand") ## Print a baseCommand section
## *****BaseCommand*****
## hisat2
printParam(cmd, position = "outputs")
## *****Outputs*****
## output1:
## type: File
## value: ./results/M1A.sam
printParam(cmd, position = "inputs", index = 1:2) ## Print by index
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
printParam(cmd, position = "inputs", index = -1:-2) ## Negative indexing printing to exclude certain indices in a position
## *****Inputs*****
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: ./data/SRR446027_1.fastq.gz
cmd2 <- subsetParam(cmd, position = "inputs", index = 1:2, trim = TRUE)
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta
cmdlist(cmd2)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta"
cmd2 <- subsetParam(cmd, position = "inputs", index = c("S", "x"), trim = TRUE)
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta
cmdlist(cmd2)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta"
cmd3 <- replaceParam(cmd, "base", index = 1, replace = list(baseCommand = "bwa"))
## Replacing baseCommand
## *****BaseCommand*****
## bwa
## *****Parsed raw command line*****
## bwa -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd3)
## $defaultid
## $defaultid$hisat2
## [1] "bwa -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
new_inputs <- new_inputs <- list(
"new_input1" = list(type = "File", preF="-b", yml ="myfile"),
"new_input2" = "-L <int: 4>"
)
cmd4 <- replaceParam(cmd, "inputs", index = 1:2, replace = new_inputs)
## Replacing inputs
## *****Inputs*****
## new_input1:
## type: File
## preF: -b
## yml: myfile
## new_input2:
## type: int
## preF: -L
## yml: 4
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: ./data/SRR446027_1.fastq.gz
## *****Parsed raw command line*****
## hisat2 -b myfile -L 4 -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd4)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -b myfile -L 4 -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
newIn <- new_inputs <- list(
"new_input1" = list(type = "File", preF="-b1", yml ="myfile1"),
"new_input2" = list(type = "File", preF="-b2", yml ="myfile2"),
"new_input3" = "-b3 <F: myfile3>"
)
cmd5 <- appendParam(cmd, "inputs", index = 1:2, append = new_inputs)
## Replacing inputs
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: ./data/SRR446027_1.fastq.gz
## new_input1:
## type: File
## preF: -b1
## yml: myfile1
## new_input2:
## type: File
## preF: -b2
## yml: myfile2
## new_input3:
## type: File
## preF: -b3
## yml: myfile3
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz -b1 myfile1 -b2 myfile2 -b3 myfile3
cmdlist(cmd5)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz -b1 myfile1 -b2 myfile2 -b3 myfile3"
cmd6 <- appendParam(cmd, "inputs", index = 1:2, after=0, append = new_inputs)
## Replacing inputs
## *****Inputs*****
## new_input1:
## type: File
## preF: -b1
## yml: myfile1
## new_input2:
## type: File
## preF: -b2
## yml: myfile2
## new_input3:
## type: File
## preF: -b3
## yml: myfile3
## S:
## type: File
## preF: -S
## yml: ./results/M1A.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: ./data/SRR446027_1.fastq.gz
## *****Parsed raw command line*****
## hisat2 -b1 myfile1 -b2 myfile2 -b3 myfile3 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
cmdlist(cmd6)
## $defaultid
## $defaultid$hisat2
## [1] "hisat2 -b1 myfile1 -b2 myfile2 -b3 myfile3 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
output
paramnew_outs <- list(
"sam_out" = "<F: $(inputs.results_path)/test.sam>"
)
cmd7 <- replaceParam(cmd, "outputs", index = 1, replace = new_outs)
## Replacing outputs
## *****Outputs*****
## sam_out:
## type: File
## value: $(inputs.results_path)/test.sam
## *****Parsed raw command line*****
## hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz
output(cmd7)
## $defaultid
## $defaultid$hisat2
## [1] "./results/test.sam"
cmd <- "
hisat2 \
-S <F, out: _SampleName_.sam> \
-x <F: ./data/tair10.fasta> \
-k <int: 1> \
-min-intronlen <int: 30> \
-max-intronlen <int: 3000> \
-threads <int: 4> \
-U <F: _FASTQ_PATH1_>
"
WF <- createParam(cmd, overwrite = TRUE, writeParamFiles = TRUE, confirm = TRUE)
## *****BaseCommand*****
## hisat2
## *****Inputs*****
## S:
## type: File
## preF: -S
## yml: _SampleName_.sam
## x:
## type: File
## preF: -x
## yml: ./data/tair10.fasta
## k:
## type: int
## preF: -k
## yml: 1
## min-intronlen:
## type: int
## preF: -min-intronlen
## yml: 30
## max-intronlen:
## type: int
## preF: -max-intronlen
## yml: 3000
## threads:
## type: int
## preF: -threads
## yml: 4
## U:
## type: File
## preF: -U
## yml: _FASTQ_PATH1_
## *****Outputs*****
## output1:
## type: File
## value: _SampleName_.sam
## *****Parsed raw command line*****
## hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U _FASTQ_PATH1_
## Written content of 'commandLine' to file:
## param/cwl/hisat2/hisat2.cwl
## Written content of 'commandLine' to file:
## param/cwl/hisat2/hisat2.yml
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
WF_test <- loadWorkflow(targets = targetspath, wf_file="hisat2.cwl",
input_file="hisat2.yml", dir_path = "param/cwl/hisat2/")
WF_test <- renderWF(WF_test, inputvars = c(FileName = "_FASTQ_PATH1_"))
WF_test
## Instance of 'SYSargs2':
## Slot names/accessors:
## targets: 18 (M1A...V12B), targetsheader: 4 (lines)
## modules: 1
## wf: 0, clt: 1, yamlinput: 9 (inputs)
## input: 18, output: 18
## cmdlist: 18
## Sub Steps:
## 1. hisat2 (rendered: TRUE)
cmdlist(WF_test)[1:2]
## $M1A
## $M1A$hisat2
## [1] "hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446027_1.fastq.gz"
##
##
## $M1B
## $M1B$hisat2
## [1] "hisat2 -S _SampleName_.sam -x ./data/tair10.fasta -k 1 -min-intronlen 30 -max-intronlen 3000 -threads 4 -U ./data/SRR446028_1.fastq.gz"
SYSargsList
steps are can be defined with two inner classes, SYSargs2
and
LineWise
. Next, more details on both classes.
SYSargs2
ClassSYSargs2
workflow control class, an S4 class, is a list-like container where
each instance stores all the input/output paths and parameter components
required for a particular data analysis step. SYSargs2
instances are
generated by two constructor functions, loadWF and renderWF, using as data
input targets or yaml files as well as two cwl parameter files (for
details see below).
In CWL, files with the extension .cwl
define the parameters of a chosen
command-line step or workflow, while files with the extension .yml
define
the input variables of command-line steps. Note, input variables provided by a
targets file can be passed on to a SYSargs2
instance via the inputvars
argument of the renderWF function.
The following imports a .cwl
file (here hisat2-mapping-se.cwl
) for
running the short read aligner HISAT2 (Kim, Langmead, and Salzberg 2015). For more details about the
file structure and how to design or customize our own software tools, please
check systemPipeR and CWL
pipeline.
The loadWF and renderWF functions render the proper command-line strings for each sample and software tool.
library(systemPipeR)
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
dir_path <- system.file("extdata/cwl", package = "systemPipeR")
WF <- loadWF(targets = targetspath, wf_file = "hisat2/hisat2-mapping-se.cwl",
input_file = "hisat2/hisat2-mapping-se.yml",
dir_path = dir_path)
WF <- renderWF(WF, inputvars = c(FileName = "_FASTQ_PATH1_",
SampleName = "_SampleName_"))
Several accessor methods are available that are named after the slot names of
the SYSargs2
object.
names(WF)
## [1] "targets" "targetsheader" "modules" "wf"
## [5] "clt" "yamlinput" "cmdlist" "input"
## [9] "output" "files" "inputvars" "cmdToCwl"
## [13] "status" "internal_outfiles"
Of particular interest is the cmdlist()
method. It constructs the system
commands for running command-line software as specified by a given .cwl
file
combined with the paths to the input samples (e.g. FASTQ files) provided by a
targets
file. The example below shows the cmdlist()
output for running
HISAT2 on the first SE read sample. Evaluating the output of cmdlist()
can
be very helpful for designing and debugging .cwl
files of new command-line
software or changing the parameter settings of existing ones.
cmdlist(WF)[1]
## $M1A
## $M1A$`hisat2-mapping-se`
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 --min-intronlen 30 --max-intronlen 3000 -U ./data/SRR446027_1.fastq.gz --threads 4"
The output components of SYSargs2
define the expected output files for each
step in the workflow; some of which are the input for the next workflow step,
here next SYSargs2
instance.
output(WF)[1]
## $M1A
## $M1A$`hisat2-mapping-se`
## [1] "./results/M1A.sam"
The targets components of SYSargs2
object can be accessed by the targets
method. Here, for single-end (SE) samples, the structure of the targets file is
defined by:
FileName
: specify the FASTQ files path;SampleName
: Unique IDs for each sample;Factor
: ID for each treatment or condition.targets(WF)[1]
## $M1A
## $M1A$FileName
## [1] "./data/SRR446027_1.fastq.gz"
##
## $M1A$SampleName
## [1] "M1A"
##
## $M1A$Factor
## [1] "M1"
##
## $M1A$SampleLong
## [1] "Mock.1h.A"
##
## $M1A$Experiment
## [1] 1
##
## $M1A$Date
## [1] "23-Mar-2012"
as(WF, "DataFrame")
## DataFrame with 18 rows and 6 columns
## FileName SampleName Factor SampleLong Experiment Date
## <character> <character> <character> <character> <character> <character>
## 1 ./data/SRR446027_1.f.. M1A M1 Mock.1h.A 1 23-Mar-2012
## 2 ./data/SRR446028_1.f.. M1B M1 Mock.1h.B 1 23-Mar-2012
## 3 ./data/SRR446029_1.f.. A1A A1 Avr.1h.A 1 23-Mar-2012
## 4 ./data/SRR446030_1.f.. A1B A1 Avr.1h.B 1 23-Mar-2012
## 5 ./data/SRR446031_1.f.. V1A V1 Vir.1h.A 1 23-Mar-2012
## ... ... ... ... ... ... ...
## 14 ./data/SRR446040_1.f.. M12B M12 Mock.12h.B 1 23-Mar-2012
## 15 ./data/SRR446041_1.f.. A12A A12 Avr.12h.A 1 23-Mar-2012
## 16 ./data/SRR446042_1.f.. A12B A12 Avr.12h.B 1 23-Mar-2012
## 17 ./data/SRR446043_1.f.. V12A V12 Vir.12h.A 1 23-Mar-2012
## 18 ./data/SRR446044_1.f.. V12B V12 Vir.12h.B 1 23-Mar-2012
Please note, to work with custom data, users need to generate a targets
file
containing the paths to their own FASTQ files and then provide under
targetspath
the path to the corresponding targets
file.
In addition, if the Environment Modules is available, it is possible to define which module should be loaded, as shown here:
modules(WF)
## module1
## "hisat2/2.1.0"
Additional information can be accessed, as the parameters files location and the
inputvars
provided to generate the object.
files(WF)
inputvars(WF)
LineWise
was designed to store all the R code chunk when an RMarkdown file is
imported as a workflow.
rmd <- system.file("extdata", "spr_simple_lw.Rmd", package = "systemPipeR")
sal_lw <- SPRproject(overwrite = TRUE)
## Recreating directory '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject'
## Creating file '/tmp/Rtmpu2iIt2/Rbuild3b77422bb7655/systemPipeR/vignettes/.SPRproject/SYSargsList.yml'
sal_lw <- importWF(sal_lw, rmd, verbose = FALSE)
## Now check if required tools are installed
## There is no commandline (SYSargs) step in this workflow, skip.
codeLine(sal_lw)
## firstStep
## mapply(function(x, y) write.csv(x, y), split(iris, factor(iris$Species)), file.path("results", paste0(names(split(iris, factor(iris$Species))), ".csv")))
## secondStep
## setosa <- read.delim("results/setosa.csv", sep = ",")
## versicolor <- read.delim("results/versicolor.csv", sep = ",")
## virginica <- read.delim("results/virginica.csv", sep = ",")
lw <- stepsWF(sal_lw)[[2]]
## Coerce
ll <- as(lw, "list")
class(ll)
## [1] "list"
lw <- as(ll, "LineWise")
lw
## Instance of 'LineWise'
## Code Chunk length: 3
length(lw)
## [1] 3
names(lw)
## [1] "codeLine" "codeChunkStart" "stepName" "dependency" "status"
## [6] "files" "runInfo"
codeLine(lw)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## versicolor <- read.delim("results/versicolor.csv", sep = ",")
## virginica <- read.delim("results/virginica.csv", sep = ",")
codeChunkStart(lw)
## integer(0)
rmdPath(lw)
## character(0)
l <- lw[2]
codeLine(l)
## versicolor <- read.delim("results/versicolor.csv", sep = ",")
l_sub <- lw[-2]
codeLine(l_sub)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## virginica <- read.delim("results/virginica.csv", sep = ",")
replaceCodeLine(lw, line = 2) <- "5+5"
codeLine(lw)
## setosa <- read.delim("results/setosa.csv", sep = ",")
## 5 + 5
## virginica <- read.delim("results/virginica.csv", sep = ",")
appendCodeLine(lw, after = 0) <- "6+7"
codeLine(lw)
## 6 + 7
## setosa <- read.delim("results/setosa.csv", sep = ",")
## 5 + 5
## virginica <- read.delim("results/virginica.csv", sep = ",")
SYSargsList
replaceCodeLine(sal_lw, step = 2, line = 2) <- LineWise(code={
"5+5"
})
codeLine(sal_lw, step = 2)
appendCodeLine(sal_lw, step = 2) <- "66+55"
codeLine(sal_lw, step = 2)
appendCodeLine(sal_lw, step = 1, after = 1) <- "66+55"
codeLine(sal_lw, step = 1)
SYSargs
: Previous versionInstances of this S4 object class are constructed by the systemArgs
function
from two simple tabular files: a targets
file and a param
file. The
latter is optional for workflow steps lacking command-line software. Typically,
a SYSargs
instance stores all sample-level inputs as well as the paths to
the corresponding outputs generated by command-line- or R-based software
generating sample-level output files, such as read preprocessors
(trimmed/filtered FASTQ files), aligners (SAM/BAM files), variant callers
(VCF/BCF files) or peak callers (BED/WIG files). Each sample level input/output
operation uses its own SYSargs
instance. The outpaths of SYSargs
usually
define the sample inputs for the next SYSargs
instance. This connectivity is
established by writing the outpaths with the writeTargetsout
function to a
new targets
file that serves as input to the next systemArgs
call.
Typically, the user has to provide only the initial targets
file. All
downstream targets
files are generated automatically. By chaining several
SYSargs
steps together one can construct complex workflows involving many
sample-level input/output file operations with any combination of command-line
or R-based software.
Current, systemPipeR provides the param
file templates for third-party
software tools. Please check the listed software tools.
Tool Name | Description | Step |
---|---|---|
<a href=’http://bio-bwa.sourceforge.net/bwa.shtml'>bwa</a>; | BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml'>Bowtie2</a>; | Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’http://hannonlab.cshl.edu/fastx_toolkit/commandline.html'>FASTX-Toolkit</a>; | FASTX-Toolkit is a collection of command line tools for Short-Reads FASTA/FASTQ files preprocessing. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(236, 119, 112, 255) !important;" >Read Preprocessing</span> |
<a href=’http://hibberdlab.com/transrate/'>TransRate</a>; | Transrate is software for de-novo transcriptome assembly quality analysis. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(217, 133, 118, 255) !important;" >Quality</span> |
<a href=’http://research-pub.gene.com/gmap/'>Gsnap</a>; | GSNAP is a genomic short-read nucleotide alignment program. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’http://www.htslib.org/doc/samtools-1.2.html'>Samtools</a>; | Samtools is a suite of programs for interacting with high-throughput sequencing data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(208, 140, 121, 255) !important;" >Post-processing</span> |
<a href=’http://www.usadellab.org/cms/?page=trimmomatic'>Trimmomatic</a>; | Trimmomatic is a flexible read trimming tool for Illumina NGS data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(236, 119, 112, 255) !important;" >Read Preprocessing</span> |
<a href=’https://bioconductor.org/packages/release/bioc/vignettes/Rsubread/inst/doc/SubreadUsersGuide.pdf'>Rsubread</a>; | Rsubread is a Bioconductor software package that provides high-performance alignment and read counting functions for RNA-seq reads. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’https://broadinstitute.github.io/picard/'>Picard</a>; | Picard is a set of command line tools for manipulating high-throughput sequencing (HTS) data and formats such as SAM/BAM/CRAM and VCF. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(180, 160, 130, 255) !important;" >Manipulating HTS data</span> |
<a href=’https://busco.ezlab.org/'>Busco</a>; | BUSCO assesses genome assembly and annotation completeness with Benchmarking Universal Single-Copy Orthologs. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(217, 133, 118, 255) !important;" >Quality</span> |
<a href=’https://ccb.jhu.edu/software/hisat2/manual.shtml'>Hisat2</a>; | HISAT2 is a fast and sensitive alignment program for mapping NGS reads (both DNA and RNA) to reference genomes. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’https://ccb.jhu.edu/software/tophat/manual.shtml'>Tophat2</a>; | TopHat is a fast splice junction mapper for RNA-Seq reads. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’https://gatk.broadinstitute.org/hc/en-us'>GATK</a>; | Variant Discovery in High-Throughput Sequencing Data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(255, 106, 106, 255) !important;" >Variant Discovery</span> |
<a href=’https://github.com/FelixKrueger/TrimGalore'>Trim_galore</a>; | Trim Galore is a wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(236, 119, 112, 255) !important;" >Read Preprocessing</span> |
<a href=’https://github.com/TransDecoder/TransDecoder/wiki'>TransDecoder</a>; | TransDecoder identifies candidate coding regions within transcript sequences. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(171, 167, 133, 255) !important;" >Find Coding Regions</span> |
<a href=’https://github.com/Trinotate/Trinotate.github.io/wiki'>Trinotate</a>; | Trinotate is a comprehensive annotation suite designed for automatic functional annotation of transcriptomes. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(245, 112, 109, 255) !important;" >Transcriptome Functional Annotation</span> |
<a href=’https://github.com/alexdobin/STAR'>STAR</a>; | STAR is an ultrafast universal RNA-seq aligner. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(143, 188, 143, 255) !important;" >Alignment</span> |
<a href=’https://github.com/trinityrnaseq/trinityrnaseq/wiki'>Trinity</a>; | Trinity assembles transcript sequences from Illumina RNA-Seq data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(161, 174, 136, 255) !important;" >denovo Transcriptome Assembly</span> |
<a href=’https://macs3-project.github.io/MACS/'>MACS2</a>; | MACS2 identifies transcription factor binding sites in ChIP-seq data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(199, 147, 124, 255) !important;" >Peak calling</span> |
<a href=’https://pachterlab.github.io/kallisto/manual'>Kallisto</a>; | kallisto is a program for quantifying abundances of transcripts from RNA-Seq data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(227, 126, 115, 255) !important;" >Read counting</span> |
<a href=’https://samtools.github.io/bcftools/howtos/index.html'>BCFtools</a>; | BCFtools is a program for variant calling and manipulating files in the Variant Call Format (VCF) and its binary counterpart BCF. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(255, 106, 106, 255) !important;" >Variant Discovery</span> |
<a href=’https://www.bioinformatics.babraham.ac.uk/projects/bismark/'>Bismark</a>; | Bismark is a program to map bisulfite treated sequencing reads to a genome of interest and perform methylation calls in a single step. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(152, 181, 139, 255) !important;" >Bisulfite mapping</span> |
<a href=’https://www.bioinformatics.babraham.ac.uk/projects/fastqc/'>Fastqc</a>; | FastQC is a quality control tool for high throughput sequence data. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(217, 133, 118, 255) !important;" >Quality</span> |
<a href=’https://www.ncbi.nlm.nih.gov/books/NBK279690/'>Blast</a>; | BLAST finds regions of similarity between biological sequences. | <span style=" font-weight: bold; color: white !important;border-radius: 4px; padding-right: 4px; padding-left: 4px; background-color: rgba(189, 153, 127, 255) !important;" >Blast</span> |
Remember, if you desire to run any of these tools, make sure to have the
respective software installed on your system and configure in the PATH
.
There are a few ways to check if the required tools/modules are installed. The easiest way is automatically performed for
users by calling the importWF
function. At the end of the import, all required tools and modules are automatically
listed and checked for users.
There are a few other methods that one could use to perform the tool validation, please read details on our website, the Before running section.
To create a Workflow within systemPipeR
, we can start by defining an empty
container and checking the directory structure:
sal <- SPRproject()
The systemPipeR
package needs to be loaded (H Backman and Girke 2016).
appendStep(sal) <- LineWise({
library(systemPipeR)
},
step_name = "load_SPR")
preprocessReads
functionThe function preprocessReads
allows to apply predefined or custom
read preprocessing functions to all FASTQ files referenced in a
SYSargsList
container, such as quality filtering or adapter trimming
routines. Internally, preprocessReads
uses the FastqStreamer
function from
the ShortRead
package to stream through large FASTQ files in a
memory-efficient manner. The following example performs adapter trimming with
the trimLRPatterns
function from the Biostrings
package.
Here, we are appending this step to the SYSargsList
object created previously.
All the parameters are defined on the preprocessReads/preprocessReads-pe.yml
file.
appendStep(sal) <- SYSargsList(
step_name = "preprocessing",
targets = "targetsPE.txt", dir = TRUE,
wf_file = "preprocessReads/preprocessReads-pe.cwl",
input_file = "preprocessReads/preprocessReads-pe.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars = c(
FileName1 = "_FASTQ_PATH1_",
FileName2 = "_FASTQ_PATH2_",
SampleName = "_SampleName_"
),
dependency = c("load_SPR"))
After the preprocessing step, the outfiles
files can be used to generate the new
targets files containing the paths to the trimmed FASTQ files. The new targets
information can be used for the next workflow step instance, e.g. running the
NGS alignments with the trimmed FASTQ files. The appendStep
function is
automatically handling this connectivity between steps. Please check the next
step for more details.
The following example shows how one can design a custom read ‘preprocessReads’
function using utilities provided by the ShortRead
package, and then run it
in batch mode with the ‘preprocessReads’ function. Here, it is possible to
replace the function used on the preprocessing
step and modify the sal
object.
Because it is a custom function, it is necessary to save the part in the R object,
and internally the preprocessReads.doc.R
is loading the custom function.
If the R object is saved with a different name (here "param/customFCT.RData"
),
please replace that accordingly in the preprocessReads.doc.R
.
Please, note that this step is not added to the workflow, here just for demonstration.
First, we defined the custom function in the workflow:
appendStep(sal) <- LineWise(
code = {
filterFct <- function(fq, cutoff = 20, Nexceptions = 0) {
qcount <- rowSums(as(quality(fq), "matrix") <= cutoff, na.rm = TRUE)
# Retains reads where Phred scores are >= cutoff with N exceptions
fq[qcount <= Nexceptions]
}
save(list = ls(), file = "param/customFCT.RData")
},
step_name = "custom_preprocessing_function",
dependency = "preprocessing"
)
After, we can edit the input parameter:
yamlinput(sal, "preprocessing")$Fct
yamlinput(sal, "preprocessing", "Fct") <- "'filterFct(fq, cutoff=20, Nexceptions=0)'"
yamlinput(sal, "preprocessing")$Fct ## check the new function
cmdlist(sal, "preprocessing", targets = 1) ## check if the command line was updated with success
TrimGalore! is a wrapper tool to consistently apply quality and adapter trimming to fastq files, with some extra functionality for removing Reduced Representation Bisulfite-Seq (RRBS) libraries.
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "trimGalore",
targets = targetspath, dir = TRUE,
wf_file = "trim_galore/trim_galore-se.cwl",
input_file = "trim_galore/trim_galore-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"),
dependency = "load_SPR",
run_step = "optional")
Trimmomatic software (Bolger, Lohse, and Usadel 2014) performs a variety of useful trimming tasks for Illumina paired-end and single ended data. Here, an example of how to perform this task using parameters template files for trimming FASTQ files.
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
appendStep(sal) <- SYSargsList(step_name = "trimmomatic",
targets = targetspath, dir = TRUE,
wf_file = "trimmomatic/trimmomatic-se.cwl",
input_file = "trimmomatic/trimmomatic-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars = c(FileName = "_FASTQ_PATH1_", SampleName = "_SampleName_"),
dependency = "load_SPR",
run_step = "optional")
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 database. The following sessions present some NGS sequence alignment software. Select the most accurate aligner and determining the optimal parameter for your custom data set project.
For all the following examples, it is necessary to install the respective software
and export the PATH
accordingly.
HISAT2
The following steps will demonstrate how to use the short read aligner Hisat2
(Kim, Langmead, and Salzberg 2015) in both interactive job submissions and batch submissions to
queuing systems of clusters using the systemPipeR's
new CWL command-line interface.
Hisat2
index.appendStep(sal) <- SYSargsList(step_name = "hisat_index",
targets = NULL, dir = FALSE,
wf_file = "hisat2/hisat2-index.cwl",
input_file = "hisat2/hisat2-index.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars = NULL,
dependency = "preprocessing")
The parameter settings of the aligner are defined in the workflow_hisat2-se.cwl
and workflow_hisat2-se.yml
files. The following shows how to construct the
corresponding SYSargsList object, and append to sal workflow.
HISAT2
and SAMtools
It possible to build an workflow with HISAT2
and SAMtools
.
appendStep(sal) <- SYSargsList(step_name = "hisat_mapping",
targets = "preprocessing", dir = TRUE,
wf_file = "workflow-hisat2/workflow_hisat2-se.cwl",
input_file = "workflow-hisat2/workflow_hisat2-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"),
dependency = c("hisat_index"),
run_session = "remote")
Tophat2
The NGS reads of this project can also be aligned against the reference genome
sequence using Bowtie2/TopHat2
(Kim et al. 2013; Langmead and Salzberg 2012).
Bowtie2
index.appendStep(sal) <- SYSargsList(step_name = "bowtie_index",
targets = NULL, dir = FALSE,
wf_file = "bowtie2/bowtie2-index.cwl",
input_file = "bowtie2/bowtie2-index.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars = NULL,
dependency = "preprocessing",
run_step = "optional")
The parameter settings of the aligner are defined in the workflow_tophat2-mapping.cwl
and tophat2-mapping-pe.yml
files. The following shows how to construct the
corresponding SYSargsList object, using the outfiles from the preprocessing
step.
appendStep(sal) <- SYSargsList(step_name = "tophat2_mapping",
targets = "preprocessing", dir = TRUE,
wf_file = "tophat2/workflow_tophat2-mapping-se.cwl",
input_file = "tophat2/tophat2-mapping-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"),
dependency = c("bowtie_index"),
run_session = "remote",
run_step = "optional")
Bowtie2
(e.g. for miRNA profiling)The following example runs Bowtie2
as a single process without submitting it to a cluster.
appendStep(sal) <- SYSargsList(step_name = "bowtie2_mapping",
targets = "preprocessing", dir = TRUE,
wf_file = "bowtie2/workflow_bowtie2-mapping-se.cwl",
input_file = "bowtie2/bowtie2-mapping-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"),
dependency = c("bowtie_index"),
run_session = "remote",
run_step = "optional")
BWA-MEM
(e.g. for VAR-Seq)The following example runs BWA-MEM as a single process without submitting it to a cluster.
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")
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
(e.g. for RNA-Seq)The following example shows how one can use within the environment the R-based aligner , allowing running from R or command-line.
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")
appendStep(sal) <- SYSargsList(step_name = "rsubread",
targets = "preprocessing", dir = TRUE,
wf_file = "rsubread/rsubread-mapping-se.cwl",
input_file = "rsubread/rsubread-mapping-se.yml",
dir_path = system.file("extdata/cwl", package = "systemPipeR"),
inputvars=c(preprocessReads_se="_FASTQ_PATH1_", SampleName="_SampleName_"),
dependency = c("rsubread_index"),
run_session = "remote",
run_step = "optional")
gsnap
(e.g. for VAR-Seq and RNA-Seq)Another R-based short read aligner is gsnap
from the gmapR
package (Wu and Nacu 2010).
The code sample below introduces how to run this aligner on multiple nodes of a compute cluster.
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")
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 Http access needs to be available in the user account (e.g. home/publichtml
) of a system. If this is not the case then the BAM files need to be moved or copied to the system where IGV runs. In the following, htmldir
defines the path to the HTML directory with http access where the symbolic links to the BAM files will be stored. The corresponding URLs will be written to a text file specified under the _urlfile
_ argument.
appendStep(sal) <- LineWise(
code = {
bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles",
column = "samtools_sort_bam")
symLink2bam(
sysargs = bampaths, htmldir = c("~/.html/", "somedir/"),
urlbase = "http://cluster.hpcc.ucr.edu/~tgirke/",
urlfile = "./results/IGVurl.txt")
},
step_name = "bam_IGV",
dependency = "hisat2_mapping",
run_step = "optional"
)
Create txdb
(needs to be done only once).
appendStep(sal) <- LineWise(code = {
library(txdbmaker)
txdb <- makeTxDbFromGFF(file="data/tair10.gff", format="gff",
dataSource="TAIR", organism="Arabidopsis thaliana")
saveDb(txdb, file="./data/tair10.sqlite")
},
step_name = "create_txdb",
dependency = "hisat_mapping")
The following performs read counting with summarizeOverlaps
in parallel mode with multiple cores.
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.
Generate a table of read and alignment counts for all samples.
appendStep(sal) <- LineWise({
read_statsDF <- alignStats(args)
write.table(read_statsDF, "results/alignStats.xls",
row.names = FALSE, quote = FALSE, sep = "\t")
},
step_name = "align_stats",
dependency = "hisat_mapping")
The following shows the first four lines of the sample alignment stats file
provided by the systemPipeR
package. For simplicity the number of PE reads
is multiplied here by 2 to approximate proper alignment frequencies where each
read in a pair is counted.
read.table(system.file("extdata", "alignStats.xls", package = "systemPipeR"), header = TRUE)[1:4,]
## FileName Nreads2x Nalign Perc_Aligned Nalign_Primary Perc_Aligned_Primary
## 1 M1A 192918 177961 92.24697 177961 92.24697
## 2 M1B 197484 159378 80.70426 159378 80.70426
## 3 A1A 189870 176055 92.72397 176055 92.72397
## 4 A1B 188854 147768 78.24457 147768 78.24457
Download miRNA
genes from miRBase
.
appendStep(sal) <- LineWise({
system("wget https://www.mirbase.org/ftp/CURRENT/genomes/ath.gff3 -P ./data/")
gff <- rtracklayer::import.gff("./data/ath.gff3")
gff <- split(gff, elementMetadata(gff)$ID)
bams <- getColumn(sal, step = "bowtie2_mapping", 'outfiles', column = 2)
bfl <- BamFileList(bams, yieldSize=50000, index=character())
countDFmiR <- summarizeOverlaps(gff, bfl, mode="Union",
ignore.strand = FALSE, inter.feature = FALSE)
countDFmiR <- assays(countDFmiR)$counts
# Note: inter.feature=FALSE important since pre and mature miRNA ranges overlap
rpkmDFmiR <- apply(countDFmiR, 2, function(x) returnRPKM(counts = x, ranges = gff))
write.table(assays(countDFmiR)$counts, "results/countDFmiR.xls",
col.names=NA, quote=FALSE, sep="\t")
write.table(rpkmDFmiR, "results/rpkmDFmiR.xls", col.names=NA, quote=FALSE, sep="\t")
},
step_name = "read_counting_mirna",
dependency = "bowtie2_mapping")
The following computes the sample-wise Spearman correlation coefficients from the rlog
(regularized-logarithm) transformed expression values generated with the DESeq2
package. After transformation to a distance matrix, hierarchical clustering is performed with the hclust
function and the result is plotted as a dendrogram (sample_tree.pdf).
appendStep(sal) <- LineWise({
library(DESeq2, warn.conflicts=FALSE, quietly=TRUE)
library(ape, warn.conflicts=FALSE)
countDFpath <- system.file("extdata", "countDFeByg.xls", package="systemPipeR")
countDF <- as.matrix(read.table(countDFpath))
colData <- data.frame(row.names = targetsWF(sal)[[2]]$SampleName,
condition=targetsWF(sal)[[2]]$Factor)
dds <- DESeqDataSetFromMatrix(countData = countDF, colData = colData,
design = ~ condition)
d <- cor(assay(rlog(dds)), method = "spearman")
hc <- hclust(dist(1-d))
plot.phylo(as.phylo(hc), type = "p", edge.col = 4, edge.width = 3,
show.node.label = TRUE, no.margin = TRUE)
},
step_name = "sample_tree_rlog",
dependency = "read_counting")
rlog
values.
edgeR
The following run_edgeR
function is a convenience wrapper for
identifying differentially expressed genes (DEGs) in batch mode with
edgeR
’s GML method (Robinson, McCarthy, and Smyth 2010) for any number of
pairwise sample comparisons specified under the cmp
argument. Users
are strongly encouraged to consult the
edgeR
vignette
for more detailed information on this topic and how to properly run edgeR
on data sets with more complex experimental designs.
appendStep(sal) <- LineWise({
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment = "#")
cmp <- readComp(file = targetspath, format = "matrix", delim = "-")
countDFeBygpath <- system.file("extdata", "countDFeByg.xls", package = "systemPipeR")
countDFeByg <- read.delim(countDFeBygpath, row.names = 1)
edgeDF <- run_edgeR(countDF = countDFeByg, targets = targets, cmp = cmp[[1]],
independent = FALSE, mdsplot = "")
DEG_list <- filterDEGs(degDF = edgeDF, filter = c(Fold = 2, FDR = 10))
},
step_name = "edger",
dependency = "read_counting")
Filter and plot DEG results for up and down-regulated genes. Because of the small size of the toy data set used by this vignette, the FDR value has been set to a relatively high threshold (here 10%). More commonly used FDR cutoffs are 1% or 5%. The definition of ‘up’ and ‘down’ is given in the corresponding help file. To open it, type ?filterDEGs
in the R console.
edgeR
.
DESeq2
The following run_DESeq2
function is a convenience wrapper for
identifying DEGs in batch mode with DESeq2
(Love, Huber, and Anders 2014) for any number of
pairwise sample comparisons specified under the cmp
argument. Users
are strongly encouraged to consult the
DESeq2
vignette
for more detailed information on this topic and how to properly run DESeq2
on data sets with more complex experimental designs.
appendStep(sal) <- LineWise({
degseqDF <- run_DESeq2(countDF=countDFeByg, targets=targets, cmp=cmp[[1]],
independent=FALSE)
DEG_list2 <- filterDEGs(degDF=degseqDF, filter=c(Fold=2, FDR=10))
},
step_name = "deseq2",
dependency = "read_counting")
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 subsection.
appendStep(sal) <- LineWise({
library("biomaRt")
listMarts() # To choose BioMart database
listMarts(host="plants.ensembl.org")
m <- useMart("plants_mart", host="https://plants.ensembl.org")
listDatasets(m)
m <- useMart("plants_mart", dataset="athaliana_eg_gene", host="https://plants.ensembl.org")
listAttributes(m) # Choose data types you want to download
go <- getBM(attributes=c("go_id", "tair_locus", "namespace_1003"), mart=m)
go <- go[go[,3]!="",]; go[,3] <- as.character(go[,3])
go[go[,3]=="molecular_function", 3] <- "F"
go[go[,3]=="biological_process", 3] <- "P"
go[go[,3]=="cellular_component", 3] <- "C"
go[1:4,]
dir.create("./data/GO")
write.table(go, "data/GO/GOannotationsBiomart_mod.txt",
quote=FALSE, row.names=FALSE, col.names=FALSE, sep="\t")
catdb <- makeCATdb(myfile="data/GO/GOannotationsBiomart_mod.txt",
lib=NULL, org="", colno=c(1,2,3), idconv=NULL)
save(catdb, file="data/GO/catdb.RData")
},
step_name = "get_go_biomart",
dependency = "edger")
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. Plotting single sample sets is achieved by subsetting the input data frame as shown in the first line of the following example.
The following example performs hierarchical clustering on the rlog
transformed expression matrix subsetted by the DEGs identified in the
above differential expression analysis. It uses a Pearson correlation-based distance measure and complete linkage for cluster join.
appendStep(sal) <- LineWise({
library(pheatmap)
geneids <- unique(as.character(unlist(DEG_list[[1]])))
y <- assay(rlog(dds))[geneids, ]
pdf("heatmap1.pdf")
pheatmap(y, scale="row", clustering_distance_rows="correlation",
clustering_distance_cols="correlation")
dev.off()
},
step_name = "hierarchical_clustering",
dependency = c("sample_tree_rlog", "edgeR"))
systemPipeR workflows instances can be visualized with the plotWF
function.
This function will make a plot of selected workflow instance and the following information is displayed on the plot:
Success
, Error
, Pending
, Warnings
;If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, etc.
plotWF(sal, show_legend = TRUE, width = "80%")
To check more details of plotWF
, visit our website.
For running the workflow, runWF
function will execute all the steps store in
the workflow container. The execution will be on a single machine without
submitting to a queuing system of a computer cluster.
sal <- runWF(sal)
Alternatively, the computation can be greatly accelerated by processing many files in parallel using several compute nodes of a cluster, where a scheduling/queuing system is used for load balancing.
The resources
list object provides the number of independent parallel cluster
processes defined under the Njobs
element in the list. The following example
will run 18 processes in parallel using each 4 CPU cores.
If the resources available on a cluster allow running all 18 processes at the
same time, then the shown sample submission will utilize in a total of 72 CPU cores.
Note, runWF
can be used with most queueing systems as it is based on utilities
from the batchtools
package, which supports the use of template files (*.tmpl
)
for defining the run parameters of different schedulers. To run the following
code, one needs to have both a conffile
(see .batchtools.conf.R
samples here)
and a template
file (see *.tmpl
samples here)
for the queueing available on a system. The following example uses the sample
conffile
and template
files for the Slurm scheduler provided by this package.
The resources can be appended when the step is generated, or it is possible to
add these resources later, as the following example using the addResources
function:
resources <- list(conffile=".batchtools.conf.R",
template="batchtools.slurm.tmpl",
Njobs=18,
walltime=120, ## minutes
ntasks=1,
ncpus=4,
memory=1024, ## Mb
partition = "short"
)
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)
To check the summary of the workflow, we can use:
sal
statusWF(sal)
systemPipeR
compiles all the workflow execution logs in one central location,
making it easier to check any standard output (stdout
) or standard error
(stderr
) for any command-line tools used on the workflow or the R code stdout.
sal <- renderLogs(sal)
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_GB
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] magrittr_2.0.3 systemPipeR_2.11.3 ShortRead_1.63.0
## [4] GenomicAlignments_1.41.0 SummarizedExperiment_1.35.1 Biobase_2.65.0
## [7] MatrixGenerics_1.17.0 matrixStats_1.3.0 BiocParallel_1.39.0
## [10] Rsamtools_2.21.0 Biostrings_2.73.1 XVector_0.45.0
## [13] GenomicRanges_1.57.1 GenomeInfoDb_1.41.1 IRanges_2.39.1
## [16] S4Vectors_0.43.1 BiocGenerics_0.51.0 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 dplyr_1.1.4 farver_2.1.2
## [5] bitops_1.0-7 fastmap_1.2.0 digest_0.6.36 lifecycle_1.0.4
## [9] pwalign_1.1.0 compiler_4.4.1 rlang_1.1.4 sass_0.4.9
## [13] tools_4.4.1 utf8_1.2.4 yaml_2.3.9 knitr_1.48
## [17] S4Arrays_1.5.3 labeling_0.4.3 htmlwidgets_1.6.4 interp_1.1-6
## [21] DelayedArray_0.31.6 xml2_1.3.6 RColorBrewer_1.1-3 abind_1.4-5
## [25] withr_3.0.0 hwriter_1.3.2.1 grid_4.4.1 fansi_1.0.6
## [29] latticeExtra_0.6-30 colorspace_2.1-0 ggplot2_3.5.1 scales_1.3.0
## [33] tinytex_0.51 cli_3.6.3 rmarkdown_2.27 crayon_1.5.3
## [37] generics_0.1.3 rstudioapi_0.16.0 httr_1.4.7 cachem_1.1.0
## [41] stringr_1.5.1 zlibbioc_1.51.1 parallel_4.4.1 BiocManager_1.30.23
## [45] vctrs_0.6.5 Matrix_1.7-0 jsonlite_1.8.8 bookdown_0.40
## [49] systemfonts_1.1.0 jpeg_0.1-10 magick_2.8.3 crosstalk_1.2.1
## [53] jquerylib_0.1.4 glue_1.7.0 codetools_0.2-20 DT_0.33
## [57] stringi_1.8.4 gtable_0.3.5 deldir_2.0-4 UCSC.utils_1.1.0
## [61] munsell_0.5.1 tibble_3.2.1 pillar_1.9.0 htmltools_0.5.8.1
## [65] GenomeInfoDbData_1.2.12 R6_2.5.1 evaluate_0.24.0 kableExtra_1.4.0
## [69] lattice_0.22-6 highr_0.11 png_0.1-8 bslib_0.7.0
## [73] Rcpp_1.0.12 svglite_2.1.3 SparseArray_1.5.17 xfun_0.45
## [77] pkgconfig_2.0.3
This project is funded by NSF award ABI-1661152.
Amstutz, Peter, Michael R Crusoe, Nebojša Tijanić, Brad Chapman, John Chilton, Michael Heuer, Andrey Kartashov, et al. 2016. “Common Workflow Language, V1.0,” July. https://doi.org/10.6084/m9.figshare.3115156.v2.
Bolger, Anthony M, Marc Lohse, and Bjoern Usadel. 2014. “Trimmomatic: A Flexible Trimmer for Illumina Sequence Data.” Bioinformatics 30 (15): 2114–20.
Crusoe, Michael R, Sanne Abeln, Alexandru Iosup, Peter Amstutz, John Chilton, Nebojša Tijanić, Hervé Ménager, Stian Soiland-Reyes, Bogdan Gavrilovic, and Carole Goble. 2021. “Methods Included: Standardizing Computational Reuse and Portability with the Common Workflow Language,” May. http://arxiv.org/abs/2105.07028.
H Backman, Tyler W, and Thomas Girke. 2016. “systemPipeR: NGS workflow and report generation environment.” BMC Bioinformatics 17 (1): 388. https://doi.org/10.1186/s12859-016-1241-0.
Kim, Daehwan, Ben Langmead, and Steven L Salzberg. 2015. “HISAT: A Fast Spliced Aligner with Low Memory Requirements.” Nat. Methods 12 (4): 357–60.
Kim, Daehwan, Geo Pertea, Cole Trapnell, Harold Pimentel, Ryan Kelley, and Steven L Salzberg. 2013. “TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions.” Genome Biol. 14 (4): R36. https://doi.org/10.1186/gb-2013-14-4-r36.
Langmead, Ben, and Steven L Salzberg. 2012. “Fast Gapped-Read Alignment with Bowtie 2.” Nat. Methods 9 (4): 357–59. https://doi.org/10.1038/nmeth.1923.
Love, Michael, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biol. 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.
Robinson, M D, D J McCarthy, and G K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.
Wu, T D, and S Nacu. 2010. “Fast and SNP-tolerant Detection of Complex Variants and Splicing in Short Reads.” Bioinformatics 26 (7): 873–81. https://doi.org/10.1093/bioinformatics/btq057.