Agentic 编程课第 14 章：科研数据分析流水线

# 任务清单

## 模块 1：导入
- 1.1 load_counts() - 读表达矩阵
- 1.2 load_samples() - 读样本分组

## 模块 2：质控
- 2.1 filter_low_expression() - 过滤低表达
- 2.2 run_pca() - PCA 检测异常样本
- 2.3 qc_summary() - 生成 QC 报告

## 模块 3：分析
- 3.1 run_deseq2() - 跑 DESeq2
- 3.2 apply_thresholds() - 应用阈值筛差异基因
- 3.3 lfc_shrinkage() - log2FC 收缩

## 模块 4：报告
- 4.1 plot_volcano() - 火山图
- 4.2 top_genes_table() - Top 基因表
- 4.3 render_report() - 整合到 HTML

## MVP
1.1 + 1.2 + 2.1 + 3.1 + 3.2 + 4.1 + 4.3 就能出最小报告。

# 差异表达分析流水线规格

## R - 背景
我是生物信息学研究员。有一批 RNA-seq 数据，
需要做"对照组 vs 处理组"的差异表达分析（DEG）。

输入：counts.csv
- 行 = 基因（gene_id, gene_name）
- 列 = 样本（control_1/2/3, treatment_1/2/3）
- 值 = read counts

输出：完整的 HTML 报告

## F - 功能

### 步骤 1：数据导入
- 读 counts.csv
- 分离 metadata（基因名）和 counts 矩阵
- 读 samples.csv（样本分组信息）

### 步骤 2：质控（QC）
- 过滤低表达基因（所有样本总和 < 10 的剔除）
- 检测异常样本（用 PCA 看聚类）
- 输出 QC 报告

### 步骤 3：差异分析
- 用 pydeseq2（Python 实现的 DESeq2）
- 设计公式：~ condition
- 显著性阈值：padj < 0.05 AND |log2FC| > 1（我的判断）
- 做 LFC shrinkage（apeglm 方法）

### 步骤 4：报告
- 火山图（volcano plot）
- Top 20 上调 / 下调基因表
- 已知管家基因列表对照（housekeeping_genes.txt）
- 全部输出到 report.html

## C - 验收标准
- 用示例数据跑通全流程
- 火山图能看出明显的上下调分布
- Top 基因里有几个我认识的真实差异基因
- 全程 < 5 分钟跑完

> 任务：实现 src/analysis.py 的 run_deseq2() 函数

需求：
1. 输入：counts_df（基因 × 样本），samples_df（含 condition 列）
2. 用 pydeseq2 包做差异分析
3. 设计公式：~ condition
4. 返回 DataFrame：gene_id / gene_name / baseMean / log2FC / lfcSE / pvalue / padj

业务约束（我的领域知识）：
- control 组样本：control_1, control_2, control_3
- treatment 组样本：treatment_1, treatment_2, treatment_3
- log2FC 方向：treatment 相对 control（treatment 升高为正）

写 tests/test_analysis.py，至少包含：
- 跑通流程
- 检查返回 DataFrame 列名
- 已知的一个上调基因 log2FC > 0
- 已知的一个下调基因 log2FC < 0

先告诉我需要 pip install 什么，再实现。

# src/analysis.py
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats
import pandas as pd

def run_deseq2(counts_df, samples_df):
    # 业务规则：保证样本顺序对齐
    samples_df = samples_df.loc[counts_df.columns]

    inference = DefaultInference(n_cpus=8)
    dds = DeseqDataSet(
        counts=counts_df,
        metadata=samples_df,
        design="~ condition",
        refit_cooks=True,
        inference=inference,
    )
    dds.deseq2()

    stats = DeseqStats(dds, inference=inference)
    stats.summary()

    result = stats.results_df.copy()
    result = result.reset_index().rename(columns={
        'log2FoldChange': 'log2FC',
        'pvalue': 'pvalue',
        'padj': 'padj'
    })
    return result

> 跑 run_deseq2，告诉我：
- 总基因数
- padj < 0.05 的有几个
- log2FC > 1 的有几个
- log2FC < -1 的有几个
- Top 5 上调基因名
- Top 5 下调基因名

> 任务：实现 src/plotting.py 的 plot_volcano() 函数

需求：
- 输入：差异分析结果 DataFrame（gene_id / gene_name / log2FC / padj）
- 输出：matplotlib Figure 对象
- 要求：
  - X 轴：log2FC，范围 [-5, 5]
  - Y 轴：-log10(padj)，范围 [0, 50]
  - 颜色规则（业务判断）：
    - 红色：padj < 0.05 AND log2FC > 1（显著上调）
    - 蓝色：padj < 0.05 AND log2FC < -1（显著下调）
    - 灰色：其他
  - 在显著基因上标注 gene_name（避免重叠，用 adjustText 库）
  - 标题："Volcano Plot: Treatment vs Control"
  - 阈值线：log2FC = ±1（虚线）、padj = 0.05（虚线）

返回 Figure，不要保存（保存由调用方负责）。

# src/plotting.py
import matplotlib.pyplot as plt
import numpy as np
from adjustText import adjust_text

def plot_volcano(deg_df):
    df = deg_df.copy()
    df['neg_log_padj'] = -np.log10(df['padj'].clip(lower=1e-50))

    # 业务规则：颜色分类
    df['color'] = 'gray'
    df.loc[(df['padj'] < 0.05) & (df['log2FC'] > 1), 'color'] = 'red'
    df.loc[(df['padj'] < 0.05) & (df['log2FC'] < -1), 'color'] = 'blue'

    fig, ax = plt.subplots(figsize=(10, 8))
    ax.scatter(df['log2FC'], df['neg_log_padj'],
               c=df['color'], s=10, alpha=0.6)

    # 阈值线
    ax.axvline(x=1, color='gray', linestyle='--', alpha=0.5)
    ax.axvline(x=-1, color='gray', linestyle='--', alpha=0.5)
    ax.axhline(y=-np.log10(0.05), color='gray', linestyle='--', alpha=0.5)

    # 标注 Top 基因
    sig = df[(df['padj'] < 0.05) & (df['log2FC'].abs() > 1)]
    top = sig.reindex(sig['padj'].abs().sort_values().head(15).index)
    texts = []
    for _, row in top.iterrows():
        texts.append(ax.text(row['log2FC'], row['neg_log_padj'],
                             row['gene_name'], fontsize=8))
    adjust_text(texts, ax=ax)

    ax.set_xlabel('log2 Fold Change')
    ax.set_ylabel('-log10(padj)')
    ax.set_title('Volcano Plot: Treatment vs Control')
    ax.set_xlim(-5, 5)
    return fig

# main.py
from pathlib import Path
from src.loader import load_counts, load_samples
from src.qc import filter_low_expression, run_pca
from src.analysis import run_deseq2, apply_thresholds, lfc_shrinkage
from src.plotting import plot_volcano, plot_heatmap
from src.report import render_report

def main():
    print("加载数据...")
    counts, gene_meta = load_counts('data/counts.csv')
    samples = load_samples('data/samples.csv')

    print("质控...")
    counts_filtered = filter_low_expression(counts, min_total=10)
    pca_fig = run_pca(counts_filtered, samples)

    print("差异分析...")
    deg = run_deseq2(counts_filtered, samples)
    deg = lfc_shrinkage(deg)

    print("生成报告...")
    volcano_fig = plot_volcano(deg)
    heatmap_fig = plot_heatmap(deg, counts_filtered)

    sig_deg = apply_thresholds(deg, padj=0.05, log2fc=1)

    render_report(
        pca_fig=pca_fig,
        volcano_fig=volcano_fig,
        heatmap_fig=heatmap_fig,
        sig_deg=sig_deg,
        output_path='output/report.html'
    )

    print(f"完成，显著差异基因：{len(sig_deg)} 个")
    print("报告：output/report.html")

if __name__ == '__main__':
    main()

python main.py

加载数据...
质控...
差异分析...
生成报告...
完成，显著差异基因：1247 个
报告：output/report.html

> 看看当前 pip 环境的所有包版本，生成 requirements.txt，
  再生成一个 environment.yml（conda 用）。
  注释说明每个包的用途。

> 把 main.py 改造成 Snakemake 工作流，
  每个步骤是一个 rule，明确输入输出。
  生成 Snakefile。

pip install dvc
dvc init
dvc add data/counts.csv
git add counts.csv.dvc .gitignore
git commit -m "数据版本控制初始化"

dvc add data/counts.csv
git commit -am "更新到 v2 数据"

> counts_df 太大，内存不够。
  改用 scipy.sparse.csr_matrix 存储，
  或者按基因分块处理。

# 业务判断：低于 10 个 read 的根本不该分析
counts_filtered = counts[counts.sum(axis=1) >= 10]

# 业务规则：counts 列名必须和 samples 索引完全一致
assert set(counts.columns) == set(samples.index), "样本名不匹配"

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']  # Mac
# plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # Windows
plt.rcParams['axes.unicode_minus'] = False