不同組織進行基因表達譜（Transcriptome）分析時，科學家驚訝地發現：健康的卵巢組織竟然高度富集（Enrichment）了原本被認為是神經系統特有的通路，如阿茲海默症（AD）、巴金森症（PD）、亨丁頓舞蹈症（HD）等通路。
最有力的證據來自 GTEx (Genotype-Tissue Expression) 計畫，這是全球最大的人體組織基因表達數據庫。 關鍵文獻： Lopes-Ramos et al. (2020), "Sex differences in gene expression and regulatory networks across 29 human tissues", Cell Reports. 證據內容： 該研究在分析各組織的調控網絡時發現，卵巢組織中的基因表達模式與大腦高度相似。
在通路分析（GSEA）中，卵巢顯著富集了 "Alzheimer's disease"、"Parkinson's disease" 與 "Huntington's disease" 等 KEGG 通路。
以下是鰻魚卵巢轉錄體以NFM的分析：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import seaborn as sns
import requests
import time
from pathlib import Path
from typing import List, Tuple, Optional, Dict
from gseapy import enrichr
import warnings
warnings.filterwarnings('ignore')
# --- 1. CNA 核心分析類別 ---
class CNA_Dual_Analyzer:
def __init__(self, species: int = 9606, max_api_genes: int = 200):
self.species = species
self.max_api_genes = max_api_genes
self.api_url = "https://string-db.org/api/json/network"
def _safe_string_api(self, genes: List[str]) -> pd.DataFrame:
if len(genes) < 2: return pd.DataFrame()
batch_genes = genes[:self.max_api_genes]
params = {
"identifiers": "%0d".join(batch_genes),
"species": self.species,
"caller_identity": "cna_v2_1"
}
try:
resp = requests.post(self.api_url, data=params, timeout=30)
resp.raise_for_status()
data = resp.json()
if not data: return pd.DataFrame()
df = pd.DataFrame(data)
return df[['preferredName_A', 'preferredName_B', 'score']].dropna()
except Exception as e:
print(f"⚠️ STRING API失敗: {e}")
return pd.DataFrame()
def compute_cna_metrics(self, cluster_data: pd.DataFrame, at_cols: List[str], mt_cols: List[str]) -> pd.DataFrame:
metrics = []
for _, row in cluster_data.iterrows():
gene = str(row.get('Gene_symbol', row.index[0])).strip().upper()
if not gene or gene in ['NAN', 'NONE', '']: continue
mean_at = row[at_cols].mean() if at_cols else 0.0
mean_mt = row[mt_cols].mean() if mt_cols else 0.0
asinh_diff = np.arcsinh(mean_mt) - np.arcsinh(mean_at)
log2fc = np.log2((mean_mt + 1e-6) / (mean_at + 1e-6))
metrics.append({
'Gene': gene, 'Mean_AT': float(mean_at), 'Mean_MT': float(mean_mt),
'asinh_Diff': float(asinh_diff), 'Log2FC': float(log2fc),
'Expr_Change': 'UP' if asinh_diff > 0 else 'DOWN'
})
return pd.DataFrame(metrics)
def build_control_network(self, genes: List[str], metrics_df: pd.DataFrame) -> Tuple[Optional[nx.Graph], pd.DataFrame]:
ppi_df = self._safe_string_api(genes)
if ppi_df.empty: return None, pd.DataFrame()
G = nx.Graph()
for _, row in ppi_df.iterrows():
G.add_edge(row['preferredName_A'], row['preferredName_B'], weight=float(row['score']) / 1000)
centrality = nx.degree_centrality(G)
results = []
for node in G.nodes():
gene_metrics = metrics_df[metrics_df['Gene'] == node]
asinh_diff = gene_metrics['asinh_Diff'].iloc[0] if not gene_metrics.empty else 0.0
ctrl_score = centrality.get(node, 0) * (abs(asinh_diff) + 1)
results.append({
'Gene': node, 'Control_Score': ctrl_score, 'Centrality': centrality.get(node, 0),
'asinh_Diff': asinh_diff, 'Log2FC': gene_metrics['Log2FC'].iloc[0] if not gene_metrics.empty else 0.0,
'Expr_Change': gene_metrics['Expr_Change'].iloc[0] if not gene_metrics.empty else 'NONE'
})
stats_df = pd.DataFrame(results).sort_values('Control_Score', ascending=False)
return G, stats_df
# --- 2. 視覺化類別 (整合 CNetplot) ---
class NetworkVisualizer:
@staticmethod
def create_cna_network(G, stats_df, output_path, cluster_id):
"""核心 PPI 控制網絡圖"""
if G is None or stats_df.empty: return
plt.figure(figsize=(12, 10), dpi=100)
pos = nx.spring_layout(G, k=0.3, iterations=50)
sizes = 100 + stats_df.set_index('Gene').reindex(G.nodes())['Control_Score'] * 2000
colors = stats_df.set_index('Gene').reindex(G.nodes())['asinh_Diff']
nx.draw_networkx_nodes(G, pos, node_size=sizes, node_color=colors, cmap='RdYlBu_r', alpha=0.8)
nx.draw_networkx_edges(G, pos, alpha=0.2)
top_genes = stats_df.head(20)['Gene'].tolist()
nx.draw_networkx_labels(G, pos, {n: n for n in G.nodes() if n in top_genes}, font_size=8)
plt.title(f"CNA Network: Cluster {cluster_id}")
plt.savefig(output_path, bbox_inches='tight')
plt.close()
@staticmethod
def create_cnet_plot(stats_df: pd.DataFrame, output_path: Path, cluster_id: str):
"""生成 CNetplot (Gene-Pathway Network)"""
print(f"🧬 正在為 Cluster {cluster_id} 生成 CNetplot...")
genes = stats_df.head(100)['Gene'].tolist() # 取控制力前100的基因
try:
enr = enrichr(gene_list=genes, gene_sets='KEGG_2021_Human', organism='Human', outdir=None)
res = enr.results.nsmallest(5, 'Adjusted P-value') # 取前5個通路
if res.empty: return
G = nx.Graph()
# 建立基因與通路的連接
gene_colors = dict(zip(stats_df['Gene'], stats_df['asinh_Diff']))
for _, row in res.iterrows():
pathway = row['Term'].split('__')[-1]
G.add_node(pathway, type='pathway')
path_genes = row['Genes'].split(';')
for g in path_genes:
if g in genes:
G.add_node(g, type='gene')
G.add_edge(pathway, g)
plt.figure(figsize=(14, 10))
pos = nx.spring_layout(G, k=0.5)
# 繪製路徑節點 (大橘色)
path_nodes = [n for n, d in G.nodes(data=True) if d['type'] == 'pathway']
nx.draw_networkx_nodes(G, pos, nodelist=path_nodes, node_color='#FF7F0E',
node_size=1200, node_shape='p', label='Pathways')
# 繪製基因節點 (顏色代表 asinh_Diff)
gene_nodes = [n for n, d in G.nodes(data=True) if d['type'] == 'gene']
g_colors = [gene_colors.get(n, 0) for n in gene_nodes]
g_sizes = [stats_df[stats_df['Gene']==n]['Control_Score'].iloc[0]*800 + 100 for n in gene_nodes]
nodes = nx.draw_networkx_nodes(G, pos, nodelist=gene_nodes, node_color=g_colors,
node_size=g_sizes, cmap='RdYlBu_r', label='Genes')
nx.draw_networkx_edges(G, pos, alpha=0.3, edge_color='gray')
nx.draw_networkx_labels(G, pos, font_size=9, font_weight='bold')
plt.colorbar(nodes, label='asinh Difference')
plt.title(f"CNetplot: Top Pathways & Control Genes (Cluster {cluster_id})")
plt.legend(scatterpoints=1)
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.close()
except Exception as e:
print(f"⚠️ CNetplot失敗: {e}")
# --- 3. 主程序 ---
def main_pipeline(excel_path: Path):
output_dir = Path.home() / "Desktop" / "CNA_Analysis_CNet"
output_dir.mkdir(exist_ok=True)
df = pd.read_excel(excel_path)
at_cols = [c for c in df.columns if 'AT' in str(c).upper()]
mt_cols = [c for c in df.columns if 'MT' in str(c).upper()]
gene_col = next(c for c in ['Gene_symbol', 'GeneID'] if c in df.columns)
analyzer = CNA_Dual_Analyzer()
visualizer = NetworkVisualizer()
clusters = sorted(df['leiden'].unique().astype(str))
for cid in clusters:
print(f"\n🚀 處理 Cluster {cid}...")
cluster_df = df[df['leiden'].astype(str) == cid]
cluster_dir = output_dir / f"Cluster_{cid}"
cluster_dir.mkdir(exist_ok=True)
# 1. 計算指標
metrics_df = analyzer.compute_cna_metrics(cluster_df, at_cols, mt_cols)
# 2. 構建 CNA 網絡
genes = metrics_df['Gene'].tolist()
G, stats_df = analyzer.build_control_network(genes, metrics_df)
if not stats_df.empty:
stats_df.to_excel(cluster_dir / "CNA_Report.xlsx", index=False)
# 3. 繪製 PPI 控制圖
visualizer.create_cna_network(G, stats_df, cluster_dir / "CNA_Network.png", cid)
# 4. 繪製 CNetplot (路徑-基因交互圖)
visualizer.create_cnet_plot(stats_df, cluster_dir / "CNetplot.png", cid)
time.sleep(1)
if __name__ == "__main__":
file_path = Path.home() / "Desktop" / "Aja 21649_LeidCom_clustered.xlsx"
main_pipeline(file_path)
import scanpy as sc
import anndata as ad
import numpy as np0
import torch
import torch.nn as nn
import torch.optim as optim
# 1. 讀取原始檔案
AData = sc.read_h5ad("6ovaries_analysis_results.h5ad")
adata = sc.read_h5ad("6ovaries_analysis_results.h5ad")
# 2. 檢查目前結構
# 如果讀進來已經是 6 x 20899，就不要執行 .T
# 如果讀進來是 20899 x 6，才執行一次 .T 把它轉正
if adata.n_obs > adata.n_vars:
adata = adata
print(f"數據維度: {adata.n_obs} 樣本, {adata.n_vars} 基因")
print("目前的 obs 欄位：", adata.obs.columns)
# 標準化與對數轉換 (深度學習對數值範圍敏感)
sc.pp.log1p(adata)
sc.pp.scale(adata, max_value=10) # 建議增加 scale 步驟
# 轉換為 Tensor
X = adata.X.toarray() if hasattr(adata.X, "toarray") else adata.X
X_tensor = torch.FloatTensor(X)
# 2. 定義自編碼器模型
class Autoencoder(nn.Module):
def __init__(self, input_dim, latent_dim=2):
super(Autoencoder, self).__init__()
# 編碼器：降低維度
self.encoder = nn.Sequential(
nn.Linear(input_dim, 128),
nn.ReLU(),
nn.Linear(128, latent_dim),
)
# 解碼器：重建原始維度
self.decoder = nn.Sequential(
nn.Linear(latent_dim, 128),
nn.ReLU(),
nn.Linear(128, input_dim),
)
def forward(self, x):
latent = self.encoder(x)
reconstructed = self.decoder(latent)
return reconstructed, latent
# 3. 初始化模型、損失函數與優化器
input_dim = X.shape[1]
latent_dim = 1000 # 只有 6 個樣本，建議將隱藏層設小 (如 2)，方便可視化
model = Autoencoder(input_dim, latent_dim)
criterion = nn.MSELoss()
optimizer = optim.Adam(model.parameters(), lr=1e-3)
# 4. 訓練模型
epochs = 1000
model.train()
for epoch in range(epochs):
# 前向傳播
reconstructed, latent = model(X_tensor)
loss = criterion(reconstructed, X_tensor)
# 反向傳播與優化
optimizer.zero_grad()
loss.backward()
optimizer.step()
if (epoch + 1) % 20 == 0:
print(f"Epoch [{epoch+1}/{epochs}], Loss: {loss.item():.4f}")
# 5. 提取隱藏層特徵 (Latent Space)
model.eval()
with torch.no_grad():
_, latent_features = model(X_tensor)
# 將訓練好的特徵存回 AnnData 的 obsm 中，方便 Scanpy 繪圖
adata.obsm['X_ae'] = latent_features.numpy()
# 6. 可視化
sc.pl.scatter(adata, basis='ae', color='treatment', title="Autoencoder Latent Space")
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns
# 1. 提取編碼器第一層的權重 (Shape: 64 x 20899)
# model.encoder[0] 對應 nn.Linear(input_dim, 64)
weights = model.encoder[0].weight.detach().numpy()
# 2. 計算每個基因（column）的綜合影響力得分
# 我們取絕對值的總和，因為正向和負向權重都代表該基因具有資訊量
gene_importance = np.sum(np.abs(weights), axis=0)
# 3. 建立 DataFrame 並與基因名稱對應
importance_df = pd.DataFrame({
'gene_symbol': adata.var_names,
'importance_score': gene_importance
})
print("-^-" * 30)
# --- 1. 計算統計基準 ---
mean_score = gene_importance.mean()
std_score = gene_importance.std()
suggested_cutoff = mean_score + 2 * std_score
print(f"\n[統計資訊]")
print(f"平均得分 (Mean): {mean_score:.4f}")
print(f"標準差 (Std): {std_score:.4f}")
print(f"建議門檻值 (Mean + 2SD): {suggested_cutoff:.4f}")
print("-" * 30)
# --- 2. 暫停並等待輸入 ---
user_input = input(f"請輸入您想要的 Cutoff 門檻值 (直接按 Enter 則使用建議值 {suggested_cutoff:.4f}): ")
# 判斷輸入邏輯
try:
if user_input.strip() == "":
cutoff = suggested_cutoff
else:
cutoff = float(user_input)
except ValueError:
print("輸入格式錯誤，將自動使用建議門檻值。")
cutoff = suggested_cutoff
print(f"\n確認使用 Cutoff = {cutoff:.4f} 進行篩選...")
plt.figure(figsize=(8, 4))
sns.histplot(importance_df['importance_score'], bins=100, kde=True)
plt.axvline(cutoff, color='red', linestyle='--', label='Mean + 2SD')
plt.title("Distribution of Gene Importance Scores")
plt.legend()
plt.show()
# --- 3. 執行篩選 ---
# 標記是否顯著
importance_df['is_significant'] = importance_df['importance_score'] > cutoff
# 挑選出符合條件的基因並排序
top_genes_df = importance_df[importance_df['is_significant']].sort_values(
by='importance_score', ascending=False
)
# --- 4. 顯示結果 ---
num_sig = len(top_genes_df)
print(f"篩選完成！共找到 {num_sig} 個顯著基因。")
print('============================'*2)
print("--- Autoencoder 識別出的最具影響力基因 --- 基因數目 ＝", len(top_genes_df))
print(top_genes_df)
import matplotlib.pyplot as plt
# 5. 可視化：將這些重要基因的得分畫出來
import seaborn as sns
plt.figure(figsize=(10, 6))
sns.barplot(data=top_genes_df, x='importance_score', y='gene_symbol', palette='viridis')
plt.title(f"Influential Genes identified by AE, n = {num_sig}")
plt.xlabel("Cumulative Absolute Weights")
plt.show()
results_dict = {}
for i in range(1, 3): # 跑兩次
print(f"\n>>> 開始第 {i} 輪訓練...")
# 這裡放入你原本的代碼：定義模型、訓練、計算 importance...
# (為了簡潔，假設你已將 AE 訓練邏輯包裝成一個函數或直接放在這裡)
# 訓練與篩選邏輯結束後：
current_genes = set(importance_df[importance_df['importance_score'] > cutoff]['gene_symbol'])
results_dict[f'run_{i}'] = current_genes
print(f"第 {i} 輪完成，篩選出 {len(current_genes)} 個基因")
# 最後提取出來
genes_run1 = results_dict['run_1']
genes_run2 = results_dict['run_2']
# 假設第一次篩選出的基因清單為 genes_run1，第二次為 genes_run2
# 這裡根據你提供的 top_genes_df 邏輯提取基因名稱
# 1. 計算交集 (Intersection)
stable_genes = genes_run1.intersection(genes_run2)
# 2. 計算聯集 (Union)
all_detected_genes = genes_run1.union(genes_run2)
# 3. 計算重合率 (Jaccard Index)
jaccard_index = len(stable_genes) / len(all_detected_genes)
print("\n" + "="*40)
print("基因穩定性分析結果")
print("="*40)
print(f"第一次運行基因數: {len(genes_run1)}")
print(f"第二次運行基因數: {len(genes_run2)}")
print(f"真正穩定的基因數 (交集): {len(stable_genes)}")
print(f"重合率 (Jaccard Index): {jaccard_index:.2%}")
# 4. 顯示前 10 個最穩定的基因
print("\n前 10 個穩定基因範例:")
print(list(stable_genes)[:10])
# 5. 存檔：將穩定基因導出
pd.DataFrame({'stable_genes': list(stable_genes)}).to_csv("stable_genes_intersection.csv", index=False)
print("\n已將穩定基因清單存至: stable_genes_intersection.csv")
import pandas as pd
import gseapy as gp
import networkx as nx
import matplotlib.pyplot as plt
from pathlib import Path
# --- 1. 路徑設定 ---
desktop = Path.home() / "Desktop"
stable_csv = Path("stable_genes_intersection.csv") # AE 產出的檔案
excel_source = desktop / "Aja 21649_LeidCom_clustered.xlsx"
output_dir = desktop / "Autoncoder_Enrichment_Results"
output_dir.mkdir(exist_ok=True)
# --- 2. 數據讀取與映射 (Mapping) ---
print("🔍 正在進行基因映射...")
try:
# 讀取 AE 篩選出的穩定基因 (假設欄位名稱為 stable_genes)
ae_genes = pd.read_csv(stable_csv)['stable_genes'].tolist()
# 讀取原始含有註解的 Excel
df_ref = pd.read_excel(excel_source)
# 進行映射：在 Excel 中找尋 Transcriptid 存在於 ae_genes 清單中的行
# 注意：請確保兩者的名稱格式一致
mapped_df = df_ref[df_ref['Transcriptid'].isin(ae_genes)].copy()
# 提取 Gene_symbol 用於富集分析 (去除空值與重複)
gene_list = mapped_df['Gene_symbol'].replace('#N/A', pd.NA).dropna().unique().tolist()
print(f"✅ 映射完成！AE 篩選 {len(ae_genes)} 個 ID，成功對應到 {len(gene_list)} 個有效基因。")
except Exception as e:
print(f"❌ 讀取或映射失敗: {e}")
exit()
# --- 3. 基因富集分析 (Enrichr) ---
print("🚀 執行 KEGG 富集分析...")
try:
enr = gp.enrichr(gene_list=gene_list,
gene_sets=['KEGG_2021_Human'], # 或選取其他適合的物種資料庫
organism='Human',
cutoff=0.5)
res = enr.results
# 儲存完整富集結果
res.to_excel(output_dir / "AE_Stable_Genes_Enrichment.xlsx", index=False)
# 篩選出顯著的通路 (Adj P < 0.05) 作為繪圖對象
top_paths = res[res['Adjusted P-value'] < 0.1].head(10) # 取前 10 名
if top_paths.empty:
print("⚠️ 無顯著富集通路，調整門檻取前 5 名展示...")
top_paths = res.nsmallest(10, 'Adjusted P-value')
except Exception as e:
print(f"❌ 富集分析失敗: {e}")
exit()
# --- 4. 繪製 Cnetplot (Gene-Concept Network) ---
print("🎨 正在生成 Cnetplot...")
try:
G = nx.Graph()
pathway_label_map = {}
for _, row in top_paths.iterrows():
pathway = row['Term']
p_val = row['Adjusted P-value']
pathway_label_map[pathway] = f"{pathway}\n(p={p_val:.2e})"
# 獲取該通路中的基因並與我們的清單取交集
path_genes = str(row['Genes']).split(';')
valid_genes = [g for g in path_genes if g in gene_list]
# 加入通路節點
G.add_node(pathway, node_type='pathway', size=len(valid_genes))
# 加入基因節點與邊 (限制顯示數量避免畫面太亂)
for gene in valid_genes[:15]:
G.add_node(gene, node_type='gene')
G.add_edge(pathway, gene)
# 繪圖佈局
plt.figure(figsize=(15, 12))
pos = nx.spring_layout(G, k=0.5, iterations=50, seed=42)
# 分別繪製通路與基因
path_nodes = [n for n, a in G.nodes(data=True) if a.get('node_type') == 'pathway']
gene_nodes = [n for n, a in G.nodes(data=True) if a.get('node_type') == 'gene']
# 繪製邊
nx.draw_networkx_edges(G, pos, edge_color='#dfe6e9', alpha=0.5)
# 繪製 Pathway 節點 (紅圈)
sizes = [G.nodes[n]['size'] * 500 for n in path_nodes]
nx.draw_networkx_nodes(G, pos, nodelist=path_nodes, node_color='#ff7675',
node_size=sizes, edgecolors='black', label='Pathways')
# 繪製 Gene 節點 (藍圈)
nx.draw_networkx_nodes(G, pos, nodelist=gene_nodes, node_color='#74b9ff',
node_size=150, alpha=0.8, label='Genes')
# 加入標籤 (Pathway)
path_labels = {n: pathway_label_map[n] for n in path_nodes}
nx.draw_networkx_labels(G, pos, labels=path_labels, font_size=10, font_weight='bold')
# 加入標籤 (Gene)
gene_labels = {n: n for n in gene_nodes}
nx.draw_networkx_labels(G, pos, labels=gene_labels, font_size=8)
plt.title("Cnetplot: AE Stable Genes & Top KEGG Pathways", fontsize=15)
plt.axis('off')
# 存檔
plt.savefig(output_dir / "AE_Stable_Cnetplot.png", dpi=300, bbox_inches='tight')
plt.show()
print(f"✅ 分析完成！結果已存至: {output_dir}")
except Exception as e:
print(f"❌ 繪圖失敗: {e}")
import pandas as pd
import gseapy as gp
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
import traceback
import sys # 修正 exit() 問題
# --- 1. 路徑與參數設定 ---
desktop = Path.home() / "Desktop"
stable_csv = Path("stable_genes_intersection.csv")
excel_source = desktop / "Aja 21649_LeidCom_clustered.xlsx"
output_dir = desktop / "AE_Stability_Network_Analysis"
output_dir.mkdir(exist_ok=True)
# 你的目標數值欄位
value_col = 'avg_diff'
# --- 2. 數據映射與預處理 ---
print("🔍 正在進行基因數據映射...")
try:
ae_ids = pd.read_csv(stable_csv)['stable_genes'].tolist()
df_ref = pd.read_excel(excel_source)
# --- 自動欄位檢查機制 ---
if value_col not in df_ref.columns:
print(f"⚠️ 找不到欄位 '{value_col}'。")
print(f"現有欄位清單: {df_ref.columns.tolist()}")
# 嘗試找尋備用欄位
alternatives = ['log2FoldChange', 'logFC', 'p_val_adj', 'cluster']
found_alt = next((alt for alt in alternatives if alt in df_ref.columns), None)
if found_alt:
value_col = found_alt
print(f"💡 自動切換至備用欄位: '{value_col}'")
else:
print("❗ 找不到任何數值欄位，將建立預設數值 (1.0) 以維持繪圖。")
df_ref['default_value'] = 1.0
value_col = 'default_value'
mapped_df = df_ref[df_ref['Transcriptid'].isin(ae_ids)].copy()
if mapped_df.empty:
print("❌ 映射後沒有任何基因！請檢查 Transcriptid 是否匹配。")
sys.exit()
# 建立對照字典
gene_value_map = mapped_df.groupby('Gene_symbol')[value_col].mean().to_dict()
gene_list = list(gene_value_map.keys())
print(f"✅ 映射完成：{len(gene_list)} 個獨特基因符號。")
except Exception as e:
print(f"❌ 映射出錯: {e}")
traceback.print_exc()
sys.exit() # 修正後的安全退出
# --- 3. 基因富集分析與長條圖繪製 ---
print("🚀 執行 KEGG 富集分析...")
try:
enr = gp.enrichr(gene_list=gene_list,
gene_sets=['KEGG_2021_Human'],
organism='Human',
cutoff=0.5)
res = enr.results
if res.empty or 'Adjusted P-value' not in res.columns:
print("⚠️ 富集分析沒有結果，可能基因清單太少或不匹配。")
else:
res['-log10_Adj_P'] = -np.log10(res['Adjusted P-value'].replace(0, 1e-10))
# 繪製 Pathway Bar Plot
# 修正：確保 n 是整數
n_top = 20
top_10 = res.nsmallest(n_top, 'Adjusted P-value').iloc[::-1]
plt.figure(figsize=(10, 8)) # 稍微增加高度以容納 20 條
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(top_10)))
bars = plt.barh(top_10['Term'], top_10['-log10_Adj_P'], color=colors)
# 加入顯著性參考線 (P < 0.05 -> -log10(P) > 1.3)
# plt.axvline(x=-np.log10(0.1), color='red', linestyle='--', alpha=0.5, label='P=0.1')
plt.xlabel('-log10(Adjusted P-value)', fontsize=12)
plt.ylabel('Pathway Term', fontsize=12)
plt.title('Top Enriched KEGG Pathways (AE Stable Genes)', fontsize=14)
plt.legend()
plt.tight_layout()
plt.savefig(output_dir / "Pathway_Enrichment_Barplot.png", dpi=300)
plt.show()
# --- 4. 繪製純基因網路圖 (Network Graph) ---
print("🎨 生成穩定基因網路圖...")
G = nx.Graph()
top_paths = res.nsmallest(5, 'Adjusted P-value')
for _, row in top_paths.iterrows():
path_genes = [g for g in str(row['Genes']).split(';') if g in gene_list]
for i in range(len(path_genes)):
for j in range(i + 1, len(path_genes)):
G.add_edge(path_genes[i], path_genes[j])
if len(G.nodes) == 0:
print("⚠️ 基因間無共同通路聯繫，跳過網路圖製作。")
else:
plt.figure(figsize=(12, 12))
pos = nx.kamada_kawai_layout(G)
node_colors = [gene_value_map.get(n, 0) for n in G.nodes()]
nx.draw_networkx_edges(G, pos, alpha=0.2, edge_color='gray')
nodes = nx.draw_networkx_nodes(G, pos, node_size=500, node_color=node_colors,
cmap=plt.cm.RdBu_r, edgecolors='black', linewidths=0.5)
nx.draw_networkx_labels(G, pos, font_size=8, font_weight='bold')
plt.colorbar(nodes, label=f'Value ({value_col})', shrink=0.7)
plt.title("Gene Interaction Network")
plt.axis('off')
plt.savefig(output_dir / "Gene_Interaction_Network.png", dpi=300)
plt.show()
print(f"✨ 分析完成！結果已儲存至: {output_dir}")
except Exception as e:
traceback.print_exc()
The late arrival of domestic cats in China via the Silk Road after 3,500 years of human-leopard cat commensalism" http://spkl.io/6187Ad0Ct
https://youtu.be/XvDfmJ4KNDM?si=09JlMhK9_tdBEy5-
全球瀕危物種保育會議今天在烏茲別克撒馬罕（Samarkand）登場，預計將就加強鯊魚保育、有限度開放犀牛角販售，以及限制鰻魚貿易等議題進行表決。
法新社報導，180多個國家的代表與保育專家現正齊聚一堂，討論數十項野生動植物貿易提案。
這場會議召集了「瀕臨絕種野生動植物國際貿易公約」（CITES）的締約國。已實施50年的CITES旨在保護全球最瀕危的動植物，同時規範逾3萬6000個物種的貿易行為。
其中最具爭議的提案之一，是加強對亞洲部分地區常見食材—鰻魚的保護。
CITES根據動植物面臨的威脅程度，分別將其列入不同的附錄進行貿易管制。被列為附錄一的物種，被認為是全球最瀕危的種類，幾乎全面禁止貿易；附錄二的物種則須經過額外審核與證明才得以交易。
歐洲鰻（Anguilla anguilla）被認為是極危物種，並於2009年被納入附錄二，但數個國家希望將其他鰻魚物種列入管制。
他們主張，各種鰻魚在幼體階段事實上難以分辨，而貿易一般集中於這一階段，導致極危的歐洲鰻常被冒充成其他受威脅程度較低的鰻魚出售。
全球最大的鰻魚消費國日本強烈反對此提案，並積極展開遊說行動。
其他備受爭議的提案，還包括允許政府販售庫存象牙與犀牛角的請求。支持者主張，各國政府應可透過出售這些庫存來資助保育工作，但許多動保團體警告，重啟販售可能會助長非法貿易。
國際動物福利基金會（International Fund for AnimalWelfare）政策資深總監柯里斯（Matt Collis）表示：「一旦開放合法貿易，非法象牙和犀牛角就會有機可乘，流入已造成嚴重問題的亞洲市場。」
會議亦將討論對7種鯊魚加強保護，包括將極危的遠洋白鰭鯊（oceanic whitetip shark）列入附錄一，實施全面貿易禁令。
專家認為，這項措施有望獲得所需的2/3多數支持，但其他修正案可能面臨更多審查，包括對因魚肝油遭捕殺的鯨鯊（whale shark）及大西洋刺鯊（gulpershark）設下貿易限制。
歐盟推動鰻魚國際貿易限制，日本要中、韓聯手反對
歐盟近日確定了將鰻魚納入（簡稱《華盛頓公約》）的提案，若該提案在11月底召開的締約國會議上獲得通過，將對東亞地區的鰻苗及加工鰻魚產品的自由貿易造成重大影響。
為阻止提案通過，日本計劃與中國和韓國合作，共同遊說其他國家，強調。
歐盟提案涵蓋全球19種鰻魚 歐盟執行機構——歐盟委員會——計劃在締約國會議上提議，將包括日本鰻在內的18種鰻魚列入《華盛頓公約》的「附錄2」。加上已列入公約的歐洲鰻，這意味著全球所有19種鰻魚將全部受到貿易限制。據悉，歐盟成員國部長理事會近期已批准該提案，標誌著歐盟的談判立場正式確定，並完成了向締約國提交提案的準備工作。
提案通過需三分之二成員國支持 接下來，，並在11月24日至12月5日於烏茲別克舉行的締約國會議上進行討論。該提案需獲得約180個成員國中三分之二以上的支持方能通過。
日本擔憂鰻魚供應鏈受阻 面對歐盟堅定的立場，日本政府將與中國和韓國——這兩個主要的鰻魚捕撈國——合作，向與會國提出「鰻魚資源充足，無滅絕風險」的主張，以爭取支持。
若提案通過，將對日本的鰻魚流通造成嚴重阻礙，因為出口國需簽發許可證。據統計，2024年日本國內鰻魚供應量為6.3萬噸，其中70%以活鰻和加工鰻魚形式從中國進口。剩餘30%在日本國內養殖，但其中一半的鰻苗仍需從中國進口，。
市場擔憂價格上漲與銷售下滑
東京豐洲市場的鰻魚批發價格在2024年維持在每公斤約4500日元，與去年持平。一位批發商表示，若限制措施實施，，可能導致價格上漲。此外，九州地區的養殖業者也擔心，隨著「」的討論升溫，零售端的銷售量可能會下降。
