Draw each reference's own genes under its own samples

CHANGELOG.md

@@ -84,6 +84,21 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Fixed
 
+- **The genome page drew one reference's genes under every sample.** The report took the GFF and the
+  masked regions of the first reference of the run, and its length, for the whole landscape, so on a
+  cohort mapped against two assemblies the second one's samples were drawn over the first one's genes,
+  and its gene table, go-to-gene and masked regions were the first one's. Every reference's GFF, masked
+  regions and FASTA index now reach the report, and the page shows one reference at a time, with its
+  own samples, length, genes and masked regions; a multi-contig reference's genes are placed contig
+  after contig, as its samples' profiles are binned.
+
+  Reading those GFFs also showed two faults of the report's GFF reader. A GFF built from a table
+  writes pandas' NaN as `Name=nan`, which it read as a name, and it merged genes by name: on one of the
+  cohort's references 3,001 of 3,989 genes became a single gene called `nan`, so the landscape knew 988
+  genes. A RefSeq GFF names every CDS after its protein, and each of H37Rv's genes came back twice,
+  once as `dnaA` and once as `NP_214515.1`. A placeholder is now read as no name, a CDS is its parent
+  gene's, and genes that share a name at different places stay two.
+
 - **The H37Rv coordinate and numbering were wrong, or missing, on any reference but H37Rv.** On the
   185-sample cohort, mapped against two assemblies, the report gave 15,653 of the second reference's
   16,296 variant sites (97%) a wrong H37Rv coordinate. Four faults:

bin/qc_report.py

@@ -27,6 +27,7 @@
 from datetime import datetime, timezone
 
 from qcreport.coverage import build_coverage, parse_depth_profiles, snp_per_kb, write_deletions
+from qcreport.genomes import build_genomes, pairs
 from qcreport.metrics import (ANC_DEF, DEF, DEFS, DIST, METRICS, build_gene_map, discover_extra_metrics,
                               flag_sample, het_frac, is_ancient, lineage_counts_parsed, lineage_fracs, robust)
 from qcreport.panels import build_dynamics, build_epistasis, build_snp_matrix
@@ -63,6 +64,14 @@ def build_parser():
                     help="TSV 'lineage<TAB>#hex' (e.g. mycolorsTB) to colour lineages with a canonical palette (optional).")
     ap.add_argument("--mask-bed", default=None,
                     help="BED of regions to mask (e.g. mtbc_mask.bed: PE/PPE, IS, DR, repeats) -> a 'mask regions' toggle.")
+    ap.add_argument("--ref-gff", nargs="*", default=[],
+                    help="REF=GFF3 for each reference of the run -> its own genes under its own samples (genome "
+                         "landscape, gene hotspots, go-to-gene). With these, --gff and --mask-bed are not used.")
+    ap.add_argument("--ref-mask", nargs="*", default=[],
+                    help="REF=BED for each reference: its repeat/exclude regions (the masked-regions toggle).")
+    ap.add_argument("--ref-fai", nargs="*", default=[],
+                    help="REF=FASTA index for each reference: its contigs' order and lengths, the axis its samples' "
+                         "profiles are binned on, where genes and masked regions are placed.")
     ap.add_argument("--gene-burden", default=None,
                     help="Cohort gene-burden TSV (collect_summary --gene-burden-out) -> the Functional gene burden panel (optional).")
     ap.add_argument("--dr-report", default=None,
@@ -370,8 +379,32 @@ def profile_of(sid):
 
     lineages = sorted({s["lineage"] for s in jsamples if s["lineage"]})
     n_ancient = sum(1 for s in jsamples if s["anc"])
-    mask_iv = parse_bed(args.mask_bed)
-    mask_bins, mask_pct = mask_profile(mask_iv, genome_len)
+    # Each reference's own genome (length, genes, masked regions), placed as its samples' profiles are
+    # binned. The top-level fields keep the reference most samples were mapped to, which is also what a
+    # single-reference run has always had there.
+    ref_gff, ref_mask, ref_fai = pairs(args.ref_gff), pairs(args.ref_mask), pairs(args.ref_fai)
+    genomes = {}
+    if ref_gff or ref_mask or ref_fai:
+        lengths = {}
+        for sid, m in summ.items():
+            v = to_float(m.get("length")) or to_float(m.get("genome_len"))
+            if v and ref_of.get(sid):
+                lengths[ref_of[sid]] = max(lengths.get(ref_of[sid], 0), int(v))
+        genomes = build_genomes(ref_gff, ref_mask, ref_fai, lengths, NBINS)
+    if genomes:
+        n_on = {}
+        for s in jsamples:
+            if s.get("ref"):
+                n_on[s["ref"]] = n_on.get(s["ref"], 0) + 1
+        main_ref = sorted(genomes, key=lambda r: (-n_on.get(r, 0), r))[0]
+        g0 = genomes[main_ref]
+        genome_len = g0["len"] or genome_len
+        genes, mask_bins, mask_pct, mask_iv = g0["genes"], g0["mask_bins"], g0["mask_pct"], g0["mask_iv"]
+        gene_map = build_gene_map(list(ref_gff.values()))
+    else:
+        genes, gene_map = parse_gff(args.gff), build_gene_map(args.gff)
+        mask_iv = parse_bed(args.mask_bed)
+        mask_bins, mask_pct = mask_profile(mask_iv, genome_len)
     provenance = {}
     for kv in (args.provenance or []):
         if "=" in kv:
@@ -407,8 +440,9 @@ def profile_of(sid):
     payload = {"generated": now, "version": clean_str(args.version) or "", "repo_url": REPO_URL,
                "counts": counts, "thresholds": thr, "dist": dist_keys,
                "genome_len": genome_len, "snp_density_ok": snp_density_ok, "defs": DEFS, "nbins": NBINS,
-               "genes": parse_gff(args.gff), "lin_colors": parse_lineage_colors(args.lineage_colors),
-               "gene_map": build_gene_map(args.gff), "aa2_label": args.aa2_label,
+               "genes": genes, "lin_colors": parse_lineage_colors(args.lineage_colors),
+               "gene_map": gene_map, "aa2_label": args.aa2_label,
+               "genomes": genomes or None,
                "ref_name": provenance.get('reference', ''),   # mapping reference name for the SNP tables' headers
                "section_info": SECTION_INFO,
                "gene_burden": parse_gene_burden(args.gene_burden) or None,

bin/qcreport/genomes.py

@@ -0,0 +1,116 @@
+"""Each reference's own genome for the report: its length, genes and masked regions, on the axis its samples'
+profiles are binned on.
+
+A cohort can be mapped against several references. Every sample's landscape profile is binned over its own
+reference, contigs end to end in the order of the reference's FASTA index (stats_to_legacy.py), so the genes
+and masked regions drawn under it have to be that reference's and placed the same way: a gene of the second
+contig sits after the whole of the first one. One genome per reference, each shaped like the single one the
+report has always carried (len, genes, mask_bins, mask_pct, mask_iv), so the page can switch between them.
+"""
+from __future__ import annotations
+
+import gzip
+import os
+
+from .parsers import NBINS, gff_features, mask_profile
+
+MASK_IV_CAP = 5000
+
+
+def _open(path):
+    return gzip.open(path, "rt", encoding="utf-8", errors="replace") if str(path).endswith(".gz") \
+        else open(path, encoding="utf-8", errors="replace")
+
+
+def _usable(path):
+    return bool(path) and os.path.exists(path) and not os.path.basename(str(path)).startswith("NO_FILE")
+
+
+def read_fai(path):
+    """[(contig, length)] of a FASTA index, in its order; [] without one."""
+    out = []
+    if not _usable(path):
+        return out
+    try:
+        with open(path, encoding="utf-8", errors="replace") as fh:
+            for line in fh:
+                c = line.rstrip("\n").split("\t")
+                if len(c) >= 2 and c[1].isdigit():
+                    out.append((c[0], int(c[1])))
+    except OSError:
+        return []
+    return out
+
+
+def bed_intervals(path):
+    """[(contig, start, end)] of a BED or of a reference's exclusion list; [] without one."""
+    out = []
+    if not _usable(path):
+        return out
+    try:
+        with _open(path) as fh:
+            for line in fh:
+                if not line.strip() or line.startswith(("#", "track", "browser")):
+                    continue
+                c = line.rstrip("\n").split("\t")
+                if len(c) >= 3 and c[1].strip().isdigit() and c[2].strip().isdigit():
+                    out.append((c[0], int(c[1]), int(c[2])))
+    except OSError:
+        return []
+    return out
+
+
+def _placer(contigs):
+    """A function placing (contig, position) on the genome-wide axis, or None where it cannot. With a single
+    contig every name is it: a GFF may call the chromosome otherwise than the FASTA does."""
+    offsets, run = {}, 0
+    for name, length in contigs:
+        offsets[name] = run
+        run += length
+    single = len(contigs) <= 1
+
+    def place(contig, pos):
+        if contig in offsets:
+            return offsets[contig] + pos
+        return pos if single else None
+    return place
+
+
+def build_genome(gff=None, mask=None, fai=None, length=0, nbins=NBINS):
+    """{len, genes, mask_bins, mask_pct, mask_iv} of one reference, positions on its genome-wide axis."""
+    contigs = read_fai(fai)
+    place = _placer(contigs)
+    genes = []
+    for contig, name, start, end in gff_features(gff):
+        s, e = place(contig, start), place(contig, end)
+        if s is not None and e is not None:
+            genes.append({"name": name, "start": s, "end": e})
+    genes.sort(key=lambda g: (g["start"], g["end"]))
+    glen = sum(n for _, n in contigs) or length or max((g["end"] for g in genes), default=0)
+    iv = []
+    for contig, start, end in bed_intervals(mask):
+        s, e = place(contig, start), place(contig, end)
+        if s is not None and e is not None:
+            iv.append((s, e))
+    bins, pct = mask_profile(iv, glen, nbins)
+    return {"len": glen, "genes": genes, "mask_bins": bins, "mask_pct": pct,
+            "mask_iv": (iv[:MASK_IV_CAP] if iv else None)}
+
+
+def build_genomes(ref_gff, ref_mask, ref_fai, lengths=None, nbins=NBINS):
+    """{reference: genome} for every reference named in any of the three {reference: path} maps. `lengths`
+    ({reference: bp}, from the samples' own summary) stands in for a reference without a FASTA index."""
+    lengths = lengths or {}
+    refs = sorted(set(ref_gff) | set(ref_mask) | set(ref_fai))
+    return {r: build_genome(ref_gff.get(r), ref_mask.get(r), ref_fai.get(r), lengths.get(r, 0), nbins)
+            for r in refs}
+
+
+def pairs(items):
+    """{reference: path} from REF=PATH arguments."""
+    out = {}
+    for it in items or []:
+        ref, sep, path = str(it).partition("=")
+        if sep and ref and path:
+            out[ref] = path
+    return out

bin/qcreport/metrics.py

@@ -11,6 +11,7 @@
 """
 from __future__ import annotations
 
+import os
 import re
 
 from .parsers import clean_str, parse_gene_locus, to_float
@@ -155,14 +156,17 @@
 _RV_LOCUS_RE = re.compile(r'^(Rv\d|MTB\d)', re.I)   # H37Rv (Mycobrowser) locus-tag schemes
 
 
-def build_gene_map(gff_path):
+def build_gene_map(gff_paths):
     """{gene: Mycobrowser (H37Rv) locus tag}. Curated resistance-gene map, overlaid with any H37Rv-style
-    locus tags from the run's GFF (so an H37Rv reference contributes every gene; a non-H37Rv reference's
-    strain-specific tags are ignored in favour of the curated H37Rv equivalents)."""
+    locus tags from the run's GFF, or from each of its references' GFFs (so an H37Rv reference contributes
+    every gene; a non-H37Rv reference's strain-specific tags are ignored in favour of the curated H37Rv
+    equivalents)."""
     m = dict(GENE_RV)
-    for gene, locus in parse_gene_locus(gff_path).items():
-        if _RV_LOCUS_RE.match(locus):
-            m[gene] = locus
+    paths = [gff_paths] if (gff_paths is None or isinstance(gff_paths, (str, os.PathLike))) else gff_paths
+    for path in paths:
+        for gene, locus in parse_gene_locus(path).items():
+            if _RV_LOCUS_RE.match(locus):
+                m[gene] = locus
     return m
 
 

bin/qcreport/parsers.py

@@ -611,23 +611,38 @@ def taxon(k, _share=share):
     return {"samples": out, "target": tgt}
 
 
+# Values that stand for "nothing here", written into a GFF as though they were data: a GFF built from a
+# table writes pandas' NaN as `nan`, and one real reference carries `Name=nan` on 3,001 of its 3,989 CDS.
+# Read as a name, every one of them was one gene called nan, so the landscape drew 988 genes of 3,989 and
+# the rest could not be found. The gene-conversion annotation already reads them as absent.
+_GFF_ABSENT = {"", ".", "-", "na", "nan", "none", "null", "n/a", "unknown"}
+
+
 def _gff_attr(attrs, key):
-    """Value of GFF3 attribute `key` from the column-9 string, or None when absent. The key is matched
-    WHOLE against each ';'-separated field (a substring search would read 'locus_tag' out of
-    'old_locus_tag' and 'gene' out of 'pseudogene', both legal GFF3 and routine in RefSeq/Prokka).
-    Only the first '=' splits, so a value that itself contains '=' survives."""
+    """Value of GFF3 attribute `key` from the column-9 string, or None when absent or a placeholder for
+    absent (nan, NA, '.'...), so a cascade of keys falls through to the next. The key is matched WHOLE
+    against each ';'-separated field (a substring search would read 'locus_tag' out of 'old_locus_tag' and
+    'gene' out of 'pseudogene', both legal GFF3 and routine in RefSeq/Prokka). Only the first '=' splits,
+    so a value that itself contains '=' survives."""
     for field in attrs.split(";"):
         k, sep, v = field.partition("=")
         if sep and k.strip() == key:
-            return v.strip()
+            v = v.strip()
+            return None if v.lower() in _GFF_ABSENT else v
     return None
 
 
-def parse_gff(path):
-    """Best-effort GFF3 parse of gene/CDS features -> [{name, start, end}] (1-based). Empty on any problem."""
+def gff_features(path):
+    """[(contig, name, start, end)] of the genes of a GFF3, 1-based inclusive; [] on any problem.
+
+    One entry per gene. A CDS stands in only where the file gives no gene for it: an annotator that writes
+    only CDS (Prokka). A CDS whose Parent is a gene of the file is that gene's, as in every RefSeq GFF, where
+    it would otherwise come back as a second gene named after its protein (NP_214515.1 beside dnaA); a gene
+    and a CDS of the same name that overlap are one gene too, the gene's span kept, while namesakes apart
+    are two genes. A name written as a placeholder (nan, NA) falls through to the locus tag or the ID."""
     if not path or not os.path.exists(path):
         return []
-    genes = {}
+    feats, gene_ids = {}, set()
     try:
         op = gzip.open if str(path).endswith(".gz") else open
         with op(path, "rt", encoding="utf-8", errors="replace") as fh:
@@ -647,12 +662,24 @@ def parse_gff(path):
                     if name is not None:
                         break
                 name = name or f"{start}-{end}"
-                # prefer a 'gene' feature over a 'CDS' of the same name
-                if name not in genes or (c[2] == "gene" and genes[name].get("type") != "gene"):
-                    genes[name] = {"name": name, "start": start, "end": end, "type": c[2]}
+                if c[2] == "gene" and _gff_attr(c[8], "ID"):
+                    gene_ids.add(_gff_attr(c[8], "ID"))
+                same = feats.setdefault((c[0], name), [])
+                hit = next((i for i, f in enumerate(same) if f[2] <= end and start <= f[3]), None)
+                if hit is None:
+                    same.append((c[0], name, start, end, c[2], _gff_attr(c[8], "Parent") or ""))
+                elif c[2] == "gene" and same[hit][4] != "gene":
+                    same[hit] = (c[0], name, start, end, "gene", "")
     except (OSError, ValueError):
         return []
-    return sorted(({"name": g["name"], "start": g["start"], "end": g["end"]} for g in genes.values()),
+    return [f[:4] for same in feats.values() for f in same
+            if not (f[4] == "CDS" and any(p in gene_ids for p in f[5].split(",") if p))]   # its gene stands for it
+
+
+def parse_gff(path):
+    """Best-effort GFF3 parse of genes -> [{name, start, end}] (1-based), in position order; see
+    gff_features for what counts as one gene. Empty on any problem."""
+    return sorted(({"name": n, "start": s, "end": e} for _, n, s, e in gff_features(path)),
                   key=lambda g: g["start"])
 
 

bin/report_assets/js/01_prelude.js

@@ -49,6 +49,6 @@ R.metrics.push({key:'het_frac',label:'Het %',kind:'pct',dir:'hi_bad'});
 R.metrics.push({key:'n_lineages',label:'# lin',kind:'int',dir:'hi_bad'});
 if(DIST.indexOf('het_frac')<0)DIST.push('het_frac');
 // ---- SNP density per callable Mb (guarded: only for a plausible genome length) ----
-if(R.snp_density_ok){var GMB=R.genome_len/1e6;R.samples.forEach(function(s){var snp=s.m.snps,cal=s.m.callable_pct;s.m.snp_density=(snp!=null&&cal!=null&&cal>0)?(snp/(cal/100*GMB)):null;});R.metrics.push({key:'snp_density',label:'SNP/Mb',kind:'float',dir:'neu'});if(DIST.indexOf('snp_density')<0)DIST.push('snp_density');}
+if(R.snp_density_ok){R.samples.forEach(function(s){var GMB=((R.genomes&&s.ref&&R.genomes[s.ref]&&R.genomes[s.ref].len)||R.genome_len)/1e6;var snp=s.m.snps,cal=s.m.callable_pct;s.m.snp_density=(snp!=null&&cal!=null&&cal>0)?(snp/(cal/100*GMB)):null;});R.metrics.push({key:'snp_density',label:'SNP/Mb',kind:'float',dir:'neu'});if(DIST.indexOf('snp_density')<0)DIST.push('snp_density');}
 // ===================== SIGNATURE ANALYSIS LAYER (QC-space PCA + robust Mahalanobis + divergence/completeness + temporal) =====================
 R.defs=R.defs||{};