Merge branch 'template_update_0624' into 'dev'

Template update

See merge request epi2melabs/workflows/wf-transcriptomes!175

.pre-commit-config.yaml

@@ -8,7 +8,7 @@ repos:
         always_run: true
         pass_filenames: false
         additional_dependencies:
-          - epi2melabs>=0.0.53
+          - epi2melabs==0.0.55
       - id: build_models
         name: build_models
         entry: datamodel-codegen --strict-nullable --base-class workflow_glue.results_schema_helpers.BaseModel --use-schema-description --disable-timestamp --input results_schema.yml --input-file-type openapi --output bin/workflow_glue/results_schema.py

README.md

@@ -38,40 +38,68 @@ ARM processor support: False
 
 ## Install and run
 
-These are instructions to install and run the workflow on command line. You can also access the workflow via the [EPI2ME application](https://labs.epi2me.io/downloads/).
 
-The workflow uses [Nextflow](https://www.nextflow.io/) to manage compute and software resources, therefore nextflow will need to be installed before attempting to run the workflow.
+These are instructions to install and run the workflow on command line.
+You can also access the workflow via the
+[EPI2ME Desktop application](https://labs.epi2me.io/downloads/).
 
-The workflow can currently be run using either [Docker](https://www.docker.com/products/docker-desktop) or
+The workflow uses [Nextflow](https://www.nextflow.io/) to manage
-[Singularity](https://docs.sylabs.io/guides/3.0/user-guide/index.html) to provide isolation of
+compute and software resources,
-the required software. Both methods are automated out-of-the-box provided
+therefore Nextflow will need to be
-either docker or singularity is installed. This is controlled by the [`-profile`](https://www.nextflow.io/docs/latest/config.html#config-profiles) parameter as exemplified below.
+installed before attempting to run the workflow.
 
-It is not required to clone or download the git repository in order to run the workflow.
+The workflow can currently be run using either
-More information on running EPI2ME workflows can be found on our [website](https://labs.epi2me.io/wfindex).
+[Docker](https://www.docker.com/products/docker-desktop
+or [Singularity](https://docs.sylabs.io/guides/3.0/user-guide/index.html)
+to provide isolation of the required software.
+Both methods are automated out-of-the-box provided
+either Docker or Singularity is installed.
+This is controlled by the
+[`-profile`](https://www.nextflow.io/docs/latest/config.html#config-profiles)
+parameter as exemplified below.
 
-The following command can be used to obtain the workflow. This will pull the repository in to the assets folder of nextflow and provide a list of all parameters available for the workflow as well as an example command:
+It is not required to clone or download the git repository
+in order to run the workflow.
+More information on running EPI2ME workflows can
+be found on our [website](https://labs.epi2me.io/wfindex).
+
+The following command can be used to obtain the workflow.
+This will pull the repository in to the assets folder of
+Nextflow and provide a list of all parameters
+available for the workflow as well as an example command:
 
 ```
-nextflow run epi2me-labs/wf-transcriptomes -–help
+nextflow run epi2me-labs/wf-transcriptomes --help
 ```
-A demo dataset is provided for testing of the workflow. It can be downloaded using:
+To update a workflow to the latest version on the command line use
+the following command:
 ```
-wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz
+nextflow pull epi2me-labs/wf-transcriptomes
-tar -xzvf differential_expression.tar.gz
 ```
-The workflow can be run with the demo data using:
+A demo dataset is provided for testing of the workflow.
+It can be downloaded and unpacked using the following commands:
+```
+wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/wf-transcriptomes-demo.tar.gz
+tar -xzvf wf-transcriptomes-demo.tar.gz
+```
+The workflow can then be run with the downloaded demo data using:
 ```
 nextflow run epi2me-labs/wf-transcriptomes \
---fastq  differential_expression/differential_expression_fastq \
+	--de_analysis \
---de_analysis --ref_genome differential_expression/hg38_chr20.fa \
+	--direct_rna \
---transcriptome-source reference-guided \
+	--fastq 'wf-transcriptomes-demo/differential_expression_fastq' \
---ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \
+	--jaffal_annotation 'genCode22' \
---direct_rna --minimap2_index_opts '-k 15'  --sample_sheet differential_expression/sample_sheet.csv \
+	--jaffal_genome 'hg38_chr20' \
---jaffal_refBase differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
+	--jaffal_refBase 'wf-transcriptomes-demo/chr20' \
--profile standard
+	--minimap2_index_opts '-k15' \
+	--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \
+	--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \
+	--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv' \
+	-profile standard
 ```
-For further information about running a workflow on the cmd line see https://labs.epi2me.io/wfquickstart/
+For further information about running a workflow on
+the command line see https://labs.epi2me.io/wfquickstart/
+
 
 
 

bin/workflow_glue/check_sample_sheet.py

@@ -38,6 +38,7 @@ def main(args):
     barcodes = []
     aliases = []
     sample_types = []
+    analysis_groups = []
     allowed_sample_types = [
         "test_sample", "positive_control", "negative_control", "no_template_control"
     ]
@@ -91,6 +92,10 @@ def main(args):
                     sample_types.append(row["type"])
                 except KeyError:
                     pass
+                try:
+                    analysis_groups.append(row["analysis_group"])
+                except KeyError:
+                    pass
     except Exception as e:
         sys.stdout.write(f"Parsing error: {e}")
         sys.exit()
@@ -136,6 +141,14 @@ def main(args):
                     sys.stdout.write(
                         f"Sample sheet requires at least 1 of {required_type}")
                     sys.exit()
+    if analysis_groups:
+        # if there was a "analysis_group" column, make sure it had values for all
+        # samples
+        if not all(analysis_groups):
+            sys.stdout.write(
+                "if an 'analysis_group' column exists, it needs values in each row"
+            )
+            sys.exit()
 
     logger.info(f"Checked sample sheet {args.sample_sheet}.")
 

bin/workflow_glue/configure_igv.py

@@ -8,16 +8,16 @@ from .util import get_named_logger, wf_parser  # noqa: ABS101
 
 
 def parse_fnames(fofn):
-    """Parse list with filenames and return them grouped as ref-, BAM-, or VCF-related.
+    """Parse list with filenames and return them grouped as ref-, XAM-, or VCF-related.
 
     :param fofn: File with list of file names (one per line)
     :return: dict of reference-related filenames (with keys 'ref', 'fai', and '.gzi' and
-        `None` as default values); lists of BAM- and VCF-related filenames
+        `None` as default values); lists of XAM- and VCF-related filenames
     """
     ref_extensions = [".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fna.gz"]
     ref_dict = {}
-    bams = []
+    xams = []
-    bam_indices = []
+    xam_indices = []
     vcfs = []
     vcf_indices = []
     with open(fofn, "r") as f:
@@ -29,10 +29,10 @@ def parse_fnames(fofn):
                 ref_dict["fai"] = fname
             elif fname.endswith(".gzi"):
                 ref_dict["gzi"] = fname
-            elif fname.endswith(".bam"):
+            elif fname.endswith(".bam") or fname.endswith(".cram"):
-                bams.append(fname)
+                xams.append(fname)
-            elif fname.endswith(".bai"):
+            elif fname.endswith(".bai") or fname.endswith(".crai"):
-                bam_indices.append(fname)
+                xam_indices.append(fname)
             elif fname.endswith(".vcf") or fname.endswith(".vcf.gz"):
                 vcfs.append(fname)
             elif fname.endswith(".csi") or fname.endswith(".tbi"):
@@ -52,20 +52,20 @@ def parse_fnames(fofn):
                 f"Found GZI reference index '{gzi}', but the reference file "
                 f"'{ref}' appears not to be compressed."
             )
-    if bam_indices:
+    if xam_indices:
-        if len(bams) != len(bam_indices):
+        if len(xams) != len(xam_indices):
-            raise ValueError("Got different number of BAM and BAM index files.")
+            raise ValueError("Got different number of XAM and XAM index files.")
     if vcf_indices:
         if len(vcfs) != len(vcf_indices):
             raise ValueError("Got different number of VCF and VCF index files.")
-    if bams and vcfs:
+    if xams and vcfs:
-        if len(bams) != len(vcfs):
+        if len(xams) != len(vcfs):
-            raise ValueError("Got different number of BAM and VCF files.")
+            raise ValueError("Got different number of XAM and VCF files.")
-    # if we got BAM or VCF indices, pair them up with their corresponding files (and
+    # if we got XAM or VCF indices, pair them up with their corresponding files (and
     # otherwise with `None`)
-    bams_with_indices = zip_longest(bams, bam_indices)
+    xams_with_indices = zip_longest(xams, xam_indices)
     vcfs_with_indices = zip_longest(vcfs, vcf_indices)
-    return ref_dict, bams_with_indices, vcfs_with_indices
+    return ref_dict, xams_with_indices, vcfs_with_indices
 
 
 def get_reference_options(ref, fai=None, gzi=None):
@@ -90,23 +90,23 @@ def get_reference_options(ref, fai=None, gzi=None):
     return ref_opts
 
 
-def get_alignment_track(bam, bai=None, extra_opts=None):
+def get_alignment_track(xam, xai=None, extra_opts=None):
     """Create dict with options for IGV alignment track.
 
-    :param bam: name of BAM file to be displayed
+    :param xam: name of XAM file to be displayed
-    :param bai: name of BAM index file
+    :param xai: name of XAM index file
     :param extra_opts: dict of extra options for the alignment track
     :return: dict with alignment track options
     """
     alignment_track_dict = {
-        "name": bam,
+        "name": xam,
         "type": "alignment",
-        "format": "bam",
+        "format": xam.split(".")[-1],
-        "url": bam,
+        "url": xam,
     }
-    # add the BAM index if present
+    # add the XAM index if present
-    if bai is not None:
+    if xai is not None:
-        alignment_track_dict["indexURL"] = bai
+        alignment_track_dict["indexURL"] = xai
     alignment_track_dict.update(extra_opts or {})
     return alignment_track_dict
 
@@ -137,7 +137,7 @@ def main(args):
     logger = get_named_logger("configIGV")
 
     # parse the FOFN
-    ref_dict, bams_with_indices, vcfs_with_indices = parse_fnames(args.fofn)
+    ref_dict, xams_with_indices, vcfs_with_indices = parse_fnames(args.fofn)
 
     # initialise the IGV options dict with the reference options
     json_dict = {"reference": get_reference_options(**ref_dict)}
@@ -145,30 +145,30 @@ def main(args):
     # if we got JSON files with extra options for the alignment / variant tracks, read
     # them
     extra_alignment_opts = {}
-    if args.extra_bam_opts is not None:
+    if args.extra_alignment_opts is not None:
-        with open(args.extra_bam_opts, "r") as f:
+        with open(args.extra_alignment_opts, "r") as f:
             extra_alignment_opts = json.load(f)
     extra_variant_opts = {}
-    if args.extra_vcf_opts is not None:
+    if args.extra_variant_opts is not None:
-        with open(args.extra_vcf_opts, "r") as f:
+        with open(args.extra_variant_opts, "r") as f:
             extra_variant_opts = json.load(f)
 
     # now add the alignment and variant tracks
     json_dict["tracks"] = []
     # we use `zip_longest` to make sure that variant and alignment tracks from the same
     # sample are added after each other
-    for (vcf, vcf_index), (bam, bam_index) in zip_longest(
+    for (vcf, vcf_index), (xam, xam_index) in zip_longest(
-        vcfs_with_indices, bams_with_indices, fillvalue=(None, None)
+        vcfs_with_indices, xams_with_indices, fillvalue=(None, None)
     ):
         if vcf is not None:
-            # add an variant track for the VCF
+            # add a variant track for the VCF
             json_dict["tracks"].append(
                 get_variant_track(vcf, vcf_index, extra_variant_opts)
             )
-        if bam is not None:
+        if xam is not None:
-            # add an alignment track for the BAM
+            # add an alignment track for the XAM
             json_dict["tracks"].append(
-                get_alignment_track(bam, bam_index, extra_alignment_opts)
+                get_alignment_track(xam, xam_index, extra_alignment_opts)
             )
 
     if args.locus is not None:
@@ -186,7 +186,7 @@ def argparser():
         "--fofn",
         required=True,
         help=(
-            "File with list of names of reference / BAM / VCF files and indices "
+            "File with list of names of reference / XAM / VCF files and indices "
             "(one filename per line)"
         ),
     )
@@ -195,11 +195,11 @@ def argparser():
         help="Locus string to set initial genomic coordinates to display in IGV",
     )
     parser.add_argument(
-        "--extra-bam-opts",
+        "--extra-alignment-opts",
         help="JSON file with extra options for alignment tracks",
     )
     parser.add_argument(
-        "--extra-vcf-opts",
+        "--extra-variant-opts",
         help="JSON file with extra options for variant tracks",
     )
     return parser

bin/workflow_glue/get_ds_records.py

@@ -1,166 +0,0 @@
-"""Extract unique values for a key from XAM RG DS headers or FASTX comments.
-
-Use pysam to read the description tag of XAM read group header(s) to
-collect values for a given key and check the expected cardinality from
-one or more XAM; or equivalently the comments of one or more FASTX to
-do the same.
-
-Use for example, to extract a single basecaller configuration name in
-order to match the input data to suitable models for downstream tools
-without troubling the user. No guarantee is made for ordering.
-"""
-
-import os
-import sys
-
-import pysam
-
-from .util import wf_parser  # noqa: ABS101
-
-
-# This is not my ideal way to raise nice user-facing errors, but embedding them in the
-#  Nextflow process itself is (a) confusing for users; as the error log includes echo
-#  commands which may or may not actually be printed and (b) a footgun for developers;
-#  who may inadvertently mishandle catching and then (re)returning a bad exit code.
-# This also rather neatly keeps intended messaging for users near the code that will
-#  cause errors to be raised.
-def get_extended_errmsg(key, expected_cardinality):
-    """Get an applicable extended error message for a given key and cardinality."""
-    if key == "basecall_model" and expected_cardinality == "zero-or-one":
-        return """
-################################################################################
-# INPUT DATA PROBLEM
-Your input data contains reads basecalled with more than one basecaller model.
-
-Our workflows automatically select appropriate configuration and models for
-downstream tools for a given basecaller model. This cannot be done reliably when
-reads with different basecaller models are mixed in the same data set.
-
-## Next steps
-To use this workflow you must separate your input files, making sure all reads
-are have been basecalled with the same basecaller model.
-################################################################################
-"""
-
-
-def path_to_lofn(input_path):
-    """Convert the input path to a list of one or more files to be checked."""
-    if os.path.isdir(input_path):
-        return [
-            os.path.join(root, f)
-            for (root, dirnames, filenames) in os.walk(input_path)
-            for f in filenames
-        ]
-    else:
-        return [input_path]
-
-
-def xam_extract_ds_key(xam_lofn, key):
-    """Extract the set of values for a given key from all RG DS tags."""
-    entries = set()
-    for xam_fn in xam_lofn:
-        with pysam.AlignmentFile(xam_fn, check_sq=False) as xam:
-            for read_group in xam.header.get("RG", []):
-                for ds_kv in read_group.get("DS", "").split():
-                    k, v = ds_kv.split("=", 1)
-                    if k == key:
-                        entries.add(v)
-    return entries
-
-
-def fastx_extract_ds_key(fastx_lofn, key, stop_after=0):
-    """Extract the set of values for a given key from all FASTQ tags."""
-    entries = set()
-    for fastx_fn in fastx_lofn:
-        with pysam.FastxFile(fastx_fn) as fastx:
-            for i, read in enumerate(fastx):
-                if stop_after > 0 and i > stop_after:
-                    break
-                for ds_kv in read.comment.split():
-                    k, v = ds_kv.split("=", 1)
-                    if k == key:
-                        entries.add(v)
-    return entries
-
-
-def check_cardinality(obs, desired_cardinality):
-    """Return whether the observed cardinality meets the desired cardinality."""
-    cardinality_lookup = {
-        0: ["zero", "zero-or-one", "zero-or-more"],
-        1: ["zero-or-one", "zero-or-more", "one", "one-or-more"],
-        2: ["zero-or-more", "one-or-more", "more-than-one"],
-    }
-    if obs > 1:
-        obs = 2
-    return desired_cardinality in cardinality_lookup[obs]
-
-
-def main(args):
-    """Script entrypoint.
-
-    Extracts values using the XAM or FASTX extractor and checks the set
-    is of the right cardinality. Prints a message to stdout (or stderr)
-    and exits appropriately.
-    """
-    if args.xam:
-        extractor = xam_extract_ds_key
-        input_path = args.xam
-    elif args.fastx:
-        extractor = fastx_extract_ds_key
-        input_path = args.fastx
-
-    lofn = path_to_lofn(input_path)
-    entries = extractor(lofn, args.key)
-    if not check_cardinality(len(entries), args.cardinality):
-        sys.stdout.write(args.sep.join(entries) + '\n')
-        extended_error = get_extended_errmsg(args.key, args.cardinality)
-        if args.explode_obviously and extended_error:
-            sys.stderr.write(extended_error)
-        else:
-            sys.stderr.write(
-                f"Required {args.cardinality} {args.key} but found {len(entries)}\n"
-            )
-        sys.exit(os.EX_DATAERR)
-    sys.stdout.write(args.sep.join(entries) + '\n')
-    sys.exit(os.EX_OK)
-
-
-def argparser():
-    """Argument parser for entrypoint."""
-    parser = wf_parser("get_ds_records")
-    parser.add_argument("--key", required=True)
-    input_arg = parser.add_mutually_exclusive_group(required=True)
-    input_arg.add_argument("--xam", help="Path to a single XAM or folder of XAM")
-    input_arg.add_argument("--fastx", help="Path to a single FASTX or folder of FASTX")
-    parser.add_argument(
-        "--sep",
-        default="\n",
-        help=(
-            "Value separator to use if more than one element is printed to stdout."
-        ),
-    )
-    parser.add_argument(
-        "--cardinality",
-        default="zero-or-more",
-        choices=[
-            "zero",
-            "zero-or-one",
-            "zero-or-more",
-            "one",
-            "one-or-more",
-            "more-than-one",
-        ],
-        help=(
-            "Expected cardinality of entries. Will exit EX_DATAERR if the wrong "
-            "number of elements are found. Defaults to zero-or-more."
-        ),
-    )
-    parser.add_argument(
-        "--explode_obviously",
-        action="store_true",
-        help=(
-            "If appropriate, print a more fulsome user facing error if the expected"
-            "cardinality has been violated."
-        )
-    )
-    return parser

docs/04_install_and_run.md

@@ -1,34 +1,61 @@
-These are instructions to install and run the workflow on command line. You can also access the workflow via the [EPI2ME application](https://labs.epi2me.io/downloads/).
 
-The workflow uses [Nextflow](https://www.nextflow.io/) to manage compute and software resources, therefore nextflow will need to be installed before attempting to run the workflow.
+These are instructions to install and run the workflow on command line.
+You can also access the workflow via the
+[EPI2ME Desktop application](https://labs.epi2me.io/downloads/).
 
-The workflow can currently be run using either [Docker](https://www.docker.com/products/docker-desktop) or
+The workflow uses [Nextflow](https://www.nextflow.io/) to manage
-[Singularity](https://docs.sylabs.io/guides/3.0/user-guide/index.html) to provide isolation of
+compute and software resources,
-the required software. Both methods are automated out-of-the-box provided
+therefore Nextflow will need to be
-either docker or singularity is installed. This is controlled by the [`-profile`](https://www.nextflow.io/docs/latest/config.html#config-profiles) parameter as exemplified below.
+installed before attempting to run the workflow.
 
-It is not required to clone or download the git repository in order to run the workflow.
+The workflow can currently be run using either
-More information on running EPI2ME workflows can be found on our [website](https://labs.epi2me.io/wfindex).
+[Docker](https://www.docker.com/products/docker-desktop
+or [Singularity](https://docs.sylabs.io/guides/3.0/user-guide/index.html)
+to provide isolation of the required software.
+Both methods are automated out-of-the-box provided
+either Docker or Singularity is installed.
+This is controlled by the
+[`-profile`](https://www.nextflow.io/docs/latest/config.html#config-profiles)
+parameter as exemplified below.
 
-The following command can be used to obtain the workflow. This will pull the repository in to the assets folder of nextflow and provide a list of all parameters available for the workflow as well as an example command:
+It is not required to clone or download the git repository
+in order to run the workflow.
+More information on running EPI2ME workflows can
+be found on our [website](https://labs.epi2me.io/wfindex).
+
+The following command can be used to obtain the workflow.
+This will pull the repository in to the assets folder of
+Nextflow and provide a list of all parameters
+available for the workflow as well as an example command:
 
 ```
-nextflow run epi2me-labs/wf-transcriptomes -–help
+nextflow run epi2me-labs/wf-transcriptomes --help
 ```
-A demo dataset is provided for testing of the workflow. It can be downloaded using:
+To update a workflow to the latest version on the command line use
+the following command:
 ```
-wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz
+nextflow pull epi2me-labs/wf-transcriptomes
-tar -xzvf differential_expression.tar.gz
 ```
-The workflow can be run with the demo data using:
+A demo dataset is provided for testing of the workflow.
+It can be downloaded and unpacked using the following commands:
+```
+wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/wf-transcriptomes-demo.tar.gz
+tar -xzvf wf-transcriptomes-demo.tar.gz
+```
+The workflow can then be run with the downloaded demo data using:
 ```
 nextflow run epi2me-labs/wf-transcriptomes \
---fastq  differential_expression/differential_expression_fastq \
+	--de_analysis \
---de_analysis --ref_genome differential_expression/hg38_chr20.fa \
+	--direct_rna \
---transcriptome-source reference-guided \
+	--fastq 'wf-transcriptomes-demo/differential_expression_fastq' \
---ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \
+	--jaffal_annotation 'genCode22' \
---direct_rna --minimap2_index_opts '-k 15'  --sample_sheet differential_expression/sample_sheet.csv \
+	--jaffal_genome 'hg38_chr20' \
---jaffal_refBase differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
+	--jaffal_refBase 'wf-transcriptomes-demo/chr20' \
--profile standard
+	--minimap2_index_opts '-k15' \
+	--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \
+	--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \
+	--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv' \
+	-profile standard
 ```
-For further information about running a workflow on the cmd line see https://labs.epi2me.io/wfquickstart/
+For further information about running a workflow on
+the command line see https://labs.epi2me.io/wfquickstart/

lib/common.nf

@@ -23,33 +23,33 @@ process configure_igv {
         // extensions
         path "file-names.txt"
         val locus_str
-        val bam_extra_opts
+        val aln_extra_opts
-        val vcf_extra_opts
+        val var_extra_opts
     output: path "igv.json"
     script:
     // the locus argument just makes sure that the initial view in IGV shows something
     // interesting
     String locus_arg = locus_str ? "--locus $locus_str" : ""
     // extra options for alignment tracks
-    def bam_opts_json_str = \
+    def aln_opts_json_str = \
-        bam_extra_opts ? new JsonBuilder(bam_extra_opts).toPrettyString() : ""
+        aln_extra_opts ? new JsonBuilder(aln_extra_opts).toPrettyString() : ""
-    String bam_extra_opts_arg = \
+    String aln_extra_opts_arg = \
-        bam_extra_opts ? "--extra-bam-opts bam-extra-opts.json" : ""
+        aln_extra_opts ? "--extra-alignment-opts extra-aln-opts.json" : ""
     // extra options for variant tracks
-    def vcf_opts_json_str = \
+    def var_opts_json_str = \
-        vcf_extra_opts ? new JsonBuilder(vcf_extra_opts).toPrettyString() : ""
+        var_extra_opts ? new JsonBuilder(var_extra_opts).toPrettyString() : ""
-    String vcf_extra_opts_arg = \
+    String var_extra_opts_arg = \
-        vcf_extra_opts ? "--extra-vcf-opts vcf-extra-opts.json" : ""
+        var_extra_opts ? "--extra-vcf-opts extra-var-opts.json" : ""
     """
     # write out JSON files with extra options for the alignment and variant tracks
-    echo '$bam_opts_json_str' > bam-extra-opts.json
+    echo '$aln_opts_json_str' > extra-aln-opts.json
-    echo '$vcf_opts_json_str' > vcf-extra-opts.json
+    echo '$var_opts_json_str' > extra-var-opts.json
 
     workflow-glue configure_igv \
         --fofn file-names.txt \
         $locus_arg \
-        $bam_extra_opts_arg \
+        $aln_extra_opts_arg \
-        $vcf_extra_opts_arg \
+        $var_extra_opts_arg \
     > igv.json
     """
 }

lib/ingress.nf

@@ -26,28 +26,40 @@ def is_target_file(Path file, List extensions) {
 /**
  * Take a channel of the shape `[meta, reads, path-to-stats-dir | null]` (or
  * `[meta, [reads, index], path-to-stats-dir | null]` in the case of XAM) and extract the
- * run IDs from the `run_ids` file in the stats directory into the metamap. If the path
+ * run IDs and basecall model, from the `run_ids` and `basecaller` files in the stats
- * to the stats dir is `null`, add an empty list.
+ * directory, into the metamap. If the path to the stats dir is `null`, add an empty list.
  *
  * @param ch: input channel of shape `[meta, reads, path-to-stats-dir | null]`
- * @return: channel with a list of run IDs added to the metamap
+ * @return: channel with lists of run IDs and basecall models added to the metamap
  */
-def add_run_IDs_to_meta(ch) {
+def add_run_IDs_and_basecall_models_to_meta(ch, boolean allow_multiple_basecall_models) {
     // HashSet for all observed run_ids
     Set<String> ingressed_run_ids = new HashSet<String>()
 
     // extract run_ids from fastcat stats / bamstats results and add to metadata as well
     // as `ingressed_run_ids`
     ch = ch | map { meta, reads, stats ->
-        ArrayList run_ids = []
         if (stats) {
             run_ids = stats.resolve("run_ids").splitText().collect { it.strip() }
             ingressed_run_ids += run_ids
+
+            basecall_models = \
+                stats.resolve("basecallers").splitText().collect { it.strip() }
+            // check if we got more than one basecall model and set reads + stats to
+            // `null` for that sample unless `allow_multiple_basecall_models`
+            if ((basecall_models.size() > 1) && !allow_multiple_basecall_models) {
+                log.warn "Found multiple basecall models for sample " + \
+                    "'$meta.alias': ${basecall_models.join(", ")}. The sample's " + \
+                    "reads were discarded."
+                reads = reads instanceof List ? [null, null] : null
+                stats = null
             }
             // `meta + [...]` returns a new map which is handy to avoid any
             // modifying-maps-in-closures weirdness
             // See https://github.com/nextflow-io/nextflow/issues/2660
-        [meta + [run_ids: run_ids], reads, stats]
+            meta = meta + [run_ids: run_ids, basecall_models: basecall_models]
+        }
+        [meta, reads, stats]
     }
     // put run_ids somewhere global for trivial access later
     // bit grim but decouples ingress metadata from workflow main.nf
@@ -125,6 +137,9 @@ def add_number_of_reads_to_meta(ch, String input_type_format) {
  *  - "required_sample_types": list of required sample types in the sample sheet
  *  - "watch_path": boolean whether to use `watchPath` and run in streaming mode
  *  - "fastq_chunk": null or a number of reads to place into chunked FASTQ files
+ *  - "allow_multiple_basecall_models": emit data of samples that had more than one
+ *     basecall model; if this is `false`, such samples will be emitted as `[meta, null,
+ *     null]`
  * @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
  *  The first element is a map with metadata, the second is the path to the
  *  `.fastq.gz` file with the (potentially concatenated) sequences and the third is
@@ -140,7 +155,6 @@ def fastq_ingress(Map arguments)
         [
             "fastcat_extra_args": "",
             "fastq_chunk": null,
-            "per_read_stats": false
         ]
     )
     margs["fastq_chunk"] ?= 0  // cant pass null through channel
@@ -196,8 +210,11 @@ def fastq_ingress(Map arguments)
             [meta + new_keys, files, stats]
         }
 
-    def ch_final = add_number_of_reads_to_meta(
+    // add number of reads, run IDs, and basecall models to meta
-        add_run_IDs_to_meta(ch_spread_result), "fastq")
+    def ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
+    ch_final = add_run_IDs_and_basecall_models_to_meta(
+        ch_final, margs.allow_multiple_basecall_models
+    )
     return ch_final
 }
 
@@ -242,7 +259,6 @@ def xam_ingress(Map arguments)
             "return_fastq": false,
             "fastcat_extra_args": "",
             "fastq_chunk": null,
-            "per_read_stats": false
         ]
     )
     margs["fastq_chunk"] ?= 0  // cant pass null through channel
@@ -276,7 +292,7 @@ def xam_ingress(Map arguments)
         [meta + [xai_fn: xai_fn], paths]
     }
     | checkBamHeaders
-    | map { meta, paths, is_unaligned_env, mixed_headers_env, is_sorted_env, ds_basecaller_env, ds_runids_env ->
+    | map { meta, paths, is_unaligned_env, mixed_headers_env, is_sorted_env ->
         // convert the env. variables from strings ('0' or '1') into bools
         boolean is_unaligned = is_unaligned_env as int as boolean
         boolean mixed_headers = mixed_headers_env as int as boolean
@@ -288,15 +304,7 @@ def xam_ingress(Map arguments)
         // add `is_unaligned` to the metamap (note the use of `+` to create a copy of
         // `meta` to avoid modifying every item in the channel;
         // https://github.com/nextflow-io/nextflow/issues/2660)
-        [
+        [meta + [is_unaligned: is_unaligned, is_sorted: is_sorted], paths]
-            meta + [
-                is_unaligned: is_unaligned,
-                is_sorted: is_sorted,
-                ds_runids: ds_runids_env.tokenize(','),
-                ds_basecall_models: ds_basecaller_env.tokenize(','),
-            ],
-            paths
-        ]
     }
     | branch { meta, paths ->
         // set `paths` to `null` for uBAM samples if unallowed (they will be added to
@@ -366,8 +374,11 @@ def xam_ingress(Map arguments)
                 [meta.findAll { it.key !in ['xai_fn', 'is_sorted'] }, path, stats]
             }
 
-        def ch_final = add_number_of_reads_to_meta(
+        // add number of reads, run IDs, and basecall models to meta
-            add_run_IDs_to_meta(ch_spread_result), "fastq")
+        def ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
+        ch_final = add_run_IDs_and_basecall_models_to_meta(
+            ch_final, margs.allow_multiple_basecall_models
+        )
         return ch_final
     }
 
@@ -457,7 +468,6 @@ def xam_ingress(Map arguments)
         | map{
             it[3] ? [it[0], [it[1], it[2]], it[3]] : it
         }
-        | add_run_IDs_to_meta
         | map{
             it.flatten()
         }
@@ -483,6 +493,7 @@ def xam_ingress(Map arguments)
     //     n_unmapped: always present, but can be `null`
     //     is_unaligned: present if there is a (u)BAM file
     // ]
+    // also, add number of reads, run IDs, and basecall models to meta
     ch_result = add_number_of_reads_to_meta(
         ch_result
             | map{
@@ -491,6 +502,9 @@ def xam_ingress(Map arguments)
             }, 
         "xam"
     )
+    ch_result = add_run_IDs_and_basecall_models_to_meta(
+        ch_result, margs.allow_multiple_basecall_models
+    )
     | map{
         it.flatten()
     }
@@ -529,7 +543,8 @@ process fastcat {
     fastcat \
         -s ${meta["alias"]} \
         -f fastcat_stats/per-file-stats.tsv \
-        -i fastcat_stats/per-file-runids.txt \
+        -i fastcat_stats/per-file-runids.tsv \
+        -l fastcat_stats/per-file-basecallers.tsv \
         --histograms histograms \
         $stats_args \
         ${fcargs["fastcat_extra_args"]} \
@@ -545,9 +560,19 @@ process fastcat {
     # get n_seqs from per-file stats - need to sum them up
     awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {c+=\$ix["n_seqs"]} END{print c}' \
         fastcat_stats/per-file-stats.tsv > fastcat_stats/n_seqs
-    # get unique run IDs
+    # get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace
-    awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {print \$ix["run_id"]}' \
+    # as field delimiter per default and thus ignores empty columns)
-        fastcat_stats/per-file-runids.txt | sort | uniq > fastcat_stats/run_ids
+    awk -F '\\t' '
+        NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
+        # only print run_id if present
+        NR>1 && \$ix["run_id"] != "" {print \$ix["run_id"]}
+    ' fastcat_stats/per-file-runids.tsv | sort | uniq > fastcat_stats/run_ids
+    # get unique basecall models
+    awk -F '\\t' '
+        NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
+        # only print basecall model if present
+        NR>1 && \$ix["basecaller"] != "" {print \$ix["basecaller"]}
+    ' fastcat_stats/per-file-basecallers.tsv | sort | uniq > fastcat_stats/basecallers
     """
 }
 
@@ -564,15 +589,11 @@ process checkBamHeaders {
             env(IS_UNALIGNED),
             env(MIXED_HEADERS),
             env(IS_SORTED),
-            env(DS_BASECALL_MODELS),
-            env(DS_RUNIDS),
         )
     script:
     """
     workflow-glue check_bam_headers_in_dir input_dir > env.vars
     source env.vars
-    DS_RUNIDS=\$(workflow-glue get_ds_records --xam input_dir --key runid --cardinality zero-or-more --sep ',')
-    DS_BASECALL_MODELS=\$(workflow-glue get_ds_records --xam input_dir --key basecall_model --cardinality zero-or-one --sep ',' --explode_obviously)
     """
 }
 
@@ -667,7 +688,8 @@ process bamstats {
     mkdir bamstats_results
     bamstats reads.bam -s $meta.alias -u \
         -f bamstats_results/bamstats.flagstat.tsv -t $bamstats_threads \
-        -i bamstats_results/bamstats.runids.txt \
+        -i bamstats_results/bamstats.runids.tsv \
+        -l bamstats_results/bamstats.basecallers.tsv \
         --histograms histograms \
     $per_read_stats_arg
     mv histograms/* bamstats_results/
@@ -675,9 +697,19 @@ process bamstats {
     # get n_seqs from flagstats - need to sum them up
     awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {c+=\$ix["total"]} END{print c}' \
         bamstats_results/bamstats.flagstat.tsv > bamstats_results/n_seqs
-    # get unique run IDs
+    # get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace
-    awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {print \$ix["run_id"]}' \
+    # as field delimiter otherwise and thus ignore empty columns)
-        bamstats_results/bamstats.runids.txt | sort | uniq > bamstats_results/run_ids
+    awk -F '\\t' '
+        NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
+        # only print run_id if present
+        NR>1 && \$ix["run_id"] != "" {print \$ix["run_id"]}
+    ' bamstats_results/bamstats.runids.tsv | sort | uniq > bamstats_results/run_ids
+    # get unique basecall models
+    awk -F '\\t' '
+        NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
+        # only print run_id if present
+        NR>1 && \$ix["basecaller"] != "" {print \$ix["basecaller"]}
+    ' bamstats_results/bamstats.basecallers.tsv | sort | uniq > bamstats_results/basecallers
     """
 }
 /**
@@ -840,7 +872,8 @@ Map parse_arguments(String func_name, Map arguments, Map extra_kwargs=[:]) {
         "stats": true,
         "required_sample_types": [],
         "watch_path": false,
-        "per_read_stats": false
+        "per_read_stats": false,
+        "allow_multiple_basecall_models": false,
     ]
     ArgumentParser parser = new ArgumentParser(
         args: required_args,
@@ -1033,6 +1066,7 @@ Map create_metamap(Map arguments) {
             "barcode": null,
             "type": "test_sample",
             "run_ids": [],
+            "basecall_models": [],
         ],
         name: "create_metamap",
     )
@@ -1077,9 +1111,31 @@ def get_sample_sheet(Path sample_sheet, ArrayList required_sample_types) {
     // concat the channel holding the path to the sample sheet to `ch_err` and call
     // `.last()` to make sure that the error-checking closure above executes before
     // emitting values from the CSV
-    return ch_err.concat(Channel.fromPath(sample_sheet)).last().splitCsv(
+    ch_sample_sheet = ch_err.concat(Channel.fromPath(sample_sheet)).last().splitCsv(
         header: true, quote: '"'
     )
+    // in case there is an 'analysis_group' column, we need to define a `groupKey` to
+    // allow for non-blocking calls of `groupTuple` later (on the values in the
+    // 'analysis_group' column); we first collect the sample sheet in a single list of
+    // maps and then count the occurrences of each group before using these to create
+    // the `groupKey` objects; note that the below doesn't do anything if there is no
+    // 'analysis_group' column
+    ch_group_counts = ch_sample_sheet
+    | collect
+    | map { rows -> rows.collect { it.analysis_group } .countBy { it } }
+
+    // now we `combine` the analysis group counts with the sample sheet channel and add
+    // the `groupKey` to the entries
+    ch_sample_sheet = ch_sample_sheet
+    | combine(ch_group_counts)
+    | map { row, group_counts ->
+        if (row.analysis_group) {
+            int counts = group_counts[row.analysis_group]
+            row = row + [analysis_group: groupKey(row.analysis_group, counts)]
+        }
+        row
+    }
+    return ch_sample_sheet
 }
 
 

nextflow.config

@@ -106,7 +106,7 @@ params {
       ]
     agent = null
     container_sha = "shae7c9f184996a384e99be68e790f0612f0c732867"
-	common_sha = "sha338caea0a2532dc0ea8f46638ccc322bb8f9af48"
+	common_sha = "sha8b5843d549bb210558cbb676fe537a153ce771d6"
     }
 }