kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Template update

Browse Source
This commit is contained in:
Neil Horner 2024-06-27 16:13:55 +01:00
parent 1fed562746
commit c662293317
9 changed files with 265 additions and 307 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.pre-commit-config.yaml
View File

@ -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

72
README.md
View File

@ -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/

13
bin/workflow_glue/check_sample_sheet.py
View File

@ -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}.")

80
bin/workflow_glue/configure_igv.py
View File

@ -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

166
bin/workflow_glue/get_ds_records.py
View File

@ -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

71
docs/04_install_and_run.md
View File

@ -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/
For further information about running a workflow on
the command line see https://labs.epi2me.io/wfquickstart/

28
lib/common.nf
View File

@ -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
"""
}

130
lib/ingress.nf
View File

@ -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<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_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,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
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
"""
}
@ -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
}

2
nextflow.config
View File

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