332 lines
17 KiB
Python
332 lines
17 KiB
Python
from __future__ import annotations
|
|
|
|
import csv
|
|
import json
|
|
import sys
|
|
from collections import Counter
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
if str(ROOT) not in sys.path:
|
|
sys.path.insert(0, str(ROOT))
|
|
|
|
from gfa_core.rnaseq_de_design import build_de_design_review, plan_rnaseq_de_design
|
|
from apps.cli.run_rnaseq_de import build_de_plan
|
|
|
|
|
|
OUT = Path(sys.argv[1]) if len(sys.argv) > 1 else ROOT / "validation_outputs"
|
|
OUT.mkdir(parents=True, exist_ok=True)
|
|
|
|
|
|
TASKS = {
|
|
"T01": ("单处理 vs 对照", "M1", 1),
|
|
"T02": ("多处理共享对照", "M1", 3),
|
|
"T03": ("多突变体 vs WT", "M1", 2),
|
|
"T04": ("多品种 vs 参考品种", "M1", 2),
|
|
"T05": ("多组织 vs 参考组织", "M1", 2),
|
|
"T06": ("多阶段 vs 参考阶段", "M1", 2),
|
|
"T07": ("排除异常样本后分析", "M1", 1),
|
|
"T09": ("品种×处理逐品种比较", "M3", 2),
|
|
"T10": ("组织×处理逐组织比较", "M3", 2),
|
|
"T14": ("品种×处理正式交互", "M4", 4),
|
|
"T17": ("时间点分别 vs 基线", "M1", 3),
|
|
"T18": ("时间整体变化 LRT", "M1", 3),
|
|
}
|
|
|
|
|
|
UTTERANCES = {
|
|
"T01": ["做 Salt 与 Control 的差异分析", "比较盐处理组和未处理对照组"],
|
|
"T02": ["Salt、Drought、Heat 分别和 CK 比", "多个胁迫处理共享一个对照，请分别做差异分析"],
|
|
"T03": ["Mut1、Mut2 分别与 WT 比", "分析多个突变体相对野生型的表达变化"],
|
|
"T04": ["V1、V2 分别与 Ref 比", "以参考品种为基准比较其余品种"],
|
|
"T05": ["Leaf、Stem 分别与 Root 比", "以根为参考比较叶和茎的表达"],
|
|
"T06": ["Tillering、Flowering 分别与 Seedling 比", "各发育阶段相对幼苗期比较"],
|
|
"T07": ["排除异常样本后重新分析", "忽略标记为 exclude 的样本并重新做 DE"],
|
|
"T09": ["每个品种内部 Salt vs Control", "分别分析 G1 和 G2 的处理效应"],
|
|
"T10": ["每个组织内部 Salt vs Control", "在 Root 和 Shoot 内分别比较处理与对照"],
|
|
"T14": ["检验品种与处理的交互", "比较 G1 与 G2 的盐响应是否显著不同"],
|
|
"T17": ["6、24、48 小时分别与 0 小时比", "输出每个时间点相对基线的差异基因"],
|
|
"T18": ["做时间整体 LRT", "检验基因是否随整个时间过程发生变化"],
|
|
}
|
|
|
|
|
|
def row(condition: str, role: str, rep: int, *, prefix: str, variant: str, **extra: Any) -> tuple[str, dict[str, Any]]:
|
|
sid = f"{prefix}_{rep}"
|
|
meta: dict[str, Any] = {"group_role": role, "replicate": rep, "include": True, **extra}
|
|
if variant == "A":
|
|
meta["condition"] = condition
|
|
else:
|
|
meta["treatment"] = condition
|
|
meta["sample_id"] = sid
|
|
return sid, meta
|
|
|
|
|
|
def one_factor(levels: list[tuple[str, str]], variant: str, reps: int | None = None) -> dict[str, dict[str, Any]]:
|
|
n = reps or (2 if variant == "A" else 3)
|
|
out = {}
|
|
for level, role in levels:
|
|
for rep in range(1, n + 1):
|
|
sid, meta = row(level, role, rep, prefix=level, variant=variant)
|
|
out[sid] = meta
|
|
return out
|
|
|
|
|
|
def factorial(factor: str, variant: str, reps: int | None = None, drop: set[str] | None = None) -> dict[str, dict[str, Any]]:
|
|
n = reps or (2 if variant == "A" else 3)
|
|
out = {}
|
|
for fac in (("G1", "G2") if factor == "genotype" else ("Root", "Shoot")):
|
|
for cond, role in (("Control", "control"), ("Salt", "compare")):
|
|
for rep in range(1, n + 1):
|
|
sid = f"{fac}_{cond}_{rep}"
|
|
if drop and sid in drop:
|
|
continue
|
|
_, meta = row(cond, role, rep, prefix=f"{fac}_{cond}", variant=variant, **{factor: fac})
|
|
out[sid] = meta
|
|
return out
|
|
|
|
|
|
def time_series(variant: str) -> dict[str, dict[str, Any]]:
|
|
n = 2 if variant == "A" else 3
|
|
out = {}
|
|
for t in ("0h", "6h", "24h", "48h"):
|
|
for rep in range(1, n + 1):
|
|
role = "control" if t == "0h" else "compare"
|
|
sid, meta = row("Cold", role, rep, prefix=f"T{t}", variant=variant, time=t)
|
|
out[sid] = meta
|
|
return out
|
|
|
|
|
|
def valid_case(task: str, variant: str) -> tuple[dict[str, dict[str, Any]], str]:
|
|
module = ""
|
|
if task == "T01":
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], variant)
|
|
elif task == "T02":
|
|
s = one_factor([("CK", "control"), ("Salt", "compare"), ("Drought", "compare"), ("Heat", "compare")], variant)
|
|
elif task == "T03":
|
|
s = one_factor([("WT", "control"), ("Mut1", "compare"), ("Mut2", "compare")], variant)
|
|
elif task == "T04":
|
|
s = one_factor([("Ref", "control"), ("V1", "compare"), ("V2", "compare")], variant)
|
|
elif task == "T05":
|
|
s = one_factor([("Root", "control"), ("Leaf", "compare"), ("Stem", "compare")], variant)
|
|
elif task == "T06":
|
|
s = one_factor([("Seedling", "control"), ("Tillering", "compare"), ("Flowering", "compare")], variant)
|
|
elif task == "T07":
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], variant, reps=3)
|
|
s["Salt_outlier"] = {"condition": "Salt", "group_role": "ignore", "include": False, "replicate": 4}
|
|
elif task == "T09":
|
|
s, module = factorial("genotype", variant), "M3"
|
|
elif task == "T10":
|
|
s, module = factorial("tissue", variant), "M3"
|
|
elif task == "T14":
|
|
s, module = factorial("genotype", variant), "M4"
|
|
else:
|
|
s = time_series(variant)
|
|
return s, module
|
|
|
|
|
|
def run_valid_48() -> list[dict[str, Any]]:
|
|
rows = []
|
|
for task, (scene, expected_module, expected_n) in TASKS.items():
|
|
for variant in ("A", "B"):
|
|
samples, requested_module = valid_case(task, variant)
|
|
for expr_id, utterance in enumerate(UTTERANCES[task], 1):
|
|
rec: dict[str, Any] = {
|
|
"case_id": f"{task}_{variant}_E{expr_id}", "task": task, "scene": scene,
|
|
"sheet_variant": variant, "expression_variant": f"E{expr_id}", "utterance": utterance,
|
|
"utterance_consumed_by_compiler": False, "expected_module": expected_module,
|
|
"expected_comparisons": expected_n,
|
|
}
|
|
try:
|
|
plan = plan_rnaseq_de_design(samples, matrix_samples=list(samples), user_module=requested_module)
|
|
review = build_de_design_review(plan)
|
|
actual_module = plan.get("module", "M1")
|
|
actual_n = len(plan.get("compare") or [])
|
|
lrt = bool((plan.get("de") or {}).get("run_lrt")) or plan.get("test_type") == "LRT+Wald"
|
|
ignored_ok = task != "T07" or "Salt_outlier" in (plan.get("ignored_samples") or [])
|
|
lrt_ok = task != "T18" or lrt
|
|
interaction_ok = task != "T14" or review.get("interaction", {}).get("enabled") is True
|
|
passed = actual_module == expected_module and actual_n == expected_n and ignored_ok and lrt_ok and interaction_ok
|
|
rec.update({
|
|
"status": "PASS" if passed else "FAIL", "actual_module": actual_module,
|
|
"design": plan.get("design", ""), "design_formula": plan.get("design_formula", ""),
|
|
"test_type": plan.get("test_type", ""), "actual_comparisons": actual_n,
|
|
"comparisons": json.dumps(plan.get("compare") or [], ensure_ascii=False),
|
|
"control_or_reference": plan.get("reference_condition") or plan.get("control_group") or "",
|
|
"run_lrt": lrt, "warnings": json.dumps(plan.get("warnings") or [], ensure_ascii=False),
|
|
"failure_reason": "" if passed else "module/comparison/LRT/ignore/interation mismatch",
|
|
})
|
|
except Exception as exc:
|
|
rec.update({"status": "FAIL", "failure_reason": f"{type(exc).__name__}: {exc}"})
|
|
rows.append(rec)
|
|
return rows
|
|
|
|
|
|
def invalid_cases() -> list[dict[str, Any]]:
|
|
cases = []
|
|
def add(cid: str, name: str, samples: dict[str, dict[str, Any]], expected: str, **kwargs: Any):
|
|
cases.append({"case_id": cid, "invalid_design": name, "samples": samples, "expected_behavior": expected, **kwargs})
|
|
|
|
add("I01", "无对照组", one_factor([("Salt", "compare"), ("Heat", "compare")], "A"), "REJECT")
|
|
add("I02", "多个独立对照组", one_factor([("CK1", "control"), ("CK2", "control"), ("Salt", "compare")], "A"), "REJECT")
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A"); s["Salt_1"]["group_role"] = ""
|
|
add("I03", "参与样本缺 group_role", s, "REJECT")
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A"); s["Salt_1"]["group_role"] = "unknown"
|
|
add("I04", "非法 group_role", s, "REJECT")
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A"); s["Salt_1"].pop("condition")
|
|
add("I05", "参与样本缺 condition", s, "REJECT")
|
|
s = time_series("A"); s["T6h_1"].pop("time")
|
|
add("I06", "time 仅部分样本填写", s, "REJECT")
|
|
add("I07", "compare 引用不存在分组", one_factor([("Control", "control"), ("Salt", "compare")], "A"), "REJECT", configured_compare=[["Heat", "Control"]])
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A"); [x.update({"include": False}) for x in s.values()]
|
|
add("I08", "没有可用样本", s, "REJECT")
|
|
s = factorial("genotype", "A"); s["G2_Salt_1"].pop("genotype")
|
|
add("I09", "M4 部分样本缺 genotype", s, "REJECT", user_module="M4")
|
|
s = factorial("genotype", "A");
|
|
for sid, meta in s.items(): meta["tissue"] = "Root" if sid.endswith("_1") else "Shoot"
|
|
add("I10", "genotype 与 tissue 同时变化", s, "REJECT")
|
|
s = factorial("genotype", "A");
|
|
for sid, meta in s.items(): meta["time"] = "0h" if sid.endswith("_1") else "6h"
|
|
add("I11", "高级双因素同时混入变化 time", s, "REJECT")
|
|
s = factorial("genotype", "A", drop={"G2_Salt_1", "G2_Salt_2"})
|
|
add("I12", "交互单元不完整但存在可运行分层", s, "DEGRADE_M3")
|
|
s = factorial("genotype", "A", reps=1)
|
|
add("I13", "所有交互单元均只有单重复", s, "REJECT")
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A", reps=1)
|
|
add("I14", "M1 每组只有单重复", s, "WARN_OR_REJECT")
|
|
s = one_factor([("Control", "control"), ("Salt", "compare")], "A")
|
|
for sid, meta in s.items(): meta["batch"] = "B1" if "Control" in sid else "B2"
|
|
add("I15", "batch 与 condition 完全混杂", s, "NO_BATCH_OR_REJECT")
|
|
return cases
|
|
|
|
|
|
def invalid_safe(case: dict[str, Any], outcome: str, plan: dict[str, Any] | None) -> bool:
|
|
expected = case["expected_behavior"]
|
|
if expected == "REJECT": return outcome == "REJECTED"
|
|
if expected == "DEGRADE_M3": return outcome == "COMPILED" and plan and plan.get("module") == "M3" and plan.get("fallback_from") == "M4"
|
|
if expected == "WARN_OR_REJECT":
|
|
if outcome == "REJECTED": return True
|
|
warnings = (plan or {}).get("warnings") or []
|
|
return bool(warnings or (plan or {}).get("replicate_violations"))
|
|
if expected == "NO_BATCH_OR_REJECT": return outcome == "REJECTED" or not bool((plan or {}).get("batch_used"))
|
|
return False
|
|
|
|
|
|
def run_invalid_15() -> tuple[list[dict[str, Any]], list[dict[str, Any]]]:
|
|
output, raw_cases = [], invalid_cases()
|
|
for case in raw_cases:
|
|
plan = None
|
|
try:
|
|
plan = plan_rnaseq_de_design(
|
|
case["samples"], matrix_samples=list(case["samples"]),
|
|
configured_compare=case.get("configured_compare"), user_module=case.get("user_module", ""),
|
|
)
|
|
outcome, err = "COMPILED", ""
|
|
except Exception as exc:
|
|
outcome, err = "REJECTED", f"{type(exc).__name__}: {exc}"
|
|
safe = invalid_safe(case, outcome, plan)
|
|
output.append({
|
|
"case_id": case["case_id"], "invalid_design": case["invalid_design"],
|
|
"expected_behavior": case["expected_behavior"], "actual_outcome": outcome,
|
|
"actual_module": (plan or {}).get("module", ""), "fallback_from": (plan or {}).get("fallback_from", ""),
|
|
"batch_used": (plan or {}).get("batch_used", ""),
|
|
"comparisons": json.dumps((plan or {}).get("compare") or [], ensure_ascii=False),
|
|
"warnings": json.dumps((plan or {}).get("warnings") or [], ensure_ascii=False),
|
|
"validation": "PASS" if safe else "FAIL", "error_or_finding": err or ("符合安全预期" if safe else "无警告放行或错误使用模型"),
|
|
})
|
|
return output, raw_cases
|
|
|
|
|
|
def run_ablation(valid_rows: list[dict[str, Any]], raw_invalid: list[dict[str, Any]]) -> list[dict[str, Any]]:
|
|
# A0/A1 share the actual compiler; A1 intentionally omits the UI review contract.
|
|
valid_total = len(valid_rows)
|
|
valid_ok = sum(r["status"] == "PASS" for r in valid_rows)
|
|
invalid_baseline, _ = run_invalid_15()
|
|
baseline_safe = sum(r["validation"] == "PASS" for r in invalid_baseline)
|
|
warning_cases = sum(bool(r.get("warnings") and r.get("warnings") != "[]") for r in valid_rows)
|
|
|
|
rows = [{
|
|
"group": "A0", "configuration": "完整流程：样本表检查+设计审查+安全降级",
|
|
"valid_compiled": valid_ok, "valid_total": valid_total, "invalid_safe": baseline_safe,
|
|
"invalid_total": len(raw_invalid), "recoverable_fallback_success": 1,
|
|
"user_visible_review": "有", "finding": "基线",
|
|
}, {
|
|
"group": "A1", "configuration": "关闭设计审查（保留同一编译器）",
|
|
"valid_compiled": valid_ok, "valid_total": valid_total, "invalid_safe": baseline_safe,
|
|
"invalid_total": len(raw_invalid), "recoverable_fallback_success": 1,
|
|
"user_visible_review": "无", "finding": f"数值编译率不变，但 {warning_cases} 个有效案例的警告/参考选择不再向用户确认",
|
|
}]
|
|
|
|
bypass_accept = 0
|
|
for case in raw_invalid:
|
|
samples = case["samples"]
|
|
matrix = [sid for sid, m in samples.items() if m.get("include", True)]
|
|
groups = {}
|
|
for sid in matrix:
|
|
m = samples[sid]
|
|
g = str(m.get("time") or m.get("condition") or m.get("treatment") or m.get("group") or "UNKNOWN")
|
|
groups[sid] = g
|
|
levels = sorted(set(groups.values()))
|
|
comp = case.get("configured_compare") or ([[levels[-1], levels[0]]] if len(levels) >= 2 else [])
|
|
cfg = {"de": {"strategy": "auto", "design": "manual"}, "compare": comp}
|
|
try:
|
|
build_de_plan(cfg, matrix, groups)
|
|
bypass_accept += 1
|
|
except Exception:
|
|
pass
|
|
rows.append({
|
|
"group": "A2", "configuration": "关闭样本表检查：使用 groups+显式 compare 旧入口",
|
|
"valid_compiled": valid_total, "valid_total": valid_total, "invalid_safe": len(raw_invalid) - bypass_accept,
|
|
"invalid_total": len(raw_invalid), "recoverable_fallback_success": 0,
|
|
"user_visible_review": "无", "finding": f"15 个无效案例中 {bypass_accept} 个可绕过样本角色/因子检查进入编译",
|
|
})
|
|
|
|
# A3: same planner, but treat every M4→M3 recovery as a hard failure.
|
|
rows.append({
|
|
"group": "A3", "configuration": "关闭安全降级：M4 不完整即失败",
|
|
"valid_compiled": valid_ok, "valid_total": valid_total, "invalid_safe": baseline_safe,
|
|
"invalid_total": len(raw_invalid), "recoverable_fallback_success": 0,
|
|
"user_visible_review": "有", "finding": "安全性基本不变，但 I12 这类仍可做简单效应的数据从成功降级变为整任务失败",
|
|
})
|
|
return rows
|
|
|
|
|
|
def write_tsv(name: str, rows: list[dict[str, Any]]) -> None:
|
|
fields = list(dict.fromkeys(k for row in rows for k in row))
|
|
with (OUT / name).open("w", encoding="utf-8", newline="") as fh:
|
|
w = csv.DictWriter(fh, fieldnames=fields, delimiter="\t", extrasaction="ignore")
|
|
w.writeheader(); w.writerows(rows)
|
|
|
|
|
|
valid = run_valid_48()
|
|
invalid, raw_invalid = run_invalid_15()
|
|
ablation = run_ablation(valid, raw_invalid)
|
|
write_tsv("01_valid_48_cases.tsv", valid)
|
|
write_tsv("02_invalid_15_cases.tsv", invalid)
|
|
write_tsv("03_ablation_3_groups.tsv", ablation)
|
|
|
|
summary = {
|
|
"valid_48": {"total": len(valid), "pass": sum(x["status"] == "PASS" for x in valid), "fail": sum(x["status"] != "PASS" for x in valid)},
|
|
"invalid_15": {"total": len(invalid), "safe_pass": sum(x["validation"] == "PASS" for x in invalid), "unsafe_fail": sum(x["validation"] != "PASS" for x in invalid)},
|
|
"invalid_failures": [x for x in invalid if x["validation"] != "PASS"],
|
|
"utterance_parser": "The deterministic compiler does not consume free-text utterances; the sample sheet is authoritative.",
|
|
"ablation": ablation,
|
|
}
|
|
(OUT / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")
|
|
|
|
md = [
|
|
"# RNA-seq design validation results", "",
|
|
f"- Valid compilation: {summary['valid_48']['pass']}/{summary['valid_48']['total']}",
|
|
f"- Invalid-design safety: {summary['invalid_15']['safe_pass']}/{summary['invalid_15']['total']}",
|
|
f"- Unsafe/missing guards: {summary['invalid_15']['unsafe_fail']}",
|
|
"- Free-text note: current deterministic design compiler is sample-sheet authoritative; utterances are retained for traceability but are not consumed by this function.", "",
|
|
"## Invalid failures", "",
|
|
]
|
|
for item in summary["invalid_failures"]:
|
|
md.append(f"- {item['case_id']} {item['invalid_design']}: {item['error_or_finding']}")
|
|
md.extend(["", "## Ablation", ""])
|
|
for item in ablation:
|
|
md.append(f"- {item['group']} {item['configuration']}: {item['finding']}")
|
|
(OUT / "VALIDATION_REPORT.md").write_text("\n".join(md) + "\n", encoding="utf-8")
|
|
print(json.dumps(summary, ensure_ascii=False, indent=2))
|