diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 26625ea..a98f0ad 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/README.md b/README.md index 6077c64..650428e 100644 --- a/README.md +++ b/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 -[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 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. -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 workflow can currently be run using either +[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 -tar -xzvf differential_expression.tar.gz +nextflow pull epi2me-labs/wf-transcriptomes ``` -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 --ref_genome differential_expression/hg38_chr20.fa \ ---transcriptome-source reference-guided \ ---ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \ ---direct_rna --minimap2_index_opts '-k 15' --sample_sheet differential_expression/sample_sheet.csv \ ---jaffal_refBase differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \ --profile standard + --de_analysis \ + --direct_rna \ + --fastq 'wf-transcriptomes-demo/differential_expression_fastq' \ + --jaffal_annotation 'genCode22' \ + --jaffal_genome 'hg38_chr20' \ + --jaffal_refBase 'wf-transcriptomes-demo/chr20' \ + --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/ + diff --git a/bin/workflow_glue/check_sample_sheet.py b/bin/workflow_glue/check_sample_sheet.py index 7fa1605..8fec3f0 100755 --- a/bin/workflow_glue/check_sample_sheet.py +++ b/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}.") diff --git a/bin/workflow_glue/configure_igv.py b/bin/workflow_glue/configure_igv.py index 45b4ebd..db158a4 100755 --- a/bin/workflow_glue/configure_igv.py +++ b/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 = [] - bam_indices = [] + xams = [] + 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"): - bams.append(fname) - elif fname.endswith(".bai"): - bam_indices.append(fname) + elif fname.endswith(".bam") or fname.endswith(".cram"): + xams.append(fname) + elif fname.endswith(".bai") or fname.endswith(".crai"): + 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 len(bams) != len(bam_indices): - raise ValueError("Got different number of BAM and BAM index files.") + if xam_indices: + if len(xams) != len(xam_indices): + 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 len(bams) != len(vcfs): - raise ValueError("Got different number of BAM and VCF files.") - # if we got BAM or VCF indices, pair them up with their corresponding files (and + if xams and vcfs: + if len(xams) != len(vcfs): + raise ValueError("Got different number of XAM and VCF files.") + # 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 bai: name of BAM index file + :param xam: name of XAM file to be displayed + :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", - "url": bam, + "format": xam.split(".")[-1], + "url": xam, } - # add the BAM index if present - if bai is not None: - alignment_track_dict["indexURL"] = bai + # add the XAM index if present + if xai is not None: + 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: - with open(args.extra_bam_opts, "r") as f: + if args.extra_alignment_opts is not None: + 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: - with open(args.extra_vcf_opts, "r") as f: + if args.extra_variant_opts is not None: + 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( - vcfs_with_indices, bams_with_indices, fillvalue=(None, None) + for (vcf, vcf_index), (xam, xam_index) in zip_longest( + 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: - # add an alignment track for the BAM + if xam is not None: + # 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 diff --git a/bin/workflow_glue/get_ds_records.py b/bin/workflow_glue/get_ds_records.py deleted file mode 100755 index 3f95336..0000000 --- a/bin/workflow_glue/get_ds_records.py +++ /dev/null @@ -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 diff --git a/docs/04_install_and_run.md b/docs/04_install_and_run.md index a8abf2b..57b482c 100644 --- a/docs/04_install_and_run.md +++ b/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 -[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 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. -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 workflow can currently be run using either +[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 -tar -xzvf differential_expression.tar.gz +nextflow pull epi2me-labs/wf-transcriptomes ``` -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 --ref_genome differential_expression/hg38_chr20.fa \ ---transcriptome-source reference-guided \ ---ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \ ---direct_rna --minimap2_index_opts '-k 15' --sample_sheet differential_expression/sample_sheet.csv \ ---jaffal_refBase differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \ --profile standard + --de_analysis \ + --direct_rna \ + --fastq 'wf-transcriptomes-demo/differential_expression_fastq' \ + --jaffal_annotation 'genCode22' \ + --jaffal_genome 'hg38_chr20' \ + --jaffal_refBase 'wf-transcriptomes-demo/chr20' \ + --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/ \ No newline at end of file +For further information about running a workflow on +the command line see https://labs.epi2me.io/wfquickstart/ diff --git a/lib/common.nf b/lib/common.nf index 9dfd5a7..2a31d49 100644 --- a/lib/common.nf +++ b/lib/common.nf @@ -23,33 +23,33 @@ process configure_igv { // extensions path "file-names.txt" val locus_str - val bam_extra_opts - val vcf_extra_opts + val aln_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 = \ - bam_extra_opts ? new JsonBuilder(bam_extra_opts).toPrettyString() : "" - String bam_extra_opts_arg = \ - bam_extra_opts ? "--extra-bam-opts bam-extra-opts.json" : "" + def aln_opts_json_str = \ + aln_extra_opts ? new JsonBuilder(aln_extra_opts).toPrettyString() : "" + String aln_extra_opts_arg = \ + aln_extra_opts ? "--extra-alignment-opts extra-aln-opts.json" : "" // extra options for variant tracks - def vcf_opts_json_str = \ - vcf_extra_opts ? new JsonBuilder(vcf_extra_opts).toPrettyString() : "" - String vcf_extra_opts_arg = \ - vcf_extra_opts ? "--extra-vcf-opts vcf-extra-opts.json" : "" + def var_opts_json_str = \ + var_extra_opts ? new JsonBuilder(var_extra_opts).toPrettyString() : "" + String var_extra_opts_arg = \ + 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 '$vcf_opts_json_str' > vcf-extra-opts.json + echo '$aln_opts_json_str' > extra-aln-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 \ - $vcf_extra_opts_arg \ + $aln_extra_opts_arg \ + $var_extra_opts_arg \ > igv.json """ } diff --git a/lib/ingress.nf b/lib/ingress.nf index b03afe7..4ba5ac9 100644 --- a/lib/ingress.nf +++ b/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 - * to the stats dir is `null`, add an empty list. + * run IDs and basecall model, from the `run_ids` and `basecaller` files in the stats + * 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 ingressed_run_ids = new HashSet() // 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 = meta + [run_ids: run_ids, basecall_models: basecall_models] } - // `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, 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_run_IDs_to_meta(ch_spread_result), "fastq") + // add number of reads, run IDs, and basecall models to meta + 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, - ds_runids: ds_runids_env.tokenize(','), - ds_basecall_models: ds_basecaller_env.tokenize(','), - ], - paths - ] + [meta + [is_unaligned: is_unaligned, is_sorted: is_sorted], 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_run_IDs_to_meta(ch_spread_result), "fastq") + // add number of reads, run IDs, and basecall models to meta + 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,10 +560,20 @@ 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 - awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {print \$ix["run_id"]}' \ - fastcat_stats/per-file-runids.txt | sort | uniq > fastcat_stats/run_ids - """ + # get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace + # as field delimiter per default and thus ignores empty columns) + 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 + """ } process checkBamHeaders { @@ -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 - awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {print \$ix["run_id"]}' \ - bamstats_results/bamstats.runids.txt | sort | uniq > bamstats_results/run_ids + # get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace + # as field delimiter otherwise and thus ignore empty columns) + 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 } diff --git a/nextflow.config b/nextflow.config index 1d37ca2..5f958d1 100644 --- a/nextflow.config +++ b/nextflow.config @@ -106,7 +106,7 @@ params { ] agent = null container_sha = "shae7c9f184996a384e99be68e790f0612f0c732867" - common_sha = "sha338caea0a2532dc0ea8f46638ccc322bb8f9af48" + common_sha = "sha8b5843d549bb210558cbb676fe537a153ce771d6" } }