-# Turning a bash script into a workflow using existing tools
+# Turning a shell script into a workflow using existing tools
In this lesson we will turn `rnaseq_analysis_on_input_file.sh` into a workflow.
For this analysis, the input parameters are the fastq file and the reference data required by STAR.
-In the original bash script, the following variables are declared:
+In the original shell script, the following variables are declared:
```
# initialize a variable with an intuitive name to store the name of the input fastq file
run: bio-cwl-tools/fastqc/fastqc_2.cwl
in:
reads_file: fq
- out: [html_file, summary_file]
+ out: [html_file]
```
5. Running alignment with STAR
qc_html:
type: File
outputSource: fastqc/html_file
- qc_summary:
- type: File
- outputSource: fastqc/summary_file
bam_sorted_indexed:
type: File
outputSource: samtools/bam_sorted_indexed
CWL input values are provided in the form of a YAML or JSON file.
Create one by right clicking on the explorer, select "New File" and
-create a called file "input.yaml".
+create a called file "main-input.yaml".
This file gives the values for parameters declared in the `inputs`
section of our workflow. Our workflow takes `fq`, `genome` and `gtf`
Type this into the terminal:
```
-cwl-runner main.cwl input.yaml
+cwl-runner main.cwl main-input.yaml
```
3. Debugging the workflow
STAR:
requirements:
ResourceRequirement:
- ramMin: 7000
+ ramMin: 8000
run: bio-cwl-tools/STAR/STAR-Align.cwl
```
"basename": "Mov10_oe_1.subset_fastqc.html",
"class": "File",
"size": 383589
- },
- "qc_summary": {
- "location": "file:///home/username/rnaseq-cwl-training-exercises/summary.txt",
- "basename": "summary.txt",
- "class": "File",
- "size": 590
}
}
```
-This has the same structure as `input.yaml`. The each output
+This has the same structure as `main-input.yaml`. The each output
parameter is listed, with the `location` field of each `File` object
indicating where the output file can be found.
It consumes some input parameters, runs a program, and produce output
values.
-Here's the original bash script
+Here is the original shell command:
```
featureCounts -T $cores -s 2 -a $gtf -o $counts $counts_input_bam
counts_input_bam: File
```
-4. The base command
+4. Specifying the program to run
-This one is easy. This is the name of program to run:
+Give the name of the program to run in `baseCommand`.
```
baseCommand: featureCounts
```
-5. The command arguments
+5. Command arguments
The easiest way to describe the command line is with an `arguments`
section. This takes a comma-separated list of command line arguments.
describes the resources allocated to running the program. Here we use
`$(runtime.cores)` to decide how many threads to request.
-File variables can also yield a partial filename, by adding
-`.nameroot`. This is the filename with the final dot-extension
-stripped off.
-
```
arguments: [-T, $(runtime.cores),
-a, $(inputs.gtf),
- -o, $(inputs.counts_input_bam.nameroot)_featurecounts.txt,
+ -o, featurecounts.tsv,
$(inputs.counts_input_bam)]
```
-6. The outputs section
+6. Outputs section
In CWL, you must explicitly identify the outputs of a program. This
-associates output parameters with specific files, and allows the
+associates output parameters with specific files, and enables the
workflow runner to know which files must be saved and which files can
be discarded.
In the previous section, we told the featureCounts program the name of
-our output files should be
-`$(inputs.counts_input_bam.nameroot)_featurecounts.txt`.
+our output files should be `featurecounts.tsv`.
We can declare an output parameter called `featurecounts` that will
have that output file as its value.
The `outputBinding` section describes how to determine the value of
the parameter. The `glob` field tells it to search for a file in the
-output directory with the
-`$(inputs.counts_input_bam.nameroot)_featurecounts.txt`
+output directory called `featurecounts.tsv`
```
outputs:
featurecounts:
type: File
outputBinding:
- glob: $(inputs.counts_input_bam.nameroot)_featurecounts.txt
+ glob: featurecounts.tsv
+```
+
+7. Running in a container
+
+In order to run the tool, it needs to be installed.
+Using software containers, a tool can be pre-installed into a
+compatible runtime environment, and that runtime environment (called a
+container image) can be downloaded and run on demand.
+
+Many bioinformatics tools are already available as containers. One
+resource is the BioContainers project. Let's find the "subread" software:
+
+ 1. Visit https://biocontainers.pro/
+ 2. Click on "Registry"
+ 3. Search for "subread"
+ 4. Click on the search result for "subread"
+ 5. Click on the tab "Packages and Containers"
+ 6. Choose a row with type "docker", then on the right side of the "Full
+Tag" column for that row, click the "copy to clipboard" button.
+
+To declare that you want to run inside a container, create a section
+called `hints` with a subsection `DockerRequirement`. Under
+`DockerRequirement`, paste the text your copied in the above step.
+Replace the text `docker pull` to `dockerPull:` and indent it so it is
+in the `DockerRequirement` section.
+
+```
+hints:
+ DockerRequirement:
+ dockerPull: quay.io/biocontainers/subread:1.5.0p3--0
+```
+
+8. Running a tool on its own
+
+When creating a tool wrapper, it is helpful to run it on its own to test it.
+
+The input to a single tool is the same kind of input parameters file
+that we used as input to a workflow in the previous lesson.
+
+featureCounts.yaml:
+
+```
+counts_input_bam:
+ class: File
+ location: Aligned.sortedByCoord.out.bam
+gtf:
+ class: File
+ location: rnaseq/reference_data/chr1-hg19_genes.gtf
```
-N.
+The invocation is also the same:
-The most portable way to run a tool is to wrap it in a Docker
-container. (Some CWL platforms, such as Arvados, require it). Many
-bioinformatics tools are already available as containers. One
-resource is the BioContainers project.
+```
+cwl-runner featureCounts.cwl featureCounts.yaml
+```
-Visit https://biocontainers.pro/
+9. Adding it to the workflow
-Click on "Registry"
+Now that we have confirmed that it works, we can add it to our workflow.
+We add it to `steps`, connecting the output of samtools to
+`counts_input_bam` and the `gtf` taking the workflow input of the same
+name.
-Search for "subread"
+```
+steps:
+ ...
+ featureCounts:
+ requirements:
+ ResourceRequirement:
+ ramMin: 500
+ run: featureCounts.cwl
+ in:
+ counts_input_bam: samtools/bam_sorted_indexed
+ gtf: gtf
+ out: [featurecounts]
+```
-Click on the search result for "subread"
+We will add the result from featurecounts to the output:
-Click on the tab "Packages and Containers"
+```
+outputs:
+ ...
+ featurecounts:
+ type: File
+ outputSource: featureCounts/featurecounts
+
+```
-Choose a row with type "docker", then click the "copy to clipboard"
-button on the right side of the"Full Tag" column for that row.
+You should now be able to re-run the workflow and it will run the
+"featureCounts" step and include "featurecounts" in the output.
--- /dev/null
+# Analyzing multiple samples
+
+Analyzing a single sample is great, but in the real world you probably
+have a batch of samples that you need to analyze and then compare.
+
+1. Subworkflows
+
+In addition to running command line tools, a workflow step can also
+execute another workflow.
+
+Let's copy "main.cwl" to "alignment.cwl".
+
+Now, edit open "main.cwl" for editing. We are going to replace the `steps` and `outputs` sections.
+
+```
+steps:
+ alignment:
+ run: alignment.cwl
+ in:
+ fq: fq
+ genome: genome
+ gtf: gtf
+ out: [qc_html, bam_sorted_indexed, featurecounts]
+```
+
+In the outputs section, all the output sources are from the alignment step:
+
+```
+outputs:
+ qc_html:
+ type: File
+ outputSource: alignment/qc_html
+ bam_sorted_indexed:
+ type: File
+ outputSource: alignment/bam_sorted_indexed
+ featurecounts:
+ type: File
+ outputSource: alignment/featurecounts
+```
+
+We also need a little boilerplate to tell the workflow runner that we want to use subworkflows:
+
+```
+requirements:
+ SubworkflowFeatureRequirement: {}
+```
+
+If you run this workflow, you will get exactly the same results as
+before, we've just wrapped the inner workflow with an outer workflow.
+
+2. Scattering
+
+The wrapper lets us do something useful. We can modify the outer
+workflow to accept a list of files, and then invoke the inner workflow
+step for every one of those files. We will need to modify the
+`inputs`, `steps`, `outputs`, and `requirements` sections.
+
+First we change the `fq` parameter to expect a list of files:
+
+```
+inputs:
+ fq: File[]
+ genome: Directory
+ gtf: File
+```
+
+Next, we add `scatter` to the alignment step. The means it will
+run `alignment.cwl` for each value in the list in the `fq` parameter.
+
+```
+steps:
+ alignment:
+ run: alignment.cwl
+ scatter: fq
+ in:
+ fq: fq
+ genome: genome
+ gtf: gtf
+ out: [qc_html, bam_sorted_indexed, featurecounts]
+```
+
+Because the scatter produces multiple outputs, each output parameter
+becomes a list as well:
+
+```
+outputs:
+ qc_html:
+ type: File[]
+ outputSource: alignment/qc_html
+ bam_sorted_indexed:
+ type: File[]
+ outputSource: alignment/bam_sorted_indexed
+ featurecounts:
+ type: File[]
+ outputSource: alignment/featurecounts
+```
+
+Finally, we need a little more boilerplate to tell the workflow runner
+that we want to use scatter:
+
+```
+requirements:
+ SubworkflowFeatureRequirement: {}
+ ScatterFeatureRequirement: {}
+```
+
+3. Running with list inputs
+
+The `fq` parameter needs to be a list. You write a list in yaml by
+starting each list item with a dash. Example `main-input.yaml`
+
+```
+fq:
+ - class: File
+ location: rnaseq/raw_fastq/Mov10_oe_1.subset.fq
+ format: http://edamontology.org/format_1930
+ - class: File
+ location: rnaseq/raw_fastq/Mov10_oe_2.subset.fq
+ format: http://edamontology.org/format_1930
+ - class: File
+ location: rnaseq/raw_fastq/Mov10_oe_3.subset.fq
+ format: http://edamontology.org/format_1930
+ - class: File
+ location: rnaseq/raw_fastq/Irrel_kd_1.subset.fq
+ format: http://edamontology.org/format_1930
+ - class: File
+ location: rnaseq/raw_fastq/Irrel_kd_2.subset.fq
+ format: http://edamontology.org/format_1930
+ - class: File
+ location: rnaseq/raw_fastq/Irrel_kd_3.subset.fq
+ format: http://edamontology.org/format_1930
+genome:
+ class: Directory
+ location: hg19-chr1-STAR-index
+gtf:
+ class: File
+ location: rnaseq/reference_data/chr1-hg19_genes.gtf
+```
+
+Now you can run the workflow the same way as in Lesson 2.
+
+4. Combining results
+
+Each instance of the alignment workflow produces its own featureCounts
+file. However, to be able to compare results easily, we need them a
+single file with all the results.
+
+The easiest way to do this is to run `featureCounts` just once at the
+end of the workflow, with all the bam files listed on the command
+line.
+
+We'll need to modify a few things.
+
+First, in `featureCounts.cwl` we need to modify it to accept either a
+single bam file or list of bam files.
+
+```
+inputs:
+ gtf: File
+ counts_input_bam:
+ - File
+ - File[]
+```
+
+Second, in `alignment.cwl` we need to remove the `featureCounts` step from alignment.cwl, as well as the `featurecounts` output parameter.
+
+Third, in `main.cwl` we need to remove `featurecounts` from the `alignment` step
+outputs, and add a new step:
+
+```
+steps:
+ alignment:
+ run: alignment.cwl
+ scatter: fq
+ in:
+ fq: fq
+ genome: genome
+ gtf: gtf
+ out: [qc_html, bam_sorted_indexed]
+ featureCounts:
+ requirements:
+ ResourceRequirement:
+ ramMin: 500
+ run: featureCounts.cwl
+ in:
+ counts_input_bam: alignment/bam_sorted_indexed
+ gtf: gtf
+ out: [featurecounts]
+```
+
+Last, we modify the `featurecounts` output parameter. Instead of a
+list of files produced by the `alignment` step, it is now a single
+file produced by the new `featureCounts` step.
+
+```
+outputs:
+ ...
+ featurecounts:
+ type: File
+ outputSource: featureCounts/featurecounts
+```
+
+Run this workflow to get a single `featurecounts.tsv` file with a column for each bam file.
--- /dev/null
+# Expressions
+
+1. Expressions on step inputs
+
+You might have noticed that the output bam files are all named
+`Aligned.sortedByCoord.out.bam`. This happens because because when we
+call STAR, it gives the output a default file name.
+
+Now, during workflow execution, this is usually not a problem. The
+workflow runner is smart enough to know that these files are different
+and keep them separate. This can even make development easier by not
+having to worry about assigning unique file names to every file.
+
+However, it is a problem for humans interpreting the output. We can
+fix this by setting the parameter `OutFileNamePrefix` on STAR.
+
+In `alignment.cwl`, we can use `valueFrom` on the `OutFileNamePrefix`
+input parameter to construct the output prefix from the input
+filename.
+
+```
+requirements:
+ StepInputExpressionRequirement: {}
+steps:
+ ...
+ STAR:
+ ...
+ in:
+ ForwardReads: fq
+ ...
+ OutFileNamePrefix: {valueFrom: "$(inputs.ForwardReads.nameroot)."}
+```
+
+The code between `$(...)` is called an "expression". It is evaluated
+when setting up the step to run, and the expression is replaced by the
+result to get the parameter value.
+
+An expression can refer to other inputs to the step that are either
+directly connected to another value, or have a default value. Here,
+we refer to the input parameter ForwardReads, which is our fastq input
+file.
+
+ForwardReads is a File object, not a plain file name, so it has some
+fields of its own. The file object has a number of fields that we can
+use. These include `basename` (the name of the file, without a
+directory), `size` (file size, in bytes), `nameext` (the last file
+extension) and `nameroot` (the name with `nameext` removed). Using
+
+Finally, our expression is embedded in a string, so after replacing
+the expression with the value of `inputs.ForwardReads.nameroot`, it
+adds the remainder of the string, which just is a dot `.`. This is to
+separate the leading part of our filename from the "Aligned.bam"
+extension that will be added by STAR.
+
+2. Organizing output files into Directories
+
+This is a more advanced example.
+
+Create `subdirs.cwl`:
+
+```
+cwlVersion: v1.2
+class: ExpressionTool
+requirements:
+ InlineJavascriptRequirement: {}
+inputs:
+ fq: File[]
+ bams: File[]
+ qc: File[]
+outputs:
+ dirs: Directory[]
+expression: |-
+ ${
+ var dirs = [];
+ for (var i = 0; i < inputs.bams.length; i++) {
+ dirs.push({
+ "class": "Directory",
+ "basename": inputs.fq[i].nameroot,
+ "listing": [inputs.bams[i], inputs.qc[i]]
+ });
+ }
+ return {"dirs": dirs};
+ }
+```
+
+Then change `main.cwl`:
+
+```
+steps:
+ ...
+ output-subdirs:
+ run: subdirs.cwl
+ in:
+ fq: fq
+ bams: alignment/bam_sorted_indexed
+ qc: alignment/qc_html
+ out: [dirs]
+outputs:
+ dirs:
+ type: Directory[]
+ outputSource: output-subdirs/dirs
+ featurecounts:
+ type: File
+ outputSource: featureCounts/featurecounts
+```
+
+3. Other uses for expressions
+
+- Creating configuration files on the fly.
projects / rnaseq-cwl-training.git / commitdiff