Snakemake动态文件列表生成及checkpoint报错问题求助

问题背景 #

首次开发Snakemake工作流，需求是从NCBI数据库下载名称未知的基因组文件，流程分为四步：

• 调用Entrez工具查询数据库，生成基因组文件链接列表
• 下载基因组压缩文件
• 解压文件
• 汇总所有已解压的基因组文件

前3步可正常运行，但汇总步骤存在问题：

• 初始通过提前读取temp目录生成输入列表的方法，仅在temp目录预先存在时有效
• 尝试用checkpoint解决动态文件名问题时，触发错误：

  InputFunctionException in line 73 (rule make_summary_table) of ~/snakemake_test/Snakefile: WildcardError: No values given for wildcard 'genome'. Wildcards: 
  

错误原因 #

原checkpoint实现中的aggregate_input函数存在两处关键错误：

1. glob_wildcards路径匹配错误：temp目录下的文件是无.fna后缀的基因组ID，原代码错误添加了该后缀
2. expand函数变量映射错误：使用了未定义的i变量，而非正确的genome变量

修正后的完整Snakefile #

import os
from snakemake.io import glob_wildcards

# 可替换为从config.yaml读取配置
DATABASE = '("Apis"[Organism] OR Apis[All Fields]) AND (latest[filter] AND "representative genome"[filter] AND all[filter] NOT anomalous[filter])'

rule all:
    input: "summary_table.txt"

checkpoint create_genome_list:
    output: directory("temp/")
    conda: "entrez_env.yaml"
    shell:
        r"""
        mkdir -p temp/
        esearch -db assembly -query '{DATABASE}' \
        | esummary \
        | xtract -pattern DocumentSummary -element FtpPath_GenBank \
        | while read -r line ; 
        do
            fname=$(echo $line | grep -o 'GCA_.*' | sed 's/$/_genomic.fna.gz/');
            wildcard=$(echo $fname | sed -e 's!.fna.gz!!');

            echo "$line/$fname" > temp/$wildcard;
        done
        """

rule download_genome:
    output: "database/{genome}/{genome}.fna.gz"
    input: "temp/{genome}"
    shell:
        r"""
        mkdir -p database/{wildcards.genome}/
        GENOME_LINK=$(cat {input})
        wget -P ./database/{wildcards.genome}/ $GENOME_LINK
        """

rule unzip_genome:
    output: "database/{genome}/{genome}.fna"
    input: "database/{genome}/{genome}.fna.gz"
    shell: "gunzip {input}"

def aggregate_input(wildcards):
    # 获取checkpoint输出的temp目录路径
    checkpoint_output = checkpoints.create_genome_list.get().output[0]
    # 匹配temp目录下的所有基因组ID文件
    genomes = glob_wildcards(os.path.join(checkpoint_output, "{genome}")).genome
    # 生成所有已解压基因组文件的路径列表
    return expand("database/{genome}/{genome}.fna", genome=genomes)

rule make_summary_table:
    output: "summary_table.txt"
    input: aggregate_input
    shell:
        """
        echo "已下载解压的基因组文件列表:" > {output}
        echo "" >> {output}
        for file in {input}; do
            echo "$file" >> {output}
        done
        """

关键修正说明 #

1. Checkpoint调用优化：checkpoints.create_genome_list.get()无需传入wildcards，因为该checkpoint未定义任何wildcard
2. 动态文件名匹配修复：glob_wildcards(os.path.join(checkpoint_output, "{genome}"))准确匹配temp目录下的基因组ID文件
3. Expand变量映射修正：expand("database/{genome}/{genome}.fna", genome=genomes)使用正确变量关联动态生成的基因组列表
4. 目录创建保障：在create_genome_list和download_genome的shell命令中添加mkdir -p，避免目录不存在导致的运行错误
5. 汇总输出优化：修改shell命令让汇总表格格式更清晰可读

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