蛋白组与转录组联合分析¶

# === 转录组-蛋白组联合分析 ===
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr

# 1. 读取两组学数据
rna = pd.read_csv("rnaseq_log2fc.csv", index_col=0)  # RNA-seq差异结果
protein = pd.read_csv("proteomics_log2fc.csv", index_col=0)  # 蛋白组差异结果

# 2. 基因ID映射和取交集
# 假设两者都用基因名（Gene Symbol）作为index
common_genes = rna.index.intersection(protein.index)  # 取交集
print(f"RNA检测基因: {len(rna)}, 蛋白检测: {len(protein)}, 交集: {len(common_genes)}")

rna_common = rna.loc[common_genes, "log2FC"]
protein_common = protein.loc[common_genes, "log2FC"]

# 3. 计算整体相关性
r_pearson, p_pearson = pearsonr(rna_common, protein_common)
r_spearman, p_spearman = spearmanr(rna_common, protein_common)
print(f"Pearson r = {r_pearson:.3f} (p={p_pearson:.2e})")
print(f"Spearman rho = {r_spearman:.3f} (p={p_spearman:.2e})")

# 4. 四象限散点图
fig, ax = plt.subplots(figsize=(10, 10))

# 定义四象限颜色
colors = []
for rna_fc, prot_fc in zip(rna_common, protein_common):
    if rna_fc > 0 and prot_fc > 0:
        colors.append("red")       # RNA↑ 蛋白↑
    elif rna_fc > 0 and prot_fc < 0:
        colors.append("blue")      # RNA↑ 蛋白↓ (翻译抑制)
    elif rna_fc < 0 and prot_fc < 0:
        colors.append("green")     # RNA↓ 蛋白↓
    else:
        colors.append("orange")    # RNA↓ 蛋白↑ (蛋白稳定)

ax.scatter(rna_common, protein_common, c=colors, alpha=0.5, s=15)
ax.axhline(0, color="gray", linestyle="--", linewidth=0.5)
ax.axvline(0, color="gray", linestyle="--", linewidth=0.5)
ax.set_xlabel("RNA log2FC", fontsize=12)
ax.set_ylabel("Protein log2FC", fontsize=12)
ax.set_title(f"Transcriptome vs Proteome (r={r_pearson:.3f})", fontsize=14)

# 标注四象限基因数
q1 = sum(1 for c in colors if c == "red")
q2 = sum(1 for c in colors if c == "blue")
q3 = sum(1 for c in colors if c == "green")
q4 = sum(1 for c in colors if c == "orange")
ax.text(0.8, 0.95, f"Q1(↑↑): {q1}", transform=ax.transAxes, color="red")
ax.text(0.8, 0.90, f"Q2(↑↓): {q2}", transform=ax.transAxes, color="blue")
ax.text(0.05, 0.05, f"Q3(↓↓): {q3}", transform=ax.transAxes, color="green")
ax.text(0.05, 0.10, f"Q4(↓↑): {q4}", transform=ax.transAxes, color="orange")

plt.savefig("rna_protein_quadrant.png", dpi=150, bbox_inches="tight")

# 5. 翻译效率分析
# 翻译效率 = 蛋白丰度 / mRNA丰度（取log后做差）
translation_eff = protein_common - rna_common  # log2(蛋白/mRNA)
top_te = translation_eff.nlargest(20)   # 翻译效率最高的20个基因
bottom_te = translation_eff.nsmallest(20)  # 翻译效率最低的20个基因
print("翻译效率最高的基因:", top_te.index.tolist()[:5])