feat(sorting-hat): per-call status entries and repeat-call calibration

analysis/calibration_truth/repeat_call_truth/accessions.tsv

@@ -0,0 +1,156 @@
+accession
+E9P9G2
+P32323
+P32781
+P32768
+P36170
+P08640
+P38894
+P39712
+P20840
+Q12140
+P40442
+Q05164
+Q6B2U8
+Q12127
+P47001
+P53301
+P29029
+P28319
+P43497
+P38248
+P23776
+P38993
+Q04433
+P22146
+P41809
+P32478
+P39005
+P36027
+P32334
+Q03178
+P54867
+P53616
+P32623
+P53832
+P32329
+P47178
+P47179
+P42835
+Q01589
+P10863
+P40552
+Q12218
+Q04739
+P22943
+Q03125
+P38082
+P40505
+Q12507
+Q5A8T4
+Q59L12
+A0A1D8PQB9
+Q5A8T7
+Q5A2Z7
+Q5A312
+A0A1D8PQ86
+G1UBC2
+P46593
+Q59PF9
+Q5AL03
+Q5AAL9
+P0CU38
+A0A1D8PIY8
+Q59XA7
+Q5A849
+Q5A7R7
+Q5A029
+Q5A1E0
+Q59XL0
+Q59XB0
+Q5A6U1
+Q5ACL7
+Q5A5M7
+Q59TP1
+Q59RR0
+Q59WH0
+Q5AMT2
+Q59WG7
+Q5AA40
+A0A1D8PMH9
+Q5AFI4
+Q59Y20
+A0A1D8PCY4
+Q5AIR7
+Q59RW5
+Q59XX2
+Q5AAN7
+P43076
+P87020
+Q59UT4
+Q5A4X3
+P0CY27
+Q5A651
+P0DJ06
+P0CY29
+Q59SU1
+Q59NP5
+Q5AJC0
+P29717
+Q59Y31
+O94072
+Q59ZB1
+A0A1D8PP43
+Q59Z29
+Q5A4F3
+P83774
+Q59U10
+P53705
+P39827
+Q9Y7W4
+Q59QH2
+A0A1D8PK00
+P53698
+A0A1D8PR83
+G1UB67
+Q59X67
+A0A1D8PQU2
+Q59XU9
+P82612
+A0A1D8PE35
+Q5A6N7
+A0A1D8PN26
+A0A1D8PKY7
+Q59SR6
+P82610
+Q00310
+P46592
+A0A1D8PD52
+A0A1D8PE87
+A0A1D8PRI0
+A0A1D8PCV9
+Q59UQ8
+O74189
+A0A1D8PIK2
+Q59XU5
+Q5AMQ6
+Q5A7M9
+Q5A287
+Q5AK51
+Q5ANF0
+Q5AQ36
+O59923
+A0A1D8PE53
+A0A1D8PJZ6
+A0A1D8PT45
+Q59Q34
+A0A1D8PTB4
+Q59UR3
+Q5ADM7
+Q5ANJ4
+A0A1D8PHU1
+P0CY34
+P13649
+Q5AH00
+Q5AP80
+Q5A2J7

analysis/calibration_truth/repeat_call_truth/build_truth.py

@@ -0,0 +1,156 @@
+#!/usr/bin/env python3
+"""Build per-species truth tables for the call `tandem_repeat_protein` ("has a tandem repeat region").
+
+Input: repeat_truth_candidates.tsv (label 1, 0 or empty per curated protein; see make_candidates.py)
+and the proteome FASTA of each run (to cluster the truth proteins at 30% identity, coverage 0.5, as
+analysis/calibration_truth/c1_truth_count.py does). Output per species: truth.<species>.tsv with the
+columns id, label, cluster (what `cellsurface_sorting_hat_calibrate truth` reads) and
+truth.<species>.annotated.tsv with the basis of each label.
+
+Labels: 1 = paper statement or at least two UniProt Repeat features. 0 = no UniProt Repeat feature
+(absence of annotation, so an ASSUMED negative: UniProt lacks features for some repeat proteins).
+Proteins with one feature, a family-inference-only claim, or no record are left out. A protein with
+several curated rows gets one row; evidence of a repeat beats an assumed negative.
+
+Usage: build_truth.py --candidates FILE --fasta Scer_S288C=PATH --fasta Calb_SC5314=PATH --out DIR
+Needs `mmseqs` on PATH (module load MMseqs2/17-b804f).
+"""
+
+import argparse
+import csv
+import subprocess
+import sys
+import tempfile
+from pathlib import Path
+
+
+def read_fasta(path):
+    seqs, name, buf = {}, None, []
+    with open(path) as f:
+        for line in f:
+            line = line.rstrip()
+            if line.startswith(">"):
+                if name:
+                    seqs[name] = "".join(buf).rstrip("*")
+                name, buf = line[1:].split()[0], []
+            else:
+                buf.append(line)
+    if name:
+        seqs[name] = "".join(buf).rstrip("*")
+    return seqs
+
+
+def cluster(seqs, threads=4):
+    """``{id: representative id}``. MMseqs2 rewrites UniProt-style IDs (sp|ACC|NAME) in its output, so
+    the sequences go in under neutral index IDs and the result is mapped back."""
+    names = list(seqs)
+    with tempfile.TemporaryDirectory() as tmp:
+        fa = Path(tmp) / "in.faa"
+        with open(fa, "w") as f:
+            for i, k in enumerate(names):
+                f.write(f">s{i}\n{seqs[k]}\n")
+        subprocess.run(
+            [
+                "mmseqs",
+                "easy-cluster",
+                str(fa),
+                f"{tmp}/clu",
+                f"{tmp}/tmp",
+                "--min-seq-id",
+                "0.3",
+                "-c",
+                "0.5",
+                "-v",
+                "1",
+                "--threads",
+                str(threads),
+            ],
+            check=True,
+            capture_output=True,
+        )
+        rep = {}
+        for line in open(f"{tmp}/clu_cluster.tsv"):
+            a, b = line.rstrip("\n").split("\t")
+            rep[names[int(b[1:])]] = names[int(a[1:])]
+    return rep
+
+
+def main():
+    ap = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+    ap.add_argument("--candidates", required=True)
+    ap.add_argument("--fasta", action="append", required=True, help="PROTEOME=PATH, repeatable")
+    ap.add_argument("--out", required=True)
+    a = ap.parse_args()
+    rows = list(csv.DictReader(open(a.candidates), delimiter="\t"))
+    out = Path(a.out)
+    out.mkdir(parents=True, exist_ok=True)
+    for spec in a.fasta:
+        prot, path = spec.split("=", 1)
+        seqs = read_fasta(path)
+        keep = [r for r in rows if r["proteome"] == prot and r["label"] in ("0", "1")]
+        # one row per protein: two curated accessions can map to one protein. Evidence of a repeat
+        # (label 1) beats an assumed negative (label 0, which is only an absent annotation).
+        by = {}
+        for r in keep:
+            by.setdefault(r["protein"], []).append(r)
+        keep = []
+        for k, v in by.items():
+            pos = [x for x in v if x["label"] == "1"]
+            if pos and len(pos) < len(v):
+                print(
+                    f"{prot}: {k}: a positive claim overrides {len(v) - len(pos)} assumed negative(s)",
+                    file=sys.stderr,
+                )
+            keep.append((pos or v)[0])
+        missing = [r["protein"] for r in keep if r["protein"] not in seqs]
+        if missing:
+            raise SystemExit(f"{prot}: {len(missing)} proteins not in the FASTA, e.g. {missing[0]}")
+        rep = cluster({r["protein"]: seqs[r["protein"]] for r in keep})
+        with open(out / f"truth.{prot}.tsv", "w", newline="") as f:
+            w = csv.writer(f, delimiter="\t", lineterminator="\n")
+            w.writerow(["id", "label", "cluster"])
+            for r in keep:
+                w.writerow([r["protein"], r["label"], rep[r["protein"]]])
+        with open(out / f"truth.{prot}.annotated.tsv", "w", newline="") as f:
+            w = csv.writer(f, delimiter="\t", lineterminator="\n")
+            w.writerow(
+                [
+                    "id",
+                    "gene",
+                    "accession",
+                    "label",
+                    "cluster",
+                    "basis",
+                    "curated_class",
+                    "evidence",
+                    "tuned_or_homolog",
+                ]
+            )
+            for r in keep:
+                w.writerow(
+                    [
+                        r["protein"],
+                        r["gene"],
+                        r["accession"],
+                        r["label"],
+                        rep[r["protein"]],
+                        r["basis"],
+                        r["curated_class"],
+                        r["evidence"],
+                        r["tuned_or_homolog"],
+                    ]
+                )
+        pos = [r for r in keep if r["label"] == "1"]
+        neg = [r for r in keep if r["label"] == "0"]
+        print(
+            f"{prot}: {len(pos)} positives in {len({rep[r['protein']] for r in pos})} clusters; "
+            f"{len(neg)} negatives in {len({rep[r['protein']] for r in neg})} clusters",
+            file=sys.stderr,
+        )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())