401_时序多组学数据分析¶

import pandas as pd  # 数据处理
import numpy as np   # 数值计算
import matplotlib.pyplot as plt  # 绘图

# ========== 1. 模拟时序多组学数据 ==========
np.random.seed(42)
timepoints = [0, 6, 12, 24, 48]  # 时间点（小时）
n_genes = 100  # 基因数量
n_samples = 3  # 每个时间点的生物学重复数

# 生成转录组时序数据
rna_data = {}
for t in timepoints:
    # 模拟不同时间点的表达量（加入时间趋势+噪声）
    expr = np.random.poisson(lam=100 + t * 2, size=(n_genes, n_samples))
    rna_data[f"T{t}h"] = pd.DataFrame(
        expr,
        index=[f"Gene_{i}" for i in range(n_genes)],
        columns=[f"Sample_{j}" for j in range(n_samples)]
    )

# ========== 2. 合并各时间点数据 ==========
# 计算每个时间点的均值表达量
mean_expr = pd.DataFrame({
    tp: df.mean(axis=1)  # 对每个基因取3个重复的均值
    for tp, df in rna_data.items()
})
print("时序表达矩阵形状:", mean_expr.shape)  # (100基因, 5时间点)

# ========== 3. 识别时序变化基因 ==========
from scipy import stats

# 对每个基因做线性回归（时间 vs 表达量）
time_numeric = np.array(timepoints)  # 将时间转为数值
results = []

for gene in mean_expr.index:
    expr_values = mean_expr.loc[gene].values  # 该基因在各时间点的表达量
    slope, intercept, r_value, p_value, std_err = stats.linregress(
        time_numeric, expr_values
    )
    results.append({
        "gene": gene,
        "slope": slope,      # 斜率：正值=随时间升高，负值=降低
        "r_squared": r_value**2,  # 决定系数：反映线性拟合程度
        "p_value": p_value   # 统计显著性
    })

results_df = pd.DataFrame(results)
# 筛选显著时序变化基因（p<0.05 且 R²>0.8）
sig_genes = results_df[(results_df["p_value"] < 0.05) & 
                        (results_df["r_squared"] > 0.8)]
print(f"显著时序变化基因数: {len(sig_genes)}")

# ========== 4. 聚类时序模式 ==========
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans

# 对表达量做Z-score归一化（消除量级差异，关注变化趋势）
scaler = StandardScaler()
expr_scaled = scaler.fit_transform(mean_expr.T).T  # 按基因归一化

# K-means聚类：找3种时序模式（上升、下降、先升后降等）
kmeans = KMeans(n_clusters=3, random_state=42)
clusters = kmeans.fit_predict(expr_scaled)

mean_expr["cluster"] = clusters  # 将聚类结果加入数据框

# ========== 5. 可视化每个聚类的时序模式 ==========
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for cluster_id in range(3):
    cluster_genes = mean_expr[mean_expr["cluster"] == cluster_id]
    ax = axes[cluster_id]
    # 画每个基因的时序曲线（浅色）
    for _, row in cluster_genes.drop("cluster", axis=1).iterrows():
        ax.plot(timepoints, row.values, alpha=0.3, color="steelblue")
    # 画聚类均值曲线（深色）
    ax.plot(timepoints, 
            cluster_genes.drop("cluster", axis=1).mean().values,
            color="red", linewidth=2, label="Mean")
    ax.set_title(f"Cluster {cluster_id} (n={len(cluster_genes)})")
    ax.set_xlabel("Time (hours)")
    ax.set_ylabel("Normalized Expression")

plt.tight_layout()
plt.savefig("temporal_clusters.png", dpi=150)
print("图片已保存")