diff --git a/CHANGELOG.md b/CHANGELOG.md index f856670..5bf3ad2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## v1.0.1dev - [date] +### `Added` + +- [#66](https://github.com/nf-core/seqsubmit/pull/66) + - Added co-assembly submission to `ASSEMBLYSUBMIT`: a samplesheet row can list several semicolon-separated run accessions (and matching read files for coverage calculation) + - When the source runs belong to different ENA samples, a virtual sample with a `composed of` attribute is registered and the co-assembly is submitted under it + - Optional `tax_id` samplesheet column to set the taxon of the virtual co-assembly sample when the source samples have different taxa + ### `Fixed` - [#117](https://github.com/nf-core/seqsubmit/pull/117) and [#119](https://github.com/nf-core/seqsubmit/pull/119) diff --git a/README.md b/README.md index ec93644..6e8a48b 100644 --- a/README.md +++ b/README.md @@ -74,8 +74,11 @@ Example: id,fasta,fastq_1,fastq_2,coverage,run_accession,assembler,assembler_version,assembly_type assembly_1,data/contigs_1.fasta.gz,data/reads_1.fastq.gz,data/reads_2.fastq.gz,,ERR011322,SPAdes,3.15.5,metagenome assembly_2,data/contigs_2.fasta.gz,,,42.7,ERR011323,MEGAHIT,1.2.9,metatranscriptome +coassembly_1,data/contigs_3.fasta.gz,data/run1_R1.fastq.gz;data/run2_R1.fastq.gz,data/run1_R2.fastq.gz;data/run2_R2.fastq.gz,,ERR011324;ERR011325,MEGAHIT,1.2.9, ``` +Co-assemblies are supported: list all source runs in `run_accession` and their read files in `fastq_1`/`fastq_2`, separated by semicolons. + See the [`metagenomic_assemblies` mode section](docs/usage.md#samplesheet-input-1) of the usage documentation for more details. ### `mags` and `bins` modes diff --git a/assets/schema_input_assembly.json b/assets/schema_input_assembly.json index a269c48..4bf1ae7 100644 --- a/assets/schema_input_assembly.json +++ b/assets/schema_input_assembly.json @@ -26,9 +26,7 @@ "anyOf": [ { "type": "string", - "format": "file-path", - "exists": true, - "pattern": "^\\S+\\.(fq|fastq)(\\.gz)?$" + "pattern": "^\\S+\\.(fq|fastq)(\\.gz)?(;\\S+\\.(fq|fastq)(\\.gz)?)*$" }, { "type": "string", @@ -37,15 +35,13 @@ ], "default": null, "errorMessage": "FASTQ file must have extension '.fq' or '.fastq' (optionally gzipped)", - "description": "Path to single-end/paired-end forward reads in FASTQ format (optionally gzipped). Required if coverage is not provided" + "description": "Path to single-end/paired-end forward reads in FASTQ format (optionally gzipped). Required if coverage is not provided. For co-assemblies, semicolon-separated paths (e.g. run1_R1.fq.gz;run2_R1.fq.gz)" }, "fastq_2": { "anyOf": [ { "type": "string", - "format": "file-path", - "exists": true, - "pattern": "^\\S+\\.(fq|fastq)(\\.gz)?$" + "pattern": "^\\S+\\.(fq|fastq)(\\.gz)?(;\\S+\\.(fq|fastq)(\\.gz)?)*$" }, { "type": "string", @@ -54,7 +50,7 @@ ], "default": null, "errorMessage": "FASTQ file for reverse reads must have extension '.fq' or '.fastq' (optionally gzipped)", - "description": "Path to reverse reads in FASTQ format for paired-end data (optionally gzipped). Leave empty for single-end reads" + "description": "Path to reverse reads in FASTQ format for paired-end data (optionally gzipped). Leave empty for single-end reads. For co-assemblies, use semicolon-separated paths (e.g. run1_R2.fq.gz;run2_R2.fq.gz)" }, "coverage": { "anyOf": [ @@ -76,7 +72,7 @@ "type": "string", "pattern": "^\\S+$", "errorMessage": "Accession must be provided and cannot contain spaces", - "description": "ENA run accession for the reads used to generate the metagenomic assembly. Reads must already be submitted to ENA", + "description": "ENA run accession(s) for the reads used to generate the metagenomic assembly. Reads must already be submitted to ENA. For co-assemblies, use semicolon-separated values (e.g. ERR000001;ERR000002)", "meta": ["run_accession"] }, "assembler": { @@ -100,6 +96,22 @@ "description": "Assembly type: metagenome or metatranscriptome. Written to the ENA manifest as ASSEMBLY_TYPE 'primary metagenome' or 'metatranscriptome'. Default: metagenome", "default": "metagenome", "meta": ["assembly_type"] + }, + "tax_id": { + "anyOf": [ + { + "type": "integer", + "exclusiveMinimum": 0 + }, + { + "type": "string", + "maxLength": 0 + } + ], + "default": null, + "errorMessage": "Taxon ID must be a positive integer or empty", + "description": "NCBI taxon ID for the virtual sample registered for a co-assembly (optional). Provide when the source samples have different taxon IDs. Ignored for assemblies that are not co-assemblies of several samples", + "meta": ["tax_id"] } }, "required": ["id", "fasta", "run_accession", "assembler", "assembler_version"], diff --git a/bin/register_coassembly_sample.py b/bin/register_coassembly_sample.py new file mode 100755 index 0000000..436a4e6 --- /dev/null +++ b/bin/register_coassembly_sample.py @@ -0,0 +1,566 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import csv +import hashlib +import logging +import os +import time +import xml.etree.ElementTree as ET +import re + +import requests +from requests.auth import HTTPBasicAuth +from retry import retry + +logger = logging.getLogger("register_coassembly_sample") + + +PORTAL_SEARCH_API = "https://www.ebi.ac.uk/ena/portal/api/search" + +CHECKLIST = "ERC000011" +CHECKLIST_XML_API = f"https://www.ebi.ac.uk/ena/browser/api/xml/{CHECKLIST}" +TAXONOMY_API = "https://www.ebi.ac.uk/ena/taxonomy/rest/tax-id" + +# ENA Portal API fields retrieved for runs and samples +RUN_SAMPLE_FIELDS = ["sample_accession", "secondary_sample_accession"] +SAMPLE_METADATA_FIELDS = ["country", "collection_date", "tax_id", "scientific_name"] + +TEST_SUBMIT_URL = "https://wwwdev.ebi.ac.uk/ena/submit/drop-box/submit/" +PROD_SUBMIT_URL = "https://www.ebi.ac.uk/ena/submit/drop-box/submit/" + +INSDC_BIOSAMPLE_ACCESSION_REGEX = re.compile( + r"^SAM(N|D|E|EA)[0-9]{7,8}$" +) +INSDC_SECONDARY_SAMPLE_ACCESSION_REGEX = re.compile( + r"^(E|D|S)RS[0-9]{6,}$" +) +EXISTING_ACCESSION_IN_ERROR_REGEX = re.compile( + r'accession:\s*"([A-Z]+[0-9]+)"', + re.IGNORECASE, +) + +DEFAULT_NA_VALUE = "not provided" + + +class CoassemblyRegistrationError(RuntimeError): + """Raised for hard validation/submission errors.""" + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Register co-assembly samples.") + parser.add_argument("--input", required=True, help="Input assembly metadata CSV") + parser.add_argument("--output", required=True, help="Output assembly metadata CSV") + parser.add_argument( + "--test", + action="store_true", + help="Register sample using the test submission endpoint", + ) + parser.add_argument( + "--private", + action="store_true", + help="Source runs/samples are private: query the ENA Portal API with Webin auth (ENA_WEBIN/ENA_WEBIN_PASSWORD)", + ) + parser.add_argument( + "--debug", + action="store_true", + help="Enable verbose debug logging", + ) + parser.add_argument("--default-country", type=str, help="Default country to use if values differ or are missing.") + parser.add_argument("--default-date", type=str, help="Default collection date to use if values differ or are missing.") + parser.add_argument("--default-taxid", type=str, help="Default taxon ID to use if values differ or are missing. The scientific name is fetched from ENA taxonomy.") + return parser.parse_args() + + +def setup_logging(debug: bool) -> None: + """Configure logging level and format based on CLI flags.""" + level = logging.DEBUG if debug else logging.INFO + logging.basicConfig(level=level, format="%(levelname)s: %(message)s", force=True) + if debug: + logging.getLogger("requests").setLevel(logging.INFO) + logging.getLogger("urllib3").setLevel(logging.INFO) + logging.getLogger("requests.packages.urllib3").setLevel(logging.INFO) + logger.debug("Debug logging enabled") + + +def get_credentials() -> tuple[str, str]: + username = os.environ.get("ENA_WEBIN", "").strip() + password = os.environ.get("ENA_WEBIN_PASSWORD", "").strip() + if not username or not password: + raise CoassemblyRegistrationError( + "Missing credentials. Set ENA_WEBIN/ENA_WEBIN_PASSWORD." + ) + return username, password + + +@retry(tries=3, delay=2, backoff=2) +def get_xml(url: str) -> ET.Element: + logger.debug(f"Requesting XML URL: {url}") + response = requests.get(url, timeout=60) + response.raise_for_status() + return ET.fromstring(response.content) + +@retry(tries=3, delay=2, backoff=2) +def get_json(url: str, params: dict[str, str], auth: HTTPBasicAuth | None = None) -> list[dict[str, str]]: + prepared = requests.Request("GET", url, params=params).prepare() + logger.debug(f"Requesting JSON URL ({'with' if auth else 'without'} Webin auth): {prepared.url}") + response = requests.get(url, params=params, auth=auth, timeout=60) + response.raise_for_status() + if not response.text.strip(): + return [] + result = response.json() + if isinstance(result, dict) and "message" in result: + raise CoassemblyRegistrationError(f"ENA Portal API error: {result['message']}") + return result + + +def search_portal(result: str, query: str, fields: list[str], private: bool = False) -> list[dict[str, str]]: + """ + Query the ENA Portal API. + + When `private` is set, the request is authenticated with the Webin account + (ENA_WEBIN/ENA_WEBIN_PASSWORD), which also returns private records owned by it. + """ + params = {"result": result, "query": query, "fields": ",".join(fields), "format": "json"} + if not private: + return get_json(PORTAL_SEARCH_API, params=params) + + username, password = get_credentials() + try: + return get_json(PORTAL_SEARCH_API, params=params, auth=HTTPBasicAuth(username, password)) + except requests.HTTPError as exc: + # ENA returns a 500 (not 401) for invalid Webin credentials, so a 500 is ambiguous + raise CoassemblyRegistrationError( + f"Authenticated ENA Portal API search failed ({exc}). This may be a temporary ENA API " + "issue, or invalid Webin credentials - check ENA_WEBIN/ENA_WEBIN_PASSWORD." + ) from exc + +def to_pretty_xml(root: ET.Element) -> str: + """Serialize an XML element to a pretty-printed UTF-8 XML string.""" + # Re-parse to avoid side effects from in-place indentation on caller-owned trees. + normalized_root = ET.fromstring(ET.tostring(root, encoding="utf-8")) + ET.indent(normalized_root, space=" ") + return ET.tostring(normalized_root, encoding="utf-8", xml_declaration=True).decode("utf-8") + +@retry(tries=3, delay=2, backoff=2) +def submit_sample_xml(auth: HTTPBasicAuth, submit_url: str, sample_xml: str) -> ET.Element: + submission_xml = """\n\n""" + files = { + "SUBMISSION": ("submission.xml", submission_xml, "application/xml"), + "SAMPLE": ("sample.xml", sample_xml, "application/xml"), + } + response = requests.post(submit_url, files=files, auth=auth, timeout=120) + response.raise_for_status() + return ET.fromstring(response.content) + +def merge_or_not_provided(values: list[str], default_value=None) -> str: + """Merge values or return default_value or DEFAULT_NA_VALUE.""" + unique_values = sorted(set(values)) + logger.info(f"Merging values: {unique_values}") + if len(unique_values) == 1 and unique_values[0]: + logger.info(f"Merge decision: using unique value '{unique_values[0]}'") + return unique_values[0] + elif default_value: + logger.info(f"Merge decision: using default value '{default_value}' to overwrite samples metadata") + return default_value + else: + logger.info(f"Merge decision: '{DEFAULT_NA_VALUE}' (multiple conflicting values or some missing)") + return DEFAULT_NA_VALUE + +def get_run_sample(run_accession: str, private: bool = False) -> str: + """ + Retrieve the sample accession linked to a run from the ENA Portal API. + Use Webin auth to retrieve private records. + + Args: + run_accession (str): The run accession to query. + private (bool): Query with Webin auth to access private data. + + Returns: + str: The secondary sample accession (e.g. ERS/SRS), or the BioSample + accession if no secondary accession is available. + + Raises: + CoassemblyRegistrationError: If the run is not found or has no linked sample. + """ + rows = search_portal( + result="read_run", + query=f'"run_accession={run_accession}"', + fields=RUN_SAMPLE_FIELDS, + private=private, + ) + if not rows: + raise CoassemblyRegistrationError( + f"Run {run_accession} not found in ENA" + ("" if private else " (use --private for private data)") + ) + + row = rows[0] + sample_accession = (row.get("secondary_sample_accession") or row.get("sample_accession") or "").strip() + if not sample_accession: + raise CoassemblyRegistrationError(f"Run {run_accession} has no linked sample accession in ENA") + + return sample_accession + + +def get_sample_metadata(sample_accession: str, private: bool = False) -> tuple[str, str, str, str]: + """ + Retrieve metadata for a given sample accession from the ENA portal. + + Args: + sample_accession (str): The sample accession to query. + private (bool): Query with Webin auth to access private data. + + Returns: + tuple[str, str, str, str]: A tuple containing tax_id, scientific_name, country, and collection_date. + """ + if INSDC_BIOSAMPLE_ACCESSION_REGEX.match(sample_accession): + query = f'"sample_accession={sample_accession}"' + elif INSDC_SECONDARY_SAMPLE_ACCESSION_REGEX.match(sample_accession): + query = f'"secondary_sample_accession={sample_accession}"' + else: + raise CoassemblyRegistrationError(f"Unrecognised sample accession format: {sample_accession}") + + rows = search_portal( + result="sample", + query=query, + fields=SAMPLE_METADATA_FIELDS, + private=private, + ) + + if not rows: + logger.warning(f"No portal metadata found for sample {sample_accession}") + return "", "", "", "" + + row = rows[0] + return ( + row.get("tax_id"), + row.get("scientific_name"), + row.get("country"), + row.get("collection_date"), + ) + +@retry(requests.RequestException, tries=3, delay=2, backoff=2) +def get_scientific_name(tax_id: str) -> str: + """ + Retrieve the scientific name for a taxon ID from the ENA taxonomy API. + + Args: + tax_id (str): The NCBI taxon ID. + + Returns: + str: The scientific name of the taxon. + + Raises: + CoassemblyRegistrationError: If the taxon is not found or cannot be used for submission. + """ + url = f"{TAXONOMY_API}/{tax_id}" + logger.debug(f"Requesting taxonomy URL: {url}") + response = requests.get(url, timeout=60) + if response.status_code == 404: + raise CoassemblyRegistrationError(f"Taxon ID {tax_id} not found in ENA taxonomy") + response.raise_for_status() + taxon = response.json() + if taxon.get("submittable") != "true": + raise CoassemblyRegistrationError(f"Taxon ID {tax_id} ({taxon.get('scientificName')}) is not submittable to ENA") + return taxon["scientificName"] + +def get_allowed_countries() -> set[str]: + """ + Retrieve the set of allowed countries from the ENA checklist XML. + + Returns: + set[str]: A set of allowed country/sea values. + """ + root = get_xml(CHECKLIST_XML_API) + values = set() + field_nodes = root.findall(".//FIELD") + for field in field_nodes: + label = field.findtext("LABEL") + name = field.findtext("NAME") + if label == "geographic location (country and/or sea)" or name == "geographic_location_country_andor_sea": + for entry in field.findall(".//TEXT_VALUE/VALUE"): + if entry.text: + values.add(entry.text.strip()) + break + return values + +def build_sample_xml( + alias: str, + title: str, + description: str, + taxon_id: str, + scientific_name: str, + composed_of: str, + collection_date: str, + geographic_location: str, +) -> str: + """ + Build the XML representation of a sample for ENA submission. + + Args: + alias (str): The sample alias. + title (str): The sample title. + description (str): The sample description. + taxon_id (str): The taxon ID of the sample. + scientific_name (str): The scientific name of the sample. + composed_of (str): The list of source samples. + collection_date (str): The collection date of the sample. + geographic_location (str): The geographic location of the sample. + + Returns: + str: The XML string representation of the sample. + """ + root = ET.Element("SAMPLE_SET") + sample = ET.SubElement(root, "SAMPLE", {"alias": alias}) + + title_node = ET.SubElement(sample, "TITLE") + title_node.text = title + + sample_name = ET.SubElement(sample, "SAMPLE_NAME") + taxon = ET.SubElement(sample_name, "TAXON_ID") + taxon.text = taxon_id + sci = ET.SubElement(sample_name, "SCIENTIFIC_NAME") + sci.text = scientific_name + + description_node = ET.SubElement(sample, "DESCRIPTION") + description_node.text = description + + attrs = ET.SubElement(sample, "SAMPLE_ATTRIBUTES") + + def add_attr(tag: str, value: str) -> None: + item = ET.SubElement(attrs, "SAMPLE_ATTRIBUTE") + tag_node = ET.SubElement(item, "TAG") + tag_node.text = tag + value_node = ET.SubElement(item, "VALUE") + value_node.text = value + + add_attr("ENA-CHECKLIST", CHECKLIST) + add_attr("organism", scientific_name) + add_attr("collection date", collection_date) + add_attr("composed of", composed_of) + add_attr("geographic location (country and/or sea)", geographic_location) + add_attr("scientific_name", scientific_name) + + return to_pretty_xml(root) + + +def build_virtual_sample_alias(source_samples: list[str], max_length: int = 50, test: bool = False) -> str: + """Build deterministic alias from source samples with an md5 suffix. + + When `test` is set, a short timestamp-based hash is appended so repeated + test submissions don't collide on alias uniqueness. + """ + sorted_samples = sorted(source_samples) + samples_hash = hashlib.md5(",".join(sorted_samples).encode("utf-8")).hexdigest()[:8] + + alias = f"coassembly_{len(sorted_samples)}_samples_{samples_hash}" + + if test: + timestamp_hash = hashlib.md5(str(time.time()).encode("utf-8")).hexdigest()[:6] + alias = f"{alias}_{timestamp_hash}" + + if len(alias) <= max_length: + return alias + + raise CoassemblyRegistrationError( + f"Hash suffix '{alias}' exceeds maximum alias length of {max_length}" + ) + + +def extract_accession_from_receipt(receipt: ET.Element) -> str: + """ + Extract the sample accession from the ENA submission receipt XML. + + Args: + receipt (ET.Element): The root element of the receipt XML. + + Returns: + str: The sample accession. + + Raises: + CoassemblyRegistrationError: If the receipt does not contain a sample accession. + """ + success = receipt.attrib.get("success", "false").lower() == "true" + if not success: + messages = [m.text for m in receipt.findall(".//MESSAGES/*") if m.text] + message_text = "; ".join(messages) if messages else "no error details in receipt" + + # If this is a duplicate submission, ENA returns the existing accession. + existing_match = EXISTING_ACCESSION_IN_ERROR_REGEX.search(message_text) + if existing_match and "already exists" in message_text.lower(): + existing_accession = existing_match.group(1) + logger.info(f"Virtual sample already exists; reusing accession: {existing_accession}") + return existing_accession + + raise CoassemblyRegistrationError( + f"ENA sample submission failed: {message_text}" + ) + + sample_node = receipt.find(".//SAMPLE") + if sample_node is not None: + ena_accession = sample_node.get("accession", "") + if ena_accession: + return ena_accession + + raise CoassemblyRegistrationError("ENA sample submission succeeded but no SAMPLE accession found in receipt XML") + + +def register_virtual_sample(run_accessions: list[str], test: bool, private: bool = False, default_country=None, default_date=None, default_taxid=None) -> str: + """ + Register a virtual sample for co-assemblies based on run accessions. + + Args: + run_accessions (list[str]): A list of run accessions for the co-assembly. + test (bool): Whether to use the test submission endpoint. + private (bool): Whether source runs/samples are private (query ENA with Webin auth). + + Returns: + str: The sample accession of the registered virtual sample. + + Raises: + CoassemblyRegistrationError: If validation or submission fails. + """ + if test: + submit_url = TEST_SUBMIT_URL + else: + submit_url = PROD_SUBMIT_URL + + sample_accessions = [] + taxon_ids = [] + scientific_names = [] + countries = [] + collection_dates = [] + + for run in run_accessions: + sample_acc = get_run_sample(run, private=private) + logger.info(f"Run {run} is linked to sample {sample_acc}") + sample_accessions.append(sample_acc) + + taxon_id, scientific_name, country, collection_date = get_sample_metadata(sample_acc, private=private) + logger.info(f"Metadata for sample {sample_acc}: tax_id={taxon_id}, scientific_name={scientific_name}, country={country}, collection_date={collection_date}") + if not taxon_id: + raise CoassemblyRegistrationError(f"Sample {sample_acc} has no tax_id in ENA portal response") + if not scientific_name: + raise CoassemblyRegistrationError(f"Sample {sample_acc} has no scientific_name in ENA portal response") + taxon_ids.append(taxon_id) + scientific_names.append(scientific_name) + countries.append(country) + collection_dates.append(collection_date) + + unique_samples = sorted(set(sample_accessions)) + if len(unique_samples) == 1: + logger.info(f"All source runs belong to one sample ({unique_samples[0]}); keeping original row unchanged") + return "" + + merged_collection_date = merge_or_not_provided(collection_dates, default_date) + merged_country = merge_or_not_provided(countries, default_country) + merged_taxid = merge_or_not_provided(taxon_ids, default_taxid) + + if merged_taxid == DEFAULT_NA_VALUE: + sample_taxa = dict(zip(sample_accessions, zip(taxon_ids, scientific_names))) + details = "\n".join( + f" {sample}: tax_id={taxid} ({name})" + for sample, (taxid, name) in sample_taxa.items() + ) + raise CoassemblyRegistrationError( + "Co-assembly source samples must share the same tax_id:\n" + f"{details}\n" + "Set the tax_id to register the co-assembly sample under with --default-taxid " + "(the 'tax_id' column of the samplesheet when running the pipeline)." + ) + + merged_tax_name = get_scientific_name(merged_taxid) + logger.info(f"Scientific name for tax_id {merged_taxid}: {merged_tax_name}") + + allowed_countries = get_allowed_countries() + if merged_country != DEFAULT_NA_VALUE and merged_country not in allowed_countries: + logger.warning( + f"Country '{merged_country}' is not in {CHECKLIST} allowed country/sea values; setting to '{DEFAULT_NA_VALUE}'" + ) + merged_country = DEFAULT_NA_VALUE + + sample_list = ",".join(unique_samples) + title = f"Combined sample from {sample_list}" + description = ( + f"This is a virtual sample representing a co-assembly of raw reads from multiple " + f"{merged_tax_name} samples. Co-assembly was performed from runs of {len(unique_samples)} samples {sample_list}." + ) + alias = build_virtual_sample_alias(unique_samples, test=test) + logger.info(f"Registering virtual sample with alias '{alias}' for co-assembly of samples: {sample_list}") + + sample_xml = build_sample_xml( + alias=alias, + title=title, + description=description, + taxon_id=merged_taxid, + scientific_name=merged_tax_name, + composed_of=sample_list, + collection_date=merged_collection_date, + geographic_location=merged_country, + ) + + with open(f"{alias}.xml", "w", encoding="utf-8") as sample_file: + sample_file.write(sample_xml) + + username, password = get_credentials() + receipt = submit_sample_xml(HTTPBasicAuth(username, password), submit_url, sample_xml) + logger.debug(f"ENA receipt XML:\n{ET.tostring(receipt, encoding='unicode')}") + virtual_sample = extract_accession_from_receipt(receipt) + logger.info(f"Registered virtual co-assembly sample: {virtual_sample}") + return virtual_sample + + +def main() -> int: + args = parse_args() + setup_logging(args.debug) + + logger.debug(f"Starting coassembly registration script: input={args.input}, output={args.output}, test={args.test}") + + if args.test: + logger.info("Running in test mode: submitting to the ENA test server. Submission will be discarded after 24 hours.") + + with open(args.input, newline="", encoding="utf-8") as handle: + reader = csv.DictReader(handle) + fieldnames = reader.fieldnames or [] + if "Runs" not in fieldnames: + raise CoassemblyRegistrationError("Input CSV must have a 'Runs' column with ENA run accessions") + rows = list(reader) + + # The script may add/update the Sample value; ensure the output schema includes it + if "Sample" not in fieldnames: + fieldnames.append("Sample") + + updated_rows: list[dict[str, str]] = [] + for row in rows: + run_accessions = row["Runs"].split(",") + + logger.info(f"Processing assembly with runs: {run_accessions}") + + # Skip non-coassembly rows by design + if len(run_accessions) <= 1: + logger.info("Single run found; skipping coassembly sample registration for this row") + updated_rows.append(row) + continue + + virtual_sample = register_virtual_sample( + run_accessions=run_accessions, + test=args.test, + private=args.private, + default_country=args.default_country, + default_date=args.default_date, + default_taxid=args.default_taxid, + ) + if virtual_sample or virtual_sample == "": + row["Sample"] = virtual_sample + updated_rows.append(row) + + logger.info(f"Writing updated assembly metadata to: {args.output}") + with open(args.output, "w", newline="", encoding="utf-8") as out_handle: + writer = csv.DictWriter(out_handle, fieldnames=fieldnames) + writer.writeheader() + writer.writerows(updated_rows) + + +if __name__ == "__main__": + main() diff --git a/conf/modules.config b/conf/modules.config index c14d7f2..c0ee3be 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -133,6 +133,7 @@ process { } withName: 'COVERM_CONTIG' { + ext.args = '--methods mean length' publishDir = [ path: { "${params.outdir}/${params.mode}/${meta.id}/coverage" }, mode: params.publish_dir_mode, @@ -182,14 +183,14 @@ process { stageInMode = 'symlink' } - withName: 'REGISTERSTUDY|CREATE_READS_MANIFEST' { + withName: 'REGISTERSTUDY|CREATE_READS_MANIFEST|REGISTER_COASSEMBLY_SAMPLE' { publishDir = [ enabled: false ] } withName: 'GENERATE_ASSEMBLY_MANIFEST' { - errorStrategy = { task.exitStatus in (1 + (130..145) + 104 + (175..177)) ? 'retry' : 'finish' } + errorStrategy = { task.exitStatus in ((130..145) + 104 + 1 + (175..177)) ? 'retry' : 'finish' } maxRetries = 3 publishDir = [ enabled: false diff --git a/conf/test_assembly_coassembly_no_coverage.config b/conf/test_assembly_coassembly_no_coverage.config new file mode 100644 index 0000000..41d1e46 --- /dev/null +++ b/conf/test_assembly_coassembly_no_coverage.config @@ -0,0 +1,35 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Nextflow config file for running minimal tests +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Defines input files and everything required to run a fast and simple pipeline test. + + Use as follows: + nextflow run nf-core/seqsubmit -profile test_assembly_coassembly_no_coverage, --outdir + +---------------------------------------------------------------------------------------- +*/ + +process { + resourceLimits = [ + cpus: 2, + memory: '8.GB', + time: '1.h' + ] +} + +params { + config_profile_name = 'Test --mode metagenomic_assemblies co-assembly without coverage profile' + config_profile_description = 'Co-assembly test with missing coverage and paired-end reads from two runs' + + // Input data + input = params.pipelines_testdata_base_path + 'seqsubmit/samplesheets/samplesheet_assembly_coassembly_no_coverage.csv' + outdir = 'test_output' + + mode = "metagenomic_assemblies" + submission_study = null + study_metadata = params.pipelines_testdata_base_path + 'seqsubmit/test_data/study_metadata/study_metadata.json' + centre_name = "TEST_CENTER" + + test_upload = true +} diff --git a/conf/test_assembly_coassembly_with_coverage.config b/conf/test_assembly_coassembly_with_coverage.config new file mode 100644 index 0000000..f688999 --- /dev/null +++ b/conf/test_assembly_coassembly_with_coverage.config @@ -0,0 +1,35 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Nextflow config file for running minimal tests +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Defines input files and everything required to run a fast and simple pipeline test. + + Use as follows: + nextflow run nf-core/seqsubmit -profile test_assembly_coassembly_with_coverage, --outdir + +---------------------------------------------------------------------------------------- +*/ + +process { + resourceLimits = [ + cpus: 2, + memory: '8.GB', + time: '1.h' + ] +} + +params { + config_profile_name = 'Test --mode metagenomic_assemblies co-assembly with coverage profile' + config_profile_description = 'Co-assembly test with pre-provided coverage and multiple run accessions' + + // Input data + input = params.pipelines_testdata_base_path + 'seqsubmit/samplesheets/samplesheet_assembly_coassembly_with_coverage.csv' + outdir = 'test_output' + + mode = "metagenomic_assemblies" + submission_study = null + study_metadata = params.pipelines_testdata_base_path + 'seqsubmit/test_data/study_metadata/study_metadata.json' + centre_name = "TEST_CENTER" + + test_upload = true +} diff --git a/docs/methods.md b/docs/methods.md index 1504fa1..38c79e8 100644 --- a/docs/methods.md +++ b/docs/methods.md @@ -12,7 +12,8 @@ The `ASSEMBLYSUBMIT` workflow: 2. Validates assembly FASTA files. 3. Reuses provided or calculates missing `coverage` values. 4. Builds the metadata table used for manifest generation. -5. Generates assembly manifests and submits the assemblies to ENA. +5. Registers a virtual ENA sample for co-assemblies whose runs come from different samples. +6. Generates assembly manifests and submits the assemblies to ENA. ### Assembly FASTA validation @@ -24,14 +25,20 @@ Only assemblies with successful validation are forwarded to coverage estimation, If the `coverage` column is already populated in the samplesheet, that value is used directly. -If `coverage` is missing, the workflow calculates coverage with `coverm contig`. +If `coverage` is missing, the workflow maps the reads of all runs listed for the assembly with `coverm contig`, which reports the mean depth and the length of each contig for each run. -`coverm contig` outputs per-contig depth. The workflow then reads this file and calculates the arithmetic mean across all contigs. +The same calculation is applied to every assembly. For an assembly from a single run ($N = 1$), it reduces to the length-weighted mean of the per-contig depths. -If the per-contig coverage values are $c_1, c_2, \ldots, c_n$, the workflow currently computes assembly coverage as an unweighted mean across contigs: +Because a co-assembly is built from the pooled reads of all runs, the depth of contig $c$ is the sum of its depths across the $N$ runs: $$ -\bar{c} = \frac{1}{n} \sum_{i=1}^{n} c_i +D_c = \sum_{i=1}^{N} D_{c,i} +$$ + +The assembly coverage is the mean of these per-contig depths, weighted by contig length $L_c$: + +$$ +D_{\text{assembly}} = \frac{\sum_{c} L_c \cdot D_c}{\sum_{c} L_c} $$ ## `GENOMESUBMIT` diff --git a/docs/usage.md b/docs/usage.md index a8ab3c5..61691bb 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -217,9 +217,11 @@ You will need to create a samplesheet with information about the metagenomic ass Example: ```csv title="samplesheet_assembly.csv" -id,fasta,fastq_1,fastq_2,coverage,run_accession,assembler,assembler_version,assembly_type -assembly_001,data/assembly_001.fasta.gz,data/assembly_001_R1.fastq.gz,data/assembly_001_R2.fastq.gz,,ERR011322,SPAdes,3.15.5,metagenome -assembly_002,data/assembly_002.fasta.gz,,,42.7,ERR011323,MEGAHIT,1.2.9,metatranscriptome +id,fasta,fastq_1,fastq_2,coverage,run_accession,assembler,assembler_version,assembly_type,tax_id +assembly_001,data/assembly_001.fasta.gz,data/assembly_001_R1.fastq.gz,data/assembly_001_R2.fastq.gz,,ERR011322,SPAdes,3.15.5,metagenome, +assembly_002,data/assembly_002.fasta.gz,,,42.7,ERR011323,MEGAHIT,1.2.9,metatranscriptome, +coassembly_001,data/coassembly_001.fasta.gz,data/run1_R1.fastq.gz;data/run2_R1.fastq.gz,data/run1_R2.fastq.gz;data/run2_R2.fastq.gz,,ERR011324;ERR011325,MEGAHIT,1.2.9,, +coassembly_002,data/coassembly_002.fasta.gz,,,35.2,ERR011326;ERR011327,MEGAHIT,1.2.9,metagenome,408172 ``` | Column | Required | Description | @@ -233,9 +235,21 @@ assembly_002,data/assembly_002.fasta.gz,,,42.7,ERR011323,MEGAHIT,1.2.9,metatrans | `assembler` | Yes | Name of the assembler software used to generate the assembly, e.g. MEGAHIT, SPAdes, metaSPAdes, Flye, etc. | | `assembler_version` | Yes | Version of the assembler software used to generate the assembly, e.g. 1.0.0. | | `assembly_type` | No | Sequencing library source: `metagenome` for metagenome assembly (would be used as _'primary metagenome'_ in manifest file following [ENA docs](https://ena-docs.readthedocs.io/en/latest/submit/assembly/metagenome/primary.html)) or `metatranscriptome` for metatranscriptome assembly (woule be used _'metatranscriptome'_ in assembly manifest following [ENA docs](https://ena-docs.readthedocs.io/en/latest/submit/assembly/metatranscriptome.html)). Default: `metagenome` | +| `tax_id` | No | NCBI taxon ID for the virtual sample registered for a co-assembly, e.g. `408172` (marine metagenome). Used only when the source samples have different taxon IDs; see [Co-assemblies](#co-assemblies). | Provide either read files (`fastq_1`, optionally `fastq_2`) or a `coverage` value for each row. If `coverage` is missing and reads are provided, the workflow calculates average coverage with CoverM. See the [Methods documentation](docs/methods.md) for more information on how coverage is calculated. +### Co-assemblies + +An assembly built from the pooled reads of several runs (a co-assembly) is described by a single samplesheet row listing all source runs in `run_accession`, separated by semicolons. If coverage has to be calculated, provide the read files of every run in `fastq_1` (and `fastq_2` for paired-end data) in the same order as the accessions; the pipeline fails if the number of read files does not match the number of run accessions. + +ENA requires each assembly to reference a single sample. The pipeline handles this as follows: + +- If all source runs belong to the same ENA sample, the assembly is submitted under that sample. +- If the runs belong to different samples, the pipeline registers a virtual ENA sample with a `composed of` attribute listing the source samples, and submits the assembly under it. The virtual sample uses the [ENA default sample checklist (`ERC000011`)](https://www.ebi.ac.uk/ena/browser/view/ERC000011), the minimal checklist, because its mandatory fields (`collection date` and `geographic location (country and/or sea)`) can be derived from any source samples, whichever checklists they were originally registered with. Its metadata is merged from the source samples: + - `collection date` and `geographic location (country and/or sea)` are kept if all source samples share the same value, otherwise they are set to `not provided`. A country that is not in the `ERC000011` list of allowed values is also set to `not provided`. + - `tax_id` must be the same for all source samples. If it differs, set the `tax_id` column for that samplesheet row to the taxon the co-assembly should be registered under (e.g. `408172` for marine metagenome when combining marine and marine sediment samples); otherwise the submission fails. The scientific name is fetched from the [ENA taxonomy](https://ena-docs.readthedocs.io/en/latest/faq/taxonomy.html) for the final `tax_id`. + ### `metagenomic_assemblies` parameters and example In addition to the [common parameters](#running-the-pipeline): @@ -274,8 +288,6 @@ nextflow run nf-core/seqsubmit \ `metagenomic_assemblies` mode does not yet support the following: -- **Co-assemblies**: Metagenomic co-assemblies are not supported. The pipeline currently expects each assembly to be derived from a single run. Tracked in [issue #61](https://github.com/nf-core/seqsubmit/issues/61), addressed in [PR #66](https://github.com/nf-core/seqsubmit/pull/66). - - **Submission of an assembly without submitting the raw reads**: Assemblies that have no associated reads accession (`run_accession` column in the samplesheet) are not supported, since the pipeline relies on that lineage to build the required metadata. Tracked in [issue #29](https://github.com/nf-core/seqsubmit/issues/29). - **Single-contig assemblies**: Assemblies that consist of only one contig cannot be uploaded through this pipeline. ENA classifies these as "chromosomal assemblies", which follow a different submission procedure with different metadata requirements. If you have a single-contig assembly to submit, please contact [ENA support](https://www.ebi.ac.uk/ena/browser/support) for guidance. diff --git a/modules/local/create_assembly_metadata_csv/main.nf b/modules/local/create_assembly_metadata_csv/main.nf index 5b25027..3e88485 100644 --- a/modules/local/create_assembly_metadata_csv/main.nf +++ b/modules/local/create_assembly_metadata_csv/main.nf @@ -19,13 +19,16 @@ process CREATE_ASSEMBLY_METADATA_CSV { script: def header = 'Runs,Coverage,Assembler,Version,Filepath,Sample,Library' + // Format run accessions: wrap in quotes if it's a list (co-assembly) or single value + def runs_joined = meta.run_accession instanceof List ? meta.run_accession.join(',') : meta.run_accession + def runs = "\"${runs_joined}\"" def row = [ - meta.run_accession, + runs, meta.coverage, meta.assembler, meta.assembler_version, fasta.name, - '', // Sample column left empty because co-assemblies are not supported + '', // Sample column is filled later for co-assemblies that require virtual sample registration meta.assembly_type ?: 'metagenome' // it would be applied by schema, but still needed as safe case for tests ].join(',') """ diff --git a/modules/local/create_assembly_metadata_csv/tests/main.nf.test.snap b/modules/local/create_assembly_metadata_csv/tests/main.nf.test.snap index af589fe..455ddaa 100644 --- a/modules/local/create_assembly_metadata_csv/tests/main.nf.test.snap +++ b/modules/local/create_assembly_metadata_csv/tests/main.nf.test.snap @@ -12,7 +12,7 @@ "assembler_version": "3.15.5", "assembly_type": "metagenome" }, - "test_sample_assembly_metadata.csv:md5,5627f0b0bc83720dd269e2ea8ac16eff" + "test_sample_assembly_metadata.csv:md5,000897a94b9f0873867fa463790583f2" ] ], "1": [ @@ -32,15 +32,15 @@ "assembler_version": "3.15.5", "assembly_type": "metagenome" }, - "test_sample_assembly_metadata.csv:md5,5627f0b0bc83720dd269e2ea8ac16eff" + "test_sample_assembly_metadata.csv:md5,000897a94b9f0873867fa463790583f2" ] ] } ], - "timestamp": "2026-10-15T11:45:01.123968", + "timestamp": "2026-10-29T15:18:32.564715", "meta": { "nf-test": "0.9.5", - "nextflow": "26.04.0" + "nextflow": "26.04.1" } }, "test_create_assembly_metadata_csv - stub": { diff --git a/modules/local/register_coassembly_sample/environment.yml b/modules/local/register_coassembly_sample/environment.yml new file mode 100644 index 0000000..66d3d46 --- /dev/null +++ b/modules/local/register_coassembly_sample/environment.yml @@ -0,0 +1,8 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - conda-forge::python=3.12.11 + - conda-forge::requests=2.32.5 + - conda-forge::retry=0.9.2 diff --git a/modules/local/register_coassembly_sample/main.nf b/modules/local/register_coassembly_sample/main.nf new file mode 100644 index 0000000..ba5f971 --- /dev/null +++ b/modules/local/register_coassembly_sample/main.nf @@ -0,0 +1,45 @@ +process REGISTER_COASSEMBLY_SAMPLE { + tag "${meta.id}" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'oras://community.wave.seqera.io/library/pip_requests_retry:d07e257656d938ab' : + 'community.wave.seqera.io/library/pip_requests_retry:eb16563fc2cb641f' }" + + input: + tuple val(meta), path(csv) + val(test) + val(is_private) + + output: + tuple val(meta), path("*_updated.csv"), emit: csv + tuple val("${task.process}"), val('python'), eval('python --version 2>&1 | sed "s/Python //g"'), topic: versions + tuple val("${task.process}"), val('requests'), eval('python -c "import requests; print(requests.__version__)"'), topic: versions + tuple val("${task.process}"), val('retry'), eval('python -c "import importlib.metadata; print(importlib.metadata.version(\'retry\'))"'), topic: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def test_arg = test ? "--test" : "" + def private_arg = is_private ? "--private" : "" + def taxid_arg = meta.tax_id ? "--default-taxid ${meta.tax_id}" : "" + """ + register_coassembly_sample.py \\ + --input ${csv} \\ + --output ${prefix}_updated.csv \\ + ${test_arg} \\ + ${private_arg} \\ + ${taxid_arg} \\ + ${args} + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + cp ${csv} ${prefix}_updated.csv + """ +} diff --git a/modules/local/register_coassembly_sample/meta.yml b/modules/local/register_coassembly_sample/meta.yml new file mode 100644 index 0000000..6f02c99 --- /dev/null +++ b/modules/local/register_coassembly_sample/meta.yml @@ -0,0 +1,39 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "register_coassembly_sample" +description: Register a virtual ENA sample for co-assemblies and write sample accession into assembly metadata CSV +keywords: + - ena + - sample + - metadata + - co-assembly +tools: + - python: + description: Python + homepage: https://www.python.org/ + - retry: + description: Retry decorator for resilient network calls + homepage: https://pypi.org/project/retry/ +input: + - - meta: + type: map + description: Groovy map containing sample metadata + - csv: + type: file + description: Assembly metadata CSV input for assembly_uploader tool + pattern: "*.csv" + - test: + type: string + description: Whether to submit to the ENA test server (`true`) instead of production. +output: + - csv: + - meta: + type: map + description: Same metadata map as input + - "*.csv": + type: file + description: Updated assembly metadata CSV with "Sample" column filled for co-assemblies + pattern: "*.csv" +authors: + - "@ochkalova" +maintainers: + - "@ochkalova" diff --git a/modules/local/register_coassembly_sample/tests/coassembly_metadata.csv b/modules/local/register_coassembly_sample/tests/coassembly_metadata.csv new file mode 100644 index 0000000..689d97b --- /dev/null +++ b/modules/local/register_coassembly_sample/tests/coassembly_metadata.csv @@ -0,0 +1,11 @@ +Runs,Coverage,Assembler,Version,Filepath +"SRR3183853,SRR3183854,SRR3183855,SRR3183856,SRR3183857,SRR3183858,SRR3183859,SRR3183860,SRR3183861,SRR3183952,SRR3183953,SRR3183954,SRR3183955,SRR3183956,SRR3183957,SRR3183958,SRR3183959,SRR3183960",10.0,MEGAHIT,1.2.9,contigs1.fasta.gz +"SRR3183850,SRR3183851,SRR3183852,SRR3183949,SRR3183950,SRR3183951",12.0,MEGAHIT,1.2.9,contigs2.fasta.gz +"SRR1631070,SRR1631071,SRR1631072,SRR1631073,SRR1631074,SRR1631075,SRR1631076,SRR1631077,SRR1631078,SRR1631079,SRR1631080,SRR1631081,SRR1631082,SRR1631083,SRR3187740,SRR3187741,SRR3187742,SRR3187743,SRR3187744,SRR3187745,SRR3187746,SRR3187747,SRR3187748,SRR3188145,SRR3188146,SRR3188147,SRR3188148,SRR3188149,SRR3188150,SRR3188151,SRR3188152,SRR3188153",15.0,MEGAHIT,1.2.9,contigs3.fasta.gz +"SRR1622826,SRR3187094,SRR3187149",8.0,MEGAHIT,1.2.9,contigs4.fasta.gz +"SRR1622910,SRR1622911",7.0,MEGAHIT,1.2.9,contigs5.fasta.gz +"SRR1624846,SRR1624847,SRR1624848,SRR1624849,SRR1624850,SRR1624851,SRR1624852,SRR1624853,SRR1624854",9.0,MEGAHIT,1.2.9,contigs6.fasta.gz +"SRR1631360,SRR1631361,SRR1631362,SRR1631363,SRR1631364,SRR1631365",6.0,MEGAHIT,1.2.9,contigs9.fasta.gz +"SRR3183597,SRR3183598,SRR3183599,SRR3183600",3.0,MEGAHIT,1.2.9,contigs12.fasta.gz +"SRR3183638,SRR3183639,SRR3183640,SRR3183641,SRR3183642,SRR3183643,SRR3183644,SRR3183645,SRR3183646,SRR3183647",2.0,MEGAHIT,1.2.9,contigs13.fasta.gz +"SRR3183844,SRR3183845,SRR3183846",1.0,MEGAHIT,1.2.9,contigs14.fasta.gz diff --git a/modules/local/register_coassembly_sample/tests/main.nf.test b/modules/local/register_coassembly_sample/tests/main.nf.test new file mode 100644 index 0000000..51df911 --- /dev/null +++ b/modules/local/register_coassembly_sample/tests/main.nf.test @@ -0,0 +1,37 @@ +nextflow_process { + name "Test Process REGISTER_COASSEMBLY_SAMPLE" + script "../main.nf" + config "./nextflow.config" + process "REGISTER_COASSEMBLY_SAMPLE" + + tag "modules" + tag "register_coassembly_sample" + + test("REGISTER_COASSEMBLY_SAMPLE completes with expected outputs") { + when { + process { + """ + input[0] = [ + [ id:'test' ], + file("${moduleDir}/tests/coassembly_metadata.csv", checkIfExists: true) + ] + input[1] = true + input[2] = false + """ + } + } + + then { + assert process.success + assertAll( + { assert process.out.csv.size() == 1 }, + { assert process.out.csv[0][1].toString().endsWith("_updated.csv") }, + { + def csvContent = path(process.out.csv[0][1]).text + assert csvContent.contains("Sample") + assert csvContent.contains("SRR3183853") + } + ) + } + } +} diff --git a/modules/local/register_coassembly_sample/tests/nextflow.config b/modules/local/register_coassembly_sample/tests/nextflow.config new file mode 100644 index 0000000..ea49989 --- /dev/null +++ b/modules/local/register_coassembly_sample/tests/nextflow.config @@ -0,0 +1,5 @@ +process { + withName: 'REGISTER_COASSEMBLY_SAMPLE' { + ext.args = '--debug' + } +} diff --git a/nextflow.config b/nextflow.config index dcfeb4c..4ec5621 100644 --- a/nextflow.config +++ b/nextflow.config @@ -255,6 +255,8 @@ profiles { test_assembly_no_coverage_single_reads { includeConfig 'conf/test_assembly_no_coverage_single_reads.config' } test_assembly_no_coverage_paired_reads { includeConfig 'conf/test_assembly_no_coverage_paired_reads.config' } test_assembly_one_contig { includeConfig 'conf/test_assembly_one_contig.config' } + test_assembly_coassembly_with_coverage { includeConfig 'conf/test_assembly_coassembly_with_coverage.config' } + test_assembly_coassembly_no_coverage { includeConfig 'conf/test_assembly_coassembly_no_coverage.config' } test_reads_paired { includeConfig 'conf/test_reads_paired.config' } test_reads_long_bam { includeConfig 'conf/test_reads_long_bam.config' } test_reads_long_bam_validate { includeConfig 'conf/test_reads_long_bam_validate.config' } diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/main.nf b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/main.nf index d4e930e..0dcba33 100644 --- a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/main.nf +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/main.nf @@ -191,6 +191,41 @@ def validateInputSamplesheet(input) { return [ metas[0], fastqs ] } +// +// Calculate co-assembly coverage from a `coverm contig --methods mean length` table. +// Rows are contigs; each read set contributes one ` Mean` and one ` Length` column. +// Per-contig depth is the sum of the `Mean` columns (reads from all runs were pooled for the +// co-assembly). Assembly depth is the length-weighted mean of the per-contig depths: +// sum_c(L_c * sum_i(D_c,i)) / sum_c(L_c) +// +def calculateCoassemblyCoverage(coverage_file) { + double weighted_depth = 0d + double total_length = 0d + + coverage_file.splitCsv(sep: '\t', header: true).each { row -> + def mean_columns = row.findAll { key, _value -> key.endsWith(' Mean') } + def length_column = row.find { key, _value -> key.endsWith(' Length') } + if (!mean_columns || !length_column) { + error("CoverM output ${coverage_file.name} must contain ' Mean' and ' Length' columns. Run 'coverm contig' with '--methods mean length'.") + } + + def depth = mean_columns.values().sum { value -> value as Double } + def length = length_column.value as Double + + weighted_depth += depth * length + total_length += length + } + + if (total_length == 0) { + error("CoverM output ${coverage_file.name} has no contigs, so assembly coverage cannot be calculated.") + } + def coverage = weighted_depth / total_length + if (coverage == 0) { + error("Calculated coverage for ${coverage_file.name} is 0: no reads mapped to the assembly. Check that the reads in 'fastq_1'/'fastq_2' belong to this assembly.") + } + return coverage +} + // // Generate methods description for MultiQC // diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_header_only.depth.tsv b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_header_only.depth.tsv new file mode 100644 index 0000000..4daac9d --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_header_only.depth.tsv @@ -0,0 +1 @@ +Contig contigs.fasta/run1_1.fastq.gz Mean contigs.fasta/run1_1.fastq.gz Length diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_no_length.depth.tsv b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_no_length.depth.tsv new file mode 100644 index 0000000..b140cc0 --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_no_length.depth.tsv @@ -0,0 +1,2 @@ +Contig contigs.fasta/run1_1.fastq.gz Mean +contig_1 10 diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_single_run.depth.tsv b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_single_run.depth.tsv new file mode 100644 index 0000000..e791392 --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_single_run.depth.tsv @@ -0,0 +1,3 @@ +Contig contigs.fasta/run1_1.fastq.gz Mean contigs.fasta/run1_1.fastq.gz Length +contig_1 10 1000 +contig_2 2 100 diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_three_runs.depth.tsv b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_three_runs.depth.tsv new file mode 100644 index 0000000..659598b --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_three_runs.depth.tsv @@ -0,0 +1,4 @@ +Contig contigs.fasta/run1_1.fastq.gz Mean contigs.fasta/run1_1.fastq.gz Length contigs.fasta/run2_1.fastq.gz Mean contigs.fasta/run2_1.fastq.gz Length contigs.fasta/run3_1.fastq.gz Mean contigs.fasta/run3_1.fastq.gz Length +contig_1 4 1000 3 1000 3 1000 +contig_2 1 3000 1 3000 0 3000 +contig_3 0 500 0 500 0 500 diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_zero_depth.depth.tsv b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_zero_depth.depth.tsv new file mode 100644 index 0000000..4ca64eb --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/fixtures/coverm_zero_depth.depth.tsv @@ -0,0 +1,3 @@ +Contig contigs.fasta/run1_1.fastq.gz Mean contigs.fasta/run1_1.fastq.gz Length +contig_1 0 1000 +contig_2 0 100 diff --git a/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/main.function.nf.test b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/main.function.nf.test new file mode 100644 index 0000000..937b342 --- /dev/null +++ b/subworkflows/local/utils_nfcore_seqsubmit_pipeline/tests/main.function.nf.test @@ -0,0 +1,102 @@ +nextflow_function { + + name "Test calculateCoassemblyCoverage" + script "../main.nf" + function "calculateCoassemblyCoverage" + tag "functions" + tag "calculateCoassemblyCoverage" + tag "mode_assembly" + + test("Single read set - length-weighted mean of contig depths") { + + when { + function { + """ + input[0] = file("${moduleDir}/tests/fixtures/coverm_single_run.depth.tsv") + """ + } + } + + then { + // (10 * 1000 + 2 * 100) / (1000 + 100) + assertAll( + { assert function.success }, + { assert Math.abs(function.result - 10200d / 1100d) < 1e-9 } + ) + } + } + + test("Several read sets - depths summed per contig, zero-depth contigs kept") { + + when { + function { + """ + input[0] = file("${moduleDir}/tests/fixtures/coverm_three_runs.depth.tsv") + """ + } + } + + then { + // contig_1: 4 + 3 + 3 = 10, contig_2: 1 + 1 + 0 = 2, contig_3: 0 + // (10 * 1000 + 2 * 3000 + 0 * 500) / (1000 + 3000 + 500) + assertAll( + { assert function.success }, + { assert Math.abs(function.result - 16000d / 4500d) < 1e-9 } + ) + } + } + + test("No contigs in CoverM output - fails with a clear error") { + + when { + function { + """ + input[0] = file("${moduleDir}/tests/fixtures/coverm_header_only.depth.tsv") + """ + } + } + + then { + assertAll( + { assert function.failed }, + { assert function.stdout.any { it.contains("has no contigs") } } + ) + } + } + + test("No reads mapped - coverage 0 fails with a clear error") { + + when { + function { + """ + input[0] = file("${moduleDir}/tests/fixtures/coverm_zero_depth.depth.tsv") + """ + } + } + + then { + assertAll( + { assert function.failed }, + { assert function.stdout.any { it.contains("no reads mapped") } } + ) + } + } + + test("CoverM run without '--methods mean length' - fails with a clear error") { + + when { + function { + """ + input[0] = file("${moduleDir}/tests/fixtures/coverm_no_length.depth.tsv") + """ + } + } + + then { + assertAll( + { assert function.failed }, + { assert function.stdout.any { it.contains("--methods mean length") } } + ) + } + } +} diff --git a/tests/assembly_coassembly_no_coverage.nf.test b/tests/assembly_coassembly_no_coverage.nf.test new file mode 100644 index 0000000..0846965 --- /dev/null +++ b/tests/assembly_coassembly_no_coverage.nf.test @@ -0,0 +1,45 @@ +nextflow_pipeline { + + name "Test assembly submission workflow - coassembly_no_coverage" + script "../main.nf" + tag "pipeline" + tag "mode_assembly" + tag "test_assembly_coassembly_no_coverage" + profile "test_assembly_coassembly_no_coverage" + + test("-profile test_assembly_coassembly_no_coverage") { + + when { + params { + outdir = "$outputDir" + } + } + + then { + // stable_name: All files + folders in ${params.outdir}/ with a stable name + def stable_name = getAllFilesFromDir(params.outdir, relative: true, includeDir: true, ignore: ['pipeline_info/*.{html,json,txt}']) + // stable_path: All files in ${params.outdir}/ with stable content + // assemblies_metadata.csv holds the co-assembly sample accession, newly registered on every --test run + def stable_path = getAllFilesFromDir(params.outdir, ignoreFile: 'tests/.nftignore', ignore: ['**/upload/assemblies_metadata.csv']) + // Early failure no need to test the rest of snapshots + assert workflow.success + def succeeded_tasks = workflow.trace.succeeded()*.name + assertAll( + // Coverage is computed with CoverM and the runs are registered as a virtual co-assembly sample + { assert succeeded_tasks.any { it.contains(':COVERM_CONTIG') } }, + { assert succeeded_tasks.any { it.contains(':REGISTER_COASSEMBLY_SAMPLE') } }, + { assert path("$outputDir/metagenomic_assemblies/upload/assemblies_metadata.csv").readLines()[1].contains(',1.263635324935214,') }, + { assert snapshot( + // Number of successful tasks + workflow.trace.succeeded().size(), + // pipeline versions.yml file for multiqc from which Nextflow version is removed because we test pipelines on multiple Nextflow versions + removeNextflowVersion("$outputDir/pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml"), + // All stable path name, with a relative path + stable_name, + // All files with stable contents + stable_path + ).match() } + ) + } + } +} diff --git a/tests/assembly_coassembly_no_coverage.nf.test.snap b/tests/assembly_coassembly_no_coverage.nf.test.snap new file mode 100644 index 0000000..063c1be --- /dev/null +++ b/tests/assembly_coassembly_no_coverage.nf.test.snap @@ -0,0 +1,88 @@ +{ + "-profile test_assembly_coassembly_no_coverage": { + "content": [ + 10, + { + "CONCAT_ACCESSIONS": { + "coreutils": 9.4, + "find": "4.6.0", + "pigz": 2.8 + }, + "CONCAT_METADATA": { + "coreutils": 9.4, + "find": "4.6.0", + "pigz": 2.8 + }, + "COVERM_CONTIG": { + "coverm": "0.7.0" + }, + "CREATE_ASSEMBLY_METADATA_CSV": { + "bash": "5.2.37(1)-release" + }, + "FALINT": { + "falint": "1.2.0" + }, + "GENERATE_ASSEMBLY_MANIFEST": { + "assembly_uploader": "1.3.6" + }, + "REGISTERSTUDY": { + "mgnify-pipelines-toolkit": "1.5.5" + }, + "REGISTER_COASSEMBLY_SAMPLE": { + "python": "3.14.4", + "requests": "2.34.1", + "retry": "0.9.2" + }, + "SUBMIT": { + "ena-webin-cli": "9.0.3", + "mgnify-pipelines-toolkit": "1.5.5", + "python": "3.12.16" + }, + "Workflow": { + "nf-core/seqsubmit": "v1.0.0" + } + }, + [ + "metagenomic_assemblies", + "metagenomic_assemblies/no_coverage_paired_reads", + "metagenomic_assemblies/no_coverage_paired_reads/coverage", + "metagenomic_assemblies/no_coverage_paired_reads/coverage/no_coverage_paired_reads.depth.tsv", + "metagenomic_assemblies/no_coverage_paired_reads/upload", + "metagenomic_assemblies/no_coverage_paired_reads/upload/webin_cli", + "metagenomic_assemblies/no_coverage_paired_reads/upload/webin_cli/no_coverage_paired_reads_accessions.tsv", + "metagenomic_assemblies/upload", + "metagenomic_assemblies/upload/assemblies_accessions.tsv", + "metagenomic_assemblies/upload/assemblies_metadata.csv", + "multiqc", + "multiqc/multiqc_data", + "multiqc/multiqc_data/llms-full.txt", + "multiqc/multiqc_data/multiqc.log", + "multiqc/multiqc_data/multiqc.parquet", + "multiqc/multiqc_data/multiqc_Assemblies assigned accessions.txt", + "multiqc/multiqc_data/multiqc_citations.txt", + "multiqc/multiqc_data/multiqc_data.json", + "multiqc/multiqc_data/multiqc_software_versions.txt", + "multiqc/multiqc_data/multiqc_sources.txt", + "multiqc/multiqc_plots", + "multiqc/multiqc_plots/pdf", + "multiqc/multiqc_plots/pdf/Assemblies assigned accessions.pdf", + "multiqc/multiqc_plots/png", + "multiqc/multiqc_plots/png/Assemblies assigned accessions.png", + "multiqc/multiqc_plots/svg", + "multiqc/multiqc_plots/svg/Assemblies assigned accessions.svg", + "multiqc/multiqc_report.html", + "pipeline_info", + "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" + ], + [ + "no_coverage_paired_reads.depth.tsv:md5,6e61f10c995d45b15f33d572d773e80e", + "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" + ] + ], + "timestamp": "2026-10-19T16:26:50.911456", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.1" + } + } +} \ No newline at end of file diff --git a/tests/assembly_coassembly_with_coverage.nf.test b/tests/assembly_coassembly_with_coverage.nf.test new file mode 100644 index 0000000..5b606e4 --- /dev/null +++ b/tests/assembly_coassembly_with_coverage.nf.test @@ -0,0 +1,43 @@ +nextflow_pipeline { + + name "Test assembly submission workflow - coassembly_with_coverage" + script "../main.nf" + tag "pipeline" + tag "mode_assembly" + tag "test_assembly_coassembly_with_coverage" + profile "test_assembly_coassembly_with_coverage" + + test("-profile test_assembly_coassembly_with_coverage") { + + when { + params { + outdir = "$outputDir" + } + } + + then { + // stable_name: All files + folders in ${params.outdir}/ with a stable name + def stable_name = getAllFilesFromDir(params.outdir, relative: true, includeDir: true, ignore: ['pipeline_info/*.{html,json,txt}']) + // stable_path: All files in ${params.outdir}/ with stable content + // assemblies_metadata.csv holds the co-assembly sample accession, newly registered on every --test run + def stable_path = getAllFilesFromDir(params.outdir, ignoreFile: 'tests/.nftignore', ignore: ['**/upload/assemblies_metadata.csv']) + // Early failure no need to test the rest of snapshots + assert workflow.success + def all_tasks = workflow.trace.tasks()*.name + assertAll( + // Coverage is provided in the samplesheet, so CoverM must not run + { assert !all_tasks.any { it.contains(':COVERM_CONTIG') } }, + { assert snapshot( + // Number of successful tasks + workflow.trace.succeeded().size(), + // pipeline versions.yml file for multiqc from which Nextflow version is removed because we test pipelines on multiple Nextflow versions + removeNextflowVersion("$outputDir/pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml"), + // All stable path name, with a relative path + stable_name, + // All files with stable contents + stable_path + ).match() } + ) + } + } +} diff --git a/tests/assembly_coassembly_with_coverage.nf.test.snap b/tests/assembly_coassembly_with_coverage.nf.test.snap new file mode 100644 index 0000000..19c51e2 --- /dev/null +++ b/tests/assembly_coassembly_with_coverage.nf.test.snap @@ -0,0 +1,82 @@ +{ + "-profile test_assembly_coassembly_with_coverage": { + "content": [ + 9, + { + "CONCAT_ACCESSIONS": { + "coreutils": 9.4, + "find": "4.6.0", + "pigz": 2.8 + }, + "CONCAT_METADATA": { + "coreutils": 9.4, + "find": "4.6.0", + "pigz": 2.8 + }, + "CREATE_ASSEMBLY_METADATA_CSV": { + "bash": "5.2.37(1)-release" + }, + "FALINT": { + "falint": "1.2.0" + }, + "GENERATE_ASSEMBLY_MANIFEST": { + "assembly_uploader": "1.3.6" + }, + "REGISTERSTUDY": { + "mgnify-pipelines-toolkit": "1.5.5" + }, + "REGISTER_COASSEMBLY_SAMPLE": { + "python": "3.14.4", + "requests": "2.34.1", + "retry": "0.9.2" + }, + "SUBMIT": { + "ena-webin-cli": "9.0.3", + "mgnify-pipelines-toolkit": "1.5.5", + "python": "3.12.16" + }, + "Workflow": { + "nf-core/seqsubmit": "v1.0.0" + } + }, + [ + "metagenomic_assemblies", + "metagenomic_assemblies/coassembly_sample", + "metagenomic_assemblies/coassembly_sample/upload", + "metagenomic_assemblies/coassembly_sample/upload/webin_cli", + "metagenomic_assemblies/coassembly_sample/upload/webin_cli/coassembly_sample_accessions.tsv", + "metagenomic_assemblies/upload", + "metagenomic_assemblies/upload/assemblies_accessions.tsv", + "metagenomic_assemblies/upload/assemblies_metadata.csv", + "multiqc", + "multiqc/multiqc_data", + "multiqc/multiqc_data/llms-full.txt", + "multiqc/multiqc_data/multiqc.log", + "multiqc/multiqc_data/multiqc.parquet", + "multiqc/multiqc_data/multiqc_Assemblies assigned accessions.txt", + "multiqc/multiqc_data/multiqc_citations.txt", + "multiqc/multiqc_data/multiqc_data.json", + "multiqc/multiqc_data/multiqc_software_versions.txt", + "multiqc/multiqc_data/multiqc_sources.txt", + "multiqc/multiqc_plots", + "multiqc/multiqc_plots/pdf", + "multiqc/multiqc_plots/pdf/Assemblies assigned accessions.pdf", + "multiqc/multiqc_plots/png", + "multiqc/multiqc_plots/png/Assemblies assigned accessions.png", + "multiqc/multiqc_plots/svg", + "multiqc/multiqc_plots/svg/Assemblies assigned accessions.svg", + "multiqc/multiqc_report.html", + "pipeline_info", + "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" + ], + [ + "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" + ] + ], + "timestamp": "2026-10-19T16:27:45.228133", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.1" + } + } +} \ No newline at end of file diff --git a/tests/assembly_complete_metadata.nf.test.snap b/tests/assembly_complete_metadata.nf.test.snap index 45ea705..daf4de4 100644 --- a/tests/assembly_complete_metadata.nf.test.snap +++ b/tests/assembly_complete_metadata.nf.test.snap @@ -69,7 +69,7 @@ "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" ], [ - "assemblies_metadata.csv:md5,4a171d2d13d999e391fd6b460f5927ec", + "assemblies_metadata.csv:md5,98b1419c74815101cfacd3daf92a435a", "multiqc_assembly_metadata.txt:md5,3187b2ddb3922a060b58f330c7609d0d", "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" ] diff --git a/tests/assembly_complete_metadata_existing_study.nf.test.snap b/tests/assembly_complete_metadata_existing_study.nf.test.snap index 9e251d5..848c6db 100644 --- a/tests/assembly_complete_metadata_existing_study.nf.test.snap +++ b/tests/assembly_complete_metadata_existing_study.nf.test.snap @@ -66,7 +66,7 @@ "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" ], [ - "assemblies_metadata.csv:md5,4a171d2d13d999e391fd6b460f5927ec", + "assemblies_metadata.csv:md5,98b1419c74815101cfacd3daf92a435a", "multiqc_assembly_metadata.txt:md5,3187b2ddb3922a060b58f330c7609d0d", "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" ] diff --git a/tests/assembly_no_coverage_paired_reads.nf.test.snap b/tests/assembly_no_coverage_paired_reads.nf.test.snap index 21cec5d..9edf9d4 100644 --- a/tests/assembly_no_coverage_paired_reads.nf.test.snap +++ b/tests/assembly_no_coverage_paired_reads.nf.test.snap @@ -74,9 +74,9 @@ "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" ], [ - "no_coverage_paired_reads.depth.tsv:md5,bb5f99e74d21df3c73e0ae0f388bcbcb", - "assemblies_metadata.csv:md5,199890227285973c56ea126e7f556017", - "multiqc_assembly_metadata.txt:md5,07b7d9fc2d6b6eaa9e779af80eb8fdde", + "no_coverage_paired_reads.depth.tsv:md5,ed807966913ebeec8b3face8e30317ed", + "assemblies_metadata.csv:md5,5a2a691d2e2ecf059a74cfa56bef5f19", + "multiqc_assembly_metadata.txt:md5,d830fad9de3ef5531c2170db07341e3a", "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" ] ], diff --git a/tests/assembly_no_coverage_single_reads.nf.test.snap b/tests/assembly_no_coverage_single_reads.nf.test.snap index 7dbf494..7fad7d9 100644 --- a/tests/assembly_no_coverage_single_reads.nf.test.snap +++ b/tests/assembly_no_coverage_single_reads.nf.test.snap @@ -74,9 +74,9 @@ "pipeline_info/nf_core_seqsubmit_software_mqc_versions.yml" ], [ - "no_coverage_single_reads.depth.tsv:md5,bd88282b17e851377b1dd223839be150", - "assemblies_metadata.csv:md5,f7506e0c624ad2909827966bf2526dfd", - "multiqc_assembly_metadata.txt:md5,e7fd323a6691296e89dc20174fcd8083", + "no_coverage_single_reads.depth.tsv:md5,7d71e53c7424f930cee76eee9fa91b04", + "assemblies_metadata.csv:md5,cfd0b7cab10c2f93d9dcff3435cff8ec", + "multiqc_assembly_metadata.txt:md5,4d08d8ca7b7d37c5db8fedeabfca82fe", "multiqc_citations.txt:md5,4c806e63a283ec1b7e78cdae3a923d4f" ] ], diff --git a/workflows/assemblysubmit.nf b/workflows/assemblysubmit.nf index bb3fe00..31a0503 100644 --- a/workflows/assemblysubmit.nf +++ b/workflows/assemblysubmit.nf @@ -6,6 +6,7 @@ include { COVERM_CONTIG } from '../modules/nf-core/coverm/contig/main' include { CREATE_ASSEMBLY_METADATA_CSV } from '../modules/local/create_assembly_metadata_csv/main' +include { REGISTER_COASSEMBLY_SAMPLE } from '../modules/local/register_coassembly_sample/main' include { GENERATE_ASSEMBLY_MANIFEST } from '../modules/local/generate_assembly_manifest/main' include { REGISTERSTUDY } from '../modules/local/registerstudy/main' include { ENA_WEBIN_CLI_WRAPPER as SUBMIT } from '../modules/local/ena_webin_cli_wrapper' @@ -20,6 +21,7 @@ include { softwareVersionsToYAML } from '../subworkflows/nf-core/ include { FASTA_VALIDATION } from '../subworkflows/local/fasta_validation/main' include { methodsDescriptionText } from '../subworkflows/local/utils_nfcore_seqsubmit_pipeline' +include { calculateCoassemblyCoverage } from '../subworkflows/local/utils_nfcore_seqsubmit_pipeline' /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -50,20 +52,49 @@ workflow ASSEMBLYSUBMIT { // --------- Create assembly channel with proper metadata structure assembly_fasta = ch_samplesheet - .map { meta, fasta, reads_1, reads_2 -> - def new_meta = meta + [single_end: !reads_2] + .map { meta, fasta, _reads_1, reads_2 -> + // support semicolon-separated values for co-assemblies (e.g. ERR000001;ERR000002) + def run_accessions = meta.run_accession.split(';')?.toList()?.unique() ?: [] + def new_meta = meta + [ + single_end: !reads_2, + run_accession: run_accessions.size() == 1 ? run_accessions[0] : run_accessions + ] [new_meta, fasta] } // --------- Create reads channel with proper metadata structure reads_fastq = ch_samplesheet .filter { row -> row[2] && row[2] != "" } // Check if fastq_1 exists and is not empty - .map { meta, fasta, reads_1, reads_2 -> - def new_meta = meta + [single_end: !reads_2] - def reads = new_meta.single_end - ? [reads_1] - : [reads_1, reads_2] - [new_meta, reads] + .map { meta, _fasta, reads_1, reads_2 -> + // support semicolon-separated values for co-assemblies (e.g. ERR000001;ERR000002) + def run_accessions = meta.run_accession.split(';')?.toList()?.unique() ?: [] + def new_meta = meta + [ + single_end: !reads_2, + run_accession: run_accessions.size() == 1 ? run_accessions[0] : run_accessions + ] + + def fastq1_list = reads_1.toString().split(';').toList() + def fastq2_list = reads_2 ? reads_2.toString().split(';').toList() : [] + + // Validation: Check that number of read files matches number of accessions + if (run_accessions.size() > 1) { + if (fastq1_list.size() != run_accessions.size()) { + error "Sample ${new_meta.id}: Number of forward read files (${fastq1_list.size()}) does not match number of run accessions (${run_accessions.size()})" + } + if (fastq2_list && fastq2_list.size() != run_accessions.size()) { + error "Sample ${new_meta.id}: Number of reverse read files (${fastq2_list.size()}) does not match number of run accessions (${run_accessions.size()})" + } + } + + // Convert paths to file objects + def fastq1_paths = fastq1_list.collect { path -> file(path.trim(), checkIfExists: true) } + def fastq2_paths = fastq2_list ? fastq2_list.collect { path -> file(path.trim(), checkIfExists: true) } : [] + + if (fastq2_paths) { + [new_meta, fastq1_paths, fastq2_paths] + } else { + [new_meta, fastq1_paths] + } } // --------- Check fasta files are properly formatted and filter out files with less than 2 contigs @@ -76,9 +107,33 @@ workflow ASSEMBLYSUBMIT { coverm_input = FASTA_VALIDATION.out.valid_fastas .filter { meta, _fasta -> !(meta.coverage) } .join(reads_fastq) - .multiMap { meta, fasta, fastq -> + .map { tuple -> + def meta = tuple[0] + def fasta = tuple[1] + def reads_data = tuple[2..-1] + + // Transform reads into flat list for CoverM + def all_reads = [] + if (meta.single_end) { + // Single-end: just flatten the R1 list + def fastq1_list = reads_data[0] + all_reads = fastq1_list + } else { + // Paired-end: interleave R1 and R2 files + // Input format: [meta, [R1_1, R1_2, ...], [R2_1, R2_2, ...]] + def fastq1_list = reads_data[0] + def fastq2_list = reads_data[1] + + // Create list as: R1_1, R2_1, R1_2, R2_2, ... + // Use collectMany to flatten the pairs + all_reads = [fastq1_list, fastq2_list].transpose().collectMany { pair -> pair } + } + + [meta, fasta, all_reads] + } + .multiMap { meta, fasta, reads -> assembly: [ meta, fasta ] - reads: [ meta, fastq ] + reads: [ meta, reads ] } COVERM_CONTIG ( @@ -89,16 +144,10 @@ workflow ASSEMBLYSUBMIT { false // enable_bam_output ) - // Calculate average coverage using splitCsv operator + // Length-weighted co-assembly coverage (see docs/methods.md) average_coverage_ch = COVERM_CONTIG.out.coverage - .splitCsv(sep: '\t', skip: 1) - .map { meta, row -> - [meta, row[1] as Double] - } - .groupTuple() - .map { meta, coverages -> - def average = coverages.sum() / coverages.size() - [meta, average] + .map { meta, coverage_tsv -> + [ meta, calculateCoassemblyCoverage(coverage_tsv) ] } // Update metadata with calculated coverage @@ -122,9 +171,26 @@ workflow ASSEMBLYSUBMIT { assemblies_with_coverage ) + // For co-assemblies (multiple run accessions), register a virtual ENA sample and fill Sample column. + // Rows with a single run accession bypass this step unchanged. + assembly_metadata_by_type = CREATE_ASSEMBLY_METADATA_CSV.out.csv + .branch { meta, _csv -> + coassembly: meta.run_accession instanceof List && meta.run_accession.size() > 1 + individual_assembly: true + } + + REGISTER_COASSEMBLY_SAMPLE( + assembly_metadata_by_type.coassembly, + test_upload, + is_private + ) + + assembly_metadata_with_sample = assembly_metadata_by_type.individual_assembly + .mix(REGISTER_COASSEMBLY_SAMPLE.out.csv) + // --------- Concatenate assembly metadata CSVs into single file to publish CONCAT_METADATA ( - CREATE_ASSEMBLY_METADATA_CSV.out.csv.map { _meta, file -> file }.collect().map { files -> [ [id: "assemblies_metadata"], files ] }, + assembly_metadata_with_sample.map { _meta, file -> file }.collect().map { files -> [ [id: "assemblies_metadata"], files ] }, 'true' // skip_header - we want to keep the header from the first file and skip it for the rest ) @@ -151,7 +217,7 @@ workflow ASSEMBLYSUBMIT { // --------- Generate assembly manifest files and submit them to ENA GENERATE_ASSEMBLY_MANIFEST( - assemblies_with_coverage.join(CREATE_ASSEMBLY_METADATA_CSV.out.csv), + assemblies_with_coverage.join(assembly_metadata_with_sample), study_accession_ch.first(), upload_tpa, test_upload,