如何在Snakemake中实现配对肿瘤-正常样本的Mutect2变异检测？

核心思路是为每一对肿瘤-正常样本单独创建一个Mutect2任务，而非将所有样本塞进同一个任务。以下是具体实现步骤：

1. 加载配置并定义样本对列表 #

在Snakefile开头加载配置文件，从配置中提取每一组配对样本：

configfile: "config.yaml"

# 生成(肿瘤样本名, 正常样本名, 样本编号)的元组列表
sample_pairs = [
    (item["tumor"], item["control"], str(item["sample"]))
    for item in config["sample_list"]
]

2. 编写Mutect2规则 #

假设肿瘤/正常样本的BAM文件路径为bams/{sample}.bam，参考基因组为ref/hg38.fasta，规则写法如下：

rule mutect2:
    input:
        tumor = "bams/{tumor_sample}.bam",
        control = "bams/{control_sample}.bam",
        ref = "ref/hg38.fasta"
    output:
        vcf = "mutect2_results/sample_{sample_id}.vcf.gz",
        idx = "mutect2_results/sample_{sample_id}.vcf.gz.tbi"
    params:
        # 可按需添加Mutect2所需其他参数，如germline资源、PON等
        germline = "resources/af-only-gnomad.hg38.vcf.gz"
    log:
        "logs/mutect2_sample_{sample_id}.log"
    shell:
        """
        gatk Mutect2 \
            -R {input.ref} \
            -I {input.tumor} -tumor {wildcards.tumor_sample} \
            -I {input.control} -normal {wildcards.control_sample} \
            --germline-resource {params.germline} \
            -O {output.vcf} \
            2> {log}
        """

3. 指定目标文件 #

在Snakefile末尾添加要生成的所有Mutect2结果文件，让Snakemake自动调度任务：

rule all:
    input:
        expand("mutect2_results/sample_{sample_id}.vcf.gz", sample_id=[pair[2] for pair in sample_pairs])

关键说明 #

• 之前用zip方法出错的原因是将所有肿瘤/正常样本作为同一个任务的输入，导致Snakemake把所有肿瘤文件合并到一个-I参数后，所有正常文件合并到另一个-I参数后，不符合Mutect2的配对样本输入逻辑。
• 上述写法通过sample_pairs为每一对样本创建独立任务，每个任务仅处理一组肿瘤-正常样本，命令中的-I参数会正确对应单组配对。

内容的提问来源于stack exchange，提问作者Franz