""" Step 1. Fetch and count what the public sequence, proteomics and metabolomics archives hold (Part 1 of the study). Everything is cached under RAW (default fermi_raw/second_origin_search/); a cached file is reused unless --refresh is given. The small derived tables go to ../data/. Database totals change daily, so every table carries the access date; the numbers in the report are the snapshot of 2026-10-12. Sources (all plain HTTPS, no key): MGnify API https://www.ebi.ac.uk/metagenomics/api/v1/ GTDB release files https://data.gtdb.ecogenomic.org/releases/latest/ NCBI E-utilities https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi (record counts only) ENA portal API https://www.ebi.ac.uk/ena/portal/api/count PRIDE, MetaboLights record counts Usage: python3 fetch_public_counts.py [--refresh] [--skip-big] """ import csv import datetime import gzip import json import os import sys import time import urllib.parse import urllib.request HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "..", "data") RAW = os.environ.get("SOS_RAW", "fermi_raw/second_origin_search") REFRESH = "--refresh" in sys.argv SKIP_BIG = "--skip-big" in sys.argv MG = "https://www.ebi.ac.uk/metagenomics/api/v1/" TODAY = datetime.date.today().isoformat() def get(url, cache, binary=False, tries=4, timeout=300): """Download url to RAW/cache (reused if present).""" path = os.path.join(RAW, cache) os.makedirs(os.path.dirname(path), exist_ok=True) if os.path.exists(path) and os.path.getsize(path) > 0 and not REFRESH: return open(path, "rb").read() if binary else open(path, encoding="utf-8", errors="replace").read() err = None for k in range(tries): try: req = urllib.request.Request(url, headers={"User-Agent": "fermi-second-origin-study/1.0"}) buf = urllib.request.urlopen(req, timeout=timeout).read() with open(path + ".part", "wb") as f: f.write(buf) os.replace(path + ".part", path) return buf if binary else buf.decode("utf-8", errors="replace") except Exception as e: # network hiccup: wait and retry err = e time.sleep(3 + 5 * k) raise RuntimeError(f"failed {url}: {err}") def jget(url, cache): return json.loads(get(url, cache)) def mg_count(endpoint, **params): q = urllib.parse.urlencode(dict(page_size=1, **params)) tag = endpoint.replace("/", "_") + "_" + "_".join(f"{k}-{v}" for k, v in sorted(params.items())).replace(":", ".").replace(" ", "") j = jget(MG + endpoint + "?" + q, f"mgnify/count_{tag}.json") return int(j["meta"]["pagination"]["count"]) def summary_dict(attrs): return {kv["key"]: kv["value"] for kv in (attrs.get("analysis-summary") or [])} def main(): os.makedirs(DATA, exist_ok=True) counts = [] # (database, quantity, value, url, access date) def add(db, what, val, url): counts.append(dict(database=db, quantity=what, value=val, source_url=url, accessed=TODAY)) print(f"{db:14s} {what:70s} {val}", flush=True) # ---------------- MGnify totals ---------------- for ep, label in [("studies", "studies"), ("samples", "samples"), ("runs", "runs"), ("assemblies", "assemblies"), ("analyses", "analyses (all pipeline versions)"), ("genomes", "species-level genomes in the genome catalogues"), ("biomes", "biome classes (GOLD ecosystem tree)")]: add("MGnify", label, mg_count(ep), MG + ep) # (Counting analyses by experiment type times out on the server for the large types, so the per-type numbers below # are the sample and run counts that the experiment-types endpoint reports, plus pipeline-5.0 analysis counts.) for et in ("assembly", "metagenomic", "metatranscriptomic"): try: add("MGnify", f"analyses, pipeline 5.0, experiment type = {et}", mg_count("analyses", experiment_type=et, pipeline_version="5.0"), MG + "analyses?pipeline_version=5.0&experiment_type=" + et) except Exception as e: print("count failed:", et, e, flush=True) j = jget(MG + "experiment-types", "mgnify/experiment_types.json") for d in j["data"]: a = d["attributes"] if a["samples-count"]: add("MGnify", f"samples with experiment type = {d['id']}", a["samples-count"], MG + "experiment-types") add("MGnify", f"runs with experiment type = {d['id']}", a["runs-count"], MG + "experiment-types") j = jget(MG + "genome-catalogues?page_size=100", "mgnify/genome_catalogues.json") cat_rows = [] for d in j["data"]: a = d["attributes"] cat_rows.append(dict(catalogue=d["id"], species_genomes=a.get("genome-count"), genomes_before_clustering=a.get("unclustered-genome-count"), version=a.get("version"), last_update=(a.get("last-update") or "")[:10])) add("MGnify", "genome catalogues", len(cat_rows), MG + "genome-catalogues") add("MGnify", "genomes before species clustering, all catalogues", sum(r["genomes_before_clustering"] or 0 for r in cat_rows), MG + "genome-catalogues") with open(os.path.join(DATA, "mgnify_genome_catalogues.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(cat_rows[0])); w.writeheader(); w.writerows(cat_rows) # ---------------- MGnify biome tree: samples per biome ---------------- rows = [] page = 1 while True: j = jget(MG + f"biomes?page_size=250&page={page}", f"mgnify/biomes_p{page}.json") for d in j["data"]: rows.append(dict(lineage=d["attributes"]["lineage"], samples=d["attributes"]["samples-count"])) if not j["links"].get("next"): break page += 1 rows.sort(key=lambda r: r["lineage"]) with open(os.path.join(DATA, "mgnify_biomes.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=["lineage", "samples"]); w.writeheader(); w.writerows(rows) add("MGnify", "biome classes with at least one sample", sum(1 for r in rows if r["samples"] > 0), MG + "biomes") # ---------------- GTDB ---------------- G = "https://data.gtdb.ecogenomic.org/releases/latest/" ver = get(G + "VERSION.txt", "gtdb/VERSION.txt").split()[0] notes = get(G + "RELEASE_NOTES.txt", "gtdb/RELEASE_NOTES.txt") add("GTDB", "release", ver, G + "VERSION.txt") tax = {} for dom, fn in (("Bacteria", "bac120_taxonomy.tsv.gz"), ("Archaea", "ar53_taxonomy.tsv.gz")): buf = get(G + fn, "gtdb/" + fn, binary=True) ranks = [set() for _ in range(7)] n = 0 for line in gzip.decompress(buf).decode().splitlines(): if not line.strip(): continue n += 1 for i, t in enumerate(line.split("\t")[1].split(";")): ranks[i].add(t) tax[dom] = dict(genomes=n, **{r: len(s) for r, s in zip(["domain", "phylum", "class", "order", "family", "genus", "species"], ranks)}) add("GTDB", f"{dom}: genomes", n, G + fn) for r in ("phylum", "class", "order", "family", "genus", "species"): add("GTDB", f"{dom}: {r} count", tax[dom][r], G + fn) add("GTDB", "genomes (both domains)", tax["Bacteria"]["genomes"] + tax["Archaea"]["genomes"], G) add("GTDB", "species clusters (both domains)", tax["Bacteria"]["species"] + tax["Archaea"]["species"], G) # how many species clusters are known only from environmental DNA (no cultured isolate genome)? needs the metadata tables if not SKIP_BIG: cat_by_sp = {} for dom, fn in (("Bacteria", "bac120_metadata.tsv.gz"), ("Archaea", "ar53_metadata.tsv.gz")): path = os.path.join(RAW, "gtdb", fn) if not (os.path.exists(path) and os.path.getsize(path) > 1e6) or REFRESH: print("downloading", fn, "(large) ...", flush=True) urllib.request.urlretrieve(G + fn, path + ".part"); os.replace(path + ".part", path) with gzip.open(path, "rt", encoding="utf-8", errors="replace") as f: hdr = f.readline().rstrip("\n").split("\t") i_tax, i_cat = hdr.index("gtdb_taxonomy"), hdr.index("ncbi_genome_category") for line in f: p = line.rstrip("\n").split("\t") sp = (dom, p[i_tax]) c = p[i_cat] env = ("metagenome" in c) or ("single cell" in c) or ("environmental" in c) st = cat_by_sp.setdefault(sp, [0, 0]) st[1 if env else 0] += 1 for dom in ("Bacteria", "Archaea"): sp = [v for (d, _), v in cat_by_sp.items() if d == dom] only_env = sum(1 for v in sp if v[0] == 0) add("GTDB", f"{dom}: species clusters with no isolate genome (metagenome- or single-cell-derived only)", only_env, G + ("bac120" if dom == "Bacteria" else "ar53") + "_metadata.tsv.gz") add("GTDB", f"{dom}: share of species clusters with no isolate genome", round(only_env / len(sp), 4), G) add("GTDB", f"{dom}: genomes derived from metagenomes or single cells", sum(v[1] for v in sp), G) # ---------------- NCBI, ENA, PRIDE, MetaboLights record counts ---------------- EU = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?retmax=0&retmode=json&" for db, term, label in [("sra", "metagenomic[Source]", "SRA experiments, library source = METAGENOMIC"), ("sra", "metatranscriptomic[Source]", "SRA experiments, library source = METATRANSCRIPTOMIC"), ("sra", "all[filter]", "SRA experiments, all"), ("assembly", "all[filter]", "genome assemblies, all"), ("assembly", "\"derived from metagenome\"[filter]", "genome assemblies derived from metagenomes"), ("biosample", "all[filter]", "BioSamples, all"), ("bioproject", "all[filter]", "BioProjects, all"), ("taxonomy", "species[Rank]", "Taxonomy: named or provisional species"), ("nuccore", "all[filter]", "nucleotide records, all"), ("protein", "all[filter]", "protein records, all")]: u = EU + urllib.parse.urlencode(dict(db=db, term=term)) try: j = jget(u, f"ncbi/{db}_{''.join(ch if ch.isalnum() else '_' for ch in term)}.json") add("NCBI", label, int(j["esearchresult"]["count"]), u) except Exception as e: print("NCBI count failed:", label, e) time.sleep(0.5) for q, label in [('library_source="METAGENOMIC"', "ENA read runs, library source = METAGENOMIC"), ('library_source="METATRANSCRIPTOMIC"', "ENA read runs, library source = METATRANSCRIPTOMIC")]: u = "https://www.ebi.ac.uk/ena/portal/api/count?" + urllib.parse.urlencode(dict(result="read_run", query=q)) try: t = get(u, f"ena/count_{label.split('=')[-1].strip()}.txt") add("ENA", label, int(t.strip().split()[-1]), u) except Exception as e: print("ENA count failed:", e) try: add("PRIDE", "proteomics projects (mass spectrometry)", int(get("https://www.ebi.ac.uk/pride/ws/archive/v2/projects/count", "pride/count.txt").strip()), "https://www.ebi.ac.uk/pride/ws/archive/v2/projects/count") except Exception as e: print("PRIDE count failed:", e) try: j = jget("https://www.ebi.ac.uk/metabolights/ws/studies", "metabolights/studies.json") add("MetaboLights", "public metabolomics studies", len(j["content"]), "https://www.ebi.ac.uk/metabolights/ws/studies") except Exception as e: print("MetaboLights count failed:", e) with open(os.path.join(DATA, "db_counts.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(counts[0])); w.writeheader(); w.writerows(counts) if SKIP_BIG: return # ---------------- MGnify per-analysis summaries (pipeline 5.0): the share of sequence with no database match ---------------- # Top-level lineages. Everything in Environmental, Engineered and Mixed; every 4th page of Host-associated (it is 70% of the archive). out = [] plan = [("assembly", "root:Environmental", 1), ("assembly", "root:Engineered", 1), ("assembly", "root:Mixed", 1), ("assembly", "root:Host-associated", 4), ("metagenomic", "root:Environmental", 1), ("metagenomic", "root:Engineered", 1), ("metagenomic", "root:Mixed", 1), ("metagenomic", "root:Host-associated", 1)] for et, lin, step in plan: page = 1 npages = None while npages is None or page <= npages: u = MG + "analyses?" + urllib.parse.urlencode(dict(experiment_type=et, pipeline_version="5.0", page_size=250, lineage=lin, page=page)) try: j = jget(u, f"mgnify/analyses/{et}_{lin.replace(':', '.')}_p{page}.json") except Exception as e: print("page failed, skipped:", et, lin, page, e, flush=True) page += step continue npages = j["meta"]["pagination"]["pages"] for d in j["data"]: s = summary_dict(d["attributes"]) rel = d["relationships"] out.append(dict(analysis=d["id"], experiment_type=et, top_biome=lin, page_step=step, study=(rel.get("study", {}).get("data") or {}).get("id", ""), sample=(rel.get("sample", {}).get("data") or {}).get("id", ""), seqs=s.get("Nucleotide sequences after undetermined bases filtering", ""), cds=s.get("Predicted CDS", ""), cds_ipr=s.get("Predicted CDS with InterProScan match", ""), contigs_cds=s.get("Contigs with predicted CDS", ""), contigs_ipr=s.get("Contigs with InterProScan match", ""), contigs_rna=s.get("Contigs with predicted RNA", ""), reads_cds=s.get("Reads with predicted CDS", ""), reads_ipr=s.get("Reads with InterProScan match", ""), reads_rna=s.get("Reads with predicted RNA", ""), ssu=s.get("Predicted SSU sequences", ""), lsu=s.get("Predicted LSU sequences", ""))) print(f"analyses {et} {lin} page {page}/{npages} rows so far {len(out)}", flush=True) page += step with open(os.path.join(DATA, "mgnify_analysis_summaries.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(out[0])); w.writeheader(); w.writerows(out) # ---------------- MGnify genome catalogue: how much of every catalogued genome matches known protein families ---------------- g = [] page = 1 npages = None while npages is None or page <= npages: try: j = jget(MG + f"genomes?page_size=250&page={page}", f"mgnify/genomes/p{page}.json") except Exception as e: print("genome page failed, skipped:", page, e, flush=True) page += 1 continue npages = j["meta"]["pagination"]["pages"] for d in j["data"]: a = d["attributes"] g.append(dict(genome=d["id"], catalogue=(d["relationships"].get("catalogue", {}).get("data") or {}).get("id", ""), type=a.get("type"), length=a.get("length"), completeness=a.get("completeness"), contamination=a.get("contamination"), num_proteins=a.get("num-proteins"), ipr_coverage=a.get("ipr-coverage"), eggnog_coverage=a.get("eggnog-coverage"), rna_16s=a.get("rna-16s"), trnas=a.get("trnas"), lineage=a.get("taxon-lineage"))) if page % 20 == 0: print(f"genomes page {page}/{npages}", flush=True) page += 1 with open(os.path.join(DATA, "mgnify_genomes.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(g[0])); w.writeheader(); w.writerows(g) print("done") if __name__ == "__main__": main()