Input data: targets file

The targets file defines the input files (e.g. FASTQ or BAM) and sample comparisons used in a data analysis workflow. It can also store any number of additional descriptive information for each sample. The following shows the first four lines of the targets file used in this workflow template.

targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4, -c(5, 6)]

##                     FileName1                   FileName2
## 1 ./data/SRR446027_1.fastq.gz ./data/SRR446027_2.fastq.gz
## 2 ./data/SRR446028_1.fastq.gz ./data/SRR446028_2.fastq.gz
## 3 ./data/SRR446029_1.fastq.gz ./data/SRR446029_2.fastq.gz
## 4 ./data/SRR446030_1.fastq.gz ./data/SRR446030_2.fastq.gz
##   SampleName Factor        Date
## 1        M1A     M1 23-Mar-2012
## 2        M1B     M1 23-Mar-2012
## 3        A1A     A1 23-Mar-2012
## 4        A1B     A1 23-Mar-2012

To work with custom data, users need to generate a targets file containing the paths to their own FASTQ files. Here is a detailed description of the structure and utility of targets files.

With preprocessReads

The preprocessReads function allows applying predefined or custom read preprocessing functions to all FASTQ files referenced in a SAL 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 uses preprocessReads to perform adapter trimming with the trimLRPatterns function from the Biostrings package. In this instance, preprocessReads is invoked through a CL interface built on docopt, that is executed from R with CWL. The parameters for running preprocessReads are specified in the corresponding cwl/yml files. It is important to point out that creating and using CL interfaces for defining R-based workflow steps is not essential in systemPipeR since LineWise offers similar capabilities while requiring less specialized knowledge from users.

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

The paths to the output files generated by the preprocessing step (here trimmed FASTQ files) are recorded in a new targets file that can be used for the next workflow step, e.g. running the NGS alignments with the trimmed FASTQ files.

The following example demonstrates how to design a custom preprocessReads function, as well as how to replace parameters in the sal object. To apply the modifications to the workflow, it needs to be saved to a file, here param/customFCT.RData which will be loaded during the workflow run by the preprocessReads.doc.R script. Please note, this step is included here solely for demonstration purposes, and thus not part of the workflow run. This is achieved by dropping spr=TRUE in the header line of the code chunk.

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 defining this step, it can be inspected and modified as follows.

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

With Trimmomatic

For demonstration purposes, this workflow uses the Trimmomatic software as an example of an external CL read trimming tool (Bolger et al. 2014). Trimmomatic offers a range of practical trimming utilities specifically designed for single- and paired-end Illumina reads.

It is important to note that while the Trimmomatic trimming step is included in this workflow, it’s not mandatory. Users can opt to use read trimming results generated by the previous preprocessReads step if preferred.

appendStep(sal) <- SYSargsList(step_name = "trimming", targets = "targetsPE.txt",
    wf_file = "trimmomatic/trimmomatic-pe.cwl", input_file = "trimmomatic/trimmomatic-pe.yml",
    dir_path = system.file("extdata/cwl", package = "systemPipeR"),
    inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"), dependency = "load_SPR",
    run_step = "optional")

FASTQ quality report

The following seeFastq and seeFastqPlot functions generate and plot a series of useful quality statistics for a set of FASTQ files, including per cycle quality box plots, base proportions, base-level quality trends, relative k-mer diversity, length, and occurrence distribution of reads, number of reads above quality cutoffs and mean quality distribution. The results can be exported to different graphics formats, such as a PNG file, here named fastqReport.png. Detailed information about the usage and visual components in the quality plots can be found in the corresponding help file (see ?seeFastq or ?seeFastqPlot).

appendStep(sal) <- LineWise(code = {
    fastq <- getColumn(sal, step = "preprocessing", "targetsWF",
        column = 1)
    fqlist <- seeFastq(fastq = fastq, batchsize = 10000, klength = 8)
    png("./results/fastqReport.png", height = 162, width = 288 *
        length(fqlist))
    seeFastqPlot(fqlist)
    dev.off()
}, step_name = "fastq_report", dependency = "preprocessing")

Read mapping with HISAT2

To use the HISAT2 short read aligner developed by Kim, Langmead, and Salzberg (2015), it is necessary to index the reference genome. HISAT2 relies on the Burrows-Wheeler index for this process.

appendStep(sal) <- SYSargsList(step_name = "hisat2_index", dir = FALSE,
    targets = NULL, wf_file = "hisat2/hisat2-index.cwl", input_file = "hisat2/hisat2-index.yml",
    dir_path = "param/cwl", dependency = "load_SPR")

HISAT2 mapping

The parameter settings of the aligner are defined in the cwl/yml files used in the following code chunk. The following shows how to construct the alignment step and append it to the SAL workflow container. Please note that the input (FASTQ) files used in this step are the output files generated by the preprocessing step (see above: step_name = "preprocessing").

appendStep(sal) <- SYSargsList(step_name = "hisat2_mapping",
    dir = TRUE, targets = "preprocessing", wf_file = "workflow-hisat2/workflow_hisat2-pe.cwl",
    input_file = "workflow-hisat2/workflow_hisat2-pe.yml", dir_path = "param/cwl",
    inputvars = c(preprocessReads_1 = "_FASTQ_PATH1_", preprocessReads_2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"), rm_targets_col = c("FileName1",
        "FileName2"), dependency = c("preprocessing", "hisat2_index"))

The cmdlist functions allows to inspect the exact CL call used for each input file (sample), here for HISAT2 alignments. Note, this step also includes the conversion of the alignment files to sorted and indexed bam files using functionalities of the SAMtools CL suite.

cmdlist(sal, step = "hisat2_mapping", targets = 1)

$hisat2_mapping
$hisat2_mapping$M1A
$hisat2_mapping$M1A$hisat2
[1] "hisat2 -S ./results/M1A.sam  -x ./data/tair10.fasta  -k 1  --min-intronlen
30  --max-intronlen 3000  -1 ./results/M1A_1.fastq_trim.gz -2 ./results/M1A_2.fa
stq_trim.gz --threads 4"

$hisat2_mapping$M1A$`samtools-view`
[1] "samtools view -bS -o ./results/M1A.bam  ./results/M1A.sam "

$hisat2_mapping$M1A$`samtools-sort`
[1] "samtools sort -o ./results/M1A.sorted.bam  ./results/M1A.bam  -@ 4"

$hisat2_mapping$M1A$`samtools-index`
[1] "samtools index -b results/M1A.sorted.bam  results/M1A.sorted.bam.bai  ./res
ults/M1A.sorted.bam "

Alignment stats

The following computes an alignment summary file (here alignStats.xls), which comprises the count of reads in each FASTQ file and the number of reads that align with the reference, presented in both total and percentage values.

appendStep(sal) <- LineWise(code = {
    fqpaths <- getColumn(sal, step = "preprocessing", "targetsWF",
        column = "FileName1")
    bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles",
        column = "samtools_sort_bam")
    read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths,
        pairEnd = TRUE)
    write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE,
        quote = FALSE, sep = "\t")
}, step_name = "align_stats", dependency = "hisat2_mapping")

The resulting alignStats.xls file can be included in the report as shown below (here restricted to the first four rows).

read.table("results/alignStats.xls", header = TRUE)[1:4, ]

##   FileName Nreads2x Nalign Perc_Aligned Nalign_Primary
## 1      M1A   115994 109977     94.81266         109977
## 2      M1B   134480 112464     83.62879         112464
## 3      A1A   127976 122427     95.66403         122427
## 4      A1B   122486 101369     82.75966         101369
##   Perc_Aligned_Primary
## 1             94.81266
## 2             83.62879
## 3             95.66403
## 4             82.75966

Gene annotation database

For efficient handling of annotation ranges obtained from GFF or GTF files, they are organized within a TxDb object. Subsequently, the object is written to a SQLite database file. It is important to note that this process only needs to be performed once for a specific version of an annotation file.

appendStep(sal) <- LineWise(code = {
    library(GenomicFeatures)
    txdb <- suppressWarnings(makeTxDbFromGFF(file = "data/tair10.gff",
        format = "gff", dataSource = "TAIR", organism = "Arabidopsis thaliana"))
    saveDb(txdb, file = "./data/tair10.sqlite")
}, step_name = "create_db", dependency = "hisat2_mapping")

Read counting with summarizeOverlaps

The provided example employs non-strand-specific read counting while disregarding overlaps between different genes. As normalization the example uses reads per kilobase per million mapped reads (RPKM). The raw read count table (countDFeByg.xls) and the corresponding RPKM table (rpkmDFeByg.xls) are written to distinct files in the project’s results directory. Parallelization across multiple CPU cores is achieved with the BiocParallel package. When supplying a BamFileList as illustrated below, the summarizeOverlaps method defaults to employing bplapply and the register interface from BiocParallel. The MulticoreParam will utilize the number of cores returned by parallel::detectCores if the number of workers is left unspecified. For further information, refer to the help documentation by typing help("summarizeOverlaps").

appendStep(sal) <- LineWise(code = {
    library(GenomicFeatures)
    library(BiocParallel)
    txdb <- loadDb("./data/tair10.sqlite")
    outpaths <- getColumn(sal, step = "hisat2_mapping", "outfiles",
        column = "samtools_sort_bam")
    eByg <- exonsBy(txdb, by = c("gene"))
    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 = FALSE,
        singleEnd = FALSE, BPPARAM = multicoreParam))
    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")
    ## Creating a SummarizedExperiment object
    colData <- data.frame(row.names = SampleName(sal, "hisat2_mapping"),
        condition = getColumn(sal, "hisat2_mapping", position = "targetsWF",
            column = "Factor"))
    colData$condition <- factor(colData$condition)
    countDF_se <- SummarizedExperiment::SummarizedExperiment(assays = countDFeByg,
        colData = colData)
    ## Add results as SummarizedExperiment to the workflow
    ## object
    SE(sal, "read_counting") <- countDF_se
}, step_name = "read_counting", dependency = "create_db")

Importantly, when conducting statistical differential expression or abundance analysis using methods like edgeR or DESeq2, the raw count values are the expected input. RPKM values should be reserved for specialized applications, such as manually inspecting expression levels across different genes or features.

Shows first 10 rows of countDFeByg.xls table.

read.delim("results/countDFeByg.xls", row.names = 1, check.names = FALSE)[1:10,
    ]

##           M1A M1B  A1A A1B V1A V1B M6A M6B A6A A6B  V6A
## AT1G01010 286 260  364 181 568 300 255 135 514 318  757
## AT1G01020 104 136  139 131 174 156 148 131 114 104  206
## AT1G01030 120 109  167  59 136 192  74  26  23  73  118
## AT1G01040 911 727 1030 627 962 918 862 618 880 639 1632
## AT1G01046  23  12   17  13  16  26  19  14  23  21   24
## AT1G01050 189 178  247 184 226 380 524 619 382 414  622
## AT1G01060  98 262   86  88  33  32   8   4   6   3    2
## AT1G01070   0   1    5   0  10   5  11  13  29   8   28
## AT1G01073   0   0    0   0   0   0   0   0   0   0    0
## AT1G01080 377 390  363 454 476 630 437 747 266 350  352
##            V6B M12A M12B A12A A12B V12A V12B
## AT1G01010  551  198  248  527  417  650  671
## AT1G01020  212   67  156  130  120   80  158
## AT1G01030  214   45   51   31   48  177  442
## AT1G01040 1552  651 1095 1324  702  671  995
## AT1G01046   36   28   23   33   13   23   20
## AT1G01050  962  666 1355  737  532  635 1004
## AT1G01060   10  220  317  501  198  164  159
## AT1G01070   14   11   65   64   39   23   24
## AT1G01073    0    0    0    0    0    0    0
## AT1G01080  765  384 1037  343  299  267  373

Sample-wise clustering

The sample-wise Spearman correlation coefficients are calculated from the rlog transformed expression values (countDF_se) generated using the DESeq2 package. These values are then converted into a distance matrix, which is subsequently used for hierarchical clustering with the hclust function. The resulting dendrogram is then saved as a PNG file named sample_tree.png.

appendStep(sal) <- LineWise(code = {
    library(DESeq2, quietly = TRUE)
    library(ape, warn.conflicts = FALSE)
    ## Extracting SummarizedExperiment object
    se <- SE(sal, "read_counting")
    dds <- DESeqDataSet(se, design = ~condition)
    d <- cor(assay(rlog(dds)), method = "spearman")
    hc <- hclust(dist(1 - d))
    png("results/sample_tree.png")
    plot.phylo(as.phylo(hc), type = "p", edge.col = "blue", edge.width = 2,
        show.node.label = TRUE, no.margin = TRUE)
    dev.off()
}, step_name = "sample_tree", dependency = "read_counting")

Run edgeR

appendStep(sal) <- LineWise(code = {
    library(edgeR)
    countDF <- read.delim("results/countDFeByg.xls", row.names = 1,
        check.names = FALSE)
    cmp <- readComp(stepsWF(sal)[["hisat2_mapping"]], format = "matrix",
        delim = "-")
    edgeDF <- run_edgeR(countDF = countDF, targets = targetsWF(sal)[["hisat2_mapping"]],
        cmp = cmp[[1]], independent = FALSE, mdsplot = "")
    write.table(edgeDF, "./results/edgeRglm_allcomp.xls", quote = FALSE,
        sep = "\t", col.names = NA)
}, step_name = "run_edger", dependency = "read_counting")

Note, to call DEGs with DESeq2 instead of edgeR, users can simply replace in the above code ‘run_edgeR’ with ‘run_DESeq2’.

Add gene descriptions

This step is optional. It appends functional descriptions obtained from BioMart to the DEG table.

appendStep(sal) <- LineWise(code = {
    library("biomaRt")
    m <- useMart("plants_mart", dataset = "athaliana_eg_gene",
        host = "https://plants.ensembl.org")
    desc <- getBM(attributes = c("tair_locus", "description"),
        mart = m)
    desc <- desc[!duplicated(desc[, 1]), ]
    descv <- as.character(desc[, 2])
    names(descv) <- as.character(desc[, 1])
    edgeDF <- data.frame(edgeDF, Desc = descv[rownames(edgeDF)],
        check.names = FALSE)
    write.table(edgeDF, "./results/edgeRglm_allcomp.xls", quote = FALSE,
        sep = "\t", col.names = NA)
}, step_name = "custom_annot", dependency = "run_edger")

Plot DEG results

Filter and plot DEG results for up and down regulated genes. The definition of up and down is given in the corresponding help file. To open it, type ?filterDEGs in the R console.

Note, due to the small number of genes in the toy dataset, the FDR cutoff in this example is set to an unreasonably large value. With real data sets this cutoff should be set to a much smaller value (often 1%, 5% or 10%).

appendStep(sal) <- LineWise(code = {
    edgeDF <- read.delim("results/edgeRglm_allcomp.xls", row.names = 1,
        check.names = FALSE)
    png("results/DEGcounts.png")
    DEG_list <- filterDEGs(degDF = edgeDF, filter = c(Fold = 2,
        FDR = 20))
    dev.off()
    write.table(DEG_list$Summary, "./results/DEGcounts.xls",
        quote = FALSE, sep = "\t", row.names = FALSE)
}, step_name = "filter_degs", dependency = "custom_annot")

Venn diagrams of DEG sets

The overLapper function can compute Venn intersects for large numbers of sample sets (up to 20 or more) and plots 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).

The overLapper function can compute Venn intersects for large numbers of sample sets (up to 20 or more) and plots 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(code = {
    vennsetup <- overLapper(DEG_list$Up[6:9], type = "vennsets")
    vennsetdown <- overLapper(DEG_list$Down[6:9], type = "vennsets")
    png("results/vennplot.png")
    vennPlot(list(vennsetup, vennsetdown), mymain = "", mysub = "",
        colmode = 2, ccol = c("blue", "red"))
    dev.off()
}, step_name = "venn_diagram", dependency = "filter_degs")

Obtain gene-to-GO mappings

The following shows how to obtain gene-to-GO mappings from biomaRt (here for A. thaliana) and how to organize them for the downstream GO term enrichment analysis. Alternatively, the gene-to-GO mappings can be obtained for many organisms from Bioconductor’s *.db genome annotation packages or GO annotation files provided by various genome databases. For each annotation this relatively slow preprocessing step needs to be performed only once. Subsequently, the preprocessed data can be loaded with the load function as shown in the next subsection.

appendStep(sal) <- LineWise(code = {
    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, ]
    if (!dir.exists("./data/GO"))
        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_annot", dependency = "filter_degs")

Batch GO term enrichment analysis

Apply the enrichment analysis to the DEG sets obtained 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 a GO slim vector for a specific organism can be obtained from BioMart.

appendStep(sal) <- LineWise(code = {
    library("biomaRt")
    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)
    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)
    write.table(BatchResultslim, "results/GOBatchSlim.xls", row.names = FALSE,
        quote = FALSE, sep = "\t")
}, step_name = "go_enrich", dependency = "get_go_annot")

Plot batch GO term results

The data.frame generated by GOCluster 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.

appendStep(sal) <- LineWise(code = {
    gos <- BatchResultslim[grep("M6-V6_up_down", BatchResultslim$CLID),
        ]
    gos <- BatchResultslim
    png("results/GOslimbarplotMF.png", height = 8, width = 10)
    goBarplot(gos, gocat = "MF")
    goBarplot(gos, gocat = "BP")
    goBarplot(gos, gocat = "CC")
    dev.off()
}, step_name = "go_plot", dependency = "go_enrich")

Parallelization on clusters

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. 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, ntasks = 1, ncpus = 4, memory = 1024,
    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 .