From cd44845d58f9bfe7dd31a37ba9dd18aea3001b62 Mon Sep 17 00:00:00 2001
From: Peter Amstutz <peter.amstutz@curii.com>
Date: Fri, 22 Jan 2021 17:35:42 -0500
Subject: [PATCH] Done with 1st draft of lesson 4, working on lesson 5

Arvados-DCO-1.1-Signed-off-by: Peter Amstutz <peter.amstutz@curii.com>
---
 lesson1/lesson1.md |   9 +-
 lesson2/lesson2.md |  14 +---
 lesson3/lesson3.md | 119 +++++++++++++++++++-------
 lesson4/lesson4.md | 203 +++++++++++++++++++++++++++++++++++++++++++++
 lesson5/lesson5.md | 109 ++++++++++++++++++++++++
 5 files changed, 410 insertions(+), 44 deletions(-)
 create mode 100644 lesson4/lesson4.md
 create mode 100644 lesson5/lesson5.md

diff --git a/lesson1/lesson1.md b/lesson1/lesson1.md
index 4846d72..b989855 100644
--- a/lesson1/lesson1.md
+++ b/lesson1/lesson1.md
@@ -1,4 +1,4 @@
-# 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.
 
@@ -54,7 +54,7 @@ series of steps, and produce output values.
 
 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
@@ -112,7 +112,7 @@ steps:
     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
@@ -186,9 +186,6 @@ outputs:
   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
diff --git a/lesson2/lesson2.md b/lesson2/lesson2.md
index 02bedfb..0df5811 100644
--- a/lesson2/lesson2.md
+++ b/lesson2/lesson2.md
@@ -4,7 +4,7 @@
 
 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`
@@ -58,7 +58,7 @@ In vscode, select "main.cwl" and then choose "Terminal -> Run task -> Run CWL wo
 Type this into the terminal:
 
 ```
-cwl-runner main.cwl input.yaml
+cwl-runner main.cwl main-input.yaml
 ```
 
 3. Debugging the workflow
@@ -103,7 +103,7 @@ request more RAM, add a "requirements" section with
   STAR:
     requirements:
       ResourceRequirement:
-        ramMin: 7000
+        ramMin: 8000
     run: bio-cwl-tools/STAR/STAR-Align.cwl
 ```
 
@@ -142,17 +142,11 @@ The CWL runner will print a results JSON object to standard output.  It will loo
         "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.
 
diff --git a/lesson3/lesson3.md b/lesson3/lesson3.md
index 93e9a55..b7e3334 100644
--- a/lesson3/lesson3.md
+++ b/lesson3/lesson3.md
@@ -22,7 +22,7 @@ A CommandLineTool describes a single invocation of a command line program.
 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
@@ -40,15 +40,15 @@ inputs:
   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.
@@ -61,60 +61,123 @@ Special variables are also available.  The runtime environment
 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.
diff --git a/lesson4/lesson4.md b/lesson4/lesson4.md
new file mode 100644
index 0000000..447726d
--- /dev/null
+++ b/lesson4/lesson4.md
@@ -0,0 +1,203 @@
+# 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.
diff --git a/lesson5/lesson5.md b/lesson5/lesson5.md
new file mode 100644
index 0000000..fe3442d
--- /dev/null
+++ b/lesson5/lesson5.md
@@ -0,0 +1,109 @@
+# 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.
-- 
2.30.2