"""
Step 2. Part 1 numbers: how much of Earth's biosphere has been sequenced, how much of the sequence goes unassigned,
and what the unassigned part turns out to be.

Inputs : ../data/db_counts.csv, ../data/mgnify_analysis_summaries.csv, ../data/mgnify_genomes.csv, ../data/mgnify_biomes.csv,
         ../data/pubmed_method_counts.csv        (all written by fetch_public_counts.py / fetch_pubmed_counts.py)
Outputs: ../data/earth_unassigned_by_biome.csv, ../data/mgnify_least_annotated_studies.csv, ../data/earth_stats.json
"""
import csv
import json
import os

import numpy as np
import pandas as pd

from fetch_public_counts import MG, jget

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "..", "data")
MIN_CONTIGS = 1000          # assemblies with fewer sequences are dropped (fragments, single genomes): the ratio is too noisy
MIN_READS = 10000
HOST_WEIGHT = 4.0           # Host-associated assemblies were fetched at every fourth API page (fetch_public_counts.py)


def q(x, p):
    return float(np.nanquantile(x, p))


def main():
    out = {}
    counts = pd.read_csv(os.path.join(DATA, "db_counts.csv"))
    out["db_counts"] = {f"{r.database}: {r.quantity}": r.value for r in counts.itertuples()}
    out["accessed"] = str(counts.accessed.iloc[0])

    d = pd.read_csv(os.path.join(DATA, "mgnify_analysis_summaries.csv"))
    for c in ["seqs", "cds", "cds_ipr", "reads_cds", "reads_ipr", "reads_rna", "ssu", "lsu"]:
        d[c] = pd.to_numeric(d[c], errors="coerce")
    d["biome"] = d.top_biome.str.replace("root:", "", regex=False)
    a = d[(d.experiment_type == "assembly")].copy()
    n_all = len(a)
    a = a[(a.seqs >= MIN_CONTIGS) & (a.cds > 0)].copy()
    a["no_match"] = 1 - a.cds_ipr / a.cds                    # share of predicted proteins with no InterPro signature
    a["has_rrna"] = (a.ssu.fillna(0) + a.lsu.fillna(0)) > 0
    m = d[(d.experiment_type == "metagenomic")].copy()
    n_all_m = len(m)
    m = m[m.seqs >= MIN_READS].copy()
    m["read_no_orf"] = 1 - m.reads_cds / m.seqs              # reads in which the gene caller finds no open reading frame
    m["read_no_match"] = 1 - m.reads_ipr / m.seqs            # reads with no InterPro-matched protein fragment
    rows = []
    for name, sub, col in (("assembled contigs: predicted proteins with no InterPro match", a, "no_match"),
                           ("raw reads: no InterPro-matched fragment", m, "read_no_match"),
                           ("raw reads: no open reading frame predicted", m, "read_no_orf")):
        for b, g in list(sub.groupby("biome")) + [("All (unweighted pool of the analyses fetched)", sub)]:
            rows.append(dict(measure=name, biome=b, analyses=len(g), studies=g.study.nunique(), samples=g["sample"].nunique(),
                             p05=q(g[col], 0.05), p25=q(g[col], 0.25), median=q(g[col], 0.5), p75=q(g[col], 0.75), p95=q(g[col], 0.95), mean=float(g[col].mean())))
    pd.DataFrame(rows).to_csv(os.path.join(DATA, "earth_unassigned_by_biome.csv"), index=False)
    out["assemblies_fetched"] = n_all; out["assemblies_used"] = len(a); out["read_analyses_fetched"] = n_all_m; out["read_analyses_used"] = len(m)
    # ---- how representative are these medians? (added in revision) ----
    # They are unweighted medians over analyses, and the analyses come from few studies. Host-associated assemblies were fetched at every
    # fourth API page, so a weight of 4 restores their share; the pooled figure weights every predicted protein equally.
    wt = np.where(a.top_biome == "root:Host-associated", HOST_WEIGHT, 1.0)
    order = np.argsort(a.no_match.values); cw = np.cumsum(wt[order]) / wt.sum()
    vc = a.study.value_counts()
    per_study = a.groupby("study").no_match.median()
    out["assembly_unmatched"] = dict(
        median_unweighted=q(a.no_match, 0.5), p05=q(a.no_match, 0.05), p95=q(a.no_match, 0.95),
        median_host_associated_weighted_x4=float(a.no_match.values[order][np.searchsorted(cw, 0.5)]),
        pooled_share_protein_weighted=float(1 - a.cds_ipr.sum() / a.cds.sum()),
        median_of_study_medians=float(per_study.median()), studies=int(a.study.nunique()), samples=int(a["sample"].nunique()),
        analyses_from_two_largest_studies=int(vc.iloc[:2].sum()), share_from_two_largest_studies=float(vc.iloc[:2].sum() / len(a)),
        analyses_by_biome={k.replace("root:", ""): int(v) for k, v in a.top_biome.value_counts().items()})
    vcm = m.study.value_counts()
    out["read_unmatched"] = dict(median_unweighted=q(m.read_no_match, 0.5), p05=q(m.read_no_match, 0.05), p95=q(m.read_no_match, 0.95),
                                 pooled_share_read_weighted=float(1 - m.reads_ipr.sum() / m.seqs.sum()), median_of_study_medians=float(m.groupby("study").read_no_match.median().median()),
                                 studies=int(m.study.nunique()), analyses_from_two_largest_studies=int(vcm.iloc[:2].sum()), share_from_two_largest_studies=float(vcm.iloc[:2].sum() / len(m)))
    # NOTE: the share of assemblies carrying rRNA genes says nothing about the unassigned fraction (any mixed community contains standard
    # cells); it is kept only as a description of the holdings and is not used as evidence in the report.
    out["assemblies_with_rrna_share"] = float(a.has_rrna.mean())
    out["assemblies_without_rrna"] = int((~a.has_rrna).sum())
    out["assemblies_over_80pct_unmatched"] = int((a.no_match > 0.8).sum())
    out["assemblies_over_95pct_unmatched"] = int((a.no_match > 0.95).sum())
    out["max_unmatched_share"] = float(a.no_match.max())
    # which studies sit in the least-annotated tail, and what are they?
    g = a.groupby("study").agg(analyses=("no_match", "size"), median_unmatched=("no_match", "median"), biome=("biome", "first")).sort_values("median_unmatched", ascending=False)
    top = g.head(15).reset_index()
    names, lineages = [], []
    for st in top.study:
        try:
            j = jget(MG + f"studies/{st}", f"mgnify/studies/{st}.json")
            names.append(j["data"]["attributes"]["study-name"])
            lineages.append("; ".join(b["id"] for b in j["data"]["relationships"]["biomes"]["data"][:2]))
        except Exception as e:
            names.append(f"(lookup failed: {e})"); lineages.append("")
    top["study_name"] = names; top["mgnify_biome"] = lineages
    top.to_csv(os.path.join(DATA, "mgnify_least_annotated_studies.csv"), index=False, quoting=csv.QUOTE_NONNUMERIC)
    out["studies_with_median_unmatched_over_60pct"] = int((g.median_unmatched > 0.6).sum()); out["studies_total"] = int(len(g))

    # the genome catalogue: every catalogued genome sits inside the known protein universe (by construction: marker genes are required)
    gpath = os.path.join(DATA, "mgnify_genomes.csv")
    if os.path.exists(gpath):
        G = pd.read_csv(gpath)
        for c in ["ipr_coverage", "eggnog_coverage", "completeness", "trnas", "rna_16s"]:
            G[c] = pd.to_numeric(G[c], errors="coerce")
        out["genome_catalogue"] = dict(genomes=int(len(G)), mags=int((G.type == "MAG").sum()), isolates=int((G.type == "Isolate").sum()),
                                       ipr_coverage_min=float(G.ipr_coverage.min()), ipr_coverage_p01=q(G.ipr_coverage, 0.01), ipr_coverage_median=q(G.ipr_coverage, 0.5),
                                       eggnog_coverage_min=float(G.eggnog_coverage.min()), eggnog_coverage_p01=q(G.eggnog_coverage, 0.01), eggnog_coverage_median=q(G.eggnog_coverage, 0.5),
                                       completeness_min=float(G.completeness.min()), genomes_with_no_trna=int((G.trnas.fillna(0) == 0).sum()),
                                       domains=G.lineage.fillna("").str.split(";").str[0].value_counts().to_dict())
    pm = pd.read_csv(os.path.join(DATA, "pubmed_method_counts.csv"))
    meth = pm[pm.method_class != "union"]
    un = pm[pm.method_class == "union"].set_index("method").pubmed_records
    seq_sum = int(meth[meth.method_class == "sequence-based"].pubmed_records.sum())
    agn_sum = int(meth[meth.method_class.isin(["activity, agnostic to the genetic polymer", "chirality", "imaging"]) & ~meth.method.str.contains("hybridisation")].pubmed_records.sum())
    seq = int(un["UNION: sequence-based (4 queries)"]); agn = int(un["UNION: polymer-agnostic (7 queries)"])
    agn_strict = int(un["UNION: polymer-agnostic, excluding nucleic-acid stable isotope probing"]); sip_na = int(un["stable isotope probing that names DNA, RNA or sequencing"])
    sip = int(meth[meth.method == "stable isotope probing"].pubmed_records.iloc[0]); expl = int(un["UNION: explicit searches (3 queries)"])
    # union queries count each record once; the summed figures (first pass) double-count overlapping records and are kept only for comparison
    out["pubmed"] = dict(sequence_based_records_union=seq, polymer_agnostic_records_union=agn, polymer_agnostic_records_union_excluding_nucleic_acid_sip=agn_strict,
                         stable_isotope_probing_records=sip, stable_isotope_probing_records_naming_dna_rna_or_sequencing=sip_na, explicit_search_records_union=expl,
                         ratio_sequence_to_agnostic_low=seq / max(agn, 1), ratio_sequence_to_agnostic_high=seq / max(agn_strict, 1),
                         first_pass_summed_sequence_based=seq_sum, first_pass_summed_polymer_agnostic=agn_sum)
    # sequencing depth of the planet: bases sequenced versus DNA in the biosphere (Landenmark et al. 2015: 5.3e31 Mb) is not computed here because
    # the archives' base totals are not exposed by a simple count endpoint; the report quotes sample and genome counts only.
    json.dump(out, open(os.path.join(DATA, "earth_stats.json"), "w"), indent=1)

    print("assemblies used:", len(a), "of", n_all, "| read-level analyses used:", len(m), "of", n_all_m)
    print(pd.DataFrame(rows)[["measure", "biome", "analyses", "studies", "p05", "median", "p95"]].to_string(index=False, float_format=lambda v: f"{v:.3f}"))
    print("assembly unmatched share:", out["assembly_unmatched"]); print("read-level unmatched share:", out["read_unmatched"])
    print("assemblies containing rRNA genes (descriptive only):", f"{out['assemblies_with_rrna_share']:.4f}", "| without:", out["assemblies_without_rrna"],
          "| >80% unmatched:", out["assemblies_over_80pct_unmatched"], "| >95%:", out["assemblies_over_95pct_unmatched"], "| max:", round(out["max_unmatched_share"], 3))
    print(top[["study", "analyses", "median_unmatched", "study_name"]].to_string(index=False, max_colwidth=90))
    if "genome_catalogue" in out:
        print("genome catalogue:", out["genome_catalogue"])
    print("PubMed:", out["pubmed"])


if __name__ == "__main__":
    main()