RcwlPipelines is a Bioconductor package that manages a collection of commonly used bioinformatics tools and pipeline based on Rcwl. These pre-built and pre-tested tools and pipelines are highly modularized with easy customization to meet different bioinformatics data analysis needs.

Rcwl and RcwlPipelines together forms a Bioconductor toolchain for use and development of reproducible bioinformatics pipelines in Common Workflow Language (CWL). The project also aims to develop a community-driven platform for open source, open development, and open review of best-practice CWL bioinformatics pipelines.

Installation

  1. Install the package from Bioconductor.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("RcwlPipelines")

The development version is also available to download from GitHub.

BiocManager::install("rworkflow/RcwlPipelines")
  1. Load the package into the R session.
library(RcwlPipelines)

Project resources

The project website https://rcwl.org/ serves as a central hub for all related resources. It provides guidance for new users and tutorials for both users and developers. Specific resources are listed below.

The R recipes and cwl scripts

The R scripts to build the CWL tools and pipelines are now residing in a dedicated GitHub repository, which is intended to be a community effort to collect and contribute Bioinformatics tools and pipelines using Rcwl and CWL.

Tutorial book

The tutorial book provides detailed instructions for developing Rcwl tools/pipelines, and also includes examples of some commonly-used tools and pipelines that covers a wide range of Bioinformatics data analysis needs.

RcwlPipelines core functions

Here we show the usage of 3 core functions: cwlUpdate, cwlSearch and cwlLoad for updating, searching, and loading the needed tools or pipelines in R.

cwlUpdate

The cwlUpdate function syncs the current Rcwl recipes and returns a cwlHub object which contains the most updated Rcwl recipes. The mcols() function returns all related information about each available tool or pipeline.

The recipes will be locally cached, so users don’t need to call cwlUpdate every time unless they want to use a tool/pipeline that is newly added to RcwlPipelines. Here we are using the recipes from Bioconductor devel version.

## For vignette use only. users don't need to do this step.
Sys.setenv(cachePath = tempdir()) 
atls <- cwlUpdate(branch = "dev") ## sync the tools/pipelines.
atls
#> cwlHub with 139 records
#> cache path:  /var/folders/7t/9l4kkf_j2sqbpn321y9g5558z96ck_/T//RtmpRbRoGf/Rcwl 
#> # last modified date:  2021-02-22 
#> # cwlSearch() to query scripts
#> # cwlLoad('title') to load the script
#> # additional mcols(): rid, rpath, Type, Container, mtime, ...
#> 
#>            title                      
#>   BFC1   | pl_alignMerge              
#>   BFC2   | pl_AnnPhaseVcf             
#>   BFC3   | pl_BaseRecal               
#>   BFC4   | pl_bwaAlign                
#>   BFC5   | pl_bwaMMRecal              
#>   ...      ...                        
#>   BFC135 | tl_VarScan2                
#>   BFC136 | tl_vcf_expression_annotator
#>   BFC137 | tl_vcf_readcount_annotator 
#>   BFC138 | tl_vep                     
#>   BFC139 | tl_vt_decompose            
#>          Command                                                            
#>   BFC1   bwaAlign+mergeBamDup                                               
#>   BFC2   VCFvep+dVCFcoverage+rVCFcoverage+VCFexpression+PhaseVcf            
#>   BFC3   BaseRecalibrator+ApplyBQSR+samtools_index+samtools_flagstat+samt...
#>   BFC4   bwa+sam2bam+sortBam+idxBam                                         
#>   BFC5   bwaAlign+mergeBamDup+BaseRecal                                     
#>   ...    ...                                                                
#>   BFC135                                                                    
#>   BFC136 vcf-expression-annotator                                           
#>   BFC137 vcf-readcount-annotator                                            
#>   BFC138 vep                                                                
#>   BFC139 vt decompose
table(mcols(atls)$Type)
#> 
#> pipeline     tool 
#>       26      113

Currently, we have integrated 113 command line tools and 26 pipelines.

cwlSearch

We can use (multiple) keywords to search for specific tools/pipelines of interest, which internally search the mcols of “rname”, “rpath”, “fpath”, “Command” and “Containers”. Here we show how to search the alignment tool bwa mem.

t1 <- cwlSearch(c("bwa", "mem"))
t1
#> cwlHub with 1 records
#> cache path:  /var/folders/7t/9l4kkf_j2sqbpn321y9g5558z96ck_/T//RtmpRbRoGf/Rcwl 
#> # last modified date:  2021-02-22 
#> # cwlSearch() to query scripts
#> # cwlLoad('title') to load the script
#> # additional mcols(): rid, rpath, Type, Container, mtime, ...
#> 
#>           title  Command
#>   BFC44 | tl_bwa bwa mem
mcols(t1)
#> DataFrame with 1 row and 14 columns
#>           rid       rname         create_time         access_time
#>   <character> <character>         <character>         <character>
#> 1       BFC44      tl_bwa 2021-02-25 04:02:11 2021-02-25 04:02:11
#>                    rpath       rtype                  fpath last_modified_time
#>              <character> <character>            <character>          <numeric>
#> 1 /var/folders/7t/9l4k..       local /var/folders/7t/9l4k..                 NA
#>          etag   expires        Type     Command              Container
#>   <character> <numeric> <character> <character>            <character>
#> 1          NA        NA        tool     bwa mem biocontainers/bwa:v0..
#>                 mtime
#>           <character>
#> 1 2021-02-22 14:45:16

cwlLoad

The last core function cwlLoad loads the Rcwl tool/pipeline into the R working environment. The code below loads the tool with a user-defined name bwa to do the read alignment.

bwa <- cwlLoad(title(t1)[1])  ## "tl_bwa"
bwa <- cwlLoad(mcols(t1)$fpath[1]) ## equivalent to the above. 
bwa
#> class: cwlProcess 
#>  cwlClass: CommandLineTool 
#>  cwlVersion: v1.0 
#>  baseCommand: bwa mem 
#> requirements:
#> - class: DockerRequirement
#>   dockerPull: biocontainers/bwa:v0.7.17-3-deb_cv1
#> inputs:
#>   threads (int): -t 
#>   RG (string): -R 
#>   Ref (File):  
#>   FQ1 (File):  
#>   FQ2 (File?):  
#> outputs:
#> sam:
#>   type: File
#>   outputBinding:
#>     glob: '*.sam'
#> stdout: bwaOutput.sam

Now the R tool of bwa is ready to use.

Customize a tool or pipeline

To fit users’ specific needs,the existing tool or pipline can be easily customized. Here we use the rnaseq_Sf pipeline to demonstrate how to access and change the arguments of a specific tool inside a pipeline. This pipeline covers RNA-seq reads quality summary by fastQC, alignment by STAR, quantification by featureCounts and quality control by RSeQC.

rnaseq_Sf <- cwlLoad("pl_rnaseq_Sf")
#> fastqc loaded
#> STAR loaded
#> sortBam loaded
#> samtools_index loaded
#> samtools_flagstat loaded
#> featureCounts loaded
#> gtfToGenePred loaded
#> genePredToBed loaded
#> read_distribution loaded
#> geneBody_coverage loaded
#> gCoverage loaded
#> STAR loaded
plotCWL(rnaseq_Sf)

There are many default arguments defined for the tool of STAR inside the pipeline. Users might want to change some of them. For example, we can change the value for --outFilterMismatchNmax argument from 2 to 5 for longer reads.

arguments(rnaseq_Sf, "STAR")[5:6]
#> [[1]]
#> [1] "--outFilterMismatchNmax"
#> 
#> [[2]]
#> [1] "2"
arguments(rnaseq_Sf, "STAR")[[6]] <- 5
arguments(rnaseq_Sf, "STAR")[5:6]
#> [[1]]
#> [1] "--outFilterMismatchNmax"
#> 
#> [[2]]
#> [1] "5"

We can also change the docker image for a specific tool (e.g., to a specific version). First, we search for all available docker images for STAR in biocontainers repository. The Source server could be quay or dockerhub.

searchContainer("STAR", repo = "biocontainers", source = "quay")
#> DataFrame with 29 rows and 5 columns
#>                  tool          repo        name          last_modified
#>           <character>   <character> <character>            <character>
#> 2.7.8a--0        STAR biocontainers   2.7.8a--0 Sun, 21 Feb 2021 12:..
#> 2.7.7a--0        STAR biocontainers   2.7.7a--0 Tue, 29 Dec 2020 13:..
#> 2.7.6a--0        STAR biocontainers   2.7.6a--0 Sun, 20 Sep 2020 09:..
#> 2.7.5c--0        STAR biocontainers   2.7.5c--0 Mon, 17 Aug 2020 09:..
#> 2.7.5b--0        STAR biocontainers   2.7.5b--0 Sat, 01 Aug 2020 17:..
#> ...               ...           ...         ...                    ...
#> 2.4.0j--0        STAR biocontainers   2.4.0j--0 Tue, 06 Mar 2018 12:..
#> 2.5.4a--0        STAR biocontainers   2.5.4a--0 Fri, 26 Jan 2018 21:..
#> 2.5.3a--0        STAR biocontainers   2.5.3a--0 Sat, 18 Mar 2017 11:..
#> 2.5.2b--0        STAR biocontainers   2.5.2b--0 Tue, 06 Sep 2016 07:..
#> 2.5.1b--0        STAR biocontainers   2.5.1b--0 Wed, 11 May 2016 08:..
#>                  size
#>           <character>
#> 2.7.8a--0     7985763
#> 2.7.7a--0     7824392
#> 2.7.6a--0     7820491
#> 2.7.5c--0     7806635
#> 2.7.5b--0     7805390
#> ...               ...
#> 2.4.0j--0     4734325
#> 2.5.4a--0     9225952
#> 2.5.3a--0     9119736
#> 2.5.2b--0     9086803
#> 2.5.1b--0    11291827

Then, we can change the STAR version into 2.7.8a (tag name: 2.7.8a–0).

requirements(rnaseq_Sf, "STAR")[[1]]
#> $class
#> [1] "DockerRequirement"
#> 
#> $dockerPull
#> [1] "quay.io/biocontainers/star:2.7.3a--0"
requirements(rnaseq_Sf, "STAR")[[1]] <- requireDocker(
    docker = "quay.io/biocontainers/star:2.7.8a--0")
requirements(rnaseq_Sf, "STAR")[[1]]
#> $class
#> [1] "DockerRequirement"
#> 
#> $dockerPull
#> [1] "quay.io/biocontainers/star:2.7.8a--0"

Run a tool or pipeline

Once the tool or pipeline is ready, we only need to assign values for each of the input parameters, and then submit using one of the functions: runCWL, runCWLBatch and cwlShiny. More detailed Usage and examples can be refer to the Rcwl vignette.

To successfully run the tool or pipeline, users either need to have all required command line tools pre-installed locally, or using the docker/singularity runtime by specifying docker = TRUE or docker = "singularity" argument inside runCWL or runCWLBatch function. Since the Bioconductor building machine doesn’t have all the tools installed, nor does it support the docker runtime, here we use some pseudo-code to demonstrate the tool/pipeline execution.

inputs(rnaseq_Sf)
rnaseq_Sf$in_seqfiles <- list("sample_R1.fq.gz",
                              "sample_R2.fq.gz")
rnaseq_Sf$in_prefix <- "sample"
rnaseq_Sf$in_genomeDir <- "genome_STAR_index_Dir"
rnaseq_Sf$in_GTFfile <- "GENCODE_version.gtf"

runCWL(rnaseq_Sf, outdir = "output/sample", docker = TRUE)

Users can also submit parallel jobs to HPC for multiple samples using runCWLBatch function. Different cluster job managers, such as “multicore”, “sge” and “slurm”, are supported using the BiocParallel::BatchtoolsParam.

library(BioParallel)
bpparam <- BatchtoolsParam(workers = 2, cluster = "sge",
                           template = batchtoolsTemplate("sge"))

inputList <- list(in_seqfiles = list(sample1 = list("sample1_R1.fq.gz",
                                                    "sample1_R2.fq.gz"),
                                     sample2 = list("sample2_R1.fq.gz",
                                                    "sample2_R2.fq.gz")),
                  in_prefix = list(sample1 = "sample1",
                                   sample2 = "sample2"))

paramList <- list(in_genomeDir = "genome_STAR_index_Dir",
                  in_GTFfile = "GENCODE_version.gtf",
                  in_runThreadN = 16)

runCWLBatch(rnaseq_Sf, outdir = "output",
            inputList, paramList,
            BPPARAM = bpparam)

SessionInfo

sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-apple-darwin13.4.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS/LAPACK: /Users/qi31566/miniconda3/envs/r-base/lib/libopenblasp-r0.3.12.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#> [1] RcwlPipelines_1.7.7  BiocFileCache_1.14.0 dbplyr_2.1.0        
#> [4] Rcwl_1.7.10          S4Vectors_0.28.1     BiocGenerics_0.36.0 
#> [7] yaml_2.2.1           BiocStyle_2.18.1    
#> 
#> loaded via a namespace (and not attached):
#>  [1] fs_1.5.0             bit64_4.0.5          filelock_1.0.2      
#>  [4] RColorBrewer_1.1-2   progress_1.2.2       httr_1.4.2          
#>  [7] rprojroot_2.0.2      tools_4.0.3          backports_1.2.1     
#> [10] bslib_0.2.4          utf8_1.1.4           R6_2.5.0            
#> [13] DBI_1.1.1            withr_2.4.1          tidyselect_1.1.0    
#> [16] prettyunits_1.1.1    bit_4.0.4            curl_4.3            
#> [19] compiler_4.0.3       git2r_0.28.0         textshaping_0.3.1   
#> [22] basilisk.utils_1.2.2 desc_1.2.0           bookdown_0.21       
#> [25] sass_0.3.1           checkmate_2.0.0      rappdirs_0.3.3      
#> [28] pkgdown_1.6.1        systemfonts_1.0.1    stringr_1.4.0       
#> [31] digest_0.6.27        rmarkdown_2.7        R.utils_2.10.1      
#> [34] basilisk_1.2.1       pkgconfig_2.0.3      htmltools_0.5.1.1   
#> [37] fastmap_1.1.0        htmlwidgets_1.5.3    rlang_0.4.10        
#> [40] rstudioapi_0.13      RSQLite_2.2.3        shiny_1.6.0         
#> [43] visNetwork_2.0.9     jquerylib_0.1.3      generics_0.1.0      
#> [46] jsonlite_1.7.2       BiocParallel_1.24.1  dplyr_1.0.4         
#> [49] R.oo_1.24.0          magrittr_2.0.1       Matrix_1.3-2        
#> [52] Rcpp_1.0.6           fansi_0.4.2          reticulate_1.18     
#> [55] lifecycle_1.0.0      R.methodsS3_1.8.1    stringi_1.5.3       
#> [58] grid_4.0.3           blob_1.2.1           promises_1.2.0.1    
#> [61] crayon_1.4.1         lattice_0.20-41      hms_1.0.0           
#> [64] batchtools_0.9.15    knitr_1.31           pillar_1.5.0        
#> [67] igraph_1.2.6         base64url_1.4        codetools_0.2-18    
#> [70] glue_1.4.2           evaluate_0.14        data.table_1.14.0   
#> [73] BiocManager_1.30.10  vctrs_0.3.6          httpuv_1.5.5        
#> [76] purrr_0.3.4          tidyr_1.1.2          assertthat_0.2.1    
#> [79] cachem_1.0.4         xfun_0.21            mime_0.10           
#> [82] xtable_1.8-4         later_1.1.0.1        ragg_1.1.0          
#> [85] tibble_3.0.6         memoise_2.0.0        DiagrammeR_1.0.6.1  
#> [88] ellipsis_0.3.1       brew_1.0-6