如何用Snakemake处理特定配对文件及Fastq样本合并问题

一、正确编写concat_fastq规则 #

你遇到的通配符推断错误，核心原因是输入的通配符无法从输出中匹配，或是expand的使用方式不对。以下是两种可靠的实现方式：

方式1：硬编码lane编号（适合固定lane数场景） #

假设原始Fastq文件存放在raw_fastq/目录下，规则可写为：

rule concat_fastq:
    output:
        "{sample_id}_{read}.fastq.gz"  # 输出格式：样本ID_R1/R2.fastq.gz
    input:
        # 动态生成对应样本+read的所有lane文件
        lambda wildcards: expand(
            "raw_fastq/{sample_id}_L{lane}_{read}_001.fastq.gz",
            lane=["001", "002"]  # 根据实际lane编号调整
        )
    shell:
        "cat {input} > {output}"

方式2：自动匹配所有lane（更灵活） #

如果lane数量不固定，用glob自动抓取所有匹配文件：

import glob

rule concat_fastq:
    output:
        "{sample_id}_{read}.fastq.gz"
    input:
        lambda wildcards: sorted(glob.glob(
            f"raw_fastq/{wildcards.sample_id}_*_{wildcards.read}_*.fastq.gz"
        ))
    shell:
        "cat {input} > {output}"

触发规则的关键：rule all #

要让Snakemake自动执行合并，需在rule all中定义所有需要的合并后文件，结合你的配置文件：

configfile: "config.yaml"

rule all:
    input:
        # 提取所有肿瘤和对照样本ID，生成合并后的文件路径
        expand(
            "{sample_id}_{read}.fastq.gz",
            sample_id=[s["tumor"] for s in config["sample_list"]] + [s["control"] for s in config["sample_list"]],
            read=["R1", "R2"]
        )

二、配置文件的合理性确认及后续分析建议 #

你的配置文件结构非常适合肿瘤-正常配对分析，每个条目明确了样本对的对应关系，后续可直接遍历config["sample_list"]处理配对样本。

后续变异检测示例规则（以GATK Mutect2为例） #

rule call_somatic_variants:
    output:
        "variants/{sample}_somatic.vcf.gz"
    input:
        tumor_bam="bam/{tumor}.bam",  # 假设已有生成BAM的比对规则
        normal_bam="bam/{control}.bam"
    params:
        ref="reference/hg38.fasta"  # 替换为你的参考基因组路径
    shell:
        """
        gatk Mutect2 \
            -R {params.ref} \
            -I {input.tumor_bam} \
            -I {input.normal_bam} \
            -tumor {wildcards.tumor} \
            -normal {wildcards.control} \
            -O {output}
        """

更新rule all加入变异检测输出：

rule all:
    input:
        # 合并后的Fastq（若后续直接用BAM可省略）
        expand(
            "{sample_id}_{read}.fastq.gz",
            sample_id=[s["tumor"] for s in config["sample_list"]] + [s["control"] for s in config["sample_list"]],
            read=["R1", "R2"]
        ),
        # 变异检测结果
        expand("variants/{sample}_somatic.vcf.gz", sample=[s["sample"] for s in config["sample_list"]])

常见问题说明 #

• 之前的expand错误：直接使用expand不结合lambda会一次性生成所有样本的文件组合，导致错误合并；用lambda结合wildcards可确保只匹配当前样本的文件。
• 通配符推断错误：输入中使用输出里没有的通配符（如{path}）时，Snakemake无法推断其值，要么在输出中包含该通配符，要么用lambda动态生成输入避免直接写通配符。

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