302 lines
11 KiB
Python
302 lines
11 KiB
Python
|
from pathlib import Path
|
|||
|
|
import subprocess
|
|||
|
|
import sys
|
|||
|
|
|
|||
|
|
import pytest
|
|||
|
|
import yaml
|
|||
|
|
|
|||
|
|
ROOT = Path(__file__).resolve().parents[2]
|
|||
|
|
if str(ROOT) not in sys.path:
|
|||
|
|
sys.path.insert(0, str(ROOT))
|
|||
|
|
|
|||
|
|
from apps.cli.run_rnaseq_de import build_de_plan, generate_modular_r_script, run_deseq2
|
|||
|
|
from apps.cli.run_rnaseq_fgsea import find_de_results
|
|||
|
|
from apps.cli.run_samples_autofill import design_from_samples
|
|||
|
|
from gfa_core.rnaseq_de_design import (
|
|||
|
|
build_de_design_review,
|
|||
|
|
plan_rnaseq_de_design,
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def _factorial(*, drop=(), time_values=None):
|
|||
|
|
samples = {}
|
|||
|
|
for genotype in ("G1", "G2"):
|
|||
|
|
for condition, role in (("Control", "control"), ("Salt", "compare")):
|
|||
|
|
for replicate in (1, 2):
|
|||
|
|
sample = f"{genotype}_{condition}_{replicate}"
|
|||
|
|
if sample in drop:
|
|||
|
|
continue
|
|||
|
|
row = {
|
|||
|
|
"condition": condition,
|
|||
|
|
"genotype": genotype,
|
|||
|
|
"group_role": role,
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
if time_values:
|
|||
|
|
row["time"] = time_values[replicate - 1]
|
|||
|
|
samples[sample] = row
|
|||
|
|
return samples
|
|||
|
|
|
|||
|
|
|
|||
|
|
def _plan(samples, module=""):
|
|||
|
|
matrix_samples = list(samples)
|
|||
|
|
groups = {sample: row["condition"] for sample, row in samples.items()}
|
|||
|
|
return build_de_plan({"samples": samples, "de": {"strategy": "auto"}}, matrix_samples, groups, user_module=module)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_complete_factorial_selects_m4_and_keeps_m3_simple_effects():
|
|||
|
|
plan = _plan(_factorial())
|
|||
|
|
assert plan["module"] == "M4"
|
|||
|
|
assert plan["advanced_factor"] == "genotype"
|
|||
|
|
assert plan["comparisons"] == [
|
|||
|
|
["G1__Salt", "G1__Control"],
|
|||
|
|
["G2__Salt", "G2__Control"],
|
|||
|
|
["G2__Control", "G1__Control"],
|
|||
|
|
["G2__Salt", "G1__Salt"],
|
|||
|
|
]
|
|||
|
|
assert plan["comparison_axis"] == "condition"
|
|||
|
|
assert plan["reference_condition"] == "Control"
|
|||
|
|
assert plan["reference_factor"] == "G1"
|
|||
|
|
assert plan["interaction_full"] == "~ genotype * condition"
|
|||
|
|
assert plan["fallback_from"] == ""
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_incomplete_factorial_degrades_to_valid_m3_strata():
|
|||
|
|
samples = _factorial(drop={"G2_Salt_1", "G2_Salt_2"})
|
|||
|
|
plan = _plan(samples)
|
|||
|
|
assert plan["module"] == "M3"
|
|||
|
|
assert plan["comparisons"] == [["G1__Salt", "G1__Control"]]
|
|||
|
|
assert plan["fallback_from"] == "M4"
|
|||
|
|
assert plan["skipped_comparisons"]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_varying_time_is_not_silently_ignored_by_m3_m4():
|
|||
|
|
with pytest.raises(ValueError, match="time"):
|
|||
|
|
_plan(_factorial(time_values=("0h", "6h")))
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_m4_script_emits_group_contrasts_and_unsigned_lrt(tmp_path):
|
|||
|
|
samples = _factorial()
|
|||
|
|
plan = _plan(samples)
|
|||
|
|
script = generate_modular_r_script(
|
|||
|
|
tmp_path / "counts.txt",
|
|||
|
|
list(samples),
|
|||
|
|
plan,
|
|||
|
|
tmp_path,
|
|||
|
|
0.05,
|
|||
|
|
1.0,
|
|||
|
|
)
|
|||
|
|
assert "dds_group" in script
|
|||
|
|
assert 'test="LRT"' in script
|
|||
|
|
assert "reduced=~ genotype + condition" in script
|
|||
|
|
assert "genotype_by_condition.interaction_LRT.DE_results.txt" in script
|
|||
|
|
assert "dds_wald" in script
|
|||
|
|
assert "interaction_Wald" in script
|
|||
|
|
assert "falling back to ~condition" not in script
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_group_role_may_define_genotype_axis_in_standard_factorial():
|
|||
|
|
samples = {}
|
|||
|
|
for genotype, role in (("WT", "control"), ("ire1dm", "compare")):
|
|||
|
|
for condition in ("DMSO", "Tm"):
|
|||
|
|
for replicate in (1, 2):
|
|||
|
|
sample = f"{genotype}_{condition}_{replicate}"
|
|||
|
|
samples[sample] = {
|
|||
|
|
"genotype": genotype,
|
|||
|
|
"condition": condition,
|
|||
|
|
"group_role": role,
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
plan = _plan(samples)
|
|||
|
|
assert plan["module"] == "M4"
|
|||
|
|
assert plan["comparison_axis"] == "genotype"
|
|||
|
|
assert plan["stratify_axis"] == "condition"
|
|||
|
|
assert plan["reference_factor"] == "WT"
|
|||
|
|
assert plan["reference_condition"] == "DMSO"
|
|||
|
|
assert plan["comparisons"] == [
|
|||
|
|
["ire1dm__DMSO", "WT__DMSO"],
|
|||
|
|
["ire1dm__Tm", "WT__Tm"],
|
|||
|
|
["WT__Tm", "WT__DMSO"],
|
|||
|
|
["ire1dm__Tm", "ire1dm__DMSO"],
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_blank_roles_use_common_reference_names_without_rejecting_m4():
|
|||
|
|
samples = {}
|
|||
|
|
for genotype in ("WT", "ire1dm"):
|
|||
|
|
for condition in ("DMSO", "Tm"):
|
|||
|
|
for replicate in (1, 2):
|
|||
|
|
sample = f"{genotype}_{condition}_{replicate}"
|
|||
|
|
samples[sample] = {
|
|||
|
|
"genotype": genotype,
|
|||
|
|
"condition": condition,
|
|||
|
|
"group_role": "",
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
plan = _plan(samples)
|
|||
|
|
assert plan["module"] == "M4"
|
|||
|
|
assert plan["reference_factor"] == "WT"
|
|||
|
|
assert plan["reference_condition"] == "DMSO"
|
|||
|
|
assert plan["factor_reference_source"] == "semantic"
|
|||
|
|
assert plan["condition_reference_source"] == "semantic"
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_fgsea_excludes_unsigned_m4_lrt(tmp_path):
|
|||
|
|
directional = tmp_path / "Os.rnaseq_de.G1__Salt_vs_G1__Control.DE_results.txt"
|
|||
|
|
interaction = tmp_path / "Os.rnaseq_de.genotype_by_condition.interaction_LRT.DE_results.txt"
|
|||
|
|
directional.write_text("gene_id\tstat\n", encoding="utf-8")
|
|||
|
|
interaction.write_text("gene_id\tstat\n", encoding="utf-8")
|
|||
|
|
assert find_de_results(tmp_path, "Os") == [directional]
|
|||
|
|
with pytest.raises(ValueError, match="not eligible"):
|
|||
|
|
find_de_results(tmp_path, "Os", "genotype_by_condition.interaction_LRT")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_frontend_review_matches_executable_m4_plan():
|
|||
|
|
samples = {}
|
|||
|
|
for genotype, role in (("WT", "control"), ("ire1dm", "compare")):
|
|||
|
|
for condition in ("DMSO", "Tm"):
|
|||
|
|
for replicate in (1, 2):
|
|||
|
|
sample = f"{genotype}_{condition}_{replicate}"
|
|||
|
|
samples[sample] = {
|
|||
|
|
"genotype": genotype,
|
|||
|
|
"condition": condition,
|
|||
|
|
"group_role": role,
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
plan = plan_rnaseq_de_design(samples, matrix_samples=list(samples))
|
|||
|
|
review = build_de_design_review(plan)
|
|||
|
|
assert review["module"] == "M4"
|
|||
|
|
assert review["comparison_axis"] == "genotype"
|
|||
|
|
assert [item["display"] for item in review["comparisons"]] == [
|
|||
|
|
"DMSO 条件下：ire1dm vs WT",
|
|||
|
|
"Tm 条件下：ire1dm vs WT",
|
|||
|
|
"WT 中：Tm vs DMSO",
|
|||
|
|
"ire1dm 中：Tm vs DMSO",
|
|||
|
|
]
|
|||
|
|
assert review["interaction"] == {
|
|||
|
|
"enabled": True,
|
|||
|
|
"display": "genotype × condition 交互效应",
|
|||
|
|
"methods": ["LRT", "Wald"],
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_saved_module_override_is_used_by_de_execution():
|
|||
|
|
samples = _factorial()
|
|||
|
|
groups = {sample: row["condition"] for sample, row in samples.items()}
|
|||
|
|
plan = build_de_plan(
|
|||
|
|
{"samples": samples, "de": {"strategy": "auto", "module": "M3"}},
|
|||
|
|
list(samples),
|
|||
|
|
groups,
|
|||
|
|
)
|
|||
|
|
assert plan["module"] == "M3"
|
|||
|
|
assert plan["fallback_from"] == "M4"
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_optional_batch_is_added_to_m1_and_frontend_review():
|
|||
|
|
samples = {}
|
|||
|
|
for condition, role in (("Control", "control"), ("Salt", "compare")):
|
|||
|
|
for replicate, batch in ((1, "B1"), (2, "B2")):
|
|||
|
|
sample = f"{condition}_{replicate}"
|
|||
|
|
samples[sample] = {
|
|||
|
|
"condition": condition,
|
|||
|
|
"group_role": role,
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"batch": batch,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
plan = plan_rnaseq_de_design(samples, matrix_samples=list(samples))
|
|||
|
|
review = build_de_design_review(plan)
|
|||
|
|
assert plan["module"] == "M1"
|
|||
|
|
assert plan["batch_requested"] is True
|
|||
|
|
assert plan["batch_used"] is True
|
|||
|
|
assert plan["design_formula"] == "~ batch + group"
|
|||
|
|
assert review["batch"]["used"] is True
|
|||
|
|
assert review["batch"]["levels"] == ["B1", "B2"]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_incomplete_batch_falls_back_without_blocking_design():
|
|||
|
|
samples = {
|
|||
|
|
"Control_1": {"condition": "Control", "group_role": "control", "batch": "B1"},
|
|||
|
|
"Control_2": {"condition": "Control", "group_role": "control", "batch": ""},
|
|||
|
|
"Salt_1": {"condition": "Salt", "group_role": "compare", "batch": "B1"},
|
|||
|
|
"Salt_2": {"condition": "Salt", "group_role": "compare", "batch": ""},
|
|||
|
|
}
|
|||
|
|
plan = plan_rnaseq_de_design(samples, matrix_samples=list(samples))
|
|||
|
|
review = build_de_design_review(plan)
|
|||
|
|
assert plan["module"] == "M1"
|
|||
|
|
assert plan["batch_requested"] is True
|
|||
|
|
assert plan["batch_used"] is False
|
|||
|
|
assert plan["batch_fallback"] is True
|
|||
|
|
assert plan["design_formula"] == "~ group"
|
|||
|
|
assert review["batch"]["fallback"] is True
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_optional_batch_is_used_by_m4_script(tmp_path):
|
|||
|
|
samples = _factorial()
|
|||
|
|
for sample, row in samples.items():
|
|||
|
|
row["batch"] = "B1" if sample.endswith("_1") else "B2"
|
|||
|
|
plan = _plan(samples)
|
|||
|
|
script = generate_modular_r_script(
|
|||
|
|
tmp_path / "counts.txt",
|
|||
|
|
list(samples),
|
|||
|
|
plan,
|
|||
|
|
tmp_path,
|
|||
|
|
0.05,
|
|||
|
|
1.0,
|
|||
|
|
)
|
|||
|
|
assert plan["batch_used"] is True
|
|||
|
|
assert plan["design_formula"] == "~ batch + genotype * condition"
|
|||
|
|
assert "batch = factor(batch_vals" in script
|
|||
|
|
assert "design=~ batch + genotype * condition" in script
|
|||
|
|
assert "reduced=~ batch + genotype + condition" in script
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test_batch_runtime_failure_retries_once_without_batch(tmp_path, monkeypatch):
|
|||
|
|
samples = {}
|
|||
|
|
for condition, role in (("Control", "control"), ("Salt", "compare")):
|
|||
|
|
for replicate, batch in ((1, "B1"), (2, "B2")):
|
|||
|
|
sample = f"{condition}_{replicate}"
|
|||
|
|
samples[sample] = {
|
|||
|
|
"condition": condition,
|
|||
|
|
"group_role": role,
|
|||
|
|
"replicate": replicate,
|
|||
|
|
"batch": batch,
|
|||
|
|
"include": True,
|
|||
|
|
}
|
|||
|
|
config = tmp_path / "samples.yaml"
|
|||
|
|
config.write_text(yaml.safe_dump({"samples": samples}), encoding="utf-8")
|
|||
|
|
count_matrix = tmp_path / "counts.txt"
|
|||
|
|
count_matrix.write_text(
|
|||
|
|
"gene_id\tControl_1\tControl_2\tSalt_1\tSalt_2\n"
|
|||
|
|
"gene1\t10\t11\t30\t31\n"
|
|||
|
|
"gene2\t20\t19\t5\t6\n",
|
|||
|
|
encoding="utf-8",
|
|||
|
|
)
|
|||
|
|
output_dir = tmp_path / "de"
|
|||
|
|
output_dir.mkdir()
|
|||
|
|
calls = []
|
|||
|
|
|
|||
|
|
def fake_rscript(*args, **kwargs):
|
|||
|
|
calls.append(args)
|
|||
|
|
if len(calls) == 1:
|
|||
|
|
raise subprocess.CalledProcessError(
|
|||
|
|
1,
|
|||
|
|
args[0],
|
|||
|
|
stderr="model matrix is not full rank",
|
|||
|
|
)
|
|||
|
|
return subprocess.CompletedProcess(args[0], 0, "", "")
|
|||
|
|
|
|||
|
|
monkeypatch.setattr("apps.cli.run_rnaseq_de.subprocess.run", fake_rscript)
|
|||
|
|
status, plan = run_deseq2(count_matrix, config, output_dir)
|
|||
|
|
assert len(calls) == 2
|
|||
|
|
assert all(status.values())
|
|||
|
|
assert plan["batch_requested"] is True
|
|||
|
|
assert plan["batch_used"] is False
|
|||
|
|
assert plan["batch_fallback"] is True
|
|||
|
|
assert (output_dir / "BATCH_FALLBACK.txt").exists()
|
|||
|
|
assert "batch +" not in (output_dir / "run_deseq2.R").read_text(encoding="utf-8")
|