如何在Snakemake规则中运行for循环避免MissingOutputException

问题背景 #

需要添加一个kraken规则，一次性处理所有样本（让Kraken2数据库驻留内存以提升效率），而非用wildcard逐个处理。但尝试用shell循环实现时，频繁触发MissingOutputException错误。

现有Snakemake核心脚本：

import os

SAMPLES = [i.replace('sample_', '') for i in os.listdir('final_raw')]

rule all:
    input:
        expand('analysis/results/unmapped/hybrid/read1_{sample}_hy.fastq', sample = SAMPLES)

# trimming、create_hybrid_genome等其余规则保持不变

尝试的错误循环代码：

for FILE in $(ls analysis/results/unmapped/hybrid/read1_*_hy.fastq | sed 's/read1_*_hy.fastq//'); do \
        kraken2 --db {input.database} \
        --memory-mapping \
        --threads 8 \
        --use-names \
        --report analysis/results/kraken2/hybrid/{{}}${{FILE}}.kraken \
        --paired read1_${{FILE}}_hy.fastq read2_${{FILE}}_hy.fastq \
        --unclassified-out /analysis/results/unclassified/hybrid/{{}}${{FILE}}_uc#.fastq; done

常规单样本处理规则（效率低，无法驻留数据库）：

rule kraken2:
    input:
        database = 'kraken2_db',
        read1 = 'analysis/results/unmapped/hybrid/read1_{sample}_hy.fastq',
        read2 = 'analysis/results/unmapped/hybrid/read2_{sample}_hy.fastq'
    output:
        'analysis/results/kraken2/hybrid/{sample}.kraken'
    shell:
        'kraken2 '
            '--db {input.database} '
            '--threads 8 '
            '--paired '
            '--use-names '
            '--unclassified-out {wildcards.sample}_uc#.fq '
            '--report {output} '
            '{input.read1} {input.read2}'

解决方案 #

核心问题：Snakemake需要明确知晓规则的所有输入输出文件，不能依赖shell动态文件查找（如ls），否则无法追踪输出是否生成，从而抛出MissingOutputException。

步骤1：更新rule all，明确所有目标输出 #

rule all:
    input:
        expand('analysis/results/unmapped/hybrid/read1_{sample}_hy.fastq', sample = SAMPLES),
        expand('analysis/results/kraken2/hybrid/{sample}.kraken', sample = SAMPLES),
        expand(['analysis/results/unclassified/hybrid/{sample}_uc_1.fastq', 
                'analysis/results/unclassified/hybrid/{sample}_uc_2.fastq'], 
               sample = SAMPLES)

步骤2：编写批量处理的kraken2_batch规则 #

rule kraken2_batch:
    input:
        database = 'kraken2_db',
        read1_list = expand('analysis/results/unmapped/hybrid/read1_{sample}_hy.fastq', sample = SAMPLES),
        read2_list = expand('analysis/results/unmapped/hybrid/read2_{sample}_hy.fastq', sample = SAMPLES)
    output:
        reports = expand('analysis/results/kraken2/hybrid/{sample}.kraken', sample = SAMPLES),
        uc_reads = expand(['analysis/results/unclassified/hybrid/{sample}_uc_1.fastq', 
                           'analysis/results/unclassified/hybrid/{sample}_uc_2.fastq'], 
                          sample = SAMPLES)
    threads: 8
    shell:
        """
        # 直接使用Snakemake定义的样本列表，避免动态查找
        for sample in {SAMPLES}; do
            kraken2 --db {input.database} \\
                --memory-mapping \\
                --threads {threads} \\
                --use-names \\
                --report analysis/results/kraken2/hybrid/${{sample}}.kraken \\
                --paired analysis/results/unmapped/hybrid/read1_${{sample}}_hy.fastq analysis/results/unmapped/hybrid/read2_${{sample}}_hy.fastq \\
                --unclassified-out analysis/results/unclassified/hybrid/${{sample}}_uc#.fastq
        done
        """

关键说明 #

1. 输入输出明确化：所有输入输出通过expand基于SAMPLES生成，Snakemake可提前知晓所有需生成的文件，避免追踪失效。
2. 循环变量一致性：直接用Snakemake定义的SAMPLES变量，确保循环样本列表与目标完全匹配。
3. 语法冲突规避：shell变量用${sample}而非{sample}，避免和Snakemake变量语法冲突。
4. 内存驻留生效：--memory-mapping参数让Kraken2加载数据库后保持内存映射，批量处理时无需重复加载，提升效率。

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