「高表達不等於作用大，且高表達基因確實往往是『功能的末端執行者』；而低表達基因則多半扮演『上游的指揮開關』。」
💡 一、高表達 = 作用大？高表達 = 作用末端？
1. 高表達 ≠ 作用大（表現量高只是因為「消耗量大」）
在細胞內，RNA 拷貝數高往往是因為該蛋白質需要大量作為結構組成或大規模參與催化。
例如：核糖體蛋白、組蛋白、細胞骨架（Actin）、DNA 複製相關蛋白。
這些基因就像工廠裡的工人與建材，雖然數量極多，但牠們只是在「執行命令」，並不是決定細胞命運的「決策者」。
2. 高表達 ≒ 作用末端 (Terminal Effectors)
你的直覺非常準確！高表達基因（如分位數 $0.75 \sim 0.95$ 的基因）通常是細胞狀態轉變後的最終表型產物（Phenotypic Endpoint）。
當細胞決定要進行有絲分裂時，它會大量轉錄 DNA 複製與細胞週期蛋白（這就是為什麼在上調高分位數中，
CELL_CYCLE的 $-\log_{10}P$ 會高達 24.92）。這些基因告訴你的是：「細胞目前正在發生什麼結果（What is happening right now）」。
🧬 二、如何界定與研究「較低表達基因」的作用？
如果高表達基因是「執行者」，那麼較低表達基因（如分位數 $0.01 \sim 0.50$）通常就是「開關與指揮官」（如轉錄因子、膜受體、激酶、細胞因子）。
細胞只需要極少量的受體或轉錄因子 mRNA，經由信號放大機制（Signal Amplification），就能驅動整個細胞產生巨大的表型變化。要界定並挖出這些低表達基因的作用，生物資訊學通常採用以下 4 種策略：
1. 分位數過濾（Quantile Stratification）—— 即我們剛才做的事
作法： 將高表達的「背景噪音/終端執行者」（如核糖體、細胞週期蛋白）剝離，單獨對低/中分位數區間進行 GO 富集或 GSEA。
作用： 讓原本被掩蓋的微環境訊號（如 GPCR 訊號轉導、離子通道、細胞間通訊）浮出水面。
2. 重視「相對變化量（Fold Change）」而非「絕對表現量」
作法： 觀察基因的 $\Delta asinh$ 或 $\log_2\text{FC}$。
作用： 一個轉錄因子的絕對表現量可能只有 $asinh = 0.5$，但若它從 $0.1 \rightarrow 0.5$（激增 5 倍），其生物學衝擊力遠大於一個從 $8.0 \rightarrow 8.5$ 的高表達結構蛋白。
3. 蛋白質互作網路（PPI Network）的「瓶頸節點（Bottleneck/Hub）」分析
作法： 將篩選出的低表達基因繪製成 PPI 網絡（如 STRING 資料庫），計算其 Betweenness Centrality（介數中心性）。
作用： 許多低表達蛋白質位於網絡的核心樞紐（Hub），連接不同的路徑。一旦刪除這個低表達基因，整個傳遞網路就會癱瘓。
4. 反推上游調控者分析（Upstream Regulator / Master Regulator Analysis）
作法： 利用如 DoRothEA、ChEA 或 SCENIC 等工具，透過「下游大量改變的高表達基因群」，反推是哪一個上游轉錄因子（TF）在操控牠們。
作用： 即使該轉錄因子本身的 mRNA 表達量極低（甚至在 RNA-seq 中抓不到），也能精準界定出它的核心控制地位。
📝 總結
高表達區間 $[0.51 \sim 0.95]$： 代表表型結果 (Phenotype)，解答「細胞演變成什麼狀態？」（如：大量增殖、癌化、DNA 損傷）。
低/中表達區間 $[0.01 \sim 0.50]$： 代表上游機制與微環境 (Upstream & Environment)，解答「是什麼訊號/受體啟動了這個變化？」（如：激素刺激、離子通道改變、免疫微環境誘導）。
import os
import sys
import json
import logging
import textwrap
import warnings
from dataclasses import dataclass, field, asdict
from pathlib import Path
from typing import Tuple, Optional, List, Dict, Any
from datetime import datetime
import numpy as np
import pandas as pd
from scipy.stats import spearmanr
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import seaborn as sns
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
try:
import gseapy as gp
GSEAPY_AVAILABLE = True
except ImportError:
gp = None
GSEAPY_AVAILABLE = False
warnings.warn("gseapy 未安裝，富集分析將使用內建 mock 結果。")
try:
from adjustText import adjust_text
ADJUSTTEXT_AVAILABLE = True
except ImportError:
ADJUSTTEXT_AVAILABLE = False
def setup_logging(level=logging.INFO, log_file=None, format_str=None) -> logging.Logger:
if format_str is None:
format_str = "%(asctime)s | %(levelname)-8s | %(name)s | %(message)s"
logger = logging.getLogger("OvaryAnalysis")
logger.setLevel(level)
logger.handlers.clear()
formatter = logging.Formatter(format_str, datefmt="%Y-%m-%d %H:%M:%S")
console_handler = logging.StreamHandler(sys.stdout)
console_handler.setLevel(level)
console_handler.setFormatter(formatter)
logger.addHandler(console_handler)
if log_file:
file_handler = logging.FileHandler(log_file, mode='w', encoding='utf-8')
file_handler.setLevel(logging.DEBUG)
file_handler.setFormatter(formatter)
logger.addHandler(file_handler)
return logger
@dataclass(frozen=True)
class AnalysisConfig:
sample_ids: Tuple[str, ...] = ('1440MT', '2003MT', '2613MT', '1440AT', '2003AT', '2613AT')
individual: Tuple[int, ...] = (0, 1, 2, 0, 1, 2)
is_treated: Tuple[int, ...] = (0, 0, 0, 1, 1, 1)
phenotype_order: Tuple[str, ...] = ('1440', '2613', '2003')
phenotype_ranks: Tuple[int, ...] = (3, 2, 1)
default_gmt_path: str = field(
default_factory=lambda: str(Path.home() / "Documents/Python/c5.go.bp.v2023.2.Hs.symbols.gmt")
)
# asinh_lower / asinh_upper 的意義取決於 asinh_use_quantile：
# asinh_use_quantile=False（預設關閉時）：兩者是絕對的 asinh 表現值門檻。
# asinh_use_quantile=True：兩者是「分位數」(0~1 之間的分數)，程式會依
# 實際資料分布動態換算成對應的 asinh 門檻值，例如 0.5 代表中位數、
# 1.0 代表最大值。這樣門檻能隨資料自動校準，不需要每次重新猜測絕對值。
asinh_lower: Optional[float] = 0.75
asinh_upper: Optional[float] = 0.95
asinh_filter_mode: str = "any"
asinh_use_quantile: bool = True
min_delta_range: float = 0.1
min_rho: float = 0.8
top_n_features: Optional[int] = None
output_dir: str = "./ovary_analysis_output"
dpi: int = 300
fig_format: str = "png"
def __post_init__(self):
valid_modes = {"mean", "all", "any"}
if self.asinh_filter_mode not in valid_modes:
raise ValueError(f"asinh_filter_mode 必須是 {valid_modes} 其中之一，收到: {self.asinh_filter_mode!r}")
if self.fig_format not in {"png", "pdf", "svg"}:
raise ValueError(f"fig_format 必須是 png/pdf/svg 其中之一，收到: {self.fig_format!r}")
if self.asinh_use_quantile:
for name, val in (("asinh_lower", self.asinh_lower), ("asinh_upper", self.asinh_upper)):
if val is not None and not (0.0 <= val <= 1.0):
raise ValueError(
f"asinh_use_quantile=True 時，{name} 必須是 0~1 之間的分位數分數，收到: {val!r}"
)
def to_dict(self) -> Dict[str, Any]:
return asdict(self)
def save(self, path: str) -> None:
with open(path, 'w', encoding='utf-8') as f:
json.dump(self.to_dict(), f, indent=2, ensure_ascii=False)
@classmethod
def load(cls, path: str) -> "AnalysisConfig":
with open(path, 'r', encoding='utf-8') as f:
data = json.load(f)
for tuple_field in ("sample_ids", "individual", "is_treated", "phenotype_order", "phenotype_ranks"):
if tuple_field in data and data[tuple_field] is not None:
data[tuple_field] = tuple(data[tuple_field])
return cls(**data)
@property
def mt_samples(self) -> List[str]:
return [sid for sid, t in zip(self.sample_ids, self.is_treated) if t == 0]
@property
def at_samples(self) -> List[str]:
return [sid for sid, t in zip(self.sample_ids, self.is_treated) if t == 1]
@property
def delta_columns(self) -> List[str]:
return [f"Delta_{p}" for p in self.phenotype_order]
@dataclass
class FeatureSelectionResult:
positive_genes: pd.DataFrame
negative_genes: pd.DataFrame
full_report: pd.DataFrame
def summary(self) -> Dict[str, Any]:
return {
"n_positive": len(self.positive_genes),
"n_negative": len(self.negative_genes),
"n_total_tested": len(self.full_report),
"top_positive": self.positive_genes.head(5).index.tolist() if not self.positive_genes.empty else [],
"top_negative": self.negative_genes.head(5).index.tolist() if not self.negative_genes.empty else [],
}
@dataclass
class EnrichmentResult:
dataframe: Optional[pd.DataFrame]
gene_list: List[str]
group_name: str
success: bool = field(init=False)
def __post_init__(self):
self.success = self.dataframe is not None and not self.dataframe.empty
if self.success:
# 依顯著性排序一次，讓後續所有取用（top_terms、繪圖、匯出）都是
# 真正依統計顯著性排序，而非 GMT/Enrichr 回傳的原始（通常是字母）順序。
sort_col = None
for candidate in ("Adjusted P-value", "P-value"):
if candidate in self.dataframe.columns:
sort_col = candidate
break
if sort_col is not None:
self.dataframe = self.dataframe.sort_values(sort_col).reset_index(drop=True)
def top_terms(self, n: int = 10) -> pd.DataFrame:
"""取得依顯著性排序後的前 N 個條目（最顯著者優先）。"""
if not self.success:
return pd.DataFrame()
return self.dataframe.head(n)
class DataLoader:
EXPECTED_COLS = "A,B,F,G,I,J,L,M"
SHEET_NAME = "data"
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
def load(self, file_path: Optional[str] = None) -> pd.DataFrame:
if file_path is None:
file_path = self._find_data_file()
if file_path and os.path.exists(file_path):
expr_matrix = self._load_and_process_data(file_path)
else:
self.logger.warning("未找到有效數據文件，生成模擬數據。")
expr_matrix = self._generate_mock_data()
self._qc_check(expr_matrix)
return expr_matrix
def _find_data_file(self) -> Optional[str]:
candidates = [
"./FC_GSEA.xlsx",
"../FC_GSEA.xlsx",
os.path.expanduser("~/Documents/Python/FC_GSEA.xlsx"),
]
for path in candidates:
if os.path.exists(path):
self.logger.info(f"找到數據文件: {path}")
return path
return None
def _load_and_process_data(self, file_path: str) -> pd.DataFrame:
self.logger.info(f"正在載入數據: {file_path}")
try:
df = pd.read_excel(file_path, sheet_name=self.SHEET_NAME, engine='openpyxl', usecols=self.EXPECTED_COLS)
except Exception as e:
self.logger.error(f"讀取 Excel 失敗: {e}")
raise
df.columns = df.columns.str.strip()
group_col = df.columns[1]
raw_samples = list(self.config.sample_ids)
missing_cols = [c for c in raw_samples if c not in df.columns]
if missing_cols:
raise ValueError(
f"Excel 檔案缺少必要的樣本欄位: {missing_cols}。"
f"現有欄位: {list(df.columns)}"
)
for col in raw_samples:
df[col] = pd.to_numeric(df[col], errors='coerce')
df_indexed = df.set_index(group_col)
n_dup = df_indexed.index.duplicated().sum()
if n_dup > 0:
self.logger.warning(f"發現 {n_dup} 個重複的基因名稱索引，將保留第一筆出現的紀錄")
df_indexed = df_indexed[~df_indexed.index.duplicated(keep='first')]
expr_matrix = df_indexed[raw_samples].copy().fillna(0.0)
self.logger.info(f"原始數據: {expr_matrix.shape[0]} 基因 x {expr_matrix.shape[1]} 樣本")
if (expr_matrix < 0).any().any():
n_neg = (expr_matrix < 0).sum().sum()
self.logger.warning(f"發現 {n_neg} 個負值，asinh 轉換仍可處理，但請確認數據來源是否正確")
expr_transformed = np.arcsinh(expr_matrix)
self.logger.info(f"asinh 轉換完成: {expr_transformed.shape[0]} 基因 x {expr_transformed.shape[1]} 樣本")
return expr_transformed
def _generate_mock_data(self) -> pd.DataFrame:
rng = np.random.default_rng(42)
n_random_genes = 300
pos_trend_genes = ["CYP19A1", "STAR", "LHCGR", "MAPK1", "CDK1"]
neg_trend_genes = ["CCNB1", "TGFBR1", "CASP3", "FOXO1"]
gene_names = pos_trend_genes + neg_trend_genes + [f"GENE_{i:04d}" for i in range(n_random_genes)]
cols = list(self.config.sample_ids)
data = np.abs(rng.normal(loc=5.0, scale=0.4, size=(len(gene_names), len(cols))))
df = pd.DataFrame(data, index=gene_names, columns=cols)
for g in pos_trend_genes:
df.loc[g, self.config.mt_samples] = 4.0
df.loc[g, "1440AT"] = 4.0 + 6.0
df.loc[g, "2613AT"] = 4.0 + 3.5
df.loc[g, "2003AT"] = 4.0 + 1.0
for g in neg_trend_genes:
df.loc[g, self.config.mt_samples] = 4.0
df.loc[g, "1440AT"] = 4.0 - 6.0
df.loc[g, "2613AT"] = 4.0 - 3.5
df.loc[g, "2003AT"] = 4.0 - 1.0
df = df.clip(lower=0.01)
self.logger.info(f"模擬數據生成: {df.shape[0]} 基因 x {df.shape[1]} 樣本")
return np.arcsinh(df)
def _qc_check(self, expr_matrix: pd.DataFrame) -> None:
missing_samples = set(self.config.sample_ids) - set(expr_matrix.columns)
if missing_samples:
self.logger.warning(f"缺少樣本: {missing_samples}")
if expr_matrix.empty:
self.logger.warning("表達矩陣為空")
return
min_val = expr_matrix.min().min()
max_val = expr_matrix.max().max()
self.logger.info(f"asinh 值域: [{min_val:.4f}, {max_val:.4f}]")
zero_genes = (expr_matrix.sum(axis=1) == 0).sum()
if zero_genes > 0:
self.logger.warning(f"發現 {zero_genes} 個全零表達基因")
class AsinhFilter:
"""asinh 數值區間篩選器。
此為全流程唯一負責 asinh 區間篩選的類別，篩選後的結果應直接往下游傳遞，
不應在其他地方（例如 PhenotypeFeatureSelector）重複套用。
"""
def __init__(self, log: logging.Logger, use_quantile: bool = False):
self.logger = log
self.use_quantile = use_quantile
def _resolve_thresholds(
self, data_df: pd.DataFrame, lower: Optional[float], upper: Optional[float], mode: str
) -> Tuple[float, float]:
"""將 lower/upper 轉換成實際用來比較的 asinh 數值門檻。
若 use_quantile=False：lower/upper 本身就是絕對 asinh 值，直接使用。
若 use_quantile=True：lower/upper 是 0~1 之間的分位數分數，會依資料
實際分布（mode="mean" 時用逐基因平均值的分布；"all"/"any" 時用整個
矩陣攤平後的分布）換算成對應的 asinh 數值，讓門檻隨資料自動校準。
"""
if not self.use_quantile:
lo = -np.inf if lower is None else lower
hi = np.inf if upper is None else upper
return lo, hi
if mode == "mean":
values = data_df.mean(axis=1).to_numpy()
else:
values = data_df.to_numpy().flatten()
values = values[~np.isnan(values)]
lo = -np.inf if lower is None else float(np.quantile(values, lower))
hi = np.inf if upper is None else float(np.quantile(values, upper))
self.logger.info(
f"分位篩選: 下分位數 {lower} -> asinh {lo:.4f}, 上分位數 {upper} -> asinh {hi:.4f}"
)
return lo, hi
def filter(self, data_df: pd.DataFrame, lower=None, upper=None, mode="any") -> pd.DataFrame:
if lower is None and upper is None:
self.logger.info("asinh 篩選: 未設置上下限，返回原始數據")
return data_df
lo, hi = self._resolve_thresholds(data_df, lower, upper, mode)
if mode == "mean":
row_stat = data_df.mean(axis=1)
mask = (row_stat >= lo) & (row_stat <= hi)
elif mode in ("all", "any"):
in_range = (data_df >= lo) & (data_df <= hi)
mask = in_range.all(axis=1) if mode == "all" else in_range.any(axis=1)
else:
raise ValueError(f"mode 必須是 'all'、'any' 或 'mean'，收到: {mode!r}")
result = data_df[mask]
n_before = data_df.shape[0]
retention = f"{result.shape[0] / n_before * 100:.1f}%" if n_before else "N/A"
label = f"分位[{lower}, {upper}] (asinh[{lo:.3f}, {hi:.3f}])" if self.use_quantile else f"[{lower}, {upper}]"
self.logger.info(
f"asinh 篩選 {label} mode={mode}: "
f"{n_before} -> {result.shape[0]} 基因 (保留率: {retention})"
)
return result
def compare_modes(self, expr_matrix: pd.DataFrame, special_genes: List[str], lower, upper) -> pd.DataFrame:
self.logger.info("展示 asinh 三種篩選模式差異...")
records = []
n_total = expr_matrix.shape[0]
for mode in ("mean", "all", "any"):
filtered = self.filter(expr_matrix, lower=lower, upper=upper, mode=mode)
kept_special = [g for g in special_genes if g in filtered.index]
records.append({
"mode": mode,
"n_genes": filtered.shape[0],
"retention_rate": f"{filtered.shape[0] / n_total * 100:.1f}%" if n_total else "N/A",
"kept_special_genes": ", ".join(kept_special) if kept_special else "None"
})
return pd.DataFrame(records)
class PhenotypeFeatureSelector:
"""表型引導特徵篩選器。
只負責「表型趨勢」篩選（Delta / Spearman rho），不重複套用 asinh 絕對數值
區間篩選 —— 該篩選已由 AsinhFilter 在 pipeline 更早的步驟完成一次。
"""
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
def select(self, data_df: pd.DataFrame, min_delta_range=None, min_rho=None, top_n=None) -> FeatureSelectionResult:
min_delta_range = self.config.min_delta_range if min_delta_range is None else min_delta_range
min_rho = self.config.min_rho if min_rho is None else min_rho
top_n = self.config.top_n_features if top_n is None else top_n
self.logger.info("開始表型引導特徵篩選...")
deltas = self._compute_deltas(data_df)
results_df = self._compute_correlations(deltas)
pos_candidates, neg_candidates = self._select_candidates(results_df, min_rho, min_delta_range, top_n)
self.logger.info(f"篩選完成: {len(pos_candidates)} 正相關, {len(neg_candidates)} 負相關")
return FeatureSelectionResult(positive_genes=pos_candidates, negative_genes=neg_candidates, full_report=results_df)
def _compute_deltas(self, data_df: pd.DataFrame) -> pd.DataFrame:
deltas = pd.DataFrame(index=data_df.index)
for pheno in self.config.phenotype_order:
at_col, mt_col = f"{pheno}AT", f"{pheno}MT"
if at_col in data_df.columns and mt_col in data_df.columns:
deltas[f"Delta_{pheno}"] = data_df[at_col] - data_df[mt_col]
else:
self.logger.warning(f"缺少 {at_col}/{mt_col}，無法計算 Delta_{pheno}")
return deltas
def _compute_correlations(self, deltas: pd.DataFrame) -> pd.DataFrame:
delta_cols = [c for c in self.config.delta_columns if c in deltas.columns]
target_ranks = list(self.config.phenotype_ranks)[:len(delta_cols)]
if len(delta_cols) != len(self.config.phenotype_ranks):
self.logger.warning("Delta 欄位數與表型等級數不匹配，已依可用欄位對齊等級")
if not delta_cols:
self.logger.error("沒有可用的 Delta 欄位，無法計算相關性")
return pd.DataFrame(columns=["Spearman_rho", "Delta_Span"]).rename_axis("Feature")
records = []
for idx, row in deltas.iterrows():
sample_deltas = [row[c] for c in delta_cols]
if len(set(sample_deltas)) == 1:
rho = 0.0
else:
rho, _ = spearmanr(sample_deltas, target_ranks)
if np.isnan(rho):
rho = 0.0
delta_span = max(sample_deltas) - min(sample_deltas)
record = {"Feature": idx, "Spearman_rho": rho, "Delta_Span": delta_span}
record.update(dict(zip(delta_cols, sample_deltas)))
records.append(record)
return pd.DataFrame(records).set_index("Feature")
def _select_candidates(self, res_df, min_rho, min_delta_range, top_n) -> Tuple[pd.DataFrame, pd.DataFrame]:
if res_df.empty:
return res_df.copy(), res_df.copy()
pos = res_df[
(res_df["Spearman_rho"] >= min_rho) & (res_df["Delta_Span"] >= min_delta_range)
].sort_values("Delta_Span", ascending=False)
neg = res_df[
(res_df["Spearman_rho"] <= -min_rho) & (res_df["Delta_Span"] >= min_delta_range)
].sort_values("Delta_Span", ascending=False)
if top_n is not None:
pos = pos.head(top_n)
neg = neg.head(top_n)
return pos, neg
class EnrichmentEngine:
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
def analyze(self, gene_list: List[str], group_name: str = "Unnamed") -> EnrichmentResult:
cleaned = self._clean_gene_list(gene_list)
if not cleaned:
self.logger.error(f"[{group_name}] 基因清單為空")
return EnrichmentResult(None, [], group_name)
self.logger.info(f"[{group_name}] 富集分析: {len(cleaned)} 基因")
df = self._try_local_gmt(cleaned)
if df is not None:
return EnrichmentResult(df, cleaned, group_name)
df = self._try_online(cleaned)
if df is not None:
return EnrichmentResult(df, cleaned, group_name)
self.logger.info(f"[{group_name}] 使用 mock 結果演示")
df = self._mock_results(cleaned)
return EnrichmentResult(df, cleaned, group_name)
def _clean_gene_list(self, gene_list: List[str]) -> List[str]:
cleaned = []
for g in gene_list:
if pd.notna(g):
g_str = str(g).strip()
if g_str and g_str.lower() != 'nan':
cleaned.append(g_str)
return sorted(set(cleaned))
def _try_local_gmt(self, gene_list: List[str]) -> Optional[pd.DataFrame]:
if not GSEAPY_AVAILABLE:
return None
gmt_path = self.config.default_gmt_path
if not gmt_path or not os.path.exists(gmt_path):
return None
try:
enr = gp.enrichr(gene_list=gene_list, gene_sets=gmt_path, outdir=None)
if enr is not None and enr.results is not None and not enr.results.empty:
self.logger.info("本地 GMT 富集分析成功")
return enr.results
except Exception as e:
self.logger.warning(f"本地 GMT 失敗: {e}")
return None
def _try_online(self, gene_list: List[str]) -> Optional[pd.DataFrame]:
if not GSEAPY_AVAILABLE:
return None
try:
enr = gp.enrichr(
gene_list=gene_list,
gene_sets=['GO_Biological_Process_2023', 'KEGG_2021_Human'],
outdir=None
)
if enr is not None and enr.results is not None and not enr.results.empty:
self.logger.info("線上 Enrichr 分析成功")
return enr.results
except Exception as e:
self.logger.warning(f"線上 API 失敗: {e}")
return None
class VisualizationSystem:
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
self._setup_style()
def _setup_style(self):
sns.set_theme(style="white", font_scale=1.0)
plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'Arial', 'PingFang TC', 'Microsoft JhengHei']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 100
def _save_fig(self, fig: plt.Figure, filename: str) -> str:
os.makedirs(self.config.output_dir, exist_ok=True)
path = os.path.join(self.config.output_dir, f"{filename}.{self.config.fig_format}")
fig.savefig(path, dpi=self.config.dpi, bbox_inches='tight', pad_inches=0.3)
self.logger.info(f"圖片已儲存: {path}")
plt.show(fig)
return path
def plot_chord_cnet(self, enrich_df, gene_directions: Dict[str, str], title_prefix="Combined",
top_n=5, max_genes=30, filename="cnetplot_chord") -> Optional[str]:
"""繪製 Chord Cnetplot（通路-基因網路圖）。
Args:
max_genes: 弧形上最多顯示的基因數。基因太多會讓標籤重疊、難以閱讀，
因此當候選基因數超過此值時，只保留「連結最多通路」的前
max_genes 個基因（即最具代表性的樞紐基因），其餘基因仍計入
通路統計但不畫在圖上。
"""
if enrich_df is None or enrich_df.empty:
self.logger.warning(f"[{title_prefix}] 富集結果為空，跳過 Chord 圖")
return None
top_pathways = enrich_df.head(top_n).copy()
pathway_genes_map: Dict[str, List[str]] = {}
gene_degree: Dict[str, int] = {}
for _, row in top_pathways.iterrows():
term = row['Term']
genes = [g.strip() for g in str(row['Genes']).split(';') if g.strip()]
pathway_genes_map[term] = genes
for g in genes:
gene_degree[g] = gene_degree.get(g, 0) + 1
pathways = list(pathway_genes_map.keys())
if not pathways or not gene_degree:
self.logger.warning("未解析出有效關係")
return None
n_all_genes = len(gene_degree)
if n_all_genes > max_genes:
# 依「連結的通路數」由多到少排序，保留最具代表性（樞紐）的基因；
# 通路數相同則依字母排序，確保結果穩定可重現。
ranked_genes = sorted(gene_degree.items(), key=lambda kv: (-kv[1], kv[0]))
genes_list = sorted(g for g, _ in ranked_genes[:max_genes])
self.logger.info(
f"[{title_prefix}] 基因數 {n_all_genes} 超過顯示上限 {max_genes}，"
f"僅保留連結最多通路的 {len(genes_list)} 個基因"
)
# 只保留有畫出來的基因的連線，避免通路節點連到不存在的基因節點
kept_set = set(genes_list)
pathway_genes_map = {
term: [g for g in genes if g in kept_set]
for term, genes in pathway_genes_map.items()
}
else:
genes_list = sorted(gene_degree.keys())
fig, ax = plt.subplots(figsize=(14, 14))
ax.set_aspect("equal")
ax.axis("off")
ax.set_xlim(-1.4, 1.4)
ax.set_ylim(-1.4, 1.55)
radius = 1.0
n_pathways, n_genes = len(pathways), len(genes_list)
# 基因數愈少，標籤與節點可以愈大，維持圖面可讀性
gene_marker_size = 160 if n_genes <= 20 else max(70, 160 - (n_genes - 20) * 3)
gene_font_size = 10 if n_genes <= 20 else max(6.5, 10 - (n_genes - 20) * 0.12)
pathway_angles = np.linspace(0.20 * np.pi, 0.80 * np.pi, n_pathways) if n_pathways > 1 else [np.pi / 2]
gene_angles = np.linspace(1.20 * np.pi, 1.80 * np.pi, n_genes) if n_genes > 1 else [1.5 * np.pi]
node_coords: Dict[str, Tuple[float, float, float]] = {}
for i, term in enumerate(pathways):
angle = pathway_angles[i]
node_coords[term] = (radius * np.cos(angle), radius * np.sin(angle), angle)
for i, gene in enumerate(genes_list):
angle = gene_angles[i]
node_coords[gene] = (radius * np.cos(angle), radius * np.sin(angle), angle)
# 使用局部別名，避免與 pathlib.Path 混淆
MPath = mpath.Path
for term, genes in pathway_genes_map.items():
px, py, _ = node_coords[term]
for gene in genes:
if gene not in node_coords:
continue
gx, gy, _ = node_coords[gene]
direction = gene_directions.get(gene, "Positive")
chord_color = "#E64B35" if direction == "Positive" else "#4DBBD5"
path_data = [
(MPath.MOVETO, (px, py)),
(MPath.CURVE3, (0.0, 0.0)),
(MPath.CURVE3, (gx, gy))
]
codes, verts = zip(*path_data)
patch = mpatches.PathPatch(
MPath(verts, codes), facecolor="none",
edgecolor=chord_color, alpha=0.30, linewidth=1.4
)
ax.add_patch(patch)
def _label(x, y, angle, text, **kwargs):
deg = np.degrees(angle) % 360
ha = "right" if 90 < deg < 270 else "left"
rot_deg = deg if ha == "left" else deg - 180
ax.text(x, y, text, ha=ha, va="center", rotation=rot_deg, rotation_mode="anchor", **kwargs)
for term in pathways:
x, y, angle = node_coords[term]
ax.scatter(x, y, s=380, color="#00A087", zorder=4, edgecolors="black", linewidth=1.2)
wrapped = "\n".join(textwrap.wrap(term.replace('_', ' '), width=26))
_label(x * 1.18, y * 1.18, angle, wrapped, fontsize=9.5, fontweight="bold")
for gene in genes_list:
x, y, angle = node_coords[gene]
direction = gene_directions.get(gene, "Positive")
color = "#E64B35" if direction == "Positive" else "#4DBBD5"
ax.scatter(x, y, s=gene_marker_size, color=color, zorder=4, edgecolors="black", linewidth=0.8)
_label(x * 1.12, y * 1.12, angle, gene, fontsize=gene_font_size, fontweight="bold", color=color)
handles = [mpatches.Patch(color="#00A087", label="Enriched Pathway / GO Term")]
if any(gene_directions.get(g) == "Positive" for g in genes_list):
handles.append(mpatches.Patch(color="#E64B35", label="Positive (Pro-development)"))
if any(gene_directions.get(g) == "Negative" for g in genes_list):
handles.append(mpatches.Patch(color="#4DBBD5", label="Negative (Inhibitor)"))
ax.legend(handles=handles, loc="lower center", bbox_to_anchor=(0.5, -0.12),
ncol=len(handles), frameon=True, fontsize=10, title="Category & Regulation")
ax.set_title(f"{title_prefix} Pathways-Genes Network", fontsize=18, fontweight="bold", pad=40, y=1.15)
return self._save_fig(fig, filename)
def plot_volcano(self, feature_result: FeatureSelectionResult, filename="volcano_plot") -> Optional[str]:
df = feature_result.full_report.copy()
if df.empty:
self.logger.warning("特徵報告為空，跳過火山圖")
return None
fig, ax = plt.subplots(figsize=(10, 8))
pos_idx = set(feature_result.positive_genes.index)
neg_idx = set(feature_result.negative_genes.index)
is_pos = df.index.to_series().isin(pos_idx)
is_neg = df.index.to_series().isin(neg_idx)
is_other = ~(is_pos | is_neg)
ax.scatter(df.loc[is_other, "Spearman_rho"], df.loc[is_other, "Delta_Span"],
c="gray", alpha=0.3, s=20, label="Non-significant")
ax.scatter(df.loc[is_pos, "Spearman_rho"], df.loc[is_pos, "Delta_Span"],
c="#E64B35", alpha=0.7, s=60, label="Positive trend")
ax.scatter(df.loc[is_neg, "Spearman_rho"], df.loc[is_neg, "Delta_Span"],
c="#4DBBD5", alpha=0.7, s=60, label="Negative trend")
ax.axvline(x=self.config.min_rho, color='red', linestyle='--', alpha=0.5)
ax.axvline(x=-self.config.min_rho, color='blue', linestyle='--', alpha=0.5)
ax.axhline(y=self.config.min_delta_range, color='green', linestyle='--', alpha=0.5)
ax.set_xlabel("Spearman rho", fontsize=12)
ax.set_ylabel("Delta Span", fontsize=12)
ax.set_title("Phenotype-guided Feature Selection", fontsize=14, fontweight='bold')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
# 只標註實際入選的正/負候選基因，並依 Delta_Span 取前 20 名，避免圖面過度擁擠
labeled = df.loc[is_pos | is_neg].sort_values("Delta_Span", ascending=False).head(20)
texts = [ax.text(row["Spearman_rho"], row["Delta_Span"], idx, fontsize=8, alpha=0.8)
for idx, row in labeled.iterrows()]
if ADJUSTTEXT_AVAILABLE and texts:
adjust_text(texts, arrowprops=dict(arrowstyle='->', color='gray', alpha=0.5))
return self._save_fig(fig, filename)
def plot_regulation_barplot(
self,
feature_result: FeatureSelectionResult,
direction: str,
top_n: int = 20,
filename: str = "regulation_barplot"
) -> Optional[str]:
"""繪製「上調」或「下調」候選基因的長條圖（只畫單一方向，兩方向各自獨立成圖）。
Args:
direction: "positive"（上調，Delta AT>MT 且與表型等級正相關）
或 "negative"（下調，負相關）。
top_n: 依 Delta_Span 由大到小，最多顯示幾個基因。
"""
if direction not in ("positive", "negative"):
raise ValueError(f"direction 必須是 'positive' 或 'negative'，收到: {direction!r}")
candidates = feature_result.positive_genes if direction == "positive" else feature_result.negative_genes
label = "Up-regulated" if direction == "positive" else "Down-regulated"
color = "#E64B35" if direction == "positive" else "#4DBBD5"
if candidates.empty:
self.logger.warning(f"[{label}] 無候選基因，跳過長條圖")
return None
plot_df = candidates.sort_values("Delta_Span", ascending=False).head(top_n)
fig, ax = plt.subplots(figsize=(10, max(6, len(plot_df) * 0.35)))
bars = ax.barh(range(len(plot_df)), plot_df["Delta_Span"], color=color, alpha=0.85,
edgecolor="black", linewidth=0.5)
ax.set_yticks(range(len(plot_df)))
ax.set_yticklabels(plot_df.index, fontsize=9)
ax.set_xlabel("Delta Span (asinh scale, |AT - MT| range across phenotypes)")
ax.set_title(
f"{label} Candidate Genes (Top {len(plot_df)} by Delta Span)",
fontsize=14, fontweight='bold'
)
ax.invert_yaxis()
ax.grid(True, axis='x', alpha=0.3)
for i, val in enumerate(plot_df["Delta_Span"]):
ax.text(val, i, f" {val:.2f}", va='center', fontsize=8)
return self._save_fig(fig, filename)
def plot_heatmap(self, expr_matrix: pd.DataFrame, feature_result: FeatureSelectionResult, filename="heatmap") -> Optional[str]:
selected_genes = list(feature_result.positive_genes.index) + list(feature_result.negative_genes.index)
if not selected_genes:
self.logger.warning("無候選基因，跳過熱力圖")
return None
if len(selected_genes) > 50:
selected_genes = selected_genes[:50]
self.logger.info("熱力圖限制顯示前 50 個基因")
missing = [g for g in selected_genes if g not in expr_matrix.index]
if missing:
self.logger.warning(f"{len(missing)} 個候選基因不在表達矩陣中（可能已被 asinh 篩選排除），將略過")
selected_genes = [g for g in selected_genes if g in expr_matrix.index]
if not selected_genes:
self.logger.warning("篩選後無可繪製之候選基因，跳過熱力圖")
return None
plot_data = expr_matrix.loc[selected_genes]
row_colors = ["#E64B35" if g in feature_result.positive_genes.index else "#4DBBD5"
for g in selected_genes]
fig, ax = plt.subplots(figsize=(10, max(6, len(selected_genes) * 0.3)))
sns.heatmap(plot_data, cmap="RdBu_r", center=plot_data.values.mean(),
annot=False, fmt=".2f", linewidths=0.5,
cbar_kws={"label": "asinh(expression)"}, ax=ax)
for i, color in enumerate(row_colors):
ax.add_patch(plt.Rectangle((-0.5, i), 0.1, 1, facecolor=color,
transform=ax.get_yaxis_transform(), clip_on=False))
ax.set_title("Candidate Gene Expression Heatmap", fontsize=14, fontweight='bold')
ax.set_xlabel("Samples")
ax.set_ylabel("Genes")
legend_elements = [
mpatches.Patch(facecolor="#E64B35", label="Positive trend"),
mpatches.Patch(facecolor="#4DBBD5", label="Negative trend")
]
ax.legend(handles=legend_elements, loc='upper left', bbox_to_anchor=(1.15, 1))
return self._save_fig(fig, filename)
def plot_pca(self, expr_matrix: pd.DataFrame, filename="pca_plot") -> Optional[str]:
if expr_matrix.shape[1] < 2:
self.logger.warning("樣本數不足以計算 PCA (需要至少 2 個樣本)，跳過 PCA 圖")
return None
data = expr_matrix.T
data.columns = data.columns.astype(str)
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)
n_components = min(2, data_scaled.shape[0] - 1, data_scaled.shape[1])
if n_components < 2:
self.logger.warning("樣本或特徵數不足以計算 2 個主成分，跳過 PCA 圖")
return None
pca = PCA(n_components=2)
pcs = pca.fit_transform(data_scaled)
sample_ids = self.config.sample_ids
treatments = ['MT' if t == 0 else 'AT' for t in self.config.is_treated]
fig, ax = plt.subplots(figsize=(10, 8))
colors = {"MT": "#4DBBD5", "AT": "#E64B35"}
markers = {0: 'o', 1: 's', 2: '^'}
for i, (pc1, pc2) in enumerate(pcs):
marker = markers.get(self.config.individual[i], 'D')
ax.scatter(pc1, pc2, c=colors[treatments[i]], marker=marker,
s=200, alpha=0.7, edgecolors='black')
ax.annotate(sample_ids[i], (pc1, pc2), xytext=(5, 5),
textcoords='offset points', fontsize=9)
ax.set_xlabel(f"PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)")
ax.set_ylabel(f"PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)")
ax.set_title("PCA of Samples", fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='gray', linestyle='-', alpha=0.3)
ax.axvline(x=0, color='gray', linestyle='-', alpha=0.3)
from matplotlib.lines import Line2D
legend_elements = [
Line2D([0], [0], marker='o', color='w', markerfacecolor="#4DBBD5", markersize=10, label='MT (Untreated)'),
Line2D([0], [0], marker='o', color='w', markerfacecolor="#E64B35", markersize=10, label='AT (Treated)')
]
ax.legend(handles=legend_elements, loc='best')
return self._save_fig(fig, filename)
def plot_enrichment_bar(self, enrich_result: EnrichmentResult, filename="enrichment_bar") -> Optional[str]:
if not enrich_result.success:
return None
df = enrich_result.top_terms(15).copy()
fig, ax = plt.subplots(figsize=(10, max(6, len(df) * 0.4)))
safe_p = df['P-value'].clip(lower=1e-300)
df['neg_log_p'] = -np.log10(safe_p)
colors = ["#E64B35" if 'GO' in gs else "#00A087" for gs in df['Gene_set']]
bars = ax.barh(range(len(df)), df['neg_log_p'], color=colors, alpha=0.8)
ax.set_yticks(range(len(df)))
ax.set_yticklabels(df['Term'], fontsize=9)
ax.set_xlabel("-log10(P-value)")
ax.set_title(f"{enrich_result.group_name} - Top Enriched Terms", fontsize=14, fontweight='bold')
ax.axvline(x=-np.log10(0.05), color='red', linestyle='--', alpha=0.5, label='p=0.05')
ax.legend()
ax.invert_yaxis()
for i, (bar, val) in enumerate(zip(bars, df['neg_log_p'])):
ax.text(val + 0.1, i, f"{val:.2f}", va='center', fontsize=8)
return self._save_fig(fig, filename)
class ReportGenerator:
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
def generate(self, feature_result: FeatureSelectionResult, enrichment_results: Dict[str, EnrichmentResult], image_paths: Dict[str, str]) -> str:
os.makedirs(self.config.output_dir, exist_ok=True)
report_path = os.path.join(self.config.output_dir, "analysis_report.html")
html = self._build_html(feature_result, enrichment_results, image_paths)
with open(report_path, 'w', encoding='utf-8') as f:
f.write(html)
self.logger.info(f"報告已生成: {report_path}")
return report_path
def _build_html(self, feature_result, enrichment_results, image_paths) -> str:
summary = feature_result.summary()
img_tags = ""
for name, path in image_paths.items():
if path and os.path.exists(path):
rel_path = os.path.relpath(path, self.config.output_dir)
img_tags += f'<h3>{name}</h3><img src="{rel_path}" style="max-width:100%;"/><br><br>'
enrich_tables = ""
for group_name, result in enrichment_results.items():
if result.success:
table_html = result.dataframe.head(10).to_html(classes='table table-striped', index=False)
enrich_tables += f"<h3>{group_name} Enrichment</h3>{table_html}<br>"
html = f"""
<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>Ovary Analysis Report</title>
<style>
body {{ font-family: Arial, sans-serif; margin: 40px; background: #f5f5f5; }}
.container {{ max-width: 1200px; margin: 0 auto; background: white; padding: 30px; box-shadow: 0 0 10px rgba(0,0,0,0.1); }}
h1 {{ color: #2c3e50; border-bottom: 3px solid #3498db; padding-bottom: 10px; }}
h2 {{ color: #34495e; margin-top: 30px; }}
h3 {{ color: #7f8c8d; }}
.summary {{ background: #ecf0f1; padding: 20px; border-radius: 8px; margin: 20px 0; }}
.summary-item {{ margin: 8px 0; }}
table {{ border-collapse: collapse; width: 100%; margin: 15px 0; }}
th, td {{ border: 1px solid #ddd; padding: 8px; text-align: left; }}
th {{ background: #3498db; color: white; }}
tr:nth-child(even) {{ background: #f9f9f9; }}
img {{ border: 1px solid #ddd; border-radius: 4px; margin: 10px 0; }}
.footer {{ margin-top: 40px; color: #95a5a6; font-size: 0.9em; text-align: center; }}
</style>
</head>
<body>
<div class="container">
<h1>Ovary Tissue Transcriptome Analysis Report</h1>
<p><strong>Generated:</strong> {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}</p>
<h2>Analysis Summary</h2>
<div class="summary">
<div class="summary-item"><strong>Positive Candidates:</strong> {summary['n_positive']}</div>
<div class="summary-item"><strong>Negative Candidates:</strong> {summary['n_negative']}</div>
<div class="summary-item"><strong>Total Genes Tested:</strong> {summary['n_total_tested']}</div>
<div class="summary-item"><strong>Top Positive:</strong> {', '.join(summary['top_positive'][:5])}</div>
<div class="summary-item"><strong>Top Negative:</strong> {', '.join(summary['top_negative'][:5])}</div>
</div>
<h2>Visualizations</h2>
{img_tags}
<h2>Enrichment Analysis</h2>
{enrich_tables}
<div class="footer">
<p>Generated by OvaryAnalysis Pipeline v3 | Configuration: {self.config.phenotype_order}</p>
</div>
</div>
</body>
</html>
"""
return html
class ResultExporter:
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
def export_feature_result(self, result: FeatureSelectionResult, prefix="features") -> Dict[str, str]:
os.makedirs(self.config.output_dir, exist_ok=True)
paths = {}
excel_path = os.path.join(self.config.output_dir, f"{prefix}.xlsx")
with pd.ExcelWriter(excel_path, engine='openpyxl') as writer:
result.positive_genes.to_excel(writer, sheet_name='Positive')
result.negative_genes.to_excel(writer, sheet_name='Negative')
result.full_report.to_excel(writer, sheet_name='Full_Report')
paths['excel'] = excel_path
self.logger.info(f"特徵結果已匯出: {excel_path}")
for name, df in [('positive', result.positive_genes), ('negative', result.negative_genes)]:
if not df.empty:
csv_path = os.path.join(self.config.output_dir, f"{prefix}_{name}.csv")
df.to_csv(csv_path)
paths[f'csv_{name}'] = csv_path
return paths
def export_enrichment_result(self, result: EnrichmentResult, prefix="enrichment") -> Optional[str]:
if not result.success:
return None
os.makedirs(self.config.output_dir, exist_ok=True)
path = os.path.join(self.config.output_dir, f"{prefix}_{result.group_name}.csv")
result.dataframe.to_csv(path, index=False)
self.logger.info(f"富集結果已匯出: {path}")
return path
class AnalysisPipeline:
def __init__(self, config: AnalysisConfig, log: logging.Logger):
self.config = config
self.logger = log
self.data_loader = DataLoader(config, log)
self.asinh_filter = AsinhFilter(log, use_quantile=config.asinh_use_quantile)
self.feature_selector = PhenotypeFeatureSelector(config, log)
self.enrichment_engine = EnrichmentEngine(config, log)
self.visualizer = VisualizationSystem(config, log)
self.report_generator = ReportGenerator(config, log)
self.exporter = ResultExporter(config, log)
self.results: Dict[str, Any] = {}
self.image_paths: Dict[str, str] = {}
def run(self, data_path: Optional[str] = None) -> Dict[str, Any]:
self.logger.info("=" * 60)
self.logger.info("開始卵巢組織轉錄組分析流程 v3")
self.logger.info("=" * 60)
try:
expr_matrix = self._step_load_data(data_path)
expr_filtered = self._step_asinh_filter(expr_matrix)
self._step_global_enrichment(expr_filtered)
feature_result = self._step_feature_selection(expr_filtered)
self._step_visualization(expr_matrix, expr_filtered, feature_result)
enrichment_results = self._step_group_enrichment(feature_result)
self._step_export(feature_result, enrichment_results)
report_path = self._step_generate_report(feature_result, enrichment_results)
self.logger.info("=" * 60)
self.logger.info("分析流程完成")
self.logger.info(f"輸出目錄: {self.config.output_dir}")
self.logger.info("=" * 60)
return {
"success": True,
"output_dir": self.config.output_dir,
"report_path": report_path,
"feature_result": feature_result,
"enrichment_results": enrichment_results,
"image_paths": self.image_paths,
}
except Exception as e:
self.logger.error(f"分析流程失敗: {e}", exc_info=True)
return {"success": False, "error": str(e)}
def _step_load_data(self, data_path: Optional[str]) -> pd.DataFrame:
self.logger.info("[Step 1/8] 載入數據...")
expr_matrix = self.data_loader.load(data_path)
self.logger.info(f"數據載入完成: {expr_matrix.shape}")
return expr_matrix
def _step_asinh_filter(self, expr_matrix: pd.DataFrame) -> pd.DataFrame:
self.logger.info("[Step 2/8] asinh 數值區間篩選...")
special_genes = ["CYP19A1", "STAR", "LHCGR", "MAPK1", "CDK1", "CCNB1", "TGFBR1", "CASP3", "FOXO1"]
comparison = self.asinh_filter.compare_modes(
expr_matrix, special_genes, self.config.asinh_lower, self.config.asinh_upper
)
print("\n========== asinh 篩選模式比較 ==========")
print(comparison.to_string(index=False))
print("=========================================\n")
expr_filtered = self.asinh_filter.filter(
expr_matrix, lower=self.config.asinh_lower,
upper=self.config.asinh_upper, mode=self.config.asinh_filter_mode
)
self.results['asinh_comparison'] = comparison
return expr_filtered
def _step_global_enrichment(self, expr_filtered: pd.DataFrame) -> None:
self.logger.info("[Step 3/8] 全域富集分析...")
gene_list = expr_filtered.index.tolist()
result = self.enrichment_engine.analyze(gene_list, "Global")
self.results['global_enrichment'] = result
if result.success:
print("\n========== 全域富集分析結果 ==========")
print(result.dataframe.head(10).to_string())
print("=======================================\n")
def _step_feature_selection(self, expr_filtered: pd.DataFrame) -> FeatureSelectionResult:
"""步驟 4: 表型引導特徵篩選。
注意：此步驟只在 asinh 篩選後的資料 (expr_filtered) 上計算 Delta 與
Spearman rho，不再依平均表現值做二次篩選，避免與 Step 2 的 asinh
篩選重複套用不同標準而造成篩選結果難以解讀。
"""
self.logger.info("[Step 4/8] 表型引導特徵篩選...")
result = self.feature_selector.select(
expr_filtered,
min_delta_range=self.config.min_delta_range,
min_rho=self.config.min_rho,
top_n=self.config.top_n_features,
)
print("\n========== 特徵篩選結果 ==========")
print(f"正相關候選: {len(result.positive_genes)}")
print(f"負相關候選: {len(result.negative_genes)}")
if not result.positive_genes.empty:
print("\nTop 5 正相關基因:")
print(result.positive_genes.head(5).to_string())
if not result.negative_genes.empty:
print("\nTop 5 負相關基因:")
print(result.negative_genes.head(5).to_string())
print("===================================\n")
return result
def _step_visualization(self, expr_matrix, expr_filtered, feature_result) -> None:
self.logger.info("[Step 5/8] 生成可視化圖表...")
path = self.visualizer.plot_pca(expr_matrix, "pca_plot")
if path:
self.image_paths['PCA'] = path
path = self.visualizer.plot_volcano(feature_result, "volcano_plot")
if path:
self.image_paths['Volcano'] = path
path = self.visualizer.plot_heatmap(expr_filtered, feature_result, "heatmap")
if path:
self.image_paths['Heatmap'] = path
# 上調 / 下調候選基因分別各自畫一張長條圖，方便單獨檢視每個方向的結果
path = self.visualizer.plot_regulation_barplot(feature_result, "positive", filename="barplot_upregulated")
if path:
self.image_paths['Barplot_Upregulated'] = path
path = self.visualizer.plot_regulation_barplot(feature_result, "negative", filename="barplot_downregulated")
if path:
self.image_paths['Barplot_Downregulated'] = path
def _step_group_enrichment(self, feature_result: FeatureSelectionResult) -> Dict[str, EnrichmentResult]:
self.logger.info("[Step 6/8] 分組富集分析...")
pos_list = feature_result.positive_genes.index.tolist()
neg_list = feature_result.negative_genes.index.tolist()
gene_directions = {g: "Positive" for g in pos_list}
gene_directions.update({g: "Negative" for g in neg_list})
enrichment_results = {}
groups = [
("Up-regulated", pos_list, "cnetplot_up"),
("Down-regulated", neg_list, "cnetplot_down"),
]
for group_name, genes, filename in groups:
if not genes:
self.logger.warning(f"[{group_name}] 無基因，跳過")
continue
result = self.enrichment_engine.analyze(genes, group_name)
enrichment_results[group_name] = result
if result.success:
print(f"\n========== {group_name} 富集分析結果 ==========")
print(result.dataframe.head(10).to_string())
print("============================================\n")
path = self.visualizer.plot_chord_cnet(
result.dataframe, gene_directions, group_name, top_n=10, filename=filename
)
if path:
self.image_paths[f'Chord_{group_name}'] = path
path = self.visualizer.plot_enrichment_bar(
result, f"enrichment_bar_{group_name.lower().replace('-', '_')}"
)
if path:
self.image_paths[f'Bar_{group_name}'] = path
combined_genes = list(set(pos_list + neg_list))
if combined_genes:
combined_result = self.enrichment_engine.analyze(combined_genes, "Combined")
enrichment_results["Combined"] = combined_result
if combined_result.success:
path = self.visualizer.plot_chord_cnet(
combined_result.dataframe, gene_directions, "Combined", top_n=10, filename="cnetplot_combined"
)
if path:
self.image_paths['Chord_Combined'] = path
return enrichment_results
def _step_export(self, feature_result: FeatureSelectionResult, enrichment_results: Dict[str, EnrichmentResult]) -> None:
self.logger.info("[Step 7/8] 匯出結果...")
self.exporter.export_feature_result(feature_result, "features")
for group_name, result in enrichment_results.items():
self.exporter.export_enrichment_result(result, "enrichment")
def _step_generate_report(self, feature_result, enrichment_results) -> Optional[str]:
self.logger.info("[Step 8/8] 生成 HTML 報告...")
return self.report_generator.generate(feature_result, enrichment_results, self.image_paths)
def main():
logger = setup_logging(level=logging.INFO)
config = AnalysisConfig()
os.makedirs(config.output_dir, exist_ok=True)
config.save(os.path.join(config.output_dir, "config.json"))
pipeline = AnalysisPipeline(config, logger)
results = pipeline.run()
if results["success"]:
print("\n" + "=" * 60)
print("分析成功完成！")
print(f"輸出目錄: {results['output_dir']}")
print(f"報告檔案: {results['report_path']}")
print("=" * 60)
else:
print(f"\n分析失敗: {results.get('error', 'Unknown error')}")
sys.exit(1)
if __name__ == "__main__":
main()
No comments:
Post a Comment