v2431/FormulaShowcase
0
1# report.py2# -*- coding: utf-8 -*-3 4from __future__ import annotations5 6from dataclasses import dataclass, asdict7from pathlib import Path8import json9import pickle10import re11import itertools12import urllib.request13from collections import OrderedDict14from typing import List, Dict, Any, Iterable, Tuple, Optional, Set15 16import numpy as np17import pandas as pd18import networkx as nx19import matplotlib.pyplot as plt20import matplotlib as mpl21import matplotlib.patches as mpatches22from matplotlib.path import Path as MPath23from matplotlib.patches import PathPatch24 25from pycirclize import Circos # type: ignore26from distinctipy import distinctipy # type: ignore27from goatools.obo_parser import GODag # type: ignore28from goatools.godag_plot import plot_gos # type: ignore29import plotly.graph_objects as go # type: ignore30 31from data_process import UniBioMapKG32from modules import FormulaSolution33from utils import EnrichmentAnalyzer34 35import warnings36warnings.filterwarnings(37 "ignore",38 message=r".*Glyph.*missing from font\(s\).*",39 category=UserWarning,40)41 42# 全局绘图风格：Arial + 字号 12 + SVG，支持中文43# 优先使用系统中支持中文的字体44plt.rcParams['font.family'] = [45 'Microsoft YaHei', # 微软雅黑46 'SimHei', # 黑体47 'SimSun', # 宋体48 'Arial', # 默认英文字体49 'DejaVu Sans',50 'sans-serif'51]52plt.rcParams['font.size'] = 1253plt.rcParams['axes.labelsize'] = 1254plt.rcParams['axes.titlesize'] = 1255plt.rcParams['xtick.labelsize'] = 1256plt.rcParams['ytick.labelsize'] = 1257plt.rcParams['legend.fontsize'] = 1258plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题59# 设置 SVG 导出时字体类型为 path，确保文本在 PDF 报告中正确显示60# 这样SVG中的文字会被转换为路径，避免了字体匹配问题61plt.rcParams['svg.fonttype'] = 'path'62# 设置PDF导出时字体类型，提高兼容性63plt.rcParams['pdf.fonttype'] = 4264# 确保SVG导出时不嵌入字体，而是将文本转换为路径65plt.rcParams['savefig.format'] = 'svg'66# 提高SVG输出的分辨率67plt.rcParams['figure.dpi'] = 15068 69 70@dataclass71class CommonResults:72 """Pipeline 核心数据容器（Step0-Step4 输出汇总）"""73 74 # KG75 kg: UniBioMapKG76 77 # 场景信息78 module_name: str79 module_nodes: Set[str]80 modules_dict: Dict[str, Set[str]] # module_name -> protein nodes81 82 # Step4: 组方方案83 presc_solutions: List[FormulaSolution]84 85 # Step2/3: 候选实体与网络86 ns_targets: Dict[str, Set[str]] # candidate_id -> protein targets87 ns_hit_modules: Dict[str, List[str]] # candidate_id -> [significant modules]88 sep_result: Dict[str, Any] # Step3 separation 结果89 90 # 可选的输入定义91 input_definition: Optional[str] = None92 93 def __post_init__(self):94 """类型转换：确保 Set 类型"""95 if not isinstance(self.module_nodes, set):96 self.module_nodes = set(self.module_nodes)97 98 if self.modules_dict:99 self.modules_dict = {k: set(v) for k, v in self.modules_dict.items()}100 101 if self.ns_targets:102 self.ns_targets = {k: set(v) for k, v in self.ns_targets.items()}103 104 105@dataclass106class FormulaSummary:107 """单个组方的关键数值指标汇总，用于写入 summary.json。"""108 109 index: int110 selected: List[str]111 size: int112 113 # Step4 基本指标114 mean_z: float115 mean_separation: float116 coverage: float117 118 # coverage 扩展119 weighted_coverage: float120 key_target_hit_ratio: float121 122 # hub bias123 hub_hit_ratio: float124 nonhub_coverage: float125 126 # group separation（pairwise 统计 + set-level 指标）127 separation_mean: float128 separation_median: float129 separation_min: float130 separation_max: float131 set_separation: float132 pos_pair_ratio: float133 134 135class FormulaReport:136 """对一批组方（FormulaSolution）生成评估报告。137 138 v2.1 新增/增强：139 - (3) GO DAG & 语义压缩（GOATOOLS）：从 GO 富集结果抽取 GO:ID，构建 enriched-term 子 DAG。140 - (5) Circos/Chord：herb–target 以及 target–(GO/KEGG) 的弦图；以及 herb–side_effect 风险结构。141 - (7) 网络打击效应：从一次性 LCC 改为 robustness curve + random / degree-matched baselines，并加入效率指标。142 143 输出结构（不破坏原有）：144 results/{module_name}/formula_{i}/145 - summary.json (原有)146 - targets.txt (原有)147 - enrichment/* (原有：由 EnrichmentAnalyzer 生成)148 - go_dag/* (新增)149 - circos/* (新增)150 - network_attack/* (新增)151 - llm_result.json (新增：聚合可读摘要)152 """153 def __init__(154 self,155 common_results: "CommonResults", # type: ignore156 *,157 out_root: str = "results",158 top_enrich: int = 15,159 verbose: bool = False,160 # --- v2.1: GO DAG ---161 go_obo_path: str = "database/go-basic.obo",162 go_dag_top_terms: int = 12,163 # --- v2.1: Circos/Chord ---164 chord_top_targets: int = 30,165 chord_top_terms: int = 12,166 chord_term_groups: Optional[List[str]] = None,167 # v2.2+: side-effects track 已弃用（保留参数以兼容旧调用，但不再绘制）168 circos_show_side_effects: bool = False,169 sankey_prefer_pdf: bool = True,170 # v2.2+: circos 中 module-track 改为 proximity z-score barplot，并支持阈值虚线171 circos_z_thres: float = -1.6,172 # --- v2.1: Network attack ---173 attack_steps: int = 20,174 attack_random_n: int = 50,175 attack_shuffle_n: int = 20,176 attack_eff_n_pairs: int = 3000,177 random_seed: int = 2025,178 ) -> None:179 """基于 CommonResults 初始化 FormulaReport。180 181 Parameters182 ----------183 common_results : CommonResults184 Pipeline 核心数据容器（Step0-Step4 输出汇总）185 out_root : str186 输出根目录187 top_enrich : int188 富集分析保留 top terms 数量189 verbose : bool190 是否打印详细日志191 ... (其他参数保持不变)192 """193 # 解包 CommonResults194 self.kg = common_results.kg195 self.module_name = common_results.module_name196 self.module_nodes = common_results.module_nodes197 self.presc_solutions = common_results.presc_solutions198 self.ns_targets = common_results.ns_targets199 self.ns_hit_modules = common_results.ns_hit_modules # 新增200 self.input_definition = common_results.input_definition # 新增201 202 self.sep_df: pd.DataFrame = common_results.sep_result.get("raw_separation")203 if not isinstance(self.sep_df, pd.DataFrame):204 raise ValueError("sep_result 中未找到有效的 'raw_separation' DataFrame")205 206 # 配置参数207 self.out_root = Path(out_root)208 self.top_enrich = int(top_enrich)209 self.verbose = bool(verbose)210 211 # 使用 ppi_gcc212 self.ppi_graph: nx.Graph = self.kg.ppi_gcc213 214 # v2.1 参数215 self.go_obo_path = Path(go_obo_path)216 self.go_dag_top_terms = int(go_dag_top_terms)217 self.chord_top_targets = int(chord_top_targets)218 self.chord_top_terms = int(chord_top_terms)219 self.chord_term_groups = chord_term_groups or ["kegg", "go"]220 # 兼容旧参数，但已不再绘制 side-effects track221 self.circos_show_side_effects = False222 self.sankey_prefer_pdf = bool(sankey_prefer_pdf)223 self.circos_z_thres = float(circos_z_thres)224 self.attack_steps = int(attack_steps)225 self.attack_random_n = int(attack_random_n)226 self.attack_shuffle_n = int(attack_shuffle_n)227 self.attack_eff_n_pairs = int(attack_eff_n_pairs)228 self.random_seed = int(random_seed)229 230 # v2.2: modules_dict（多模块信息）与配色231 self._prepare_modules_dict(common_results.modules_dict)232 233 # 预计算：模块诱导子图上的 degree234 self._prepare_module_centrality()235 236 # 富集分析器237 self.enrich_libraries = [238 "kegg", "go_bp", "go_cc", "go_mf",239 "hpo", "tf", "kinase", "tissue", "hallmark",240 ]241 self.enrichment_analyzer = EnrichmentAnalyzer(242 enrich_libraries=self.enrich_libraries,243 top_enrich=top_enrich,244 verbose=verbose,245 )246 247 # # v2.2+: 用于 circos barplot 的 proximity z-score 矩阵（如存在就读取）248 # self.prox_z_df: Optional[pd.DataFrame] = self._load_proximity_z_scores()249 250 # 预处理 herb 显示数据251 self._prepare_herb_display_data()252 253 # v2.2+: 尝试加载 Step2 输出的 z-score 矩阵（用于 circos module-track barplot）254 self.prox_z_df: Optional[pd.DataFrame] = self._try_load_proximity_z_scores()255 256 257 # ------------------------------------------------------------------258 # 公共入口259 # ------------------------------------------------------------------260 def run(self) -> List[FormulaSummary]:261 """对所有组方生成报告，返回数值 summary 列表。"""262 base_dir = self.out_root / self.module_name / "report"263 base_dir.mkdir(parents=True, exist_ok=True)264 265 summaries: List[FormulaSummary] = []266 267 for idx, sol in enumerate(self.presc_solutions, start=1):268 presc_dir = base_dir / f"formula_{idx}"269 presc_dir.mkdir(parents=True, exist_ok=True)270 271 summary, llm_result = self._analyze_single_formula(272 index=idx, sol=sol, presc_dir=presc_dir273 )274 summaries.append(summary)275 276 with open(presc_dir / "llm_result.json", "w", encoding="utf-8") as f:277 json.dump(llm_result, f, ensure_ascii=False, indent=2)278 279 if self.verbose:280 print(281 f"[REPORT] {self.module_name} - formula_{idx}: "282 f"size={summary.size}, coverage={summary.coverage:.3f}, "283 f"weighted_cov={summary.weighted_coverage:.3f}, "284 f"hub_hit_ratio={summary.hub_hit_ratio:.3f}, "285 f"set_sep={summary.set_separation:.3f}"286 )287 288 # 输出总览表289 summary_df = pd.DataFrame([asdict(s) for s in summaries])290 summary_df.to_csv(base_dir / "formulas_summary.tsv", sep="\t", index=False)291 292 return summaries293 294 # ------------------------------------------------------------------295 # 内部预处理：centrality / hub / 关键靶标296 # ------------------------------------------------------------------297 def _prepare_module_centrality(self) -> None:298 """在模块诱导子图上计算 degree，并定义 hub / key targets。"""299 module_nodes_in_graph = self.module_nodes & set(self.ppi_graph.nodes())300 self.module_subgraph = self.ppi_graph.subgraph(module_nodes_in_graph).copy()301 302 deg_dict = dict(self.module_subgraph.degree())303 self.module_degree: Dict[str, float] = deg_dict304 305 if len(deg_dict) == 0:306 self.hub_nodes = set()307 self.degree_sum_all = 1.0308 self.key_targets = set()309 return310 311 degrees = np.array(list(deg_dict.values()), dtype=float)312 313 # hub 定义为度数 >= 95 百分位314 thresh = np.percentile(degrees, 95)315 self.hub_nodes = {n for n, d in deg_dict.items() if d >= thresh}316 317 # key targets：中心性 top 10% 的节点318 key_thresh = np.percentile(degrees, 90)319 self.key_targets = {n for n, d in deg_dict.items() if d >= key_thresh}320 321 self.degree_sum_all = float(degrees.sum()) if degrees.sum() > 0 else 1.0322 323 # ------------------------------------------------------------------324 # 单个组方分析325 # ------------------------------------------------------------------326 def _analyze_single_formula(327 self,328 index: int,329 sol: FormulaSolution,330 presc_dir: Path,331 ) -> Tuple[FormulaSummary, Dict[str, Any]]:332 selected = list(sol.selected)333 size = len(selected)334 335 # 1) 靶标集合336 herb_targets = self._collect_targets_by_herb(selected)337 targets = set().union(*herb_targets.values()) if herb_targets else set()338 module_targets = targets & self.module_nodes339 340 # 保存靶标列表341 (presc_dir / "targets.txt").write_text("\n".join(sorted(targets)), encoding="utf-8")342 343 # 2) coverage / weighted coverage / key target 命中率344 coverage, weighted_cov, key_hit_ratio = self._compute_coverage_metrics(module_targets)345 346 # 3) hub bias347 hub_hit_ratio, nonhub_cov = self._compute_hub_bias(module_targets)348 349 # 4) separation（原有）350 (351 sep_mean,352 sep_median,353 sep_min,354 sep_max,355 set_sep,356 pos_ratio,357 ) = self._compute_group_separation_metrics(selected)358 359 # 4.1) separation heatmap（新增可视化）360 self._plot_separation_heatmap(selected, presc_dir / "separation_heatmap.svg")361 362 # 5) 富集分析（原有：由 EnrichmentAnalyzer 生成）363 enrich_dir = presc_dir / "enrichment"364 self.enrichment_analyzer.run_enrichment(365 proteins=targets,366 output_dir=enrich_dir,367 prefix="",368 )369 370 # 6) (3) GO DAG & 语义压缩（新增）371 go_dag_dir = presc_dir / "go_dag"372 go_dag_result = self._build_go_dag_from_enrichment(373 enrich_dir=enrich_dir,374 out_dir=go_dag_dir,375 top_terms=self.go_dag_top_terms,376 )377 378 # 7) (5) Circos/Chord（新增）：herb-target、target-term、herb-side_effect379 circos_dir = presc_dir / "circos"380 circos_result = self._build_circos_chords(381 selected=selected,382 herb_targets=herb_targets,383 enrich_dir=enrich_dir,384 out_dir=circos_dir,385 )386 387 # 8) (7) 网络打击效应（重写）：robustness curve + baselines388 attack_dir = presc_dir / "network_attack"389 attack_result = self._network_attack_robustness(390 hit_nodes=module_targets,391 out_dir=attack_dir,392 )393 394 # 9) 写 summary.json（保持原有字段不变）395 summary = FormulaSummary(396 index=index,397 selected=selected,398 size=size,399 mean_z=float(sol.mean_z),400 mean_separation=float(sol.mean_separation),401 coverage=float(coverage),402 weighted_coverage=float(weighted_cov),403 key_target_hit_ratio=float(key_hit_ratio),404 hub_hit_ratio=float(hub_hit_ratio),405 nonhub_coverage=float(nonhub_cov),406 separation_mean=float(sep_mean),407 separation_median=float(sep_median),408 separation_min=float(sep_min),409 separation_max=float(sep_max),410 set_separation=float(set_sep),411 pos_pair_ratio=float(pos_ratio),412 )413 414 with open(presc_dir / "summary.json", "w", encoding="utf-8") as f:415 json.dump(asdict(summary), f, ensure_ascii=False, indent=2)416 417 # 10) LLM 友好结果（尽量少数值，多结构）418 # separation：把 separation>0 的 pair 直接列出来，更符合“互补”语义419 complementary_pairs = self._get_complementary_pairs(selected, thres=0.0)420 421 # enrichment：将 enrich_libraries 的显著条目汇总进 json（只保留 Term + Adjusted P-value）422 enrichment_summary = self._collect_enrichment_summary(enrich_dir, topk=self.top_enrich)423 424 llm_result = {425 "module_name": self.module_name,426 "formula_index": index,427 "selected": selected,428 "targets": {429 "n_total": int(len(targets)),430 "n_in_module": int(len(module_targets)),431 "multi_herb_targets_top": self._top_multi_support_targets(herb_targets, topk=30),432 },433 "separation": {434 "complementary_pairs": complementary_pairs,435 },436 "enrichment": enrichment_summary,437 "go_dag": go_dag_result,438 "circos": circos_result,439 "network_attack": attack_result.get("summary_llm", {}),440 }441 442 # 添加input_definition（如果存在）443 if self.input_definition:444 llm_result["input_definition"] = self.input_definition445 446 return summary, llm_result447 448 # ------------------------------------------------------------------449 # 靶标集合450 # ------------------------------------------------------------------451 def _collect_targets_by_herb(self, selected: List[str]) -> Dict[str, Set[str]]:452 """返回：herb -> set(UniProt)。"""453 return {h: set(self.ns_targets.get(h, [])) for h in selected}454 455 # ------------------------------------------------------------------456 # coverage & weighted coverage & key targets457 # ------------------------------------------------------------------458 def _compute_coverage_metrics(self, module_targets: set) -> Tuple[float, float, float]:459 """返回：coverage, weighted_coverage, key_target_hit_ratio"""460 module_size = len(self.module_nodes)461 if module_size == 0:462 return 0.0, 0.0, 0.0463 464 coverage = len(module_targets) / float(module_size)465 466 if len(self.module_degree) == 0:467 weighted_cov = coverage468 else:469 num = float(sum(self.module_degree.get(g, 0.0) for g in module_targets))470 weighted_cov = num / float(self.degree_sum_all) if self.degree_sum_all > 0 else 0.0471 472 if len(self.key_targets) == 0:473 key_hit_ratio = 0.0474 else:475 key_hits = len(module_targets & self.key_targets)476 key_hit_ratio = key_hits / float(len(self.key_targets))477 478 return coverage, weighted_cov, key_hit_ratio479 480 # ------------------------------------------------------------------481 # hub bias482 # ------------------------------------------------------------------483 def _compute_hub_bias(self, module_targets: set) -> Tuple[float, float]:484 """返回：hub_hit_ratio, nonhub_coverage"""485 if len(self.module_degree) == 0:486 return 0.0, 0.0487 488 if len(module_targets) == 0:489 hub_hit_ratio = 0.0490 else:491 hub_hits = len(module_targets & self.hub_nodes)492 hub_hit_ratio = hub_hits / float(len(module_targets))493 494 nonhub_nodes = self.module_nodes - self.hub_nodes495 if len(nonhub_nodes) == 0:496 nonhub_cov = 0.0497 else:498 nonhub_hits = len(module_targets & nonhub_nodes)499 nonhub_cov = nonhub_hits / float(len(nonhub_nodes))500 501 return hub_hit_ratio, nonhub_cov502 503 # ------------------------------------------------------------------504 # 组方整体分离度：S_set + pairwise 统计505 # ------------------------------------------------------------------506 def _compute_group_separation_metrics(507 self, selected: List[str]508 ) -> Tuple[float, float, float, float, float, float]:509 """使用 pairwise separation 矩阵，计算组方内部的统计量。"""510 idx = [h for h in selected if h in self.sep_df.index and h in self.sep_df.columns]511 if len(idx) < 2:512 nan = float("nan")513 return nan, nan, nan, nan, nan, nan514 515 sub = self.sep_df.loc[idx, idx].astype(float)516 mat = sub.values.copy()517 n = mat.shape[0]518 519 mask = ~np.eye(n, dtype=bool)520 vals = mat[mask]521 vals = vals[~np.isnan(vals)]522 if vals.size == 0:523 nan = float("nan")524 return nan, nan, nan, nan, nan, nan525 526 sep_mean = float(np.mean(vals))527 sep_median = float(np.median(vals))528 sep_min = float(np.min(vals))529 sep_max = float(np.max(vals))530 531 nn_mins: List[float] = []532 for i in range(n):533 row = mat[i, :].astype(float)534 row[i] = np.nan535 row = row[~np.isnan(row)]536 if row.size == 0:537 continue538 nn_mins.append(float(np.min(row)))539 set_sep = float(np.mean(nn_mins)) if nn_mins else float("nan")540 541 pos = vals[vals > 0]542 pos_pair_ratio = float(len(pos) / len(vals)) if len(vals) > 0 else float("nan")543 544 return sep_mean, sep_median, sep_min, sep_max, set_sep, pos_pair_ratio545 546 # ------------------------------------------------------------------547 # (3) GO DAG & 语义压缩（GOATOOLS）548 # ------------------------------------------------------------------549 550 def _build_go_dag_from_enrichment(551 self,552 enrich_dir: Path,553 out_dir: Path,554 top_terms: int,555 ) -> Dict[str, Any]:556 """从 GO 富集结果构建 GO DAG，并产出“更适合读”的 GO 总结与官方风格图。557 558 这版做两件事：559 1) GO lineage/DAG：优先用 GOATOOLS 的绘图（plot_gos），失败再回退到 networkx 简易布局。560 2) GO 语义压缩：保留 leaf enriched term + GO-slim 映射汇总（更像论文里“主题级”GO总结）。561 """562 out_dir.mkdir(parents=True, exist_ok=True)563 564 # 读取 GO 富集结果565 go_files = {566 "go_bp": enrich_dir / "go_bp_enrichment.tsv",567 "go_cc": enrich_dir / "go_cc_enrichment.tsv",568 "go_mf": enrich_dir / "go_mf_enrichment.tsv",569 }570 571 # 1) 加载 go-basic.obo572 go_obo = self.go_obo_path573 if not go_obo.exists():574 self._download_go_obo(go_obo)575 576 godag = GODag(str(go_obo))577 578 # 2) 加载 goslim（可选）579 goslim_obo = self.go_obo_path.parent / "goslim_generic.obo"580 goslim_dag = None581 try:582 if not goslim_obo.exists():583 self._download_goslim_obo(goslim_obo)584 goslim_dag = GODag(str(goslim_obo))585 goslim_terms = set(goslim_dag.keys())586 except Exception as e:587 print(f"[GO-SLIM] Load failed -> skip slim mapping: {e}")588 goslim_terms = set()589 590 def _sort_df(df: pd.DataFrame) -> pd.DataFrame:591 # EnrichmentAnalyzer 固定输出 "Adjusted P-value" 列592 if "Adjusted P-value" in df.columns:593 return df.sort_values("Adjusted P-value", ascending=True)594 return df595 596 results: Dict[str, Any] = {}597 for tag, fp in go_files.items():598 if not fp.exists():599 continue600 df = pd.read_csv(fp, sep="\t")601 if df.empty or "Term" not in df.columns:602 continue603 604 df = _sort_df(df)605 606 # enriched GO IDs（按显著性排序，截断 top_terms）607 enriched: List[str] = []608 for term in df["Term"].tolist():609 go_id = self._extract_go_id(str(term))610 if go_id and go_id in godag:611 enriched.append(go_id)612 enriched = list(dict.fromkeys(enriched))[: int(top_terms)]613 614 if not enriched:615 continue616 617 # 语义压缩：保留“更具体”的 enriched term（剔除 enriched term 之间的祖先项）618 leaf_terms = self._semantic_compress_terms(enriched, godag)619 620 # GO-slim 映射：把每个 enriched term 投到上位 slim term（用于汇总主题）621 slim_counter: Dict[str, int] = {}622 if goslim_terms:623 for go_id in enriched:624 ancestors = set(godag[go_id].get_all_parents()) | {go_id}625 hits = ancestors & goslim_terms626 for sid in hits:627 slim_counter[sid] = slim_counter.get(sid, 0) + 1628 629 # 构建 DAG（只画 enriched + 其祖先上下文；避免全 DAG 爆炸）630 nodes_set: Set[str] = set()631 edges: List[Tuple[str, str]] = []632 for go_id in enriched:633 nodes_set.add(go_id)634 nodes_set.update(godag[go_id].get_all_parents())635 636 for go_id in list(nodes_set):637 term_obj = godag.get(go_id)638 if term_obj is None:639 continue640 for parent in term_obj.parents:641 pid = parent.id642 if pid in nodes_set:643 edges.append((pid, go_id))644 645 nodes_df = pd.DataFrame(646 [647 {648 "go_id": nid,649 "name": godag[nid].name if nid in godag else "",650 "namespace": godag[nid].namespace if nid in godag else "",651 "level": int(godag[nid].level) if nid in godag and godag[nid].level is not None else -1,652 "is_enriched": int(nid in enriched),653 "is_leaf_enriched": int(nid in leaf_terms),654 "is_goslim": int(nid in goslim_terms) if goslim_terms else 0,655 }656 for nid in sorted(nodes_set)657 ]658 )659 edges_df = pd.DataFrame(edges, columns=["parent", "child"])660 661 nodes_df.to_csv(out_dir / f"{tag}_dag_nodes.tsv", sep="\t", index=False)662 edges_df.to_csv(out_dir / f"{tag}_dag_edges.tsv", sep="\t", index=False)663 664 # 画图：优先 goatools.godag_plot.plot_gos，失败回退 networkx665 fig_path = out_dir / f"{tag}_go_dag.png"666 fig_fallback = out_dir / f"{tag}_go_dag.svg"667 668 try:669 plot_gos(str(fig_path), leaf_terms[: min(len(leaf_terms), int(top_terms))], godag, title=f"{tag.upper()} GO lineage")670 plotted = fig_path.exists()671 except Exception:672 plotted = False673 674 if not plotted:675 self._plot_go_dag_networkx(nodes_df, edges_df, fig_fallback)676 677 # plot_gos(str(fig_path), leaf_terms[: min(len(leaf_terms), int(top_terms))], godag, title=f"{tag.upper()} GO lineage")678 679 680 # GO-slim 汇总表（生物学家更易读）681 slim_tbl_path = None682 top_slim = []683 if slim_counter:684 items = sorted(slim_counter.items(), key=lambda x: x[1], reverse=True)685 top_slim = items[: min(20, len(items))]686 slim_rows = []687 for sid, cnt in top_slim:688 slim_rows.append(689 {690 "goslim_id": sid,691 "goslim_name": goslim_dag[sid].name if goslim_dag and sid in goslim_dag else (godag[sid].name if sid in godag else ""),692 "count_enriched_terms_mapped": int(cnt),693 }694 )695 slim_df = pd.DataFrame(slim_rows)696 slim_tbl_path = out_dir / f"{tag}_goslim_summary.tsv"697 slim_df.to_csv(slim_tbl_path, sep="\t", index=False)698 699 # LLM 友好：用 GO 英文名称替代 GO:ID（避免纯 ID 不可读）700 def _to_name(go_id: str) -> str:701 return str(godag[go_id].name) if go_id in godag else str(go_id)702 703 results[tag] = {704 "available": True,705 "n_enriched": int(len(enriched)),706 "top_enriched_terms": [_to_name(x) for x in enriched[: min(len(enriched), int(top_terms))]],707 "compressed_leaf_terms": [_to_name(x) for x in leaf_terms],708 "goslim_top": [709 {710 "goslim_term": (goslim_dag[sid].name if goslim_dag and sid in goslim_dag else _to_name(sid)),711 "count": int(cnt),712 }713 for sid, cnt in top_slim714 ],715 }716 717 return {"available": True, "ontologies": results}718 719 720 def _download_go_obo(self, out_path: Path) -> None:721 """下载 go-basic.obo（只需一次）。"""722 out_path.parent.mkdir(parents=True, exist_ok=True)723 url = "https://purl.obolibrary.org/obo/go/go-basic.obo"724 725 print(f"[GO-DAG] Downloading go-basic.obo -> {out_path}")726 try:727 urllib.request.urlretrieve(url, str(out_path))728 except Exception as e:729 raise RuntimeError(730 f"Failed to download go-basic.obo: {e}. "731 f"Please manually download from {url} and place it at {out_path}."732 )733 734 def _download_goslim_obo(self, out_path: Path) -> None:735 """下载 goslim_generic.obo（用于 GO 语义压缩/分组）。736 737 GOATOOLS 文档建议从 current.geneontology.org 获取 goslim_generic.obo。738 """739 out_path.parent.mkdir(parents=True, exist_ok=True)740 url = "http://current.geneontology.org/ontology/subsets/goslim_generic.obo"741 742 print(f"[GO-SLIM] Downloading goslim_generic.obo -> {out_path}")743 try:744 urllib.request.urlretrieve(url, str(out_path))745 except Exception as e:746 raise RuntimeError(747 f"Failed to download goslim_generic.obo: {e}. "748 f"Please manually download from {url} and place it at {out_path}."749 )750 751 @staticmethod752 def _extract_go_id(term: str) -> Optional[str]:753 """从 enrichr 的 Term 字段中提取 GO:xxxxxxx。"""754 m = re.search(r"(GO:\d{7})", str(term))755 return m.group(1) if m else None756 757 @staticmethod758 def _semantic_compress_terms(enriched: List[str], godag) -> List[str]:759 """保留更具体（非祖先）的 enriched terms 作为压缩后的代表。"""760 enriched_set = set(enriched)761 ancestor_set = set()762 for go_id in enriched:763 ancestor_set.update(godag[go_id].get_all_parents())764 # 如果某 enriched term 是另一个 enriched term 的祖先，则剔除765 leaf = [go_id for go_id in enriched if go_id not in (ancestor_set & enriched_set)]766 return leaf if leaf else enriched767 768 @staticmethod769 def _plot_go_dag_networkx(nodes_df: pd.DataFrame, edges_df: pd.DataFrame, out_path: Path) -> None:770 """用 networkx 画一个简易 DAG（层级布局），避免 graphviz 依赖。"""771 if nodes_df.empty or edges_df.empty:772 return773 774 G = nx.DiGraph()775 for _, r in nodes_df.iterrows():776 G.add_node(r["go_id"], level=int(r["level"]), enriched=int(r["is_enriched"]), leaf=int(r["is_leaf_enriched"]))777 for _, r in edges_df.iterrows():778 G.add_edge(r["parent"], r["child"])779 780 # 按 level 分层；level=-1 的放到最底层781 levels = {n: (G.nodes[n].get("level", -1)) for n in G.nodes()}782 uniq_levels = sorted(set(levels.values()))783 level_to_nodes: Dict[int, List[str]] = {lv: [] for lv in uniq_levels}784 for n, lv in levels.items():785 level_to_nodes[lv].append(n)786 for lv in level_to_nodes:787 level_to_nodes[lv] = sorted(level_to_nodes[lv])788 789 # 布局：同一层水平排列，层级越深 y 越低790 pos: Dict[str, Tuple[float, float]] = {}791 for yi, lv in enumerate(uniq_levels):792 nodes = level_to_nodes[lv]793 if not nodes:794 continue795 xs = np.linspace(0, 1, num=len(nodes)) if len(nodes) > 1 else np.array([0.5])796 y = -yi797 for x, n in zip(xs, nodes):798 pos[n] = (float(x), float(y))799 800 # 绘制801 plt.figure(figsize=(12, 8))802 # 边803 nx.draw_networkx_edges(G, pos, alpha=0.35, arrows=False)804 805 # 节点：enriched / leaf_enriched 强调806 base_nodes = [n for n in G.nodes() if G.nodes[n].get("enriched", 0) == 0]807 enr_nodes = [n for n in G.nodes() if G.nodes[n].get("enriched", 0) == 1]808 leaf_nodes = [n for n in G.nodes() if G.nodes[n].get("leaf", 0) == 1]809 810 nx.draw_networkx_nodes(G, pos, nodelist=base_nodes, node_size=40, alpha=0.35)811 nx.draw_networkx_nodes(G, pos, nodelist=enr_nodes, node_size=120, alpha=0.9)812 nx.draw_networkx_nodes(G, pos, nodelist=leaf_nodes, node_size=200, alpha=0.95)813 814 # 标签：只给 enriched term 打标签，避免图爆炸815 labels = {n: n for n in enr_nodes}816 nx.draw_networkx_labels(G, pos, labels=labels, font_size=8)817 818 plt.title("GO DAG (enriched terms in context)")819 plt.axis("off")820 plt.tight_layout()821 out_path.parent.mkdir(parents=True, exist_ok=True)822 plt.savefig(out_path, format="svg", bbox_inches="tight", dpi=150)823 # 同时保存为PNG格式，用于PDF报告（避免SVG字体解析问题）824 png_path = out_path.with_suffix(".png")825 plt.savefig(png_path, format="png", bbox_inches="tight", dpi=300)826 plt.close()827 828 # ------------------------------------------------------------------829 # (5) Circos/Chord830 # ------------------------------------------------------------------831 832 def _build_circos_chords(833 self,834 selected: List[str],835 herb_targets: Dict[str, Set[str]],836 enrich_dir: Path,837 out_dir: Path,838 ) -> Dict[str, Any]:839 """生成 Circos/Chord 图（v2.3）：840 841 本次改动：842 - 不再绘制 chord_herb_side_effect.svg 与 chord_herb_target.svg（信息量较低）。843 - chord_target_term 拆分为多幅图：按 term 类型分别绘制（默认 ["kegg", "go"]）。844 - circos / sankey 中 component 名称仅用于绘图时做简化（用 '|' 分隔取第一段）。845 - sankey 优先导出 PDF（Plotly + Kaleido 支持 PDF/PNG/SVG 等静态导出）。846 """847 out_dir.mkdir(parents=True, exist_ok=True)848 849 # --------------------------------------------------------------850 # 1) chord：Target–Term（按 term 类型分别绘制）851 # --------------------------------------------------------------852 top_targets = self._select_top_targets(herb_targets, topk=self.chord_top_targets)853 854 chord_term_files: Dict[str, str] = {}855 top_terms_by_group: Dict[str, List[str]] = {}856 857 for grp in (self.chord_term_groups or []):858 terms = self._select_top_terms_from_enrichment_group(enrich_dir, group=grp, topk=self.chord_top_terms)859 if not terms:860 continue861 top_terms_by_group[grp] = terms862 edges_tt = self._build_target_term_edges(863 enrich_dir,864 top_targets=top_targets,865 terms=terms,866 term_group=grp,867 )868 if not edges_tt:869 continue870 fig_tt = out_dir / f"chord_target_term_{grp}.svg"871 self._plot_chord_auto(edges_tt, fig_tt, title=f"Target–{grp.upper()} chord (most significant terms)")872 if fig_tt.exists():873 chord_term_files[grp] = fig_tt.name874 875 # --------------------------------------------------------------876 # 2) Circos：modules（side-effects track 已移除）877 # --------------------------------------------------------------878 fig_circos = out_dir / "circos_modules.svg"879 circos_ok = False880 881 if self.modules_dict and len(self.modules_dict) >= 2:882 self._plot_circos_modules_zbar(883 selected=selected,884 herb_targets=herb_targets,885 out_path=fig_circos,886 title="Source-target circos",887 z_thres=self.circos_z_thres,888 )889 circos_ok = fig_circos.exists()890 else:891 print("[CIRCOS] modules_dict not provided or <2 modules -> skip circos_modules.svg")892 893 # --------------------------------------------------------------894 # 3) Sankey：Herb → Protein → Module → Term（优先 PDF）895 # --------------------------------------------------------------896 sankey_ext = ".pdf" if self.sankey_prefer_pdf else ".html"897 sankey_go_req = out_dir / f"sankey_herb_protein_module_go{sankey_ext}"898 sankey_pw_req = out_dir / f"sankey_herb_protein_module_pathway{sankey_ext}"899 900 sankey_go_out = None901 sankey_pw_out = None902 903 if self.modules_dict and len(self.modules_dict) >= 1:904 sankey_go_out = self._plot_sankey_formula(905 selected=selected,906 herb_targets=herb_targets,907 enrich_dir=enrich_dir,908 out_path=sankey_go_req,909 mode="go",910 )911 sankey_pw_out = self._plot_sankey_formula(912 selected=selected,913 herb_targets=herb_targets,914 enrich_dir=enrich_dir,915 out_path=sankey_pw_req,916 mode="pathway",917 )918 else:919 print("[SANKEY] modules_dict not provided -> skip sankey outputs")920 921 # 注：具体图片文件路径不写入 llm_result.json（避免喂给 LLM 无关信息）。922 return {923 "available": True,924 "top_targets": [self._to_gene(t) for t in top_targets],925 "top_terms": top_terms_by_group,926 "has_circos": bool(circos_ok),927 "has_sankey": bool((sankey_go_out is not None and sankey_go_out.exists()) or (sankey_pw_out is not None and sankey_pw_out.exists())),928 "has_target_term_chord": bool(chord_term_files),929 }930 931 def _plot_chord_auto(self, edges: List[Tuple[str, str]], out_path: Path, title: str) -> None:932 """使用 pyCirclize 画 chord diagram。933 934 - term 扇区直接使用分配的颜色，连线颜色使用 term 颜色935 """936 if not edges:937 return938 939 left = [a for a, _ in edges]940 right = [b for _, b in edges]941 left_names = sorted(set(left))942 right_names = sorted(set(right))943 944 # 清理 term 名称：删除括号及内容945 def clean_term_name(term: str) -> str:946 """remove parentheses and content: 'some function (GO:12345)' -> 'some function'"""947 cleaned = re.sub(r'\s*\([^)]*\)\s*', '', str(term)).strip()948 return cleaned if cleaned else str(term)949 950 # 为 term 生成显示 label（最多20字符）和全名映射951 term_full_names: Dict[str, str] = {}952 term_display_names: Dict[str, str] = {}953 for term in right_names:954 cleaned = clean_term_name(term)955 term_full_names[term] = cleaned956 if len(cleaned) > 20:957 term_display_names[term] = cleaned[:20] + '...'958 else:959 term_display_names[term] = cleaned960 961 # 构造矩阵（使用原始term名称作为列，保证唯一性）962 mat = pd.DataFrame(0, index=left_names, columns=right_names, dtype=int)963 for a, b in edges:964 if a in mat.index and b in mat.columns:965 mat.loc[a, b] += 1966 967 # 为每个原始term生成颜色（使用 distinctipy）968 term_colors_list = self._get_distinct_colors(len(right_names), pastel_factor=0.3, colorblind_type="Deuteranomaly", rng=114514)969 term_colors: Dict[str, str] = {name: term_colors_list[i] for i, name in enumerate(right_names)}970 971 # 手动构建 circos（不使用 chord_diagram 自动模式，以便自定义颜色）972 # sectors 使用原始term名称作为key，保证数据正确性973 sectors = OrderedDict()974 975 # 添加 left (proteins) 扇区 - 使用 .values 确保获取标量976 for name in left_names:977 row_sum = mat.loc[name, :].values.sum() # 使用 .values.sum() 获取标量978 sectors[name] = int(row_sum)979 980 # 添加 right (terms) 扇区 - 使用原始term名称，使用 .values 确保获取标量981 for name in right_names:982 col_sum = mat.loc[:, name].values.sum() # 使用 .values.sum() 获取标量983 sectors[name] = int(col_sum)984 985 circos = Circos(sectors, space=2)986 987 # 绘制扇区和标签988 for sector in circos.sectors:989 name = sector.name990 991 # term 扇区：使用分配的颜色作为背景，标签使用简化名称992 if name in term_colors:993 track = sector.add_track((92, 100))994 track.axis(fc=term_colors[name], ec="white", lw=0.5, alpha=0.8)995 # 使用简化的显示名称作为标签996 display_label = term_display_names.get(name, name)997 sector.text(display_label, r=105, size=10, orientation="vertical", adjust_rotation=True)998 # protein 扇区：使用默认浅色999 else:1000 track = sector.add_track((92, 100))1001 track.axis(fc="#e0e0e0", ec="white", lw=0.5)1002 sector.text(name, r=105, size=9, orientation="vertical", adjust_rotation=True)1003 1004 # 绘制连线：使用 term 颜色1005 for (source, target), count in mat.stack().items():1006 if count > 0:1007 target_color = term_colors.get(target, "#888888")1008 # 使用 .values.sum() 获取标量值1009 source_sum = int(mat.loc[source, :].values.sum())1010 target_sum = int(mat.loc[:, target].values.sum())1011 circos.link(1012 (source, 0, source_sum),1013 (target, 0, target_sum),1014 color=target_color,1015 alpha=0.4,1016 lw=0.0,1017 )1018 1019 fig = circos.plotfig()1020 fig.suptitle(title, y=0.98)1021 1022 # 添加 legend 展示 term 全名1023 legend_handles = []1024 for term_orig in right_names:1025 full_name = term_full_names[term_orig]1026 color = term_colors[term_orig] # 使用原始term名称获取颜色1027 # legend 显示全名，最多50字符1028 label_text = full_name[:50] + '...' if len(full_name) > 50 else full_name1029 legend_handles.append(mpatches.Patch(color=color, label=label_text))1030 1031 if legend_handles:1032 fig.legend(1033 handles=legend_handles,1034 title="Terms (full names)",1035 loc="center left",1036 bbox_to_anchor=(1.0, 0.5),1037 fontsize=7,1038 title_fontsize=8,1039 frameon=True,1040 )1041 1042 fig.savefig(out_path, format="svg", dpi=300, bbox_inches="tight")1043 # 同时保存为PNG格式，用于PDF报告（避免SVG字体解析问题）1044 png_path = out_path.with_suffix(".png")1045 fig.savefig(png_path, format="png", dpi=300, bbox_inches="tight")1046 plt.close(fig)1047 # ------------------------------------------------------------------1048 # v2.2: modules_dict（多模块信息）与配色1049 # ------------------------------------------------------------------1050 1051 def _prepare_modules_dict(self, modules_dict: Optional[Dict[str, Iterable[str]]]) -> None:1052 """保存 modules_dict（module_name -> protein nodes），并预生成 module 配色。1053 1054 注意：1055 - modules_dict 的节点应与 KG 的 PPI 节点同 ID 空间（通常 UniProt）。1056 - 报告内部展示时可能用 Gene Symbol 作为 label，但 membership 判断仍用原始节点 ID。1057 """1058 self.modules_dict: Dict[str, Set[str]] = {}1059 if modules_dict:1060 self.modules_dict = {k: set(v) for k, v in modules_dict.items()}1061 1062 # module 顺序：保持输入顺序（dict 在 py3.7+ 有序），否则按 key 排1063 self.module_names: List[str] = list(self.modules_dict.keys()) if self.modules_dict else []1064 self.module_count: int = len(self.module_names)1065 1066 # module 配色：使用 distinctipy 生成高区分度颜色（更适合 module 数量较多的场景）1067 # - module: 高饱和度（pastel_factor 较低）1068 mod_colors = self._get_distinct_colors(self.module_count,1069 pastel_factor=0.2,1070 colorblind_type="Deuteranomaly",1071 rng=666)1072 self.module_colors: Dict[str, str] = {m: mod_colors[i] for i, m in enumerate(self.module_names)}1073 1074 def _prepare_herb_display_data(self) -> None:1075 """预处理 herb 显示数据和颜色映射，用于 circos/sankey 绘图。1076 1077 在所有组方生成前预先准备好这些数据，避免在每次绘图时重复计算：1078 - herb 显示标签映射（简化且唯一）1079 - herb 颜色映射（固定色板）1080 - herb 排序（按覆盖蛋白数量降序）1081 """1082 # 收集所有组方中的 selected1083 all_selected: Set[str] = set()1084 for sol in self.presc_solutions:1085 all_selected.update(sol.selected)1086 1087 # 为每个 herb 生成简化显示标签（保证唯一性）1088 self.herb_label_map: Dict[str, str] = self._make_unique_display_labels(list(all_selected))1089 self.disp2herb: Dict[str, str] = {v: k for k, v in self.herb_label_map.items()}1090 1091 # 生成固定的 herb 颜色映射（离散色板，易区分）1092 self.herb_palette = [1093 "#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854",1094 "#ffd92f", "#e5c494", "#b3b3b3", "#1b9e77", "#d95f02",1095 ]1096 1097 def _protein_to_modules(self, protein_id: str) -> List[str]:1098 """返回 protein 属于哪些 module（按 module_names 顺序）。"""1099 if not self.modules_dict:1100 return []1101 out = []1102 for m in self.module_names:1103 if protein_id in self.modules_dict.get(m, set()):1104 out.append(m)1105 return out1106 1107 def _herb_to_modules(self, herb: str, herb_targets: Dict[str, Set[str]]) -> Set[str]:1108 """返回 herb 显著命中的 module（基于 Step2 proximity z-score < -1.6）。1109 1110 Parameters1111 ----------1112 herb : str1113 候选实体 ID1114 herb_targets : Dict[str, Set[str]]1115 备用：若 ns_hit_modules 中无该 herb，回退到靶标覆盖判断1116 1117 Returns1118 -------1119 Set[str]1120 该 herb 显著命中的 module 集合1121 """1122 # 优先使用 Step2 的显著性结果1123 if herb in self.ns_hit_modules:1124 return set(self.ns_hit_modules[herb])1125 1126 # 回退：基于靶标覆盖（任一靶标落入 module 即认为覆盖，但不保证显著性）1127 if not self.modules_dict:1128 return set()1129 1130 ts = herb_targets.get(herb, set()) or set()1131 hit = set()1132 for m, nodes in self.modules_dict.items():1133 if ts & nodes:1134 hit.add(m)1135 1136 return hit1137 1138 # ------------------------------------------------------------------1139 # v2.2: CircOS（modules + side-effect 方格风格）1140 # ------------------------------------------------------------------1141 1142 def _plot_circos_modules_zbar(1143 self,1144 selected: List[str],1145 herb_targets: Dict[str, Set[str]],1146 out_path: Path,1147 title: str,1148 z_thres: float,1149 max_proteins: int = 60,1150 ) -> None:1151 """Circos：Component–Protein links + module-track（z-score barplot）。1152 1153 改动点：1154 - 移除 side-effects track1155 - component 的 module-track 把 proximity 的 z-score 画成条形图1156 - 在条形图上加一条 z_thres 的虚线作为显著阈值1157 - 条形颜色保持与 module 的配色一致1158 """1159 out_path.parent.mkdir(parents=True, exist_ok=True)1160 if not selected or not herb_targets or not self.modules_dict:1161 return1162 1163 # --- 1) 选择 proteins（避免节点爆炸） ---1164 top_targets = self._select_top_targets(herb_targets, topk=int(max_proteins))1165 top_targets_set = set(top_targets)1166 1167 herb2t: Dict[str, List[str]] = {}1168 for h in selected:1169 ts = sorted(set(herb_targets.get(h, set()) & top_targets_set))1170 if ts:1171 herb2t[h] = ts1172 if not herb2t:1173 return1174 1175 herb2t_disp: Dict[str, List[str]] = {self.herb_label_map.get(h, h): ts for h, ts in herb2t.items()}1176 1177 def _uniq_labels(uniprots: List[str]) -> Tuple[List[str], Dict[str, str]]:1178 label2id: Dict[str, str] = {}1179 labels: List[str] = []1180 seen: Set[str] = set()1181 for uid in uniprots:1182 g = self._to_gene(uid)1183 lab = g1184 if lab in seen:1185 lab = f"{g}|{uid}"1186 seen.add(lab)1187 labels.append(lab)1188 label2id[lab] = uid1189 return labels, label2id1190 1191 protein_labels, label2uid = _uniq_labels(list(dict.fromkeys(top_targets)))1192 1193 # --- 2) 扇区 size 设计：herb 大弧度，protein 小弧度 ---1194 M = max(int(self.module_count), 1)1195 sectors = OrderedDict()1196 1197 herb_sorted = sorted(list(herb2t_disp.keys()), key=lambda x: len(herb2t_disp.get(x, [])), reverse=True)1198 for h in herb_sorted:1199 n_link = len(herb2t_disp.get(h, []))1200 sectors[h] = int(max(12, n_link, M * 2))