diff --git a/docs/src/en/ref.md b/docs/src/en/ref.md index 0d799ab64..1f7ce6975 100644 --- a/docs/src/en/ref.md +++ b/docs/src/en/ref.md @@ -10,7 +10,8 @@ Return format: dictionary/JSON. Species for which the FTPs will be fetched in the format genus_species, e.g. homo_sapiens. Supports all available vertebrate and invertebrate (plants, fungi, protists, and invertebrate metazoa) genomes from Ensembl, except bacteria. Note: Not required when using flags `--list_species` or `--list_iv_species`. -Supported shortcuts: 'human', 'mouse', 'human_grch37' (accesses the GRCh37 genome assembly) +Supported shortcuts: 'human', 'mouse', 'human_grch37' (accesses the GRCh37 genome assembly) +When using the `--assembly_report` flag, this is instead an NCBI assembly accession, e.g. GCF_000001405.40. The version suffix is optional; if omitted (e.g. GCF_000001405), the latest available version is used. **Optional arguments** `-w` `--which` @@ -43,6 +44,22 @@ Lists all available invertebrate species. (Python: combine with `species=None`.) `-ftp` `--ftp` Returns only the requested FTP links. +`-ar` `--assembly_report` +Returns the [NCBI assembly report](https://www.ncbi.nlm.nih.gov/datasets/docs/v2/troubleshooting/faq/#what-is-an-assembly-report) for the NCBI assembly accession passed as the positional argument (instead of Ensembl FTP links). The report maps sequence/chromosome names across the Ensembl/short, GenBank, RefSeq, and UCSC naming conventions, which is useful for translating chromosome names between databases. +Python: returns a `pandas` DataFrame (use `assembly_report=True`). + +`-tx` `--taxon` +Only used with `--assembly_report`: interpret the positional argument as an organism/taxon name (e.g. `"homo sapiens"`) and resolve it to that taxon's NCBI reference assembly before fetching the report. +Python: use `taxon=True`. + +`-la` `--list_assemblies` +Only used with `--assembly_report`: interpret the positional argument as an organism/taxon name and list all of that taxon's NCBI assemblies (accession, assembly_name, refseq_category, assembly_level, organism), reference/representative first, instead of a report. Use it to find a specific accession to fetch. +Python: use `list_assemblies=True`. + +`-csv` `--csv` +Command-line only. Only used with `--assembly_report`: returns the report in csv format instead of json. +Python: Use `json=True` to return a list of dictionaries instead of a DataFrame. + `-d` `--download` Command-line only. Downloads the requested FTPs to the directory specified by `out_dir` (requires [curl](https://curl.se/docs/) to be installed). @@ -98,6 +115,47 @@ gget.ref(species=None, list_species=True, release=103)

+**Get the NCBI assembly report to translate chromosome names between conventions:** +```bash +gget ref GCF_000001405.40 --assembly_report +``` +```python +# Python +gget.ref("GCF_000001405.40", assembly_report=True) +``` +→ Returns the NCBI assembly report for the human GRCh38.p14 assembly, mapping each sequence across the Ensembl/short (`Sequence-Name`), GenBank (`GenBank-Accn`), RefSeq (`RefSeq-Accn`), and UCSC (`UCSC-style-name`) naming conventions: + +| Sequence-Name | Sequence-Role | Assigned-Molecule | ... | GenBank-Accn | RefSeq-Accn | ... | UCSC-style-name | +| --- | --- | --- | --- | --- | --- | --- | --- | +| 1 | assembled-molecule | 1 | ... | CM000663.2 | NC_000001.11 | ... | chr1 | +| 2 | assembled-molecule | 2 | ... | CM000664.2 | NC_000002.12 | ... | chr2 | + +

+ +**Get the assembly report by organism name instead of accession (resolves to the reference assembly):** +```bash +gget ref "homo sapiens" --assembly_report --taxon +``` +```python +# Python +gget.ref("homo sapiens", assembly_report=True, taxon=True) +``` +→ Resolves *homo sapiens* to its NCBI reference assembly (`GCF_000001405.40`) and returns the same report as above. (If you already know the specific assembly you want, pass its accession directly without `--taxon`.) + +

+ +**List all assemblies for an organism (to pick a specific non-reference assembly):** +```bash +gget ref "homo sapiens" --assembly_report --list_assemblies +``` +```python +# Python +gget.ref("homo sapiens", assembly_report=True, list_assemblies=True) +``` +→ Returns a table of all NCBI assemblies for *homo sapiens* (reference/representative first), e.g. `GCF_000001405.40` (GRCh38.p14), `GCF_009914755.1` (T2T-CHM13v2.0), .... Pick the accession you want and pass it back to `--assembly_report` to fetch its report. + +

+ **Use `gget ref` in combination with [kallisto | bustools](https://www.kallistobus.tools/kb_usage/kb_ref/) to build a reference index:** ```bash kb ref \ diff --git a/docs/src/en/updates.md b/docs/src/en/updates.md index 89be22281..49148812b 100644 --- a/docs/src/en/updates.md +++ b/docs/src/en/updates.md @@ -5,6 +5,12 @@ #### *gget* officially became part of [*scverse*](https://scverse.org/) on June 9, 2026. 🥳🥳🥳 **Version ≥ 0.30.9** (XXX XX, 2026): +- [`gget ref`](ref.md): Added support for fetching [NCBI assembly reports](https://www.ncbi.nlm.nih.gov/datasets/docs/v2/troubleshooting/faq/#what-is-an-assembly-report) (fixes [issue 179](https://github.com/scverse/gget/issues/179)). + - New `assembly_report=True` option (command line: `--assembly_report`/`-ar`) interprets the positional argument as an NCBI assembly accession (e.g. `GCF_000001405.40`) and returns the assembly report, which maps sequence/chromosome names across the Ensembl/short, GenBank, RefSeq, and UCSC naming conventions (e.g. `chr1` ↔ `1` ↔ `CM000663.2` ↔ `NC_000001.11`). + - The accession version suffix is optional; if omitted (e.g. `GCF_000001405`), the latest available version is used. + - New `taxon=True` option (command line: `--taxon`/`-tx`) interprets the input as an organism/taxon name (e.g. `"homo sapiens"`) and resolves it to that taxon's NCBI reference assembly before fetching the report. + - New `list_assemblies=True` option (command line: `--list_assemblies`/`-la`) lists all NCBI assemblies for an organism/taxon name (reference/representative first) so a specific non-reference assembly can be selected. + - Returns a `pandas` DataFrame in Python (pass `json=True` for a list of dictionaries). On the command line the report prints as JSON by default; use `--csv` for CSV. The new `gget.assembly_report()` function is also exposed for direct use. **Version ≥ 0.30.8** (Jun 28, 2026): - [`gget g2p`](g2p.md): Either `gene` or `--uniprot_id` is now sufficient — whichever is missing is resolved via UniProt and cached. Gene→UniProt picks the canonical reviewed human Swiss-Prot entry; the resolution and its limitations are logged. The canonical pair is **always** prepended to the result as `gene_name` / `uniprot_id` columns (and stored on `df.attrs`), so the output schema is invariant regardless of input mode. Existing call sites continue to work. diff --git a/gget/__init__.py b/gget/__init__.py index f56cbcddc..b5e70883e 100644 --- a/gget/__init__.py +++ b/gget/__init__.py @@ -22,7 +22,7 @@ from .gget_mutate import mutate from .gget_opentargets import opentargets from .gget_pdb import pdb -from .gget_ref import ref +from .gget_ref import assembly_report, ref from .gget_search import search from .gget_seq import seq from .gget_setup import setup diff --git a/gget/constants.py b/gget/constants.py index 38f90119f..f3ba21fb2 100644 --- a/gget/constants.py +++ b/gget/constants.py @@ -17,6 +17,9 @@ # NCBI URL for gget info NCBI_URL = "https://www.ncbi.nlm.nih.gov" +# NCBI genomes FTP root for gget ref assembly reports +NCBI_FTP_GENOMES_URL = "https://ftp.ncbi.nlm.nih.gov/genomes/all/" + # NCBI VIRUS REST API URL for gget virus - Version 2 API endpoint NCBI_API_BASE = "https://api.ncbi.nlm.nih.gov/datasets/v2" diff --git a/gget/gget_ref.py b/gget/gget_ref.py index b16391c41..14f11ba82 100644 --- a/gget/gget_ref.py +++ b/gget/gget_ref.py @@ -1,8 +1,11 @@ from __future__ import annotations -import json -from typing import Any +import json as json_package +import re +import subprocess +from typing import Any, cast +import pandas as pd import requests from bs4 import BeautifulSoup @@ -21,9 +24,271 @@ ENSEMBL_FTP_URL, ENSEMBL_FTP_URL_GRCH37, ENSEMBL_FTP_URL_NV, + NCBI_FTP_GENOMES_URL, ) +def _resolve_taxon_to_accession(taxon_name: str, verbose: bool = True) -> str: + """Resolve an organism / taxon name to its NCBI reference assembly accession. + + Uses the bundled NCBI `datasets` CLI (`datasets summary genome taxon --reference`) + and returns the reference assembly's accession, e.g. "homo sapiens" -> "GCF_000001405.40". + """ + # Imported lazily so that plain `gget ref` / accession-mode assembly_report do not pull + # in the heavier gget_virus module (only taxon-name resolution needs the datasets CLI). + from .gget_virus import _get_datasets_path + + datasets_path = _get_datasets_path() + try: + result = subprocess.run( + [datasets_path, "summary", "genome", "taxon", taxon_name, "--reference", "--as-json-lines"], + capture_output=True, + text=True, + timeout=120, + ) + except (subprocess.SubprocessError, OSError) as e: + raise RuntimeError(f"Failed to run the NCBI datasets CLI to resolve taxon '{taxon_name}': {e}\n") from e + + if result.returncode != 0: + raise RuntimeError( + f"The NCBI datasets CLI failed while resolving taxon '{taxon_name}' (exit {result.returncode}). " + "This is often a transient network issue; please try again.\n" + ) + + for line in result.stdout.splitlines(): + line = line.strip() + if not line.startswith("{"): + continue # skip non-JSON notices (e.g. update banners) + try: + record = json_package.loads(line) + except ValueError: + continue + accession = record.get("accession") + if accession: + if verbose: + logger.info(f"Resolved taxon '{taxon_name}' to reference assembly {accession}.") + return accession + + raise ValueError( + f"Could not find a reference assembly for taxon '{taxon_name}'. Please check the name, " + "or pass an NCBI assembly accession (e.g. 'GCF_000001405.40') directly.\n" + ) + + +def _list_taxon_assemblies( + taxon_name: str, json: bool = False, save: bool = False, verbose: bool = True +) -> pd.DataFrame | list[dict[str, Any]]: + """List all NCBI genome assemblies available for an organism/taxon name. + + Returns a DataFrame (or list of dicts) with columns accession, assembly_name, + refseq_category, assembly_level and organism, with the reference / representative + assemblies listed first. Use it to discover a specific accession to then pass to + assembly_report(). + """ + from .gget_virus import _get_datasets_path + + datasets_path = _get_datasets_path() + try: + result = subprocess.run( + [datasets_path, "summary", "genome", "taxon", taxon_name, "--as-json-lines"], + capture_output=True, + text=True, + timeout=120, + ) + except (subprocess.SubprocessError, OSError) as e: + raise RuntimeError(f"Failed to run the NCBI datasets CLI to list assemblies for '{taxon_name}': {e}\n") from e + + if result.returncode != 0: + raise RuntimeError( + f"The NCBI datasets CLI failed while listing assemblies for '{taxon_name}' (exit {result.returncode}). " + "This is often a transient network issue; please try again.\n" + ) + + rows: list[dict[str, Any]] = [] + for line in result.stdout.splitlines(): + line = line.strip() + if not line.startswith("{"): + continue # skip non-JSON notices (e.g. update banners) + try: + record = json_package.loads(line) + except ValueError: + continue + info = record.get("assembly_info", {}) + rows.append( + { + "accession": record.get("accession", ""), + "assembly_name": info.get("assembly_name", ""), + "refseq_category": info.get("refseq_category") or "na", + "assembly_level": info.get("assembly_level", ""), + "organism": record.get("organism", {}).get("organism_name", ""), + } + ) + + if not rows: + raise ValueError(f"No assemblies found for taxon '{taxon_name}'. Please check the name.\n") + + # Surface the reference / representative assemblies first, then by accession. + category_rank = {"reference genome": 0, "representative genome": 1} + rows.sort(key=lambda r: (category_rank.get(r["refseq_category"], 2), r["accession"])) + df = pd.DataFrame(rows) + + if verbose: + logger.info(f"Found {len(df)} assemblies for taxon '{taxon_name}'.") + + file_stem = taxon_name.replace(" ", "_") + if json: + result_list = cast("list[dict[str, Any]]", df.to_dict(orient="records")) + if save: + with open(f"{file_stem}_assemblies.json", "w", encoding="utf-8") as f: + json_package.dump(result_list, f, ensure_ascii=False, indent=4) + return result_list + + if save: + df.to_csv(f"{file_stem}_assemblies.csv", index=False) + + return df + + +def assembly_report( + accession: str, + json: bool = False, + save: bool = False, + taxon: bool = False, + list_assemblies: bool = False, + verbose: bool = True, +) -> pd.DataFrame | list[dict[str, Any]]: + """Fetch the NCBI assembly report for a genome assembly accession. + + The assembly report maps sequence (e.g. chromosome) names across the + different naming conventions (Ensembl/short names, GenBank accessions, + RefSeq accessions, UCSC-style names), which is useful for translating + chromosome names between databases. + + Args: + - accession NCBI assembly accession, e.g. "GCF_000001405.40" (RefSeq) or + "GCA_000001405.29" (GenBank). The version suffix is optional; if omitted + (e.g. "GCF_000001405"), the latest available version is used. + When taxon=True, this is instead an organism/taxon name (e.g. "homo sapiens"). + - json If True, returns the report as a list of dictionaries instead of + a pandas DataFrame. Default: False. + - save If True, saves the report to '{accession}_assembly_report.csv' + (or .json if json=True) in the current working directory. Default: False. + - taxon If True, `accession` is interpreted as an organism/taxon name and resolved to + that taxon's NCBI reference assembly accession (via the bundled datasets CLI) + before fetching the report. Default: False. + - list_assemblies If True, `accession` is interpreted as an organism/taxon name and, instead of a + report, a table of all NCBI assemblies for that taxon is returned (accession, + assembly_name, refseq_category, assembly_level, organism), reference/representative + first. Use it to find a specific accession to pass back in. Default: False. + - verbose True/False whether to print progress information. Default: True. + + Returns a pandas DataFrame (or list of dictionaries if json=True) with one row per + sequence and the columns provided by the NCBI assembly report + (Sequence-Name, Sequence-Role, Assigned-Molecule, Assigned-Molecule-Location/Type, + GenBank-Accn, Relationship, RefSeq-Accn, Assembly-Unit, Sequence-Length, UCSC-style-name). + """ + # When list_assemblies=True, return the taxon's assembly catalogue instead of a single report. + if list_assemblies: + return _list_taxon_assemblies(accession, json=json, save=save, verbose=verbose) + + # When taxon=True, resolve the organism/taxon name to its reference assembly accession first. + if taxon: + accession = _resolve_taxon_to_accession(accession, verbose=verbose) + + # Validate the accession format. The version suffix is optional: when omitted + # (e.g. "GCF_000001405"), the latest available version is used. + accession = accession.strip() if isinstance(accession, str) else accession + match = re.fullmatch(r"(GC[AF]_\d{9})(\.\d+)?", accession) if isinstance(accession, str) else None + if not match: + raise ValueError( + f"Invalid NCBI assembly accession: '{accession}'. " + "Expected format 'GCA_########.#' or 'GCF_########.#'; the version suffix is optional, " + "e.g. 'GCF_000001405.40' or 'GCF_000001405'.\n" + ) + base_accession = match.group(1) # accession without version, e.g. GCF_000001405 + has_version = match.group(2) is not None + + # Build the path to the assembly's FTP directory from the accession digits + # (the version does not affect this sharded path). + prefix = base_accession[:3] # GCA or GCF + digits = base_accession.split("_")[1] # e.g. 000001405 + parent_url = f"{NCBI_FTP_GENOMES_URL}{prefix}/{digits[0:3]}/{digits[3:6]}/{digits[6:9]}/" + + # List the parent directory to find the full assembly folder name (accession.version + assembly name) + parent_html = requests.get(parent_url, timeout=DEFAULT_REQUESTS_TIMEOUT) + if parent_html.status_code != 200: + raise RuntimeError( + f"NCBI assembly directory for accession '{accession}' returned status code " + f"{parent_html.status_code}. Please double-check the accession.\n" + ) + + soup = BeautifulSoup(parent_html.text, "html.parser") + hrefs = [href.rstrip("/") for href in (a.get("href", "") for a in soup.find_all("a"))] + if has_version: + # Match on the trailing underscore so that e.g. 'GCF_000001405.4' does not + # prefix-match the folder for 'GCF_000001405.40'. + folders = [h for h in hrefs if h.startswith(accession + "_")] + else: + # No version given: match any version of this accession and pick the latest. + versioned = re.compile(rf"{re.escape(base_accession)}\.(\d+)_") + candidates = [(int(m.group(1)), h) for h in hrefs if (m := versioned.match(h))] + folders = [max(candidates)[1]] if candidates else [] + if not folders: + raise RuntimeError( + f"No assembly folder found for accession '{accession}' at {parent_url}. " + "Please double-check the accession.\n" + ) + + folder = folders[0] + # Recover the accession with its resolved version from the folder name (e.g. GCF_000001405.40), + # so saved files and logs always carry the exact version even when the input omitted it. + resolved_accession = "_".join(folder.split("_")[:2]) + report_url = f"{parent_url}{folder}/{folder}_assembly_report.txt" + + if verbose: + logger.info(f"Fetching NCBI assembly report for {resolved_accession} from {report_url}.") + + report_html = requests.get(report_url, timeout=DEFAULT_REQUESTS_TIMEOUT) + if report_html.status_code != 200: + raise RuntimeError( + f"Assembly report for accession '{accession}' returned status code " + f"{report_html.status_code} ({report_url}).\n" + ) + + # Parse the tab-delimited report: comment lines start with '#'; the column + # header is the comment line that starts with '# Sequence-Name'. + columns: list[str] | None = None + rows: list[list[str]] = [] + for line in report_html.text.splitlines(): + if line.startswith("#"): + if line.startswith("# Sequence-Name"): + columns = line.lstrip("#").strip().split("\t") + continue + if not line.strip(): + continue + rows.append(line.split("\t")) + + if columns is None: + raise RuntimeError( + f"Could not parse the assembly report for accession '{accession}' " + f"(no column header found at {report_url}).\n" + ) + + df = pd.DataFrame(rows, columns=columns) + + if json: + result = cast("list[dict[str, Any]]", df.to_dict(orient="records")) + if save: + with open(f"{resolved_accession}_assembly_report.json", "w", encoding="utf-8") as f: + json_package.dump(result, f, ensure_ascii=False, indent=4) + return result + + if save: + df.to_csv(f"{resolved_accession}_assembly_report.csv", index=False) + + return df + + def find_FTP_link(url: str, link_substring: str) -> tuple[str | None, str | None, str | None]: """Helper function for gget ref to find an FTP link, its release date and size. @@ -58,6 +323,11 @@ def find_FTP_link(url: str, link_substring: str) -> tuple[str | None, str | None return link_str, date_str, size_str +# Module-level alias so ref() can delegate without its boolean ``assembly_report`` +# parameter shadowing the function of the same name. +_assembly_report_fn = assembly_report + + def ref( species: str | None, which: str | list[str] = "all", @@ -66,6 +336,10 @@ def ref( save: bool = False, list_species: bool = False, list_iv_species: bool = False, + assembly_report: bool = False, + json: bool = False, + taxon: bool = False, + list_assemblies: bool = False, verbose: bool = True, ) -> Any: """Fetch FTPs for reference genomes and annotations by species from Ensembl. @@ -74,6 +348,9 @@ def ref( - species Defines the species for which the reference should be fetched in the format "_", e.g. species = "homo_sapiens". Supported shortcuts: "human", "mouse", "human_grch37" (accesses the GRCh37 genome assembly) + When `assembly_report=True`, this is instead an NCBI assembly accession, + e.g. "GCF_000001405.40" (the version suffix is optional; if omitted, the + latest available version is used). - which Defines which results to return. Default: 'all' -> Returns all available results. Possible entries are one or a combination (as a list of strings) of the following: @@ -91,11 +368,34 @@ def ref( (Can be combined with the `release` argument to get the available species from a specific Ensembl release.) - list_iv_species If True and `species=None`, returns a list of all available INVERTEBRATE species from the Ensembl database (default: False). (Can be combined with the `release` argument to get the available species from a specific Ensembl release.) + - assembly_report If True, `species` is interpreted as an NCBI assembly accession (e.g. "GCF_000001405.40") + and the NCBI assembly report is returned instead of Ensembl FTP links. Useful for translating + sequence/chromosome names between Ensembl, GenBank, RefSeq and UCSC conventions (default: False). + - json Only used when `assembly_report=True`: if True, returns the report as a list of + dictionaries instead of a pandas DataFrame (default: False). + - taxon Only used when `assembly_report=True`: if True, `species` is interpreted as an + organism/taxon name (e.g. "homo sapiens") and resolved to that taxon's NCBI + reference assembly accession before fetching the report (default: False). + - list_assemblies Only used when `assembly_report=True`: if True, `species` is interpreted as an + organism/taxon name and a table of all NCBI assemblies for that taxon is + returned instead of a report (default: False). - verbose True/False whether to print progress information (default: True). Returns a dictionary containing the requested URLs with their respective Ensembl version and release date and time. (If FTP=True, returns a list containing only the URLs.) + (If assembly_report=True, returns the NCBI assembly report as a pandas DataFrame, or a list of dictionaries if json=True.) """ + # Return the NCBI assembly report instead of Ensembl FTP links + if assembly_report: + if species is None: + raise ValueError( + "An NCBI assembly accession (e.g. 'GCF_000001405.40') must be provided as the " + "`species` argument when `assembly_report=True`.\n" + ) + return _assembly_report_fn( + species, json=json, save=save, taxon=taxon, list_assemblies=list_assemblies, verbose=verbose + ) + # Return list of all available species if list_species: if release is None: @@ -481,7 +781,7 @@ def ref( if save: with open("gget_ref_results.json", "w", encoding="utf-8") as file: - json.dump(ref_dict, file, ensure_ascii=False, indent=4) + json_package.dump(ref_dict, file, ensure_ascii=False, indent=4) if verbose: logger.info(f"Fetching reference information for {species} from Ensembl release: {ENS_rel}.") return ref_dict diff --git a/gget/main.py b/gget/main.py index 7a2944b09..c3553e223 100644 --- a/gget/main.py +++ b/gget/main.py @@ -199,6 +199,51 @@ def main() -> None: required=False, help="Return only the FTP link(s).", ) + parser_ref.add_argument( + "-ar", + "--assembly_report", + default=False, + action="store_true", + required=False, + help=( + "Return the NCBI assembly report instead of Ensembl FTP links.\n" + "When set, the positional argument is interpreted as an NCBI assembly accession,\n" + "e.g. 'gget ref GCF_000001405.40 --assembly_report'.\n" + "Useful for translating sequence/chromosome names between Ensembl, GenBank, RefSeq and UCSC." + ), + ) + parser_ref.add_argument( + "-tx", + "--taxon", + default=False, + action="store_true", + required=False, + help=( + "Only used with --assembly_report: interpret the positional argument as an organism/taxon\n" + "name (e.g. 'gget ref \"homo sapiens\" --assembly_report --taxon') and resolve it to that\n" + "taxon's NCBI reference assembly before fetching the report." + ), + ) + parser_ref.add_argument( + "-la", + "--list_assemblies", + default=False, + action="store_true", + required=False, + help=( + "Only used with --assembly_report: interpret the positional argument as an organism/taxon\n" + "name and list all of that taxon's NCBI assemblies (accession, name, category, level) instead\n" + "of a report, so you can pick a specific accession to fetch." + ), + ) + parser_ref.add_argument( + "-csv", + "--csv", + default=True, + action="store_false", + required=False, + help="Only used with --assembly_report: return the report in csv format instead of json.", + ) parser_ref.add_argument( "-d", "--download", @@ -3323,6 +3368,42 @@ def main() -> None: ## ref return if args.command == "ref": + # Return the NCBI assembly report instead of Ensembl FTP links + if args.assembly_report: + if args.species is None: + parser_ref.error( + "the following argument is required with --assembly_report: an NCBI assembly accession\n" + "Example: 'gget ref GCF_000001405.40 --assembly_report'" + ) + + report_results = ref( + species=args.species, + assembly_report=True, + json=args.csv, + taxon=args.taxon, + list_assemblies=args.list_assemblies, + verbose=args.quiet, + ) + + # Save in specified file if -o specified + if args.out: + directory = "/".join(args.out.split("/")[:-1]) + if directory != "": + os.makedirs(directory, exist_ok=True) + if args.csv: + with open(args.out, "w", encoding="utf-8") as f: + json.dump(report_results, f, ensure_ascii=False, indent=4) + else: + report_results.to_csv(args.out, index=False) + # Otherwise print to standard out + else: + if args.csv: + print(json.dumps(report_results, ensure_ascii=False, indent=4)) + else: + report_results.to_csv(sys.stdout, index=False) + + return + # Return all vertebrate available species if args.list_species: species_list = ref(species=None, release=args.release, list_species=args.list_species) diff --git a/tests/fixtures/test_ref.json b/tests/fixtures/test_ref.json index c88a364db..b022d2953 100644 --- a/tests/fixtures/test_ref.json +++ b/tests/fixtures/test_ref.json @@ -648,5 +648,14 @@ "ftp": false }, "expected_result": "RuntimeError" + }, + "test_ref_assembly_report_bad_accession": { + "type": "error", + "args": { + "species": "banana", + "assembly_report": true + }, + "expected_result": "ValueError", + "expected_msg": "Invalid NCBI assembly accession: 'banana'. Expected format 'GCA_########.#' or 'GCF_########.#'; the version suffix is optional, e.g. 'GCF_000001405.40' or 'GCF_000001405'.\n" } } \ No newline at end of file diff --git a/tests/test_ref.py b/tests/test_ref.py index c6ab3e0e8..d07e177ed 100644 --- a/tests/test_ref.py +++ b/tests/test_ref.py @@ -1,7 +1,12 @@ import json +import os +import tempfile import unittest +from unittest.mock import Mock, patch -from gget.gget_ref import ref +import gget.gget_ref as gget_ref +import requests +from gget.gget_ref import assembly_report, ref from .from_json import from_json @@ -12,3 +17,224 @@ class TestRef(unittest.TestCase, metaclass=from_json(ref_dict, ref)): pass # all tests are loaded from json + + +class _FakeResp: + """Minimal stand-in for a requests.Response used to test assembly_report offline.""" + + def __init__(self, text="", status_code=200): + self.text = text + self.status_code = status_code + + +# A parent-directory listing containing the assembly's folder, and a minimal +# tab-delimited assembly report (with a comment header and a blank line). +_PARENT_HTML = 'GCF_000001405.40_GRCh38.p14/' +_REPORT_TEXT = ( + "# Assembly name: GRCh38.p14\n" + "# Sequence-Name\tSequence-Role\tGenBank-Accn\tRefSeq-Accn\tUCSC-style-name\n" + "1\tassembled-molecule\tCM000663.2\tNC_000001.11\tchr1\n" + "\n" +) + + +class TestAssemblyReportOffline(unittest.TestCase): + """Network-free tests of assembly_report and the ref() delegation (issue #179).""" + + def test_invalid_accession_raises(self): + with self.assertRaises(ValueError): + assembly_report("not-an-accession", verbose=False) + + @patch.object(gget_ref.requests, "get") + def test_parse_verbose_json_and_save(self, mock_get): + # Happy path: verbose log, blank-line skip, parsing, then json/save/csv branches. + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp(_REPORT_TEXT)] + df = assembly_report("GCF_000001405.40", verbose=True) + self.assertEqual(list(df["Sequence-Name"]), ["1"]) + self.assertEqual(df.iloc[0]["UCSC-style-name"], "chr1") + + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp(_REPORT_TEXT)] + result = assembly_report("GCF_000001405.40", json=True, verbose=False) + self.assertIsInstance(result, list) + self.assertEqual(result[0]["RefSeq-Accn"], "NC_000001.11") + + with tempfile.TemporaryDirectory() as tmp: + cwd = os.getcwd() + os.chdir(tmp) + try: + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp(_REPORT_TEXT)] + assembly_report("GCF_000001405.40", save=True, verbose=False) + self.assertTrue(os.path.exists("GCF_000001405.40_assembly_report.csv")) + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp(_REPORT_TEXT)] + assembly_report("GCF_000001405.40", json=True, save=True, verbose=False) + self.assertTrue(os.path.exists("GCF_000001405.40_assembly_report.json")) + finally: + os.chdir(cwd) + + @patch.object(gget_ref.requests, "get") + def test_parent_dir_non_200_raises(self, mock_get): + mock_get.return_value = _FakeResp("", status_code=404) + with self.assertRaises(RuntimeError): + assembly_report("GCF_000001405.40", verbose=False) + + @patch.object(gget_ref.requests, "get") + def test_no_folder_found_raises(self, mock_get): + mock_get.return_value = _FakeResp("no matching links", status_code=200) + with self.assertRaises(RuntimeError): + assembly_report("GCF_000001405.40", verbose=False) + + @patch.object(gget_ref.requests, "get") + def test_report_non_200_raises(self, mock_get): + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp("", status_code=500)] + with self.assertRaises(RuntimeError): + assembly_report("GCF_000001405.40", verbose=False) + + @patch.object(gget_ref.requests, "get") + def test_report_missing_header_raises(self, mock_get): + # No "# Sequence-Name" header line -> columns stay None -> RuntimeError. + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp("# only comments\nfoo\tbar\n")] + with self.assertRaises(RuntimeError): + assembly_report("GCF_000001405.40", verbose=False) + + def test_ref_assembly_report_requires_species(self): + with self.assertRaises(ValueError): + ref(None, assembly_report=True, verbose=False) + + @patch.object(gget_ref, "_assembly_report_fn") + def test_ref_delegates_to_assembly_report(self, mock_fn): + mock_fn.return_value = "DELEGATED" + self.assertEqual(ref("GCF_000001405.40", assembly_report=True, verbose=False), "DELEGATED") + mock_fn.assert_called_once() + + @patch.object(gget_ref.requests, "get") + def test_version_prefix_no_collision(self, mock_get): + # Only version .40 exists; querying .4 must NOT prefix-match the .40 folder. + mock_get.return_value = _FakeResp('x') + with self.assertRaises(RuntimeError): + assembly_report("GCF_000001405.4", verbose=False) + + @patch.object(gget_ref.requests, "get") + def test_versionless_resolves_latest(self, mock_get): + # No version given -> resolve to the latest version folder (.40, not .39). + parent = 'ab' + with tempfile.TemporaryDirectory() as tmp: + cwd = os.getcwd() + os.chdir(tmp) + try: + mock_get.side_effect = [_FakeResp(parent), _FakeResp(_REPORT_TEXT)] + assembly_report("GCF_000001405", save=True, verbose=False) + # The resolved (latest) version is reflected in the saved filename. + self.assertTrue(os.path.exists("GCF_000001405.40_assembly_report.csv")) + finally: + os.chdir(cwd) + + @patch.object(gget_ref.requests, "get") + @patch.object(gget_ref, "_resolve_taxon_to_accession", return_value="GCF_000001405.40") + def test_taxon_resolves_then_fetches(self, mock_resolve, mock_get): + # taxon=True routes the input through the resolver, then fetches the report normally. + mock_get.side_effect = [_FakeResp(_PARENT_HTML), _FakeResp(_REPORT_TEXT)] + df = assembly_report("homo sapiens", taxon=True, verbose=False) + mock_resolve.assert_called_once_with("homo sapiens", verbose=False) + self.assertEqual(df.iloc[0]["UCSC-style-name"], "chr1") + + @patch("gget.gget_virus._get_datasets_path", return_value="datasets") + @patch.object(gget_ref.subprocess, "run") + def test_resolve_taxon_parses_accession(self, mock_run, _mock_path): + mock_run.return_value = Mock( + returncode=0, + stdout='new version notice\n{"accession": "GCF_000001405.40", "organism": {}}\n', + ) + self.assertEqual(gget_ref._resolve_taxon_to_accession("homo sapiens", verbose=False), "GCF_000001405.40") + + @patch("gget.gget_virus._get_datasets_path", return_value="datasets") + @patch.object(gget_ref.subprocess, "run") + def test_resolve_taxon_not_found_raises(self, mock_run, _mock_path): + mock_run.return_value = Mock(returncode=0, stdout="") # no records -> unknown taxon + with self.assertRaises(ValueError): + gget_ref._resolve_taxon_to_accession("not-a-species", verbose=False) + + @patch("gget.gget_virus._get_datasets_path", return_value="datasets") + @patch.object(gget_ref.subprocess, "run") + def test_resolve_taxon_datasets_failure_raises(self, mock_run, _mock_path): + mock_run.return_value = Mock(returncode=1, stdout="") # datasets CLI error (often transient) + with self.assertRaises(RuntimeError): + gget_ref._resolve_taxon_to_accession("homo sapiens", verbose=False) + + @patch("gget.gget_virus._get_datasets_path", return_value="datasets") + @patch.object(gget_ref.subprocess, "run") + def test_list_taxon_assemblies_columns_and_order(self, mock_run, _mock_path): + # Two records; the "reference genome" must sort ahead of the "na" one. + mock_run.return_value = Mock( + returncode=0, + stdout=( + '{"accession": "GCF_000002125.1", "assembly_info": {"assembly_name": "HuRef", ' + '"assembly_level": "Chromosome"}, "organism": {"organism_name": "Homo sapiens"}}\n' + '{"accession": "GCF_000001405.40", "assembly_info": {"assembly_name": "GRCh38.p14", ' + '"refseq_category": "reference genome", "assembly_level": "Chromosome"}, ' + '"organism": {"organism_name": "Homo sapiens"}}\n' + ), + ) + df = gget_ref._list_taxon_assemblies("homo sapiens", verbose=False) + self.assertEqual( + list(df.columns), ["accession", "assembly_name", "refseq_category", "assembly_level", "organism"] + ) + self.assertEqual(df.iloc[0]["accession"], "GCF_000001405.40") # reference first + self.assertEqual(df.iloc[1]["refseq_category"], "na") # missing category -> "na" + + @patch("gget.gget_virus._get_datasets_path", return_value="datasets") + @patch.object(gget_ref.subprocess, "run") + def test_list_taxon_assemblies_empty_raises(self, mock_run, _mock_path): + mock_run.return_value = Mock(returncode=0, stdout="") + with self.assertRaises(ValueError): + gget_ref._list_taxon_assemblies("not-a-species", verbose=False) + + @patch.object(gget_ref, "_list_taxon_assemblies") + def test_assembly_report_list_assemblies_delegates(self, mock_list): + mock_list.return_value = "CATALOGUE" + self.assertEqual(assembly_report("homo sapiens", list_assemblies=True, verbose=False), "CATALOGUE") + mock_list.assert_called_once() + + +class TestAssemblyReportLive(unittest.TestCase): + """Live test hitting the real NCBI FTP (issue #179). + + Anchored to SARS-CoV-2 (GCF_009858895.2), a frozen single-sequence reference, so + the identity columns are stable. Skips (rather than fails) on a network error or a + transient non-200 instead of reddening CI. + """ + + def test_sars_cov2_report_key_columns(self): + try: + df = assembly_report("GCF_009858895.2", verbose=False) + except requests.RequestException as e: + self.skipTest(f"Network error reaching NCBI FTP: {e}") + except RuntimeError as e: + self.skipTest(f"NCBI FTP transient error: {e}") + # Anchor on stable identity columns, not the full row. + self.assertEqual(len(df), 1) + row = df.iloc[0] + self.assertEqual(row["RefSeq-Accn"], "NC_045512.2") + self.assertEqual(row["GenBank-Accn"], "MN908947.3") + self.assertEqual(row["Sequence-Length"], "29903") + + def test_taxon_human_resolves_to_reference(self): + # taxon=True resolves "homo sapiens" (via the bundled datasets CLI) to the human + # reference assembly and returns its report, which must contain chromosome 1. + try: + df = assembly_report("homo sapiens", taxon=True, verbose=False) + except requests.RequestException as e: + self.skipTest(f"Network error reaching NCBI: {e}") + except RuntimeError as e: + self.skipTest(f"NCBI datasets/FTP transient error: {e}") + self.assertIn("NC_000001.11", set(df["RefSeq-Accn"])) + + def test_list_assemblies_human_contains_reference(self): + # list_assemblies=True returns the taxon's assembly catalogue; the human reference must be in it. + try: + df = assembly_report("homo sapiens", list_assemblies=True, verbose=False) + except requests.RequestException as e: + self.skipTest(f"Network error reaching NCBI: {e}") + except RuntimeError as e: + self.skipTest(f"NCBI datasets transient error: {e}") + self.assertIn("GCF_000001405.40", set(df["accession"])) + self.assertEqual(df.iloc[0]["refseq_category"], "reference genome") # reference sorted first