Read each sample against the cohort: deletions, relatedness, series and minority variants

CHANGELOG.md

@@ -8,6 +8,57 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Added
 
+- **How low an allele frequency can be trusted.** A *Minority variants* panel counts each
+  sample's calls below fixation, their spread of allele fractions and how many rest on three
+  alternate reads or fewer, where an error and a real minority look the same. When the
+  samplesheet names each sample's DNA extract (a column such as `dna_id`, `extract` or
+  `biosample`), libraries of the same DNA are compared: a real minority is in the DNA and another
+  library calls it too, an error is not reproduced. Pairs that share a FASTQ file (a merged sample
+  and its runs) or were mapped against different references are left out, and a call only counts,
+  reproduced or not, where the other library could have called it: above `--consensus_min_dp`
+  reads, with at least five alternate reads expected at the call's fraction.
+
+- **What each series gained since its first time point.** For every samplesheet group followed
+  over time (a patient, a passage line), the report lists the SNPs each later sample carries fixed
+  that the group's first time point did not: *new* where that first time point was read at the
+  site without the allele, *risen* where it held the allele as a minority, *unknown* where it was
+  not read there, and *lost* for those fixed at the start and read without the allele later. The
+  depth comes from the SNP matrix, which the report now reads, so a site the start never read is
+  not counted as new. The first time point is the earliest with a sample the QC keeps and places
+  in the series, and a time that is a date is read as one (`2021-03-01`, `15/01/2020`), with
+  `baseline` or `pre` first. A chart follows each series through the median of each time point,
+  coloured by treatment, and leaves out the samples the QC fails or places outside their series.
+  On the *Resistance* page, each series' grade 1-2 mutations are laid out per time point with
+  those acquired since the start set apart: on the cohort that motivated it, the three bedaquiline
+  lines of one lineage acquire atpE I66M at passages 18 and 19, two of them with E61D.
+
+- **How close the samples are to each other.** A new step, `SNP_DISTANCES`, counts the SNPs
+  between every two consensus sequences over the positions both called, so a no-call, a masked
+  position or a mixed site is never a difference: the distances a tree built on the same
+  sequences sees. They are written as `<samplesheet>_snp_distances.tsv` and, one row per pair
+  with how many positions each comparison rests on, `_pairs.tsv`. The report gains a
+  *Relatedness* page: a heatmap in the order that keeps every single-linkage cluster together,
+  the clusters at a threshold that can be moved live (`--snp_cluster_threshold`, 12 by default,
+  the usual *M. tuberculosis* cut), and the samples far from the rest of their samplesheet group.
+  Those are flagged `GROUP_MISMATCH` (WARN) when their group is otherwise tight: on the cohort
+  that motivated it, with distances approximated from its SNP matrix, 18 of the 20 samples it
+  flags are the ones typed as another lineage than their line. A sample mapped against two
+  references is compared on each, as `sample@reference`.
+
+- **Deletions, and SNPs per callable kb along the genome.** Each sample's all-positions VCF is
+  reduced to its depth per window, the stretches no read covers and a per-gene table
+  (`DEPTH_PROFILE`, published under `stats/`). The report compares those stretches across the
+  samples of each reference: one a sample lacks while the others read it is a deletion, private or
+  shared; one nearly every sample lacks (90% or more) is a repeat or a part of the reference none
+  of these genomes has, and is listed apart instead of being counted as anyone's. Two samples'
+  stretches are one deletion when each covers at least half of the other. A *Deletions* panel
+  lists them with the genes they remove and zooms the landscape onto each, and
+  `<samplesheet>_deletions.tsv` keeps them with every carrier's own coordinates. The landscape
+  gains a *Deletions* track and a *SNPs / kb* track that divides each bin's SNPs by the positions
+  deep enough for the consensus to call a base (`--allpos_min_cov` reads), so a bin half of which
+  was not read no longer looks half as variable. Samples under `--report_depth_min` are not
+  assessed: stretches without reads turn up there by chance.
+
 - **Two guards against a sample mapped to a genome it does not belong to.** In a 185-sample
   cohort, 22 samples annotated as one lineage were really another, were routed to that lineage's
   reference, carried about 2,000 SNPs where their line mates carried 5, and produced 507 of the

assets/NO_FILE.README.md

@@ -11,6 +11,8 @@ stages one of these empty files instead, and the process recognises it by the
 | `NO_FILE_DR` | drug-resistance calls (`QC_REPORT`) | `--run_pathotypr false` |
 | `NO_FILE_H37RV` | canonically annotated VCFs (`QC_REPORT`) | `--annotate_canonical false` |
 | `NO_FILE_LIFTOVER` | variant coordinate map (`QC_REPORT`) | `--variant_liftover false` |
+| `NO_FILE_DISTANCES` | pairwise SNP distances (`QC_REPORT`) | `--make_snp_distances false`, or no consensus sequences |
+| `NO_FILE_MATRIX` | the master SNP matrix (`QC_REPORT`) | `--make_snp_matrix false` |
 
 Two properties matter, and both are load-bearing:
 

bin/depth_profile.py

@@ -0,0 +1,256 @@
+#!/usr/bin/env python3
+"""What one sample's reads cover, from its all-positions VCF.
+
+The all-positions VCF has a record for every reference position, so it holds the depth of the
+whole genome. This reduces it to three small tables for the cohort report:
+
+  <prefix>.depth_windows.tsv  one row per window: mean depth, the share of positions no read
+                              covers, and the share deep enough for the consensus to call a
+                              base (depth of --callable-dp or more)
+  <prefix>.zero_depth.tsv     every stretch of at least --min-run positions no read covers
+  <prefix>.gene_depth.tsv     per gene of the GFF: the share of its positions read at all, the
+                              share the consensus can call, its mean depth, and that depth
+                              against the sample's genome-wide median (with --gff)
+
+A stretch no read covers is a deletion only if the rest of the cohort reads it; otherwise it is
+a part of the reference none of these genomes has, or one no read can be placed on. That
+comparison needs the cohort, so it is the report's to make, not this script's.
+
+Each table opens with '#' lines naming the sample, the reference and the sample's genome-wide
+depth, so a table can be read without its file name.
+"""
+from __future__ import annotations
+
+import argparse
+import gzip
+import os
+import sys
+
+import numpy as np
+
+from qcreport.parsers import _gff_attr
+
+
+def _open(p):
+    return gzip.open(p, 'rt', encoding='utf-8', errors='replace') if str(p).endswith('.gz') \
+        else open(p, encoding='utf-8', errors='replace')
+
+
+def read_depths(path):
+    """(sample, {contig: depth array}, {contig: read-at-all mask}) from an all-positions VCF.
+
+    Arrays are sized from the ##contig lengths and grown when a record lies beyond them. A
+    position without a record is not a position of depth 0, so the mask keeps the two apart.
+    """
+    lengths, order, depth, known = {}, [], {}, {}
+    sample = None
+    dp_index = {}                  # FORMAT string -> index of DP in it (-1 when it has none)
+
+    def arrays(contig, need):
+        if contig not in depth:
+            n = max(lengths.get(contig, 0), need)
+            depth[contig] = np.zeros(n, dtype=np.int32)
+            known[contig] = np.zeros(n, dtype=bool)
+            order.append(contig)
+        elif need > len(depth[contig]):
+            grow = max(need, 2 * len(depth[contig]))
+            depth[contig] = np.concatenate([depth[contig], np.zeros(grow - len(depth[contig]), np.int32)])
+            known[contig] = np.concatenate([known[contig], np.zeros(grow - len(known[contig]), bool)])
+        return depth[contig], known[contig]
+
+    last = (None, None, None)      # (contig, depth array, mask): most lines stay on one contig
+    with _open(path) as fh:
+        for line in fh:
+            if line[0] == '#':
+                if line.startswith('##contig='):
+                    body = line.strip()[len('##contig=<'):-1]
+                    fields = dict(f.split('=', 1) for f in body.split(',') if '=' in f)
+                    try:
+                        lengths[fields['ID']] = int(fields.get('length', 0))
+                    except (KeyError, ValueError):
+                        pass
+                elif line.startswith('#CHROM'):
+                    cols = line.rstrip('\n').split('\t')
+                    sample = cols[9] if len(cols) > 9 else None
+                continue
+            c = line.rstrip('\n').split('\t')
+            if len(c) < 10:
+                continue
+            k = dp_index.get(c[8])
+            if k is None:
+                keys = c[8].split(':')
+                k = dp_index[c[8]] = keys.index('DP') if 'DP' in keys else -1
+            try:
+                pos = int(c[1])
+            except ValueError:
+                continue
+            vals = c[9].split(':')
+            raw = vals[k] if 0 <= k < len(vals) else '.'
+            try:
+                dp = int(raw)
+            except ValueError:
+                continue                   # a depth the record does not state is not a depth of 0
+            contig = c[0]
+            if last[0] != contig or pos > len(last[1]):
+                d, m = arrays(contig, pos)
+                last = (contig, d, m)
+            last[1][pos - 1] = dp
+            last[2][pos - 1] = True
+    # Trim what growing over-allocated past the last record of a contig with no stated length.
+    for contig in order:
+        if not lengths.get(contig):
+            idx = np.flatnonzero(known[contig])
+            n = int(idx[-1]) + 1 if idx.size else 0
+            depth[contig], known[contig] = depth[contig][:n], known[contig][:n]
+    return sample, {c: depth[c] for c in order}, {c: known[c] for c in order}
+
+
+def windows(depth, known, size, callable_dp):
+    """[(contig, start, end, mean depth, zero fraction, callable fraction)], 1-based inclusive.
+
+    Fractions are of the positions with a record; a window with none has no values.
+    """
+    out = []
+    for contig, d in depth.items():
+        m = known[contig]
+        for s in range(0, len(d), size):
+            e = min(len(d), s + size)
+            dm, km = d[s:e], m[s:e]
+            n = int(km.sum())
+            if n == 0:
+                out.append((contig, s + 1, e, None, None, None))
+                continue
+            dk = dm[km]
+            out.append((contig, s + 1, e, float(dk.mean()), float((dk == 0).mean()),
+                        float((dk >= callable_dp).mean())))
+    return out
+
+
+def zero_runs(depth, known, min_run):
+    """[(contig, start, end)] of the stretches of at least `min_run` read positions at depth 0."""
+    out = []
+    for contig, d in depth.items():
+        z = ((d == 0) & known[contig]).astype(np.int8)
+        edges = np.diff(np.concatenate([[0], z, [0]]))
+        starts, ends = np.flatnonzero(edges == 1), np.flatnonzero(edges == -1)
+        for s, e in zip(starts, ends):
+            if e - s >= min_run:
+                out.append((contig, int(s) + 1, int(e)))
+    return out
+
+
+def read_genes(path):
+    """[(contig, start, end, name, locus_tag)] of the gene features of a GFF3, CDS where there
+    are no genes (a GFF from an annotator that writes only CDS). Empty without a readable file."""
+    if not path or not os.path.exists(path) or os.path.basename(path).startswith('NO_FILE'):
+        return []
+    genes, cds = [], []
+    try:
+        with _open(path) as fh:
+            for line in fh:
+                if line.startswith('#') or '\t' not in line:
+                    continue
+                c = line.rstrip('\n').split('\t')
+                if len(c) < 9 or c[2] not in ('gene', 'CDS'):
+                    continue
+                try:
+                    start, end = int(c[3]), int(c[4])
+                except ValueError:
+                    continue
+                name = None
+                for key in ('Name', 'gene', 'locus_tag', 'ID'):
+                    name = _gff_attr(c[8], key)
+                    if name:
+                        break
+                row = (c[0], start, end, name or f'{start}-{end}', _gff_attr(c[8], 'locus_tag') or '')
+                (genes if c[2] == 'gene' else cds).append(row)
+    except OSError:
+        return []
+    return sorted(genes or cds, key=lambda g: (g[0], g[1]))
+
+
+def gene_depths(depth, known, genes, callable_dp, median):
+    """[(contig, start, end, name, locus_tag, length, breadth, callable, mean depth, relative)].
+
+    breadth is the share of the gene's read positions at depth 1 or more, callable the share at
+    `callable_dp` or more, relative the mean depth against the genome-wide median (None at
+    median 0).
+    """
+    out = []
+    for contig, start, end, name, locus in genes:
+        d = depth.get(contig)
+        if d is None:
+            continue
+        s, e = max(0, start - 1), min(len(d), end)
+        if e <= s:
+            continue
+        km = known[contig][s:e]
+        dk = d[s:e][km]
+        if dk.size == 0:
+            out.append((contig, start, end, name, locus, end - start + 1, None, None, None, None))
+            continue
+        mean = float(dk.mean())
+        out.append((contig, start, end, name, locus, end - start + 1, float((dk > 0).mean()),
+                    float((dk >= callable_dp).mean()), mean, (mean / median) if median else None))
+    return out
+
+
+def _fmt(v, nd=4):
+    return '' if v is None else (f'{v:.{nd}f}' if isinstance(v, float) else str(v))
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument('--vcf', required=True, help="the sample's all-positions VCF")
+    ap.add_argument('--gff', default=None, help="the reference's GFF3, for the per-gene table")
+    ap.add_argument('--reference', default='', help="the reference id, written into each table")
+    ap.add_argument('--window', type=int, default=1000, help="window size in bp")
+    # The consensus writes a reference base only where the backbone calls a confident reference,
+    # which it does from --allpos_min_cov reads (30); between --consensus_min_dp and that it
+    # writes N. A position read at 20x is not callable, and counting it so would make every
+    # SNP-per-callable-kb density read low.
+    ap.add_argument('--callable-dp', type=int, default=30,
+                    help="depth from which the consensus can call a base (--allpos_min_cov)")
+    ap.add_argument('--min-run', type=int, default=50,
+                    help="shortest stretch without reads that is written out")
+    ap.add_argument('--out-prefix', required=True)
+    a = ap.parse_args(argv)
+
+    sample, depth, known = read_depths(a.vcf)
+    if sample is None:
+        sys.stderr.write(f"[depth_profile] {a.vcf} names no sample\n")
+        return 1
+    read = np.concatenate([depth[c][known[c]] for c in depth]) if depth else np.zeros(0, np.int32)
+    median = float(np.median(read)) if read.size else 0.0
+    mean = float(read.mean()) if read.size else 0.0
+    head = [f"# sample={sample}", f"# reference={a.reference}", f"# genome_len={int(read.size)}",
+            f"# median_dp={median:g}", f"# mean_dp={mean:.2f}",
+            f"# zero_frac={float((read == 0).mean()) if read.size else 0:.4f}",
+            f"# callable_frac={float((read >= a.callable_dp).mean()) if read.size else 0:.4f}",
+            f"# callable_dp={a.callable_dp}", f"# window={a.window}", f"# min_run={a.min_run}",
+            "# contigs=" + ",".join(f"{c}:{len(depth[c])}" for c in depth)]
+
+    with open(f"{a.out_prefix}.depth_windows.tsv", 'w') as fh:
+        fh.write("\n".join(head) + "\ncontig\tstart\tend\tmean_dp\tzero_frac\tcallable_frac\n")
+        for row in windows(depth, known, a.window, a.callable_dp):
+            fh.write("\t".join(_fmt(v, 2 if i == 3 else 4) for i, v in enumerate(row)) + "\n")
+    runs = zero_runs(depth, known, a.min_run)
+    with open(f"{a.out_prefix}.zero_depth.tsv", 'w') as fh:
+        fh.write("\n".join(head) + "\ncontig\tstart\tend\tlength\n")
+        for contig, s, e in runs:
+            fh.write(f"{contig}\t{s}\t{e}\t{e - s + 1}\n")
+    genes = read_genes(a.gff)
+    if genes:
+        with open(f"{a.out_prefix}.gene_depth.tsv", 'w') as fh:
+            fh.write("\n".join(head) + "\ncontig\tstart\tend\tgene\tlocus_tag\tlength\tbreadth\t"
+                                       "callable\tmean_dp\trel_depth\n")
+            for row in gene_depths(depth, known, genes, a.callable_dp, median):
+                fh.write("\t".join(_fmt(v, 2 if i == 8 else 4) for i, v in enumerate(row)) + "\n")
+    sys.stderr.write(f"[depth_profile] {sample}: {int(read.size)} positions, median depth {median:g}, "
+                     f"{len(runs)} stretch(es) of {a.min_run}+ bp without reads"
+                     + (f", {len(genes)} gene(s)" if genes else "") + "\n")
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())