62. AI自动化生信Pipeline¶

用户目标："帮我做RNA-seq差异表达分析"
    ↓
Agent大脑（LLM）思考：
  → 需要做哪些步骤？质控→比对→定量→差异分析
  → 数据是什么格式？看看输入文件…是FASTQ
  → 用什么工具？fastp质控、HISAT2比对、featureCounts定量、DESeq2差异
    ↓
Agent手脚（Tool Use）执行：
  → 调用fastp工具，设定参数
  → 检查输出质量，决定是否需要调参重跑
  → 继续下一步…

# 白话：给AI一个"工具箱"，AI根据需要自己挑工具用
tools = [
    {
        "name": "run_fastp",           # 工具名
        "description": "RNA-seq原始数据质控，去接头序列、低质量碱基", # AI看这个描述来决定什么时候用
        "parameters": {                # AI需要填的参数
            "input_r1": "正向reads文件路径",
            "input_r2": "反向reads文件路径",
            "qualified_quality": "碱基质量阈值，默认20",
            "length_required": "最短reads长度，默认50"
        }
    },
    {
        "name": "run_hisat2",
        "description": "将reads比对到参考基因组",
        "parameters": {
            "index": "HISAT2索引路径",
            "input_r1": "质控后的正向reads",
            "input_r2": "质控后的反向reads"
        }
    }
]

# ============================================================
# 用LangGraph构建RNA-seq自动分析Agent
# 环境：pip install langgraph langchain-openai
# ============================================================

import subprocess                    # 用于调用命令行工具
import os                            # 文件路径操作
from typing import TypedDict, Literal  # 类型注解
from langgraph.graph import StateGraph, END  # LangGraph核心：状态图
from langchain_openai import ChatOpenAI      # OpenAI接口（也可用本地模型）

# ---------- 第1步：定义Agent的"记忆"（State） ----------
class PipelineState(TypedDict):
    """Agent在整个分析过程中需要记住的信息"""
    task: str                    # 用户的分析任务描述
    input_files: list[str]       # 输入的FASTQ文件列表
    reference_genome: str        # 参考基因组路径
    current_step: str            # 当前执行到哪一步
    qc_passed: bool              # 质控是否通过
    aligned_bam: str             # 比对后的BAM文件路径
    counts_matrix: str           # 定量后的表达矩阵路径
    de_results: str              # 差异分析结果路径
    error_log: list[str]         # 错误记录
    quality_metrics: dict        # 各步骤的质量指标

# ---------- 第2步：定义每个分析节点（函数） ----------
def quality_control(state: PipelineState) -> PipelineState:
    """节点1：用fastp做质控"""
    print("🔬 执行质控：fastp...")

    r1 = state["input_files"][0]  # 正向reads
    r2 = state["input_files"][1]  # 反向reads（如果是双端）

    # 构造fastp命令
    cmd = [
        "fastp",
        "-i", r1,                      # 输入正向
        "-I", r2,                      # 输入反向
        "-o", "clean_R1.fastq.gz",     # 输出正向
        "-O", "clean_R2.fastq.gz",     # 输出反向
        "--qualified_quality_phred", "20",  # Q20质量过滤
        "--length_required", "50",     # 最短50bp
        "--json", "fastp_report.json", # JSON报告
        "--thread", "8"                # 8线程加速
    ]

    # 执行命令
    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode == 0:
        # 解析fastp报告判断质量
        import json
        with open("fastp_report.json") as f:
            report = json.load(f)

        # 提取关键质量指标
        q30_rate = report["summary"]["after_filtering"]["q30_rate"]

        state["qc_passed"] = q30_rate > 0.85  # Q30>85%算通过
        state["quality_metrics"]["qc_q30"] = q30_rate
        state["current_step"] = "qc_done"
    else:
        state["error_log"].append(f"fastp失败: {result.stderr}")
        state["qc_passed"] = False

    return state

def alignment(state: PipelineState) -> PipelineState:
    """节点2：用HISAT2做比对"""
    print("🧬 执行比对：HISAT2...")

    cmd = [
        "hisat2",
        "-x", state["reference_genome"],   # 参考基因组索引
        "-1", "clean_R1.fastq.gz",         # 质控后正向reads
        "-2", "clean_R2.fastq.gz",         # 质控后反向reads
        "--dta",                           # 下游用StringTie时加这个
        "-p", "8",                         # 8线程
        "-S", "aligned.sam"                # 输出SAM文件
    ]

    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode == 0:
        # SAM转BAM并排序
        subprocess.run(["samtools", "sort", "-@", "8", 
                       "-o", "aligned_sorted.bam", "aligned.sam"])
        subprocess.run(["samtools", "index", "aligned_sorted.bam"])

        state["aligned_bam"] = "aligned_sorted.bam"
        state["current_step"] = "alignment_done"

        # 提取比对率
        for line in result.stderr.split("\n"):
            if "overall alignment rate" in line:
                rate = float(line.split("%")[0].strip())
                state["quality_metrics"]["alignment_rate"] = rate
    else:
        state["error_log"].append(f"HISAT2失败: {result.stderr}")

    return state

def quantification(state: PipelineState) -> PipelineState:
    """节点3：用featureCounts做定量"""
    print("📊 执行定量：featureCounts...")

    cmd = [
        "featureCounts",
        "-a", "annotation.gtf",        # 基因注释文件
        "-o", "counts.txt",            # 输出计数矩阵
        "-T", "8",                     # 8线程
        "-p",                          # 双端测序
        "--countReadPairs",            # 按read pair计数
        state["aligned_bam"]           # 输入BAM文件
    ]

    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode == 0:
        state["counts_matrix"] = "counts.txt"
        state["current_step"] = "quantification_done"
    else:
        state["error_log"].append(f"featureCounts失败: {result.stderr}")

    return state

def differential_expression(state: PipelineState) -> PipelineState:
    """节点4：用DESeq2做差异表达分析"""
    print("📈 执行差异分析：DESeq2...")

    # 调用R脚本执行DESeq2
    r_script = """
    library(DESeq2)
    counts <- read.table("counts.txt", header=TRUE, row.names=1, skip=1)
    counts <- counts[, 6:ncol(counts)]  # 去掉前5列注释信息

    # 构建实验设计（示例：对照 vs 处理）
    condition <- factor(c(rep("control", 3), rep("treatment", 3)))
    colData <- data.frame(condition=condition)

    dds <- DESeqDataSetFromMatrix(countData=counts, colData=colData, design=~condition)
    dds <- DESeq(dds)
    res <- results(dds, alpha=0.05)

    write.csv(as.data.frame(res), file="DE_results.csv")
    """

    with open("run_deseq2.R", "w") as f:
        f.write(r_script)

    result = subprocess.run(["Rscript", "run_deseq2.R"], capture_output=True, text=True)

    if result.returncode == 0:
        state["de_results"] = "DE_results.csv"
        state["current_step"] = "de_done"
    else:
        state["error_log"].append(f"DESeq2失败: {result.stderr}")

    return state

# ---------- 第3步：定义路由逻辑（AI决策） ----------
def should_continue_after_qc(state: PipelineState) -> Literal["alignment", "retry_qc"]:
    """质控后AI决定：通过→继续比对，不通过→重试"""
    if state["qc_passed"]:
        return "alignment"
    else:
        return "retry_qc"  # 可以调参重跑

def check_alignment_quality(state: PipelineState) -> Literal["quantification", "retry_alignment"]:
    """检查比对率决定下一步"""
    rate = state["quality_metrics"].get("alignment_rate", 0)
    if rate > 70:  # 比对率>70%算合格
        return "quantification"
    else:
        return "retry_alignment"

# ---------- 第4步：构建LangGraph状态图 ----------
workflow = StateGraph(PipelineState)

# 添加节点
workflow.add_node("quality_control", quality_control)
workflow.add_node("alignment", alignment)
workflow.add_node("quantification", quantification)
workflow.add_node("differential_expression", differential_expression)

# 添加边（定义流转逻辑）
workflow.set_entry_point("quality_control")          # 入口：质控
workflow.add_conditional_edges(                       # 质控后条件判断
    "quality_control",
    should_continue_after_qc,
    {"alignment": "alignment", "retry_qc": "quality_control"}
)
workflow.add_conditional_edges(                       # 比对后条件判断
    "alignment",
    check_alignment_quality,
    {"quantification": "quantification", "retry_alignment": "alignment"}
)
workflow.add_edge("quantification", "differential_expression")  # 定量→差异分析
workflow.add_edge("differential_expression", END)               # 差异分析→结束

# 编译图
app = workflow.compile()

# ---------- 第5步：运行Agent ----------
if __name__ == "__main__":
    # 初始化状态
    initial_state = PipelineState(
        task="RNA-seq差异表达分析",
        input_files=["sample_R1.fastq.gz", "sample_R2.fastq.gz"],
        reference_genome="/ref/hisat2_index/hg38",
        current_step="start",
        qc_passed=False,
        aligned_bam="",
        counts_matrix="",
        de_results="",
        error_log=[],
        quality_metrics={}
    )

    # 执行整个pipeline
    final_state = app.invoke(initial_state)

    # 输出结果摘要
    print(f"\n{'='*50}")
    print(f"分析完成！")
    print(f"差异表达结果：{final_state['de_results']}")
    print(f"质量指标：{final_state['quality_metrics']}")
    if final_state['error_log']:
        print(f"警告日志：{final_state['error_log']}")

┌─────────────────────────────────────────────┐
│              LangGraph StateGraph             │
│                                              │
│  [quality_control] ──条件──→ [alignment]     │
│        ↑ 不通过则重试            │            │
│        └────────────────────────┘            │
│                                   ↓条件判断   │
│                            [quantification]  │
│                                   ↓          │
│                      [differential_expression]│
│                                   ↓          │
│                                 [END]        │
└─────────────────────────────────────────────┘

# 用Hermes模型（通过Ollama本地运行）做生信Agent
# 安装：ollama pull hermes3:8b

from langchain_ollama import ChatOllama  # 连接本地Ollama

# 初始化本地Hermes模型
llm = ChatOllama(
    model="hermes3:8b",        # 本地运行的Hermes 3模型
    temperature=0,             # 生信分析要确定性输出
)

# 定义生信工具集
bioinformatics_tools = [
    {"name": "blast_search", "desc": "序列相似性搜索"},
    {"name": "run_fastp", "desc": "测序数据质控"},
    {"name": "run_bwa", "desc": "短reads比对"},
    {"name": "variant_calling", "desc": "变异检测"},
    {"name": "gene_annotation", "desc": "基因功能注释"},
]

# Hermes会自动根据用户需求选择合适的工具

def validate_de_results(results_file: str) -> dict:
    """验证差异表达结果的合理性"""
    import pandas as pd

    df = pd.read_csv(results_file)

    checks = {
        "total_genes": len(df),                              # 总基因数
        "deg_count": len(df[df["padj"] < 0.05]),            # 显著差异基因数
        "deg_ratio": len(df[df["padj"] < 0.05]) / len(df), # DEG占比
        "max_log2fc": df["log2FoldChange"].abs().max(),     # 最大fold change
        "na_ratio": df["padj"].isna().sum() / len(df),      # NA值占比
    }

    # 合理性判断
    warnings = []
    if checks["deg_ratio"] > 0.3:
        warnings.append("⚠️ DEG占比超过30%，可能有批次效应")
    if checks["max_log2fc"] > 10:
        warnings.append("⚠️ 存在极端fold change，检查是否有异常样本")
    if checks["total_genes"] < 10000:
        warnings.append("⚠️ 检出基因数偏少，检查比对和定量步骤")

    checks["warnings"] = warnings
    return checks

# 安装LangGraph
pip install langgraph langchain-openai langchain-ollama

# 本地运行Hermes模型
ollama pull hermes3:8b
ollama run hermes3:8b

# 可视化LangGraph（生成流程图）
from IPython.display import Image
Image(app.get_graph().draw_mermaid_png())