WMEAN 的降維機制： WMEAN 利用「先驗知識（Prior Knowledge）」（基因集定義，如
KEGG 或 GO 資料庫），把這 20,000 個基因進行「歸類與加權權衡」。
透過矩陣運算，它把原本 [樣本 × 20,000 個基因] 的龐大矩陣，直接壓縮、轉換成 [樣本
× 5 條核心通路] 的極精簡矩陣。
當要計算 AMPK_SIGNALING_PATHWAY 的活性時，WMEAN
就像雷達一樣，自動從兩萬個基因中定位出這 5 個核心成員：PRKAA1, PRKAA2, STK11,
CAMKK2,
PPARGC1A。其他不相關的一萬九千多個基因，在此通路的計算中直接被「無視」（權重設為
0）。
SPIA (Signaling Pathway Impact Analysis) 是一個非常經典且強大的演算法。
它與一般 GSEA 或 WMEAN
最大的不同在於：它不僅看基因的表現量變化，還把通路內部的「網絡拓撲結構（Topology）」納入矩陣運算。
要計算並解讀 SPIA (Signaling Pathway Impact Analysis) 的網絡擾動分析
(Perturbation
Network)，其核心邏輯是將「基因表現的變異」轉化為「網路上的訊號傳遞與擾動」。
1.
核心發現：神經系統發育的強烈擾動 目前的結果顯示，你的處理對 「神經生成
(NEUROGENESIS)」 具有壓倒性的影響力。 圧倒性的第一名：NEUROGENESIS 的
Total_Abs_Perturbation 高達 85.86，遠高於第二名的 SYNAPTIC_FUNCTION (30.22)。
解讀：這意味著在你的實驗處理下，從幹細胞轉化為神經元的訊號級聯反應，在數學模型中產生了極大的偏離（擾動）。這不僅是個別基因表現改變，而是整個神經發育網路的狀態發生了劇烈位移。
2. 網路階層的連鎖反應 你的前三名通路分別是：
NEUROGENESIS (85.86)
SYNAPTIC_FUNCTION (30.22)
AXON_GUIDANCE (25.47)
這三個通路在生物學上高度相關，形成了一個「功能軸線」。這暗示你的處理可能正影響著神經發育的早期決定階段（Neurogenesis），進而連鎖引發了突觸功能的成熟（Synaptic
Function）以及軸突的生長導向（Axon
Guidance）變化。這是一個非常具備一致性的生理調控模型。 3. 低影響通路的排除
排在最後的 HIPPO_SIGNALING_PATHWAY (0.71) 與 PI3K_AKT_MTOR_SIGNALING (6.79)
相對穩定。這是一個重要的負向控制結果，顯示你的實驗處理並非「全面性地擾亂細胞」，而是具備高度的「細胞功能專一性」，主要針對神經發育相關機制，而非基礎代謝或器官大小調控機制。
下一步深入分析建議
為了將這些發現推向報告或論文層次，建議你進行以下兩項深度分析： 開啟
SPIA_NEUROGENESIS_Impact.csv：找出該通路中 B 值（Perturbation
B）絕對值最大的那幾個基因。這些基因就是你這項處理的「主控因子 (Master
Regulators)」。
觀察網路拓撲圖：打開
SPIA_NEUROGENESIS_Network.png。觀察訊號是如何從特定的轉錄因子（如 SOX2,
PAX6）流向下游，並觀察這些節點是如何連接的。
後續觀察： 在這份 NEUROGENESIS
的詳細 CSV 表格中，B
值最高的基因是哪一個？如果該基因是轉錄因子，那它很可能就是你實驗處理導致神經發育表現改變的最核心調控者。
數值最大的是「腦源性神經營養因子（BDNF）」不僅在神經系統中發揮作用，在女性生殖系統中也扮演著至關重要的調控角色。
在卵巢中，BDNF
主要透過與其高親和力受體 TrkB（酪氨酸激酶受體
B）結合，參與卵泡的生長與卵母細胞的成熟過程。
以下歸納 BDNF
在卵巢生理發育中的主要功能：
1. 促進卵泡發育與生長 (Follicular Development)
2.
調控卵母細胞成熟 (Oocyte Maturation)
3. 參與排卵與類固醇激素合成
TrkB（酪氨酸激酶受體 B）信息路徑總結為： 生存路線
(PI3K/AKT)：決定神經元「活不活」。 分化與記憶路線
(MAPK/ERK)：決定神經元「怎麼變」與「記住什麼」。 突觸功能路線
(PLC$\gamma$)：決定神經元「如何溝通」
import argparse
from pathlib import Path
from typing import Dict, List, Sequence, Tuple, Optional
import networkx as nx
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# =====================================================================
# 1. 系統設定與日誌
# =====================================================================
sns.set_theme(style="ticks", context="talk")
plt.rcParams["font.family"] = "sans-serif"
plt.rcParams["axes.unicode_minus"] = False
def log(msg: str) -> None:
"""標準化輸出"""
print(msg)
from typing import Dict, List, Tuple
# =====================================================================
# 2. 基因集與拓撲資料庫 (Pathway Database)
# =====================================================================
PATHWAY_GENESETS: Dict[str, List[str]] = {
# 既有路徑
"AMPK_SIGNALING_PATHWAY": ["PRKAA1", "PRKAA2", "STK11", "CAMKK2", "PPARGC1A", "PRKAB1", "PRKAG1", "TSC2", "EEF2K", "PRKAG2", "PRKAB2", "FOXO3", "SIRT1", "ULK1", "ATG13", "RPTOR", "SLC2A4", "ACACA"],
"AXON_GUIDANCE": ["SEMA3A", "EPHA2", "SLIT2", "ROBO1", "PLXNA1", "NETRIN1", "DCC", "TRIO", "RAC1", "SRGAP1", "EPHB2", "EFNB2", "CDC42", "PAK1", "ENAH", "ABLIM1", "DPYSL2", "RHOA"],
"HIPPO_SIGNALING_PATHWAY": ["MST1", "LATS1", "LATS2", "YAP1", "TAZ", "SAV1", "MOB1A", "TEAD1", "BIRC5", "NF2", "FRMD6", "AMOTL2", "CCN2", "ANKRD1", "TEAD2", "TEAD4", "MOB1B"],
"ESTROGEN_RESPONSE": ["ESR1", "ESR2", "PGR", "MYC", "GREB1", "NCOA1", "NCOA2", "FOS", "JUN", "NCOA3", "NRIP1", "GATA3", "FOXA1", "TFF1", "XBP1", "KRT19"],
"PI3K_AKT_MTOR_SIGNALING": ["PIK3CA", "AKT1", "MTOR", "PTEN", "RPS6KB1", "EGFR", "PIK3R1", "AKT2", "GSK3B", "EIF4EBP1", "AKT3", "PDPK1", "RICTOR", "TSC1", "RHEB", "FOXO1", "MAPKAP1", "IGF1R", "PIK3CB"],
"WNT_BETA_CATENIN_SIGNALING": ["WNT1", "WNT3A", "FZD1", "CTNNB1", "AXIN1", "APC", "GSK3B", "TCF7L2", "CCND1", "LRP5", "LRP6", "DVL1", "AXIN2", "CSNK1A1", "LEF1", "BTRC", "SFRP1"],
"APOPTOSIS_PATHWAY": ["FAS", "BAX", "BAK1", "BCL2", "BCL2L1", "CYCS", "CASP9", "CASP3", "DIABLO", "CASP8", "CASP7", "APAF1", "TP53", "BID", "BAD", "XIAP", "TNFRSF10A"],
"NOTCH_SIGNALING_PATHWAY": ["NOTCH1", "NOTCH2", "DLL1", "DLL4", "JAG1", "RBPJ", "HES1", "HEY1", "NOTCH3", "NOTCH4", "JAG2", "MAML1", "HEY2", "DTX1", "PSEN1", "ADAM17"],
"CELL_CYCLE_G1_S": ["CCND1", "CCNE1", "CDK4", "CDK6", "CDK2", "RB1", "E2F1", "CDKN1A", "CDKN2A", "MYC"],
"TGF_BETA_SIGNALING": ["TGFB1", "TGFBR1", "TGFBR2", "SMAD2", "SMAD3", "SMAD4", "SMAD7", "SKI", "ZEB1"],
"JAK_STAT_SIGNALING": ["IL6R", "IL6ST", "JAK1", "JAK2", "STAT3", "STAT1", "SOCS3", "PIAS1", "IL6"],
# 🟢 新增路徑
"ANGIOGENESIS": ["VEGFA", "VEGFB", "VEGFC", "KDR", "FLT1", "ANGPT1", "ANGPT2", "TEK", "HIF1A", "NOS3", "CDH5", "PECAM1", "FGF2", "FGFR1", "MMP9", "MMP2", "EPAS1"],
"NEUROGENESIS": ["SOX2", "PAX6", "NEUROD1", "NES", "DCX", "TUBB3", "BDNF", "NTRK2", "NGF", "NTRK1", "CREB1", "ASCL1", "HES1", "MAP2", "NCAM1", "FOXG1", "MASH1"],
"SYNAPTIC_FUNCTION": ["SYP", "SYN1", "SNAP25", "VAMP2", "STX1A", "GRIA1", "GRIA2", "GRIN1", "GRIN2A", "GRIN2B", "CAMK2A", "DLG4", "HOMER1", "SYT1", "SLC17A7", "GABRA1"]
}
EdgeType = Tuple[str, str, str]
PATHWAY_TOPOLOGIES: Dict[str, List[EdgeType]] = {
# 既有拓撲
"HIPPO_SIGNALING_PATHWAY": [("MST1", "LATS1", "activation"), ("LATS1", "YAP1", "inhibition"), ("LATS1", "TAZ", "inhibition"), ("YAP1", "BIRC5", "activation"), ("TAZ", "BIRC5", "activation")],
"AXON_GUIDANCE": [("NETRIN1", "DCC", "activation"), ("DCC", "TRIO", "activation"), ("TRIO", "RAC1", "activation"), ("RAC1", "ACTIN_POLYMERIZATION", "activation"), ("SLIT2", "ROBO1", "activation"), ("ROBO1", "SRGAP1", "activation"), ("SRGAP1", "RAC1", "inhibition")],
"AMPK_SIGNALING_PATHWAY": [("STK11", "PRKAA1", "activation"), ("CAMKK2", "PRKAA1", "activation"), ("PRKAA1", "TSC2", "activation"), ("TSC2", "RPTOR", "inhibition"), ("PRKAA1", "EEF2K", "activation"), ("PRKAA1", "PPARGC1A", "activation")],
"PI3K_AKT_MTOR_SIGNALING": [("EGFR", "PIK3CA", "activation"), ("PIK3CA", "AKT1", "activation"), ("PTEN", "AKT1", "inhibition"), ("AKT1", "TSC2", "inhibition"), ("TSC2", "RHEB", "inhibition"), ("RHEB", "MTOR", "activation"), ("MTOR", "RPS6KB1", "activation"), ("RPS6KB1", "EIF4EBP1", "activation")],
"WNT_BETA_CATENIN_SIGNALING": [("WNT3A", "FZD1", "activation"), ("FZD1", "DVL1", "activation"), ("DVL1", "AXIN1", "inhibition"), ("AXIN1", "CTNNB1", "inhibition"), ("APC", "CTNNB1", "inhibition"), ("CTNNB1", "TCF7L2", "activation"), ("TCF7L2", "CCND1", "activation")],
"APOPTOSIS_PATHWAY": [("FAS", "CASP8", "activation"), ("CASP8", "BID", "activation"), ("BID", "BAX", "activation"), ("BCL2", "BAX", "inhibition"), ("BAX", "CYCS", "activation"), ("CYCS", "APAF1", "activation"), ("APAF1", "CASP9", "activation"), ("CASP9", "CASP3", "activation"), ("XIAP", "CASP3", "inhibition")],
"NOTCH_SIGNALING_PATHWAY": [("DLL1", "NOTCH1", "activation"), ("JAG1", "NOTCH1", "activation"), ("ADAM17", "NOTCH1", "activation"), ("PSEN1", "NOTCH1", "activation"), ("NOTCH1", "RBPJ", "activation"), ("RBPJ", "HES1", "activation"), ("HES1", "HEY1", "activation")],
# 🟢 新增拓撲
"ANGIOGENESIS": [
("HIF1A", "VEGFA", "activation"), # 缺氧誘導 VEGF 表現
("VEGFA", "KDR", "activation"), # VEGF 結合並活化 VEGFR2 (KDR)
("KDR", "NOS3", "activation"), # VEGFR2 活化內皮一氧化氮合酶
("ANGPT1", "TEK", "activation"), # Angiopoietin 1 活化 Tie2 (TEK) 受體促進血管穩定
("ANGPT2", "TEK", "inhibition"), # Angiopoietin 2 拮抗 Tie2
("FGF2", "FGFR1", "activation") # FGF 訊號傳遞
],
"NEUROGENESIS": [
("BDNF", "NTRK2", "activation"), # BDNF 結合 TrkB (NTRK2) 促進神經存活分化
("NTRK2", "CREB1", "activation"), # TrkB 下游活化轉錄因子 CREB1
("NGF", "NTRK1", "activation"), # NGF 活化 TrkA (NTRK1)
("SOX2", "PAX6", "activation"), # 神經幹細胞核心轉錄因子調控
("PAX6", "NEUROD1", "activation"), # 促進神經元分化
("HES1", "ASCL1", "inhibition") # Notch 下游 HES1 抑制促神經分化因子 ASCL1
],
"SYNAPTIC_FUNCTION": [
("GRIN1", "CAMK2A", "activation"), # NMDA 受體活化引發鈣離子內流，活化 CaMKII
("CAMK2A", "GRIA1", "activation"), # CaMKII 磷酸化 AMPA 受體 (GRIA1) 促進突觸可塑性 (LTP)
("CAMK2A", "DLG4", "activation"), # CaMKII 與 PSD-95 (DLG4) 交互作用
("VAMP2", "SNAP25", "activation"), # SNARE 複合體形成 (突觸小泡釋放)
("STX1A", "SNAP25", "activation"), # SNARE 複合體核心結合
("BDNF", "NTRK2", "activation") # BDNF 在突觸端的局部活化作用
],
"MITOCHONDRIAL_SIGNALING": [
"MT-CO1", "MT-ND1", "MT-ATP6", # 呼吸鏈核心組件
"TFAM", "POLG", # 粒線體 DNA 複製與轉錄
"OPA1", "DNM1L", # 粒線體融合與分裂
"PINK1", "PRKN", # 粒線體自噬 (Mitophagy)
"CYCS", "VDAC1", # 細胞凋亡關鍵因子
"PPARGC1A", "NRF1", # 粒線體生物合成調控
"UCP2", "SOD2" # 氧化壓力與解偶聯
]
}
PATHWAY_TOPOLOGIES.update({
"ANGIOGENESIS": [
("HIF1A", "VEGFA", "activation"),
("VEGFA", "KDR", "activation"),
("KDR", "NOS3", "activation"),
("ANGPT1", "TEK", "activation"),
("ANGPT2", "TEK", "inhibition"),
("FGF2", "FGFR1", "activation"),
# --- 擴充部分 ---
("KDR", "PIK3CA", "activation"), # 血管生成訊號與代謝通路連結
("HIF1A", "NOS3", "activation"), # 缺氧直接促進 NO 生成
("TEK", "PTEN", "inhibition"), # 血管穩定性調控
("NOS3", "VEGFA", "activation") # 血管生成的正回饋迴路
],
"NEUROGENESIS": [
("BDNF", "NTRK2", "activation"),
("NTRK2", "CREB1", "activation"),
("NGF", "NTRK1", "activation"),
("SOX2", "PAX6", "activation"),
("PAX6", "NEUROD1", "activation"),
("HES1", "ASCL1", "inhibition"),
# --- 擴充部分 ---
("CREB1", "BDNF", "activation"), # BDNF 自體正回饋，維持神經發育活性
("NTRK2", "MAPK1", "activation"), # 擴充 MAPK 傳遞路徑
("ASCL1", "NEUROD1", "activation"),# 分化級聯
("HES1", "PAX6", "inhibition"), # Notch 抑制早期神經幹細胞標記
("NEUROD1", "MAP2", "activation") # 終末分化標記
],
"SYNAPTIC_FUNCTION": [
("GRIN1", "CAMK2A", "activation"),
("CAMK2A", "GRIA1", "activation"),
("CAMK2A", "DLG4", "activation"),
("VAMP2", "SNAP25", "activation"),
("STX1A", "SNAP25", "activation"),
("BDNF", "NTRK2", "activation"),
# --- 擴充部分 ---
("CAMK2A", "CREB1", "activation"), # 突觸活動影響轉錄
("GRIN2B", "CAMK2A", "activation"),# NMDA 亞基特異性活化
("SNAP25", "SYP", "activation"), # 突觸小泡釋放機制的連鎖
("DLG4", "GRIA1", "activation"), # 突觸後緻密區 (PSD) 的結構穩定
("NTRK2", "VAMP2", "activation") # 神經營養因子促進小泡轉運
],
"MITOCHONDRIAL_SIGNALING": [
("PPARGC1A", "NRF1", "activation"), # PGC-1a 活化 NRF1 促進粒線體生物合成
("NRF1", "TFAM", "activation"), # NRF1 促進 TFAM 表達
("PINK1", "PRKN", "activation"), # PINK1 招募 PRKN 啟動粒線體自噬
("PRKN", "DNM1L", "inhibition"), # 抑制分裂，維持粒線體網狀結構
("DNM1L", "OPA1", "inhibition"), # 分裂與融合的拮抗平衡
("VDAC1", "CYCS", "activation"), # 粒線體外膜通透性增加釋放細胞色素 c
("SOD2", "ROS", "inhibition"), # SOD2 清除粒線體產生的活性氧
("MT-CO1", "ATP_PRODUCTION", "activation")
]
})
# =====================================================================
# 3. 核心工具層 (Common Utils)
# =====================================================================
def parse_cols(df: pd.DataFrame, cols_str: str) -> List[str]:
"""解析 Excel 欄位（支援數字索引、字母、字串名稱）"""
if not cols_str: return list(df.columns)
resolved, seen = [], set()
for p in (x.strip() for x in cols_str.split(",")):
if p.isdigit() and int(p) < len(df.columns):
resolved.append(df.columns[int(p)])
elif p.isalpha() and p.isupper() and len(p) <= 2:
idx = sum((ord(c) - 64) * (26 ** i) for i, c in enumerate(p[::-1])) - 1
if idx < len(df.columns): resolved.append(df.columns[idx])
elif p in df.columns:
resolved.append(p)
return [c for c in resolved if not (c in seen or seen.add(c))]
def load_expression_matrix(path: Path, sheet: str, usecol: str) -> Tuple[pd.DataFrame, List[str]]:
"""讀取 Excel、標準化並進行基因去重與平均"""
if not path.exists(): raise FileNotFoundError(f"❌ 找不到檔案: {path}")
df_raw = pd.read_excel(path, sheet_name=sheet)
df_raw.columns = df_raw.columns.astype(str).str.strip()
df_raw.rename(columns={df_raw.columns[0]: "Transcriptid", df_raw.columns[1]: "Gene_symbol"}, inplace=True)
target_cols = parse_cols(df_raw, usecol)
df_filt = df_raw[target_cols].copy()
df_filt["Gene_symbol"] = df_filt["Gene_symbol"].astype(str).str.strip().str.upper()
df_filt = df_filt[~df_filt["Gene_symbol"].isin(["NAN", "#N/A"])]
sample_cols = [c for c in target_cols if c not in ("Transcriptid", "Gene_symbol")]
return df_filt.groupby("Gene_symbol")[sample_cols].mean(), sample_cols
def compute_log2fc(df_pivot: pd.DataFrame, epsilon: float = 1e-5) -> pd.DataFrame:
"""計算 Log2FC，自動識別對照(MT)與實驗(AT)樣本"""
cols = list(df_pivot.columns)
# 智慧分組：根據樣本名稱中的 MT / AT 進行分組並排序
ctrl = sorted([c for c in cols if "MT" in str(c)])
treat = sorted([c for c in cols if "AT" in str(c)])
# 容錯機制：如果找不到 MT/AT，則沿用舊的「對半分割」邏輯
if not ctrl or not treat:
half = len(cols) // 2
ctrl, treat = cols[:half], cols[half:]
log(f" -> 基準對照組樣本: {ctrl}\n -> 處理實驗組樣本: {treat}")
mean_c, mean_t = df_pivot[ctrl].mean(axis=1), df_pivot[treat].mean(axis=1)
return pd.DataFrame({"Log2FC": np.log2((mean_t + epsilon) / (mean_c + epsilon))}, index=df_pivot.index)
# =====================================================================
# 4. 視覺化模組 (Visualization)
# =====================================================================
def plot_results(df_scores: pd.DataFrame, out_dir: Path) -> None:
"""WMEAN 結果視覺化 (維持原有功能)"""
if df_scores.empty: return
# Heatmap
df_zscore = df_scores.apply(lambda x: (x - x.mean()) / x.std() if x.std() else x, axis=0)
plt.figure(figsize=(10, 8))
sns.heatmap(df_zscore.T, cmap="RdBu_r", center=0, robust=True, annot=True, fmt=".2f", linewidths=0.5)
plt.title("Pathway Activity Profiles (WMEAN Z-score)", pad=20, weight="bold")
plt.tight_layout()
plt.savefig(out_dir / "WMEAN_Pathway_Heatmap.png", dpi=300)
plt.show()
# Paired Dynamics
at_samples = sorted([s for s in df_scores.index if "AT" in str(s)])
mt_samples = sorted([s for s in df_scores.index if "MT" in str(s)])
if len(at_samples) == len(mt_samples) and at_samples:
rows = []
for at_s, mt_s in zip(at_samples, mt_samples):
pair_id = at_s.replace("AT", "").replace("MT", "")
for pw in df_scores.columns:
rows.extend([
{"Pair": pair_id, "Condition": "Control (MT)", "Score": np.log2(df_scores.loc[mt_s, pw] + 1e-5), "Pathway": pw},
{"Pair": pair_id, "Condition": "Experimental (AT)", "Score": np.log2(df_scores.loc[at_s, pw] + 1e-5), "Pathway": pw}
])
g = sns.catplot(data=pd.DataFrame(rows), x="Condition", y="Score", hue="Condition", col="Pathway",
col_wrap=3, kind="box", palette=["#3498db", "#e74c3c"], height=5, aspect=0.8, legend=False)
for ax, (_, sub_df) in zip(g.axes.flat, pd.DataFrame(rows).groupby("Pathway")):
for pid in sub_df["Pair"].unique():
pair_sub = sub_df[sub_df["Pair"] == pid].sort_values("Condition")
ax.plot(pair_sub["Condition"], pair_sub["Score"], color="#7f8c8d", alpha=0.5, ls="--")
sns.stripplot(data=sub_df, x="Condition", y="Score", color="#2c3e50", alpha=0.6, ax=ax)
g.fig.suptitle("Pathway Activation Dynamics", weight="bold", y=1.02)
g.savefig(out_dir / "WMEAN_Pathway_Dynamics_Log.png", dpi=300)
plt.show("all")
def plot_spia_network(G: nx.DiGraph, res_df: pd.DataFrame, pw_name: str, out_dir: Path) -> None:
"""SPIA 新增: 繪製 Pathway 拓撲網絡與節點擾動熱力映射"""
plt.figure(figsize=(10, 8))
pos = nx.spring_layout(G, seed=42) # 使用彈簧佈局
# 獲取節點擾動值作為顏色 (若無數據則設為 0)
node_colors = [res_df.loc[n, "SPIA_Perturbation_B"] if n in res_df.index else 0.0 for n in G.nodes]
vmax = max(abs(min(node_colors, default=0)), abs(max(node_colors, default=0)), 1e-5)
# 定義邊緣顏色 (紅 = 活化, 藍 = 抑制)
edge_colors = ["#e74c3c" if d.get("type") == "activation" else "#3498db" for u, v, d in G.edges(data=True)]
# 繪製網路
nodes = nx.draw_networkx_nodes(G, pos, node_color=node_colors, cmap="RdBu_r",
vmin=-vmax, vmax=vmax, node_size=800, edgecolors="gray")
nx.draw_networkx_edges(G, pos, edge_color=edge_colors, arrowsize=15,
connectionstyle="arc3,rad=0.1", width=1.5, alpha=0.7)
nx.draw_networkx_labels(G, pos, font_size=9, font_weight="bold")
plt.colorbar(nodes, label="Perturbation Level (B)")
plt.title(f"Topology Perturbation Network:\n{pw_name}", weight="bold", pad=15)
plt.axis("off")
plt.tight_layout()
plt.savefig(out_dir / f"SPIA_{pw_name}_Network.png", dpi=300)
plt.show()
def plot_spia_summary(df_summary: pd.DataFrame, out_dir: Path) -> None:
"""SPIA 新增: 繪製所有 Pathway 的總擾動排名"""
if df_summary.empty: return
plt.figure(figsize=(10, 6))
df_sorted = df_summary.sort_values("Total_Abs_Perturbation", ascending=False)
sns.barplot(data=df_sorted, x="Total_Abs_Perturbation", y="Pathway", palette="viridis")
plt.title("SPIA: Total Accumulative Perturbation Rank", weight="bold")
plt.xlabel("Total Absolute Perturbation (|B|)")
plt.ylabel("")
plt.tight_layout()
plt.savefig(out_dir / "SPIA_Summary_Ranking.png", dpi=300)
plt.show()
# =====================================================================
# 5. 分析主體 (WMEAN & exact-SPIA)
# =====================================================================
def run_wmean(df_expr: pd.DataFrame, out_dir: Path) -> pd.DataFrame:
"""WMEAN 計算"""
results = {}
for pw, genes in PATHWAY_GENESETS.items():
valid_genes = [g for g in genes if g in df_expr.columns]
if valid_genes:
results[pw] = df_expr[valid_genes].mean(axis=1)
df_scores = pd.DataFrame(results)
df_scores.to_csv(out_dir / "WMEAN_Pathway_Scores.csv")
plot_results(df_scores, out_dir)
return df_scores
def run_spia(df_degs: pd.DataFrame, out_dir: Path) -> None:
"""SPIA 精確解計算與視覺化整合"""
summary = []
for pw_name, edges in PATHWAY_TOPOLOGIES.items():
# 建構有向圖
G_nx = nx.DiGraph()
G_nx.add_edges_from([(u, v, {"type": t}) for u, v, t in edges])
nodes = list(G_nx.nodes)
n = len(nodes)
node_idx = {node: i for i, node in enumerate(nodes)}
# 建立 G 向量與 W 權重矩陣
G_vec, W_mat = np.zeros(n), np.zeros((n, n))
fc_col = df_degs.columns[0]
for node in nodes:
if node in df_degs.index:
G_vec[node_idx[node]] = df_degs.loc[node, fc_col]
for u, v, data in G_nx.edges(data=True):
sign = 1.0 if data.get("type") == "activation" else -1.0
W_mat[node_idx[v], node_idx[u]] = sign / G_nx.out_degree(u)
# 解矩陣
I = np.eye(n)
try:
B_exact = np.linalg.solve(I - W_mat, G_vec)
except np.linalg.LinAlgError:
B_exact = np.dot(np.linalg.pinv(I - W_mat), G_vec)
res_df = pd.DataFrame({
"Log2FC_G": G_vec,
"SPIA_Perturbation_B": B_exact,
"Accumulated_Impact": B_exact - G_vec
}, index=nodes)
# 輸出 CSV
res_df.to_csv(out_dir / f"SPIA_{pw_name}_Impact.csv")
# 🟢 新增：繪製此 Pathway 的拓撲網路圖
plot_spia_network(G_nx, res_df, pw_name, out_dir)
# 記錄總結
summary.append({
"Pathway": pw_name,
"N_Nodes": n,
"Total_Abs_Perturbation": res_df["SPIA_Perturbation_B"].abs().sum()
})
# 🟢 新增：輸出總結 CSV 與繪製整體排名長條圖
df_summary = pd.DataFrame(summary).sort_values("Total_Abs_Perturbation", ascending=False)
df_summary.to_csv(out_dir / "SPIA_Summary.csv", index=False)
plot_spia_summary(df_summary, out_dir)
# =====================================================================
# 6. 主程式入口
# =====================================================================
def main():
parser = argparse.ArgumentParser()
parser.add_argument("--excel", default="/Users/yshuang/Documents/Python/FC_GSEA.xlsx")
parser.add_argument("--sheet", default="data")
parser.add_argument("--wmean-cols", default="A, B, F, G, I, J, L, M")
# 根據你實際 Excel 檔案中 1440MT, 2003MT, 2613MT, 1440AT, 2003AT, 2613AT 所在的欄位調整字母
#parser.add_argument("--spia-cols", default="A, B, C, D, E, F, G, H")
parser.add_argument("--spia-cols", default="A, B, F, G, I, J, L, M")
parser.add_argument("--outdir", default="/Users/yshuang/Documents/Python/wMEAN_SPIA_Results")
args = parser.parse_args()
out_dir = Path(args.outdir)
out_dir.mkdir(parents=True, exist_ok=True)
log("🚀 啟動多體學通路分析流程...\n")
# 共用資料載入
df_pivot_wmean, _ = load_expression_matrix(Path(args.excel), args.sheet, args.wmean_cols)
df_pivot_spia, _ = load_expression_matrix(Path(args.excel), args.sheet, args.spia_cols)
log("📊 [階段 1] 執行 WMEAN 活性推算...")
run_wmean(df_pivot_wmean.T, out_dir)
log("\n📈 [階段 2] 執行 SPIA 拓撲擾動分析...")
run_spia(compute_log2fc(df_pivot_spia), out_dir)
log(f"\n✅ 全部分析完成！結果已輸出至: {out_dir}")
if __name__ == "__main__":
main()
No comments:
Post a Comment