一、這段程式在做什麼（原理）
它其實是一個 「先驗拓撲 × 資料權重」的雙層網路建模，四個階段：
| 階段 | 做的事 | 輸出 |
|---|---|---|
| ① 先驗層 | 從 KEGG KGML 抽出 hsa04360 的 entry / relation | 「哪些基因之間可能有交互作用」 |
| ② 資料層 | 讀表達矩陣 → 正規化 → 算基因間 Pearson r | 「這些交互作用在你的樣本裡共變得多強」 |
| ③ 分群 | 只保留「KEGG 有邊 且 |r| ≥ 門檻」的邊，權重 = |r|，跑加權 Leiden | 模組 (module) |
| ④ 詮釋 | degree / betweenness / log2FC + 五張圖 + 多工作表 Excel | 可發表的結果 |
核心邏輯是 交集式約束：資料驅動的共表達網路（如 WGCNA）容易產生大量無生物意義的邊；純 KEGG 拓撲則是靜態的、與你的樣本無關。這裡讓一條邊必須同時通過「文獻先驗」與「你的數據」兩道關卡，是合理且常見的策略（概念上接近 pathway-constrained co-expression network）。
BUGFIX 清單裡真正有價值的修正（不是在湊數）：
- [1] AnnData 方向 — 這是最致命的。
normalize_total是**沿著 obs（列）**做總量對齊。原本 obs = 基因，等於「把每個基因跨樣本的總和拉成一樣」，直接抹平樣本間差異，後續相關性算的是雜訊。修正正確。 - [2]
ModularityVertexPartition不吃resolution_parameter— 事實如此，會 TypeError；改RBConfigurationVertexPartition是對的（res=1 時等價於標準 modularity）。 - [4] group entry 展開成 component — KEGG 的
group是蛋白複合體，沒有graphics/@name，原程式會生出一個叫undefined的幽靈節點並吸走大量邊。這個修正很關鍵。 - [3] 排序決定論、[6] 過濾 maplink、[7] 移除孤立節點 — 都是真問題，修得對。
所以：作為工程品質的重構，這份程式明顯優於原始版本。
二、有沒有解決「有向定量」的問題？
直說：「定量」解決了一半，「有向」完全沒有解決 —— 而且是被主動放棄的。
1. 有向性在解析階段就被丟掉了
G = nx.Graph() # ← parse_kgml 裡是無向圖第一段程式碼原本用的是 nx.DiGraph()，重構版反而降級成 nx.Graph()。KGML 裡 entry1 → entry2 的方向、以及 subtype 的 activation / inhibition，正是「有向」資訊的全部來源。程式雖然把 subtypes 存成邊屬性，但從頭到尾沒有任何一行用它——不影響權重、不影響分群、不影響圖上的箭頭。等於抄下來當註解。
2. Pearson r 在數學上不可能有方向
corr.loc[u,v] == corr.loc[v,u]，相關係數是對稱的。無論怎麼包裝，用 r 當權重就永遠得不到有向網路。而且 r 也不是因果量——KEGG 說「A 磷酸化 B」，你算的卻是 A、B 的 mRNA 共變，兩者在生物層次上根本不是同一件事。
3. 這點對 hsa04360 特別致命
Axon guidance 是訊息傳導通路，主要事件是磷酸化、GTP/GDP 交換、複合體組裝——全部發生在蛋白質與轉譯後修飾層次。RHOA → ROCK1 這條邊的強弱跟兩者 mRNA 是否共變幾乎無關。用共表達當這類邊的「定量」，理論基礎是薄弱的（用在轉錄調控網路 GErel 上會合理得多）。
4. n=6 讓「定量」本身也站不住
3 對 AT/MT = 6 個樣本，df=4。Pearson r 的雙尾 p=0.05 臨界值是 r=0.811。
--min-abs-r 0.30→ p ≈ 0.56，等於沒有過濾，留下的邊有一半以上是純雜訊。- 就算調到 r=0.8，p ≈ 0.056，仍然只是勉強擦邊，且未做多重檢定校正。
更麻煩的是，這 6 欄的變異主要來自「病人 1440 / 2003 / 2613」的個體差異，不是 AT/MT 的生物學。若目標是腫瘤相關的調控，至少該先做 病人內置中（within-patient centering） 再算相關，否則算出來的是「病人指紋」。
5. 其他幾個實質 bug（不在 BUGFIX 清單裡）
btw = nx.betweenness_centrality(W, weight=None) # 辛苦建的權重完全沒用上而且不能只是改成 weight="weight" —— networkx 把 weight 當成距離，會導致「相關性越強 = 距離越遠」，語義完全相反。正確做法是另加一個 distance = 1 - |r| 屬性再傳進去。目前 top_hub 是純拓撲樞紐，與表達量無關。
"p_paired": fc_tbl.loc[src, "p_paired"].mean() # p 值取平均沒有統計意義多個 transcript 對應同一 symbol 時，應該用 Fisher / Stouffer 合併，或選代表性 transcript。另外全流程沒有 FDR 校正。
id2keggids 抓了 KEGG 的 hsa:xxxx（Entrez ID）卻從未使用。你的表格既然有 Transcriptid，用 Entrez/Ensembl 做映射會比手刻 alias 字典穩健一個數量級——後者永遠補不完。
還有：檔名是 FC_GSEA.xlsx。若那些欄位已經是正規化過的表達量（甚至是 log ratio），再跑一次 normalize_total + log1p 是重複正規化；若含負值，log1p 會直接壞掉。程式沒有任何檢查。
三、要真正做到「有向定量」，該怎麼改
最低成本的一步（改 3 行，馬上有方向）：
g = ig.Graph(n=len(nodes), directed=True) # leidenalg 的 RBConfiguration 支援有向圖igraph 會改用 directed modularity。但這只是形式上有向，權重仍是對稱的 r，實質意義有限。
真正對的路線，依投入程度排序：
- 符號一致性權重（最務實） — 從 subtype 取期望符號（activation/expression = +1，inhibition/repression = −1），權重 = 「觀測到的 r 符號是否符合 KEGG 預期」。這是有向的、有生物學意義的，而且直接回答「這條通路在你的樣本裡是否照教科書運作」。
- 以先驗為約束的迴歸 — 對每個節點，只用它在 KEGG 中的上游 parent 做 ridge/lasso 迴歸，係數就是有向且定量的邊權重。n=6 仍然吃緊，但至少方向是真的。
- 既有的有向通路方法 — 你要的東西已經有成熟工具：
- SPIA — 就是用 KGML 的有向帶符號拓撲做擾動傳播
- CARNIVAL / CausalR — 有向帶符號的因果推論（配 DoRothEA / PROGENy）
- Infomap — 分群方面，它是基於流（flow）的，天生處理有向圖，對「訊號沿通路流動」的語義比 modularity 貼切得多
- 穩健性檢驗 — 不論用哪種權重，n=6 都必須加 bootstrap 或 permutation 的邊穩定度評估，只保留重複出現的邊。或者用 TCGA 之類的大 cohort 估相關結構，再把你的 3 對樣本投影上去。
四、結論
| 面向 | 評價 |
|---|---|
| 工程品質（快取、重試、可重現、繪圖） | 明顯改善，BUGFIX 清單多數是真 bug |
| 「定量」 | 形式上有（|r| 當權重），但 n=6 使其在統計上不可辯護 |
| 「有向」 | 完全未解決，且比第一版的 DiGraph 更退步；subtype 被解析後閒置 |
| 生物學適切性 | 對訊號傳導通路用 mRNA 共表達當邊權重，前提有問題 |
它是一份寫得很好的無向、共表達加權模組偵測程式，不是有向定量通路模型。BUGFIX 註解的自信程度（「9 個問題全修好」）反而容易讓人以為核心方法學問題也一併解決了。
import xml.etree.ElementTree as ET
import networkx as nx
import requests
# 1. 從 KEGG REST API 下載 hsa04390 (Hippo 通路) 的 KGML (XML) 拓撲圖
url = "https://rest.kegg.jp/get/hsa04360/kgml"
print("🌐 正在從 KEGG 伺服器下載 hsa04360 (axon guidance)...")
response = requests.get(url)
if response.status_code != 200:
raise RuntimeError(f"下載失敗，HTTP 狀態碼: {response.status_code}")
root = ET.fromstring(response.text)
# 2. 解析基因節點 (Entries)
entry_map = {}
G = nx.DiGraph()
for entry in root.findall("entry"):
entry_id = entry.get("id")
entry_type = entry.get("type")
if entry_type in ["gene", "group"]:
graphics = entry.find("graphics")
if graphics is not None and "name" in graphics.attrib:
# 取得簡短基因名稱 (如 YAP1, MST1...)
gene_name = (
graphics.get("name").split(",")[0].replace("...", "").strip()
)
else:
gene_name = entry.get("name", entry_id)
entry_map[entry_id] = gene_name
G.add_node(gene_name, entry_id=entry_id, type=entry_type)
# 3. 解析交互作用邊 (Relations)
for relation in root.findall("relation"):
e1 = relation.get("entry1")
e2 = relation.get("entry2")
if e1 in entry_map and e2 in entry_map:
node1 = entry_map[e1]
node2 = entry_map[e2]
G.add_edge(node1, node2)
# 4. 輸出結果
print("✅ 解析成功！")
print(f"📊 Axon guidance 通路基因節點數: {G.number_of_nodes()}")
print(f"📊 基因間調控關係數: {G.number_of_edges()}")
print("前 10 個解析到的基因:", list(G.nodes())[:10])
import xml.etree.ElementTree as ET
from pathlib import Path
import anndata as ad
import igraph as ig
import leidenalg
import networkx as nx
import numpy as np
import pandas as pd
import requests
import scanpy as sc
# ------------------------------------------------------------------
# 1. 讀取數據與建構 AnnData (基因為列，樣本為欄)
# ------------------------------------------------------------------
input_folder = Path("/Users/yshuang/Documents/Python")
file_name = "FC_GSEA.xlsx"
full_path = input_folder / file_name
if not full_path.exists():
raise FileNotFoundError(f"找不到輸入檔案: {full_path}")
xls = pd.ExcelFile(full_path)
df_raw = pd.read_excel(xls, sheet_name="data")
# 自動偵測 Symbol 與 ID 欄位
symbol_col = next(
(
c
for c in ["gene_symb", "gene_symbol", "Gene_symbol", "Gene_Symbol"]
if c in df_raw.columns
),
None,
)
id_col = next(
(
c
for c in ["Transcriptid", "Transcript_ID", "Transcript ID"]
if c in df_raw.columns
),
None,
)
expression_cols = [
c
for c in ["1440AT", "1440MT", "2003AT", "2003MT", "2613AT", "2613MT"]
if c in df_raw.columns
]
# 資料清理與數值化
df_raw = df_raw.dropna(subset=[symbol_col if symbol_col else id_col])
numeric_data = (
df_raw[expression_cols].apply(pd.to_numeric, errors="coerce").fillna(0)
)
# 建立以 Gene Symbol 為主鍵的 AnnData
gene_names = (
df_raw[symbol_col].astype(str).str.strip() if symbol_col else df_raw[id_col]
)
adata = ad.AnnData(
X=numeric_data.values,
obs=pd.DataFrame(index=gene_names),
var=pd.DataFrame(index=expression_cols),
)
# 去除重複基因並進行標準化 (Normalize + Log1p)
adata = adata[~adata.obs_names.duplicated(keep="first")].copy()
sc.pp.filter_cells(adata, min_genes=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
# 擷取標準化表達量矩陣 (DataFrame 格式，用於計算相關性)
df_expr = pd.DataFrame(
adata.X.copy(), index=adata.obs_names, columns=adata.var_names
)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
================================================================================
KEGG Hippo (hsa04390) 拓撲 × 表達量權重 × Leiden 分群 × 可視化 —— 整合修正版
================================================================================
修正原始程式的 9 個問題（詳見 README 區塊），並加入：
* KEGG KGML 本地快取 + 重試 + 離線備援網路
* 「全別名」比對（KGML graphics/@name 的每一個同義字都納入對照表）
* 正確的 scanpy 資料方向（obs = 樣本、var = 基因）
* 可重現的 Leiden 分群（RBConfiguration + resolution + 固定順序）
* log2FC（MT vs AT 配對）與模組層級統計
* 5 張出版級圖表 + 多工作表 Excel
執行：
python hippo_kegg_leiden.py # 完整線上流程
python hippo_kegg_leiden.py --offline # 不連 KEGG，用內建拓撲
python hippo_kegg_leiden.py --demo # 自動產生模擬資料跑通全流程
python hippo_kegg_leiden.py --resolution 1.2 --min-abs-r 0.4
依賴：
pip install networkx python-igraph leidenalg scanpy anndata \
pandas numpy scipy matplotlib openpyxl requests
================================================================================
"""
import argparse
import re
import sys
import time
import warnings
import xml.etree.ElementTree as ET
from dataclasses import dataclass, field
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
from matplotlib.colors import LinearSegmentedColormap, Normalize, TwoSlopeNorm
from matplotlib.lines import Line2D
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)
# ==============================================================================
# ❖ 原始程式的問題清單（修正對照）
# ==============================================================================
BUGFIX_NOTES = """
[1] scanpy 資料方向錯誤 —— 原本 obs=基因、var=樣本，導致 normalize_total()
是「對每個基因跨樣本」正規化，等於把樣本間的表達差異抹平。
修正：AnnData 建成 obs=樣本、var=基因（scanpy 的標準語意），
正規化後再轉置成基因×樣本供相關性計算。
[2] ModularityVertexPartition 不接受 resolution_parameter —— 會直接 TypeError。
修正：改用 RBConfigurationVertexPartition（真正支援 resolution）。
[3] list(set(...).intersection(...)) 順序隨 PYTHONHASHSEED 改變 ——
seed=42 也救不了，每次跑出的模組編號與切分都可能不同。
修正：一律 sorted()，並固定 n_iterations=-1 跑到收斂。
[4] type="group" 的 entry 沒有 graphics/@name，fallback 會抓到 "undefined"，
製造出一個假基因節點並吸走大量邊。
修正：group 展開成其 <component> 成員，不自成節點。
[5] 只取 graphics/@name 的第一個別名，丟掉其餘同義字 ——
資料若寫 CTGF、MST1、TAZ，與 KEGG 的 CCN2、STK4、WWTR1 對不上，
交集會莫名其妙變得很少。
修正：建立「全別名 → 官方 symbol」對照表再比對。
[6] common_genes.index(u) 在迴圈內是 O(n) 查找；且 relations 未過濾 type，
把 maplink（通路對通路的連結）也當成基因交互作用。
修正：改用 dict 索引；只保留 PPrel/GErel/PCrel 且 subtype 有意義的關係。
[7] weights 可能為空 list —— g_ig.es["weight"] = [] 在有邊時會炸掉；
孤立節點也會變成一堆單基因「模組」污染結果。
修正：先過濾零變異基因、套用 |r| 門檻，再移除孤立節點並明確回報。
[8] corr 用 .abs() 後就丟掉了方向 —— 正/負相關被當成一樣。
修正：權重用 |r|（Leiden 需要非負），但把 r 的正負號保留在邊屬性與圖上。
[9] 每次執行都重打 KEGG 伺服器、且沒有錯誤重試。
修正：本地 XML 快取 + 3 次退避重試 + 離線備援拓撲。
"""
# ==============================================================================
# ❖ 設定
# ==============================================================================
@dataclass
class Config:
input_folder: Path = Path("/Users/yshuang/Documents/Python")
file_name: str = "FC_GSEA.xlsx"
sheet_name: str = "data"
pathway_id: str = "hsa04360"
sample_cols: list[str] = field(default_factory=lambda: [
"1440AT", "1440MT", "2003AT", "2003MT", "2613AT", "2613MT"])
symbol_candidates: tuple = ("gene_symb", "gene_symbol", "Gene_symbol",
"Gene_Symbol", "Symbol", "SYMBOL")
id_candidates: tuple = ("Transcriptid", "Transcript_ID", "Transcript ID")
resolution: float = 1.0 # 越大 → 模組越多越細
min_abs_r: float = 0.30 # 邊的 |Pearson r| 門檻
corr_method: str = "pearson" # pearson / spearman
seed: int = 42
dpi: int = 300
@property
def kgml_url(self) -> str:
return f"https://rest.kegg.jp/get/{self.pathway_id}/kgml"
@property
def out_dir(self) -> Path:
return self.input_folder / f"{self.pathway_id}_leiden_output"
# ------------------------------------------------------------------ 調色盤
# 來源：已驗證的類別色板。網路圖屬於「全配對比較」情境，
# 前 4 槽通過全配對 CVD/正常視覺門檻（第 5 槽起會失敗），
# 因此超過 4 個模組時折疊為 Other，並額外輸出「模組分面圖」保證身份可辨。
CAT4 = ["#2a78d6", "#eb6834", "#1baf7a", "#4a3aa7"] # blue / orange / aqua / violet
MARKERS = ["o", "s", "^", "D", "v", "P", "X", "h"] # 次要編碼：形狀
GRAY_OTHER = "#898781"
INK, INK2, MUTED = "#0b0b0b", "#52514e", "#898781"
GRID, SURFACE = "#e1e0d9", "#fcfcfb"
DIVERGING = LinearSegmentedColormap.from_list(
"blue_gray_red", ["#184f95", "#2a78d6", "#f0efec", "#e34948", "#a02020"])
plt.rcParams.update({
"font.family": ["DejaVu Sans"], "font.size": 9,
"axes.edgecolor": "#c3c2b7", "axes.labelcolor": INK2,
"xtick.color": MUTED, "ytick.color": MUTED,
"figure.facecolor": SURFACE, "axes.facecolor": SURFACE,
"savefig.facecolor": SURFACE, "axes.titlecolor": INK,
})
# ==============================================================================
# ❖ 1. KEGG KGML —— 下載（含快取/重試）、解析、離線備援
# ==============================================================================
MEANINGFUL_SUBTYPES = {
"activation", "inhibition", "expression", "repression",
"phosphorylation", "dephosphorylation", "ubiquitination",
"binding/association", "dissociation", "indirect effect", "state change",
}
def fetch_kgml(cfg: Config, retries: int = 3) -> str:
"""下載 KGML，帶本地快取與指數退避重試。"""
cache = cfg.out_dir / f"{cfg.pathway_id}.xml"
cache.parent.mkdir(parents=True, exist_ok=True)
if cache.exists() and cache.stat().st_size > 1000:
print(f"📁 使用本地快取: {cache.name}")
return cache.read_text(encoding="utf-8")
import requests
for attempt in range(1, retries + 1):
try:
print(f"🌐 下載 KEGG {cfg.pathway_id} KGML (第 {attempt} 次)...")
r = requests.get(cfg.kgml_url, timeout=30)
r.raise_for_status()
if "<pathway" not in r.text:
raise RuntimeError("回應不是合法 KGML")
cache.write_text(r.text, encoding="utf-8")
return r.text
except Exception as exc:
if attempt == retries:
raise
print(f" ⚠️ 失敗（{exc}），{2 ** attempt}s 後重試...")
time.sleep(2 ** attempt)
raise RuntimeError("unreachable")
def parse_kgml(xml_text: str) -> tuple[nx.Graph, dict[str, str]]:
"""
回傳 (無向圖, alias2symbol)
* group entry 展開成 component，不自成節點 → 修正 [4]
* graphics/@name 的每個別名都進 alias2symbol → 修正 [5]
* 只保留有意義的 relation type / subtype → 修正 [6]
"""
root = ET.fromstring(xml_text)
id2symbol: dict[str, str] = {}
id2group: dict[str, list[str]] = {}
alias2symbol: dict[str, str] = {}
id2keggids: dict[str, str] = {}
for entry in root.findall("entry"):
eid, etype = entry.get("id"), entry.get("type")
if etype == "gene":
g = entry.find("graphics")
if g is None or not g.get("name"):
continue
aliases = [a.strip().rstrip(".").strip()
for a in g.get("name").split(",")]
aliases = [a for a in aliases if a]
if not aliases:
continue
symbol = aliases[0]
id2symbol[eid] = symbol
id2keggids[eid] = entry.get("name", "")
for a in aliases: # 全別名對照
alias2symbol.setdefault(a.upper(), symbol)
elif etype == "group":
id2group[eid] = [c.get("id") for c in entry.findall("component")]
def expand(eid: str) -> list[str]:
if eid in id2symbol:
return [id2symbol[eid]]
return [id2symbol[m] for m in id2group.get(eid, []) if m in id2symbol]
G = nx.Graph()
for eid, sym in id2symbol.items():
G.add_node(sym, kegg_ids=id2keggids.get(eid, ""))
for rel in root.findall("relation"):
if rel.get("type") not in ("PPrel", "GErel", "PCrel"):
continue
subs = sorted({s.get("name") for s in rel.findall("subtype")}
& MEANINGFUL_SUBTYPES)
if not subs:
continue
for a in expand(rel.get("entry1")):
for b in expand(rel.get("entry2")):
if a == b:
continue
if G.has_edge(a, b):
G[a][b]["subtypes"] = sorted(set(G[a][b]["subtypes"]) | set(subs))
else:
G.add_edge(a, b, subtypes=subs)
return G, alias2symbol
import networkx as nx
# --------------------------------------------------- 離線備援拓撲（60 基因）
OFFLINE_AXON_GUIDANCE_EDGES = [
# Netrin 配體與受體 (Attraction / Repulsion)
("NTN1", "DCC"), ("NTN1", "NEO1"), ("NTN1", "UNC5A"), ("NTN1", "UNC5B"), ("NTN1", "UNC5C"), ("NTN1", "UNC5D"),
("DCC", "UNC5A"), ("DCC", "UNC5B"),
# DCC 下游吸引信號 (Cytoskeletal reorganisation)
("DCC", "PTK2"), ("DCC", "FYN"), ("DCC", "SRC"), ("DCC", "NCK1"),
("PTK2", "RAC1"), ("PTK2", "CDC42"), ("FYN", "ABL1"), ("NCK1", "PAK1"),
# UNC5 下游排斥信號
("UNC5B", "ABL1"), ("UNC5B", "RHOA"), ("UNC5B", "PTPN11"),
# Slit / Robo 配體受體與排斥信號
("SLIT1", "ROBO1"), ("SLIT2", "ROBO1"), ("SLIT2", "ROBO2"), ("SLIT3", "ROBO1"), ("SLIT3", "ROBO2"),
("ROBO1", "SRGAP1"), ("ROBO1", "ABL1"), ("ROBO1", "ENAH"),
("SRGAP1", "CDC42"), ("SRGAP1", "RAC1"), ("ABL1", "CFL1"), ("ENAH", "ACTR2"),
# Ephrin-A / EphA 信號 (Forward & Reverse)
("EFNA1", "EPHA2"), ("EFNA1", "EPHA4"), ("EFNA5", "EPHA4"), ("EFNA5", "EPHA7"),
("EPHA4", "ADAM10"), ("EPHA4", "VAV2"), ("EPHA2", "RHOA"), ("EPHA4", "FGFR1"),
("VAV2", "RAC1"), ("VAV2", "RHOA"),
# Ephrin-B / EphB 信號
("EFNB1", "EPHB1"), ("EFNB2", "EPHB1"), ("EFNB2", "EPHB2"), ("EFNB3", "EPHB3"),
("EPHB2", "KALRN"), ("EPHB2", "GRIP1"), ("EPHB2", "SRC"),
("KALRN", "RAC1"), ("GRIP1", "GRIA1"),
# Semaphorin / Neuropilin & Plexin 複合體
("SEMA3A", "NRP1"), ("NRP1", "PLXNA1"), ("SEMA3A", "PLXNA1"), ("SEMA3A", "PLXNA2"),
("SEMA4D", "PLXNB1"), ("SEMA6D", "PLXNA1"),
("PLXNA1", "FARP2"), ("PLXNA1", "RND1"), ("PLXNA1", "FES"),
("PLXNB1", "ARHGEF11"), ("PLXNB1", "MET"),
("FARP2", "RAC1"), ("RND1", "PLXNA1"), ("ARHGEF11", "RHOA"),
# 下游 Rho GTPases 與細胞骨架動力學核心 (Downstream Convergence)
("RHOA", "ROCK1"), ("RHOA", "ROCK2"), ("ROCK1", "LIMK1"), ("ROCK2", "LIMK2"),
("RAC1", "PAK1"), ("CDC42", "PAK1"), ("PAK1", "LIMK1"),
("LIMK1", "CFL1"), ("LIMK2", "CFL1"),
("RAC1", "WASF1"), ("CDC42", "WASL"),
("WASL", "ACTR2"), ("WASL", "ACTR3"), ("WASF1", "ACTR2"),
# 軸突生長與微管調控 (Microtubule & Microfilament dynamics)
("AKT1", "GSK3B"), ("GSK3B", "DPYSL2"), ("GSK3B", "MAPT"),
("DPYSL2", "TUBA1A"), ("MAPT", "TUBB"),
]
# 軸突引導通路常見別名/舊名映射至 KEGG 官方 HGNC Symbol
OFFLINE_AXON_GUIDANCE_ALIASES = {
# 常用別名 / 舊稱
"FAK": "PTK2", "PTK2": "PTK2",
"CRMP2": "DPYSL2", "CRMP-2": "DPYSL2", "DPYSL2": "DPYSL2",
"MENA": "ENAH", "ENAH": "ENAH",
"NETRIN1": "NTN1", "NETRIN-1": "NTN1", "NTN1": "NTN1",
"NETRIN2": "NTN2", "NTN2": "NTN2",
"PDZ-RHOGEF": "ARHGEF11", "ARHGEF11": "ARHGEF11",
"N-WASP": "WASL", "WASL": "WASL",
"WASP": "WAS", "WAS": "WAS",
"WAVE1": "WASF1", "WASF1": "WASF1",
"ARP2": "ACTR2", "ACTR2": "ACTR2",
"ARP3": "ACTR3", "ACTR3": "ACTR3",
"COFILIN": "CFL1", "CFL1": "CFL1",
"TAU": "MAPT", "MAPT": "MAPT",
"SHP2": "PTPN11", "PTPN11": "PTPN11",
"EPHRIN-A1": "EFNA1", "EFNA1": "EFNA1",
"EPHRIN-B1": "EFNB1", "EFNB1": "EFNB1",
}
def build_kegg_network(cfg, offline: bool = False) -> tuple[nx.Graph, dict]:
"""建立 KEGG Axon Guidance 網路圖與別名對照字典。"""
if offline:
print("📦 離線模式：使用內建 Axon Guidance 拓撲")
G = nx.Graph()
G.add_edges_from((a, b, {"subtypes": ["curated"]})
for a, b in OFFLINE_AXON_GUIDANCE_EDGES)
alias = {g.upper(): g for g in G.nodes()}
alias.update(OFFLINE_AXON_GUIDANCE_ALIASES)
return G, alias
# 線上模式解析 KGML
G, alias = parse_kgml(fetch_kgml(cfg))
for k, v in OFFLINE_AXON_GUIDANCE_ALIASES.items(): # 補上 KGML 可能沒收錄的舊名
if v in G:
alias.setdefault(k, v)
return G, alias
# ==============================================================================
# ❖ 2. 表達資料 —— 正確方向的 AnnData 正規化 → 修正 [1]
# ==============================================================================
def load_expression(cfg: Config) -> tuple[pd.DataFrame, pd.DataFrame]:
"""回傳 (log 正規化後的 基因×樣本, 原始 基因×樣本)"""
path = cfg.input_folder / cfg.file_name
if not path.exists():
raise FileNotFoundError(f"找不到輸入檔案: {path}")
df_raw = pd.read_excel(pd.ExcelFile(path), sheet_name=cfg.sheet_name)
symbol_col = next((c for c in cfg.symbol_candidates if c in df_raw.columns), None)
id_col = next((c for c in cfg.id_candidates if c in df_raw.columns), None)
if symbol_col is None and id_col is None:
raise ValueError(f"找不到基因名稱欄位，現有欄位: {list(df_raw.columns)}")
key_col = symbol_col or id_col
cols = [c for c in cfg.sample_cols if c in df_raw.columns]
if len(cols) < 3:
raise ValueError(f"可用樣本欄位不足（{cols}），至少需要 3 個")
df_raw = df_raw.dropna(subset=[key_col]).copy()
df_raw["__gene__"] = df_raw[key_col].astype(str).str.strip().str.upper()
mat = df_raw[cols].apply(pd.to_numeric, errors="coerce").fillna(0.0)
mat.index = df_raw["__gene__"].values
# 重複 symbol：取平均（比 keep-first 穩健）
n_dup = mat.index.duplicated().sum()
if n_dup:
print(f"🔁 合併 {n_dup} 個重複 gene symbol（取平均）")
mat = mat.groupby(level=0).mean()
mat = mat.loc[mat.sum(axis=1) > 0] # 移除全零基因
# AnnData：obs = 樣本、var = 基因（scanpy 標準語意）
import anndata as ad
import scanpy as sc
adata = ad.AnnData(
X=mat.T.values.astype(np.float64),
obs=pd.DataFrame(index=pd.Index(cols, name="sample")),
var=pd.DataFrame(index=pd.Index(mat.index, name="gene")),
)
sc.pp.normalize_total(adata, target_sum=1e4) # 每個「樣本」總量對齊 ✔
sc.pp.log1p(adata)
df_norm = pd.DataFrame(np.asarray(adata.X).T, index=mat.index, columns=cols)
print(f"📥 表達矩陣: {df_norm.shape[0]} 基因 × {df_norm.shape[1]} 樣本")
return df_norm, mat
def compute_log2fc(mat_raw: pd.DataFrame) -> pd.DataFrame:
"""依欄名 <case><AT|MT> 自動配對，計算 MT vs AT 的 log2FC 與配對 p 值。"""
pat = re.compile(r"^(?P<case>.+?)(?P<grp>AT|MT)$", re.IGNORECASE)
at, mt = [], []
for c in mat_raw.columns:
m = pat.match(str(c))
if m:
(mt if m.group("grp").upper() == "MT" else at).append(c)
if not at or not mt:
return pd.DataFrame(index=mat_raw.index,
columns=["log2FC", "p_paired"], dtype=float)
cpm = mat_raw.div(mat_raw.sum(axis=0).replace(0, np.nan), axis=1) * 1e6
lg = np.log2(cpm + 1)
log2fc = lg[mt].mean(axis=1) - lg[at].mean(axis=1)
p = pd.Series(np.nan, index=mat_raw.index)
if len(at) == len(mt) and len(at) >= 2:
try:
from scipy import stats
t, pv = stats.ttest_rel(lg[mt].values, lg[at].values,
axis=1, nan_policy="omit")
p = pd.Series(np.asarray(pv, dtype=float), index=mat_raw.index)
except Exception:
pass
return pd.DataFrame({"log2FC": log2fc, "p_paired": p})
# ==============================================================================
# ❖ 3. 比對 + 定量權重 → 修正 [3][5][7][8]
# ==============================================================================
def match_genes(kegg_net: nx.Graph, alias2symbol: dict,
df_expr: pd.DataFrame) -> tuple[list[str], pd.DataFrame, list[str]]:
"""用全別名表把資料的基因名映射到 KEGG 官方 symbol。"""
rows, hits = {}, []
for g in df_expr.index: # index 已是大寫
sym = alias2symbol.get(g)
if sym and sym in kegg_net:
hits.append((g, sym))
rows.setdefault(sym, []).append(g)
if not rows:
return [], pd.DataFrame(), sorted(kegg_net.nodes())
mapped = pd.DataFrame(
{sym: df_expr.loc[src].mean(axis=0) for sym, src in rows.items()}).T
mapped.index.name = "gene"
common = sorted(mapped.index) # 固定順序 → 可重現 ✔
missing = sorted(set(kegg_net.nodes()) - set(common))
alias_used = [f"{s}→{t}" for s, t in hits if s != t]
if alias_used:
print(f"🔗 別名比對救回 {len(alias_used)} 個基因，例如: "
f"{', '.join(alias_used[:6])}")
return common, mapped.loc[common], missing
def build_weighted_graph(kegg_net: nx.Graph, expr: pd.DataFrame,
cfg: Config) -> tuple[nx.Graph, pd.DataFrame]:
sub = kegg_net.subgraph(expr.index).copy()
var_ok = expr.loc[expr.std(axis=1) > 1e-9].index
dropped = sorted(set(expr.index) - set(var_ok))
if dropped:
print(f"⚠️ 移除 {len(dropped)} 個零變異基因（無法算相關性）")
corr = expr.loc[var_ok].T.corr(method=cfg.corr_method)
W = nx.Graph()
W.add_nodes_from(var_ok)
kept = []
for u, v, d in sub.edges(data=True):
if u not in corr.index or v not in corr.index:
continue
r = float(corr.loc[u, v])
if not np.isfinite(r) or abs(r) < cfg.min_abs_r:
continue
W.add_edge(u, v, weight=abs(r), r=r,
sign="+" if r >= 0 else "-",
subtypes=";".join(d.get("subtypes", [])))
kept.append({"source": u, "target": v, "pearson_r": r,
"abs_r": abs(r), "kegg_subtypes": ";".join(d.get("subtypes", []))})
iso = [n for n in W.nodes() if W.degree(n) == 0]
W.remove_nodes_from(iso) # 避免一堆單基因假模組 ✔
print(f"🔗 KEGG 子圖 {sub.number_of_edges()} 邊 → 過門檻 "
f"(|r|≥{cfg.min_abs_r}) 保留 {W.number_of_edges()} 邊；"
f"移除 {len(iso)} 個孤立節點")
if W.number_of_edges() == 0:
raise ValueError("🚨 沒有任何邊通過 |r| 門檻，請調低 --min-abs-r")
return W, pd.DataFrame(kept)
# ==============================================================================
# ❖ 4. 加權 Leiden → 修正 [2][3]
# ==============================================================================
def weighted_leiden(W: nx.Graph, cfg: Config) -> tuple[dict[str, int], float]:
import igraph as ig
import leidenalg
nodes = sorted(W.nodes()) # 固定順序 ✔
idx = {g: i for i, g in enumerate(nodes)} # O(1) 查找 ✔
g = ig.Graph(n=len(nodes))
g.vs["name"] = nodes
edges = [(idx[u], idx[v]) for u, v in W.edges()]
g.add_edges(edges)
g.es["weight"] = [W[u][v]["weight"] for u, v in W.edges()]
part = leidenalg.find_partition(
g,
leidenalg.RBConfigurationVertexPartition, # 支援 resolution ✔
weights="weight",
resolution_parameter=cfg.resolution,
n_iterations=-1, # 跑到收斂
seed=cfg.seed,
)
membership = {nodes[i]: int(c) for i, c in enumerate(part.membership)}
# 依模組大小重新編號，讓 Module 0 永遠是最大的模組（跨次執行可讀）
order = pd.Series(membership).value_counts().index.tolist()
remap = {old: new for new, old in enumerate(order)}
membership = {g_: remap[c] for g_, c in membership.items()}
q = nx.community.modularity(
W, [{g_ for g_, c in membership.items() if c == k}
for k in sorted(set(membership.values()))], weight="weight")
print(f"🧬 Leiden(resolution={cfg.resolution}) → "
f"{len(set(membership.values()))} 個模組，加權模組度 Q={q:.4f}")
return membership, q
# ==============================================================================
# ❖ 5. 可視化
# ==============================================================================
def module_color(m: int) -> str:
"""顏色 × 形狀成對編碼：4 色 × 8 形 = 32 組唯一組合，且色數不超過全配對安全上限。"""
return CAT4[m % len(CAT4)]
def module_marker(m: int) -> str:
return MARKERS[(m // len(CAT4)) % len(MARKERS)]
def _layout(W: nx.Graph, seed: int) -> dict:
"""
連通分量各自排版後再打包 —— 直接對整張圖跑 spring_layout 時，
分離的小分量會被推到極遠處，害主體擠成一團。
"""
comps = sorted(nx.connected_components(W), key=len, reverse=True)
boxes, sizes = [], []
for c in comps:
sub = W.subgraph(c)
if len(c) == 1:
p = {next(iter(c)): np.array([0.0, 0.0])}
elif len(c) <= 12:
p = nx.kamada_kawai_layout(sub)
else:
p = nx.spring_layout(sub, seed=seed, iterations=600,
k=2.2 / np.sqrt(len(c)), weight="weight")
arr = np.array([p[n] for n in sub.nodes()], dtype=float)
if len(arr) > 1:
arr -= arr.mean(0)
span = np.abs(arr).max() or 1.0
arr /= span
boxes.append(dict(zip(sub.nodes(), arr)))
sizes.append(np.sqrt(len(c)))
# 主分量占左側大格，其餘依序排入右側欄
pos, main_w = {}, 1.0
for n, xy in boxes[0].items():
pos[n] = xy * main_w
y_cursor, x_off = 1.0, main_w + 0.55
col_w = 0.30
for box, s in zip(boxes[1:], sizes[1:]):
h = max(0.16, 0.10 * s)
for n, xy in box.items():
pos[n] = np.array([x_off + xy[0] * col_w, y_cursor - h + xy[1] * h])
y_cursor -= 2.4 * h + 0.16
if y_cursor < -1.05:
y_cursor, x_off = 1.0, x_off + 2 * col_w + 0.30
return pos
def _ring(r: float, n: int = 8, start: float = 90.0):
return [(r * np.cos(np.deg2rad(start + 360 * i / n)),
r * np.sin(np.deg2rad(start + 360 * i / n))) for i in range(n)]
# 由近而遠的候選位置；超過 18pt 時自動補一條細引線
_LABEL_OFFSETS = ([(0, 11), (0, -12)] + _ring(13) + _ring(20, start=67.5)
+ _ring(28) + _ring(37, start=67.5) + _ring(47))
def _place_labels(ax, fig, pos, names, fontsize=6.4, color=INK,
node_radius_pt=7.0):
"""
貪婪式標籤避讓：逐一嘗試候選偏移，選第一個不與「已放置標籤或節點」重疊的位置。
注意 fontsize / offset 都是 point，碰撞判斷在 pixel，必須換算（dpi/72）。
"""
fig.canvas.draw()
k = fig.dpi / 72.0
placed: list[tuple[float, float, float, float]] = []
for n in pos: # 節點本身先佔位
x, y = ax.transData.transform(pos[n])
r = node_radius_pt * k
placed.append((x - r, y - r, x + r, y + r))
for n in sorted(names, key=lambda s: -len(s)):
x, y = ax.transData.transform(pos[n])
w = len(n) * fontsize * 0.62 * k
h = fontsize * 1.25 * k
for dx, dy in _LABEL_OFFSETS:
cx, cy = x + dx * k, y + dy * k
box = (cx - w / 2, cy - h / 2, cx + w / 2, cy + h / 2)
if all(box[2] < p[0] or box[0] > p[2] or
box[3] < p[1] or box[1] > p[3] for p in placed):
placed.append(box)
far = (dx ** 2 + dy ** 2) ** 0.5 > 18
ax.annotate(n, pos[n], fontsize=fontsize, color=color,
ha="center", va="center", zorder=6,
xytext=(dx, dy), textcoords="offset points",
arrowprops=dict(arrowstyle="-", lw=0.45,
color="#c3c2b7", shrinkA=1,
shrinkB=5) if far else None)
break
def _limits(pos, pad=0.10):
xs = [p[0] for p in pos.values()]
ys = [p[1] for p in pos.values()]
dx = (max(xs) - min(xs)) or 1.0
dy = (max(ys) - min(ys)) or 1.0
return ((min(xs) - pad * dx, max(xs) + pad * dx),
(min(ys) - pad * dy, max(ys) + pad * dy))
def _frame(ax, pos, pad=0.10):
(x0, x1), (y0, y1) = _limits(pos, pad)
ax.set_xlim(x0, x1); ax.set_ylim(y0, y1)
ax.set_aspect("equal")
ax.axis("off")
def _network_canvas(pos, ax_w=10.5, right_in=2.9, top_in=1.05, pad=0.10):
"""依版面實際長寬比開圖，避免 equal aspect 造成上下大片留白。"""
(x0, x1), (y0, y1) = _limits(pos, pad)
ax_h = max(3.0, ax_w * (y1 - y0) / (x1 - x0))
fig = plt.figure(figsize=(ax_w + right_in, ax_h + top_in))
fig.subplots_adjust(left=0.02 / (ax_w + right_in),
right=(ax_w) / (ax_w + right_in),
bottom=0.04 / (ax_h + top_in),
top=(ax_h) / (ax_h + top_in))
ax = fig.add_subplot(111)
ax.set_position([0.015, 0.02, ax_w / (ax_w + right_in),
ax_h / (ax_h + top_in)])
return fig, ax
def _titles(fig, title, subtitle):
h = fig.get_size_inches()[1]
fig.text(0.012, 1 - 0.30 / h, title, fontsize=13.5, color=INK,
ha="left", va="top")
fig.text(0.012, 1 - 0.62 / h, subtitle, fontsize=8.8, color=INK2,
ha="left", va="top")
def fig_network_modules(W, pos, membership, cfg, out: Path):
"""圖1：模組著色網路（顏色+形狀+直接標籤三重編碼）"""
fig, ax = _network_canvas(pos, ax_w=10.5, right_in=2.7)
for u, v, d in W.edges(data=True):
ax.plot(*zip(pos[u], pos[v]), lw=0.7 + 2.8 * d["weight"],
color="#e34948" if d["r"] < 0 else "#c3c2b7",
alpha=0.7 if d["r"] < 0 else 0.6, zorder=1,
solid_capstyle="round")
deg = dict(W.degree())
mods = sorted(set(membership.values()))
for m in mods:
ns = [n for n in W.nodes() if membership[n] == m]
ax.scatter([pos[n][0] for n in ns], [pos[n][1] for n in ns],
s=[54 + 30 * deg[n] for n in ns],
c=module_color(m), marker=module_marker(m),
edgecolors=SURFACE, linewidths=1.8, zorder=3,
label=f"Module {m} n={len(ns)}")
# 模組編號直接標在質心 → 身份不靠顏色單獨承擔
cx, cy = np.mean([pos[n] for n in ns], axis=0)
ax.annotate(f"M{m}", (cx, cy), fontsize=12, fontweight="bold",
color=module_color(m), ha="center", va="center", zorder=2,
alpha=0.30)
_frame(ax, pos)
_place_labels(ax, fig, pos, list(W.nodes()))
_titles(fig, f"Axon guidance pathway (KEGG {cfg.pathway_id}) — "
"expression-weighted Leiden modules",
"node size = degree · edge width = |Pearson r| · "
"red edge = negative correlation · marker shape + M-label = module identity")
ax.legend(loc="upper left", bbox_to_anchor=(1.01, 1.0), frameon=False,
fontsize=9, labelcolor=INK2, handletextpad=0.7, labelspacing=0.85)
fig.savefig(out, dpi=cfg.dpi, bbox_inches="tight")
plt.close(fig)
def fig_module_facets(W, pos, membership, cfg, out: Path):
"""圖2：模組分面小倍數 —— 模組數 >4 時保證身份不靠顏色"""
mods = sorted(set(membership.values()))
ncol = min(4, len(mods))
nrow = int(np.ceil(len(mods) / ncol))
fig, axes = plt.subplots(nrow, ncol, figsize=(3.7 * ncol, 3.7 * nrow),
squeeze=False)
for ax in axes.ravel():
ax.axis("off")
for k, m in enumerate(mods):
ax = axes[k // ncol][k % ncol]
ns = [n for n in W.nodes() if membership[n] == m]
for u, v, d in W.edges(data=True):
hot = membership[u] == m and membership[v] == m
ax.plot(*zip(pos[u], pos[v]),
lw=2.0 * d["weight"] if hot else 0.35,
color=CAT4[0] if hot else "#e1e0d9",
alpha=0.8 if hot else 0.55, zorder=2 if hot else 1)
others = [n for n in W.nodes() if membership[n] != m]
ax.scatter([pos[n][0] for n in others], [pos[n][1] for n in others],
s=11, c="#e1e0d9", zorder=2)
ax.scatter([pos[n][0] for n in ns], [pos[n][1] for n in ns],
s=46, c=CAT4[0], edgecolors=SURFACE, linewidths=1.2, zorder=3)
_frame(ax, pos, pad=0.10)
_place_labels(ax, fig, {n: pos[n] for n in ns}, ns, fontsize=5.6)
ax.set_title(f"Module {m} · {len(ns)} genes", fontsize=10,
color=INK, loc="left", pad=6)
fig.suptitle("Leiden modules, one panel each — identity without relying on hue",
fontsize=12.5, x=0.012, ha="left", color=INK)
fig.tight_layout(rect=(0, 0, 1, 0.955))
fig.savefig(out, dpi=cfg.dpi, bbox_inches="tight")
plt.close(fig)
def fig_network_log2fc(W, pos, membership, node_tbl, cfg, out: Path):
"""圖3：同版面，節點改以 log2FC 發散配色（極性）"""
fc = node_tbl.set_index("Gene_Symbol")["log2FC"].reindex(list(W.nodes()))
if fc.isna().all():
return False
lim = float(np.nanpercentile(np.abs(fc.values), 97)) or 1.0
norm = TwoSlopeNorm(vmin=-lim, vcenter=0.0, vmax=lim)
fig, ax = _network_canvas(pos, ax_w=10.5, right_in=1.9)
for u, v, d in W.edges(data=True):
ax.plot(*zip(pos[u], pos[v]), lw=0.5 + 2.2 * d["weight"],
color="#e1e0d9", zorder=1, solid_capstyle="round")
ns = list(W.nodes())
sc_ = ax.scatter([pos[n][0] for n in ns], [pos[n][1] for n in ns],
s=[54 + 30 * W.degree(n) for n in ns],
c=fc.values, cmap=DIVERGING, norm=norm,
edgecolors=SURFACE, linewidths=1.6, zorder=3)
_frame(ax, pos)
_place_labels(ax, fig, pos, ns)
cb = fig.colorbar(sc_, ax=ax, fraction=0.022, pad=0.015,
shrink=0.52, aspect=26)
cb.set_label("log2FC (MT vs AT)", color=INK2, fontsize=9)
cb.outline.set_visible(False)
_titles(fig, "Differential expression mapped onto the Axon guidance topology",
"blue = down in tumour · gray ≈ unchanged · red = up in tumour · "
"node size = degree · exact values in the Genes_Modules sheet")
fig.savefig(out, dpi=cfg.dpi, bbox_inches="tight")
plt.close(fig)
return True
def fig_corr_heatmap(expr, membership, cfg, out: Path):
"""圖4：依模組排序的基因-基因相關性熱圖（發散）"""
genes = sorted(membership, key=lambda g: (membership[g], g))
corr = expr.loc[genes].T.corr(method=cfg.corr_method)
n = len(genes)
fig, ax = plt.subplots(figsize=(max(6.5, n * 0.20), max(5.6, n * 0.20)))
im = ax.imshow(corr.values, cmap=DIVERGING, vmin=-1, vmax=1)
ax.set_xticks(range(n)); ax.set_yticks(range(n))
ax.set_xticklabels(genes, rotation=90, fontsize=5.6)
ax.set_yticklabels(genes, fontsize=5.6)
bounds, prev = [], None
for i, g in enumerate(genes):
if membership[g] != prev:
bounds.append(i); prev = membership[g]
for b in bounds[1:]:
ax.axhline(b - .5, color=INK, lw=1.1)
ax.axvline(b - .5, color=INK, lw=1.1)
cb = fig.colorbar(im, ax=ax, fraction=0.030, pad=0.02)
cb.set_label(f"{cfg.corr_method.title()} r", color=INK2, fontsize=9)
cb.outline.set_visible(False)
ax.set_title("Gene–gene correlation, ordered by Leiden module",
fontsize=12, pad=12, loc="left")
for s in ax.spines.values():
s.set_visible(False)
ax.tick_params(length=0)
fig.tight_layout()
fig.savefig(out, dpi=cfg.dpi, bbox_inches="tight")
plt.close(fig)
def fig_hubs_and_sizes(node_tbl, membership, cfg, out: Path):
"""圖5：模組規模 + 樞紐基因排名（單一序列，直接標值，無圖例）"""
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.2, 4.9),
gridspec_kw={"width_ratios": [1, 1.35]})
sizes = pd.Series(membership).value_counts().sort_index()
ax1.bar(sizes.index.astype(str), sizes.values, color=CAT4[0], width=0.62)
for x, v in zip(range(len(sizes)), sizes.values):
ax1.annotate(str(v), (x, v), ha="center", va="bottom",
fontsize=9, color=INK2, xytext=(0, 3),
textcoords="offset points")
ax1.set_title("Genes per Leiden module", fontsize=11.5, loc="left", pad=10)
ax1.set_xlabel("module"); ax1.set_ylabel("genes")
ax1.spines[["top", "right"]].set_visible(False)
ax1.grid(axis="y", color=GRID, lw=0.7); ax1.set_axisbelow(True)
top = node_tbl.nlargest(12, "betweenness").iloc[::-1]
ax2.barh(top["Gene_Symbol"], top["betweenness"], color=CAT4[0], height=0.66)
for y, v in enumerate(top["betweenness"]):
ax2.annotate(f"{v:.3f}", (v, y), va="center", fontsize=8.4,
color=INK2, xytext=(4, 0), textcoords="offset points")
ax2.set_title("Top bottleneck genes (betweenness centrality)",
fontsize=11.5, loc="left", pad=10)
ax2.set_xlabel("betweenness")
ax2.spines[["top", "right"]].set_visible(False)
ax2.grid(axis="x", color=GRID, lw=0.7); ax2.set_axisbelow(True)
ax2.margins(x=0.16)
fig.tight_layout()
fig.savefig(out, dpi=cfg.dpi, bbox_inches="tight")
plt.close(fig)
# ==============================================================================
# ❖ 6. 主流程
# ==============================================================================
def make_demo_data(cfg: Config) -> None:
"""--demo：合成一份帶模組結構的表達量 Excel，用來驗證整條流程。"""
rng = np.random.default_rng(cfg.seed)
genes = sorted({g for e in OFFLINE_AXON_GUIDANCE_EDGES for g in e})
genes = [{"CCN2": "CTGF", "CCN1": "CYR61", "STK4": "MST1",
"WWTR1": "TAZ"}.get(g, g) for g in genes] # 故意混入舊名
n = len(genes)
base = rng.lognormal(4.2, 1.05, n)
blocks = rng.integers(0, 4, n) # 4 個潛在模組
prog = rng.normal(0, 1, (4, 6))
X = np.empty((n, 6))
for i in range(n):
X[i] = base[i] * (1 + 0.55 * prog[blocks[i]] + rng.normal(0, .18, 6))
X = np.clip(X, 0, None)
X[:, [1, 3, 5]] *= rng.lognormal(0.25, 0.30, (n, 1)) # MT 擾動
df = pd.DataFrame(np.round(X, 2), columns=cfg.sample_cols)
df.insert(0, "gene_symb", genes)
df.insert(1, "Transcriptid", [f"ENST{i:011d}" for i in range(n)])
cfg.input_folder.mkdir(parents=True, exist_ok=True)
with pd.ExcelWriter(cfg.input_folder / cfg.file_name) as w:
df.to_excel(w, sheet_name=cfg.sheet_name, index=False)
print(f"🧪 已產生示範資料: {cfg.input_folder / cfg.file_name} ({n} 基因)")
def main(cfg: Config, offline: bool) -> None:
cfg.out_dir.mkdir(parents=True, exist_ok=True)
# 1) KEGG 拓撲
kegg_net, alias2symbol = build_kegg_network(cfg, offline)
print(f"📊 KEGG 網路: {kegg_net.number_of_nodes()} 節點 / "
f"{kegg_net.number_of_edges()} 邊")
# 2) 表達量
df_expr, mat_raw = load_expression(cfg)
fc_tbl = compute_log2fc(mat_raw)
# 3) 比對
common, expr, missing = match_genes(kegg_net, alias2symbol, df_expr)
print(f"✅ 配對成功的 Axon guidance 基因: {len(common)} / {kegg_net.number_of_nodes()}")
if len(common) < 3:
raise ValueError("🚨 配對基因過少，請確認 Excel 的 Gene Symbol 格式")
# 4) 加權子圖 + Leiden
W, edge_tbl = build_weighted_graph(kegg_net, expr, cfg)
membership, q = weighted_leiden(W, cfg)
# 5) 節點指標
deg_c = nx.degree_centrality(W)
btw = nx.betweenness_centrality(W, weight=None)
fc_map = {}
for sym in W.nodes():
src = [s for s in df_expr.index if alias2symbol.get(s) == sym]
fc_map[sym] = (fc_tbl.loc[src, "log2FC"].mean() if src else np.nan,
fc_tbl.loc[src, "p_paired"].mean() if src else np.nan)
node_tbl = pd.DataFrame({
"Gene_Symbol": list(W.nodes()),
"KEGG_Leiden_Module": [membership[n] for n in W.nodes()],
"degree": [W.degree(n) for n in W.nodes()],
"degree_centrality": [deg_c[n] for n in W.nodes()],
"betweenness": [btw[n] for n in W.nodes()],
"log2FC": [fc_map[n][0] for n in W.nodes()],
"p_paired": [fc_map[n][1] for n in W.nodes()],
}).sort_values(["KEGG_Leiden_Module", "betweenness"],
ascending=[True, False]).reset_index(drop=True)
node_tbl = node_tbl.merge(expr, left_on="Gene_Symbol",
right_index=True, how="left")
mod_tbl = (node_tbl.groupby("KEGG_Leiden_Module")
.agg(n_genes=("Gene_Symbol", "size"),
mean_log2FC=("log2FC", "mean"),
top_hub=("Gene_Symbol", "first"),
genes=("Gene_Symbol", lambda s: ", ".join(sorted(s))))
.reset_index())
edge_tbl["module_source"] = edge_tbl["source"].map(membership)
edge_tbl["module_target"] = edge_tbl["target"].map(membership)
edge_tbl["is_intra_module"] = (edge_tbl["module_source"]
== edge_tbl["module_target"])
# 6) 圖表
print("🎨 產生圖表...")
pos = _layout(W, cfg.seed)
figs = []
p = cfg.out_dir / "fig1_network_modules.png"
fig_network_modules(W, pos, membership, cfg, p); figs.append(p)
p = cfg.out_dir / "fig2_module_facets.png"
fig_module_facets(W, pos, membership, cfg, p); figs.append(p)
p = cfg.out_dir / "fig3_network_log2fc.png"
if fig_network_log2fc(W, pos, membership, node_tbl, cfg, p):
figs.append(p)
p = cfg.out_dir / "fig4_corr_heatmap.png"
fig_corr_heatmap(expr, membership, cfg, p); figs.append(p)
p = cfg.out_dir / "fig5_modules_and_hubs.png"
fig_hubs_and_sizes(node_tbl, membership, cfg, p); figs.append(p)
# 7) Excel（多工作表 = 圖表的表格版，滿足無障礙「表格檢視」）
out_xlsx = cfg.out_dir / f"Axon guidance_{cfg.pathway_id}_Weighted_Leiden.xlsx"
summary = pd.DataFrame({
"item": ["pathway", "KEGG nodes", "KEGG edges", "matched genes",
"graph nodes", "graph edges", "modules", "modularity Q",
"resolution", "min |r|", "corr method", "seed"],
"value": [cfg.pathway_id, kegg_net.number_of_nodes(),
kegg_net.number_of_edges(), len(common),
W.number_of_nodes(), W.number_of_edges(),
len(set(membership.values())), round(q, 4),
cfg.resolution, cfg.min_abs_r, cfg.corr_method, cfg.seed],
})
with pd.ExcelWriter(out_xlsx, engine="openpyxl") as w:
summary.to_excel(w, sheet_name="Summary", index=False)
node_tbl.to_excel(w, sheet_name="Genes_Modules", index=False)
mod_tbl.to_excel(w, sheet_name="Module_Summary", index=False)
edge_tbl.to_excel(w, sheet_name="Edge_List", index=False)
pd.DataFrame({"unmatched_KEGG_gene": missing}).to_excel(
w, sheet_name="Unmatched", index=False)
print("\n🎉 完成！")
print(f"💾 Excel : {out_xlsx}")
for f in figs:
print(f"🖼️ Figure: {f}")
print("\n各模組基因數：")
print(mod_tbl[["KEGG_Leiden_Module", "n_genes", "top_hub"]].to_string(index=False))
if __name__ == "__main__":
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--offline", action="store_true", help="不連 KEGG，用內建拓撲")
ap.add_argument("--demo", action="store_true", help="產生模擬資料並跑通全流程")
ap.add_argument("--input-folder", type=Path)
ap.add_argument("--file-name", type=str)
ap.add_argument("--resolution", type=float, default=1.0)
ap.add_argument("--min-abs-r", type=float, default=0.30)
ap.add_argument("--corr-method", choices=["pearson", "spearman"], default="pearson")
ap.add_argument("--show-bugfixes", action="store_true", help="列出修正清單後結束")
a = ap.parse_args()
if a.show_bugfixes:
print(BUGFIX_NOTES)
sys.exit(0)
cfg = Config(resolution=a.resolution, min_abs_r=a.min_abs_r,
corr_method=a.corr_method)
if a.input_folder:
cfg.input_folder = a.input_folder
if a.file_name:
cfg.file_name = a.file_name
if a.demo:
make_demo_data(cfg)
try:
main(cfg, offline=a.offline or a.demo)
except Exception as exc:
print(f"\n❌ 執行失敗: {type(exc).__name__}: {exc}", file=sys.stderr)
sys.exit(1)
No comments:
Post a Comment