📚 课程目录第4章批量样本并行生成框架

批量样本并行生成框架：基于 Snakemake / Nextflow 的流程设计

好，咱们进入正题。

前几章聊了 GVCF 的原理，也讲了单样本怎么跑。但说实话，真正干活的时候，你面对的不是一个样本，而是几百上千个。一个个手动跑？那估计得跑到明年。

这一章，我就带你看看怎么用工作流引擎把这套流程自动化、并行化。我个人比较常用的是 Snakemake 和 Nextflow，各有各的脾气。但核心思想就一个：把样本拆开，同时跑，最后汇总。

为什么要用工作流框架？

说白了，就是管理依赖和并行。你想想看：

依赖管理：每个样本必须先完成比对、排序、去重，然后才能跑 HaplotypeCaller。框架会自动判断哪个先跑、哪个后跑。
并行调度：几百个样本，互不依赖，框架会根据你的集群资源（比如 SLURM、SGE、本地 CPU）同时扔出去跑。
断点续跑：跑一半集群挂了？没事，框架会记录哪些任务完成了，重启后只跑剩下的。嗯，这点我深有体会……

核心思想：把每个样本当作一个独立任务，框架自动调度。你只需要定义输入、输出和规则。

Snakemake 实现方案

Snakemake 是基于 Python 的，规则用类似 Makefile 的语法写。我个人习惯用它来做中小规模的流程，灵活，上手快。

一个典型的 GVCF 并行流程，核心就是 rule 的定义。比如下面这个 HaplotypeCaller 规则：

# Snakefile — GVCF 并行流程示例
configfile: "config.yaml"

rule all:
    input:
        expand("results/gvcf/{sample}.g.vcf.gz", 
               sample=SAMPLES)

rule haplotype_caller:
    input:
        bam = "results/bam/{sample}.sorted.dedup.bam",
        bai = "results/bam/{sample}.sorted.dedup.bam.bai",
        ref = config["reference"]
    output:
        gvcf = "results/gvcf/{sample}.g.vcf.gz"
    params:
        interval = config.get("interval", None)
    threads: 4
    run:
        cmd = f"gatk HaplotypeCaller \
                -R {input.ref} \
                -I {input.bam} \
                -O {output.gvcf} \
                -ERC GVCF \
                --native-pair-hmm-threads {threads}"
        if params.interval:
            cmd += f" -L {params.interval}"
        shell(cmd)

然后你有一个 config.yaml 文件，里面定义样本列表：

# config.yaml
samples:
  - sample_001
  - sample_002
  - sample_003
  # ... 几百个样本
reference: /path/to/reference/Homo_sapiens_assembly38.fasta
interval: /path/to/intervals/regions.interval_list

运行起来也很简单：

snakemake -s Snakefile --cores 40 --use-conda

如果你是用集群，比如 SLURM，就加个 --cluster 参数：

snakemake -s Snakefile --cluster "sbatch -c {threads} --mem=8G -t 2:00:00" \
          --jobs 20 --latency-wait 60

Nextflow 实现方案

Nextflow 是用 Groovy 写的，但不用特意去学 Groovy 语法，基本流程用起来很简单。它最大的优势是容器化和云原生支持——Docker/Singularity 一挂，环境问题少很多。

一个 Nextflow 流程大概长这样：

// main.nf — GVCF 并行流程
#!/usr/bin/env nextflow
nextflow.enable.dsl=2

params.samples = "samples.csv"
params.reference = "/path/to/reference.fasta"

process HAPLOTYPE_CALLER {
    tag "${sample_id}"
    cpus 4
    memory '8 GB'

    input:
    tuple val(sample_id), path(bam), path(bai)

    output:
    path("${sample_id}.g.vcf.gz"), emit: gvcf

    script:
    """
    gatk HaplotypeCaller \
        -R $params.reference \
        -I $bam \
        -O ${sample_id}.g.vcf.gz \
        -ERC GVCF \
        --native-pair-hmm-threads $task.cpus
    """
}

workflow {
    samples_ch = Channel.fromPath(params.samples)
        .splitCsv(header: true)
        .map { row -> [row.sample_id, file(row.bam), file(row.bai)] }

    HAPLOTYPE_CALLER(samples_ch)
    
    HAPLOTYPE_CALLER.out.gvcf.collect().set { all_gvcf }

    // 接下来可以接 CombineGVCFs 或 GenomicsDBImport
}

对应的 samples.csv 长这样：

sample_id,bam,bai
sample_001,results/bam/sample_001.sorted.dedup.bam,results/bam/sample_001.sorted.dedup.bam.bai
sample_002,results/bam/sample_002.sorted.dedup.bam,results/bam/sample_002.sorted.dedup.bam.bai

运行命令：

nextflow run main.nf -resume -with-docker

我的一点建议：如果你的团队已经有 Docker/Singularity 环境，Nextflow 的 -with-docker 选项实在太香了。环境不一致踩过的坑，我之前在 Snakemake 里遇到过好多次。

Snakemake vs Nextflow：我该怎么选？

这两者都在用，但适用场景有点不一样。我列个表格，你看着选：

维度	Snakemake	Nextflow
语言风格	Python，规则像 Makefile	Groovy，流程像 DSL
上手难度	低，Python 用户友好	中等，需要适应 Groovy 语法
容器支持	支持，但需要 conda 或 singularity	原生支持，docker/singularity 一键集成
断点续跑	优秀，`--rerun-incomplete`	优秀，`-resume`
集群调度	通过 `--cluster` 参数	通过 `executor` 配置
云原生	有限	强，AWS Batch / Google Cloud 直接支持
个人偏好	小团队、本地集群	需要容器化、云环境、大型项目

一句话总结：你团队 Python 多、环境简单，选 Snakemake；你追求标准化、容器化、云上部署，选 Nextflow。

实战中的坑与避坑指南

坑一：文件路径不一致。我曾经在 Snakemake 里因为 input 和 output 的路径写错一个斜杠，导致整个流程跑完才发现所有 gVCF 都是空的……建议用 Path 对象或者绝对路径。
坑二：资源分配不合理。HaplotypeCaller 默认会用多线程，但你给每个任务 8 个核，集群却只有 64 核，一次跑 8 个任务，反而比跑 16 个任务（每任务 4 核）慢。嗯，这里要默默计算一下资源。
坑三：中间文件清理。几百个样本的中间 bam 文件很占空间。Snakemake 的 --delete-all-output 或者 Nextflow 的 clean 命令，跑完后记得清理。

重要提醒：不要直接在生产环境跑几百个样本的 --cores 1 串行模式。那样会跑吐血的。并行粒度一定要和集群资源匹配，我建议每个任务 4~8 核，同时跑 CPU 总数 / 每任务核数 个任务。

一个完整的并行流程配置示例（Snakemake）

这里我给一个稍微完整点的例子，包含 config.yaml、Snakefile 和 cluster.json：

# config.yaml
samples: ["sample_001", "sample_002", "sample_003"]
reference: "/data/reference/Homo_sapiens_assembly38.fasta"
interval: "/data/interval/wgs.interval_list"
outdir: "/data/results/gvcf/"

# cluster.json
__default__:
    time: "2:00:00"
    mem: "8G"
    cpus: 4
    partition: "batch"

运行：

snakemake -s Snakefile \
          --cluster-config cluster.json \
          --cluster "sbatch --time {cluster.time} --mem {cluster.mem} --cpus-per-task {cluster.cpus} --partition {cluster.partition}" \
          --jobs 20 \
          --latency-wait 60
          --rerun-incomplete

这样，20 个样本同时跑，每个任务用 4 核 8G 内存，跑 2 小时。完美。

好，这一章就到这。你心里应该有个底了：用什么框架、怎么配置、怎么避免踩坑。