差异表达分析

import pandas as pd
import numpy as np
from pydeseq2 import DESeq2
from pydeseq2.preprocessing import deseq2_norm

# 加载计数矩阵
counts_df = pd.read_csv('count_matrix.csv', index_col=0)

# 加载样本信息
sample_info = pd.DataFrame({
    'condition': ['control']*5 + ['treated']*5  # 对照组(n=5) + 处理组(n=5)
}, index=counts_df.columns)

# 查看数据维度
print(f"基因数: {counts_df.shape[0]}")
print(f"样本数: {counts_df.shape[1]}")

# 预过滤：去除低表达基因
keep = (counts_df >= 10).sum(axis=1) >= 3
counts_filtered = counts_df[keep]

print(f"过滤后基因数: {counts_filtered.shape[0]}")

# 创建DESeq2对象
dds = DESeq2(
    counts=counts_filtered,
    design_factors=sample_info,
    design_formula="~ condition"
)

# 运行DESeq2流程
# 包括：size factor估计、离散度估计、Wald检验
dds.run_deseq()

# 获取结果
res = dds.results_df

# 筛选显著差异基因
sig_genes = res[
    (res['padj'] < 0.05) &      # BH校正p值 < 0.05
    (abs(res['log2FoldChange']) > 1)  # |log2FC| > 1
]

print(f"发现 {len(sig_genes)} 个显著差异基因")
print(sig_genes.head(10))

import matplotlib.pyplot as plt
import seaborn as sns

# 1. 火山图
fig, ax = plt.subplots(figsize=(10, 6))

# 标记显著基因
res['significant'] = (res['padj'] < 0.05) & (abs(res['log2FoldChange']) > 1)

# 绘制散点
sns.scatterplot(
    data=res,
    x='log2FoldChange',
    y='-np.log10(padj)',
    hue='significant',
    palette={False: 'gray', True: '#8B5CF6'},
    ax=ax
)

# 标注top 10基因
top_genes = sig_genes.nsmallest(10, 'padj')
for _, gene in top_genes.iterrows():
    ax.annotate(gene.name, (gene['log2FoldChange'], -np.log10(gene['padj'])))

ax.set_xlabel('Log2 Fold Change')
ax.set_ylabel('-Log10(Adjusted P-value)')
ax.set_title('Volcano Plot')
plt.savefig('volcano_plot.png', dpi=300, bbox_inches='tight')

# 2. MA图
fig, ax = plt.subplots(figsize=(10, 6))
sns.scatterplot(
    data=res,
    x='baseMean',
    y='log2FoldChange',
    hue='significant',
    palette={False: 'gray', True: '#8B5CF6'},
    ax=ax
)
ax.axhline(y=0, color='red', linestyle='--')
ax.set_xlabel('Mean of Normalized Counts')
ax.set_ylabel('Log2 Fold Change')
ax.set_title('MA Plot')
plt.savefig('ma_plot.png', dpi=300, bbox_inches='tight')

# 3. PCA图（样本聚类）
normalized_counts = dds.normalized_counts
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
pca_result = pca.fit_transform(normalized_counts.T)

fig, ax = plt.subplots(figsize=(8, 6))
sns.scatterplot(
    x=pca_result[:, 0],
    y=pca_result[:, 1],
    hue=sample_info['condition'],
    palette={'control': '#3B82F6', 'treated': '#8B5CF6'},
    ax=ax
)
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_title('PCA Plot')
plt.savefig('pca_plot.png', dpi=300, bbox_inches='tight')

print("可视化图表已保存！")

import seaborn as sns

# 选择top 50差异基因
top50_genes = sig_genes.nsmallest(50, 'padj').index

# 提取标准化表达值
heatmap_data = dds.normalized_counts.loc[top50_genes]

# z-score标准化
heatmap_data_z = (heatmap_data - heatmap_data.mean(axis=1).values.reshape(-1)) / heatmap_data.std(axis=1).values.reshape(-1)

# 绘制热图
g = sns.clustermap(
    heatmap_data_z,
    cmap='coolwarm',
    center=0,
    figsize=(12, 10),
    row_cluster=True,
    col_cluster=True,
    cbar_kws={'label': 'Z-score'}
)

# 添加样本分组注释
col_colors = sample_info['condition'].map({'control': '#3B82F6', 'treated': '#8B5CF6'})
g.ax_heatmap.set_xlabel('Samples')
g.ax_heatmap.set_ylabel('Genes')
g.ax_col_dendrogram.set_visible(False)

plt.savefig('heatmap_top50.png', dpi=300, bbox_inches='tight')
print("热图已保存！")

import gget

# 选择top 20差异基因进行注释
top20_genes = sig_genes.nsmallest(20, 'padj').index.tolist()

# 批量查询基因信息
gene_annotations = []
for gene in top20_genes:
    try:
        info = gget.info(gene, "ensembl")
        gene_annotations.append({
            'gene': gene,
            'full_name': info.get('name', 'N/A'),
            'description': info.get('description', 'N/A'),
            'log2FC': sig_genes.loc[gene, 'log2FoldChange'],
            'padj': sig_genes.loc[gene, 'padj']
        })
    except:
        gene_annotations.append({
            'gene': gene,
            'full_name': 'N/A',
            'description': 'N/A',
            'log2FC': sig_genes.loc[gene, 'log2FoldChange'],
            'padj': sig_genes.loc[gene, 'padj']
        })

# 创建注释数据框
annotation_df = pd.DataFrame(gene_annotations)

# 保存注释结果
annotation_df.to_csv('gene_annotations.csv', index=False)

print("Top 20差异基因注释：")
print(annotation_df.to_string())

# 查询GO注释（以某个基因为例）
top_gene = top20_genes[0]
try:
    go_result = gget.go(top_gene, "ensembl")
    print(f"\n{top_gene} GO注释:")
    print(go_result.head(10))
except:
    print("GO注释查询失败")

# 导出所有基因的差异分析结果
res.to_csv('deseq2_results_all.csv', index=True)

# 导出显著差异基因
sig_genes.to_csv('deseq2_results_significant.csv', index=True)

# 导出标准化表达矩阵
dds.normalized_counts.to_csv('normalized_counts.csv')

# 生成分析报告
report = f"""
差异表达分析报告
================

样本信息:
- 对照组: {sum(sample_info['condition']=='control')} 样本
- 处理组: {sum(sample_info['condition']=='treated')} 样本

分析结果:
- 总基因数: {len(res)}
- 显著差异基因数: {len(sig_genes)}
  * 上调基因: {sum(sig_genes['log2FoldChange'] > 0)}
  * 下调基因: {sum(sig_genes['log2FoldChange'] < 0)}

筛选标准:
- |log2FC| > 1
- padj < 0.05

生成文件:
- volcano_plot.png (火山图)
- ma_plot.png (MA图)
- pca_plot.png (PCA图)
- heatmap_top50.png (热图)
- deseq2_results_all.csv (所有基因)
- deseq2_results_significant.csv (显著基因)
- gene_annotations.csv (基因注释)
"""

with open('deseq2_report.txt', 'w') as f:
    f.write(report)

print("分析完成！报告已保存。")
print(report)