Snakemake动态cluster目录wildcard正则匹配失效问题求助

问题场景与代码 #

我在Snakemake中尝试为ReferenceDatabase规则设置Wildcard，匹配动态生成的Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}格式目录（dirname对应cluster1、cluster2这类运行前无法预知数量的目录），编写的Snakefile如下：

import glob

# Need sample name and dirname
SAMPLES, = glob_wildcards("Campylobacter/core_genome/core/{sample}.fa.align")

dirnames, = glob_wildcards("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}", "Campylobacter/Gene_Flow/DatabaseQuery/{dirname}/{dirname}")

wildcard_constraints:
    dirname="cluster[0-9]+"

rule all:
    input:
        distmat_out = "Campylobacter/ANI_results/ani/ani.distmat",
        parse_distances_out = "Campylobacter/ANI_results/genome_pairs.csv",
        cluster_genomes_out = "Campylobacter/ANI_results/cluster_genomes.csv",
        liste_genomes = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt", dirname=dirnames),
        core_genome_within_species = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/concat.fa", dirname=dirnames),
        distances_between_genomes_r = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/distances.dist", dirname=dirnames)

rule define_ANI_species:
    input:
        fasta = "Campylobacter/core_genome/concat.fa",
        dir = "Campylobacter"
    output:
        distmat = "Campylobacter/ANI_results/ani/ani.distmat",
        parse_distances = "Campylobacter/ANI_results/genome_pairs.csv",
        cluster_genomes = "Campylobacter/ANI_results/cluster_genomes.csv",
    shell:
        """
        mkdir -p Campylobacter/ANI_results/ani
        distmat  -sequence {input.fasta} -nucmethod 0 -outfile {output.distmat}
        python pipelines/ANI/parse_distances.py {input.dir}
        python pipelines/ANI/cluster_genomes.py {input.dir}
        """

rule ReferenceDatabase:
    input:
        cluster_genomes = "Campylobacter/ANI_results/cluster_genomes.csv",
        dir = "Campylobacter"
    output:
        liste = "Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt"
    shell:
        "python pipelines/ConSpecifix/create_Refdb.py {input.dir}"

rule core_genome_within_species:
    input:
        dir = "Campylobacter/genomes",
        liste = "Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt"
    output:
        fasta = "Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/concat.fa",
        family = "Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/families_core.txt"
    params:
        dir = directory("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome")
    shell:
        "python pipelines/CoreCruncher/corecruncher_master.py -in {input.dir} -out {params.dir} -list {input.liste} -freq 85 -prog usearch -ext .fa -length 80 -score 70 -align mafft"

运行错误 #

运行后出现如下错误：

rule ReferenceDatabase:
input: Campylobacter/ANI_results/genome_clusters.csv, Campylobacter
output: Campylobacter/Gene_Flow/ReferenceDatabase/cluster[0-9]+/path_to_genome_list.txt
jobid: 18
wildcards: dirname=cluster[0-9]+
Waiting at most 5 seconds for missing files.
MissingOutputException in line 171 of /Users/home//Bioinformatic_tool/Snakefile:
Job completed successfully, but some output files are missing. Missing files after 5 seconds:
Campylobacter/Gene_Flow/ReferenceDatabase/cluster[0-9]+/path_to_genome_list.txt
This might be due to filesystem latency. If that is the case, consider to increase the wait time with --latency-wait.

用户疑问：看起来Snakemake未识别用于dirname的正则表达式"[0-9]+"，是否存在可匹配cluster1、cluster2等带整数后缀目录的Wildcard？

问题根源 #

1. 提前读取不存在的目录：运行初期Campylobacter/Gene_Flow/ReferenceDatabase/下还没有生成任何cluster目录，glob_wildcards无法捕获到有效dirname值，导致后续expand无法生成具体文件路径，Snakemake只能将Wildcard约束的正则表达式cluster[0-9]+当作实际值去查找文件，自然找不到。
2. 规则输出与实际脚本不匹配：ReferenceDatabase规则声明了带{dirname}的输出路径，但shell脚本并没有根据cluster信息生成对应目录和文件，相当于规则只承诺了生成匹配正则的路径，但实际没有输出具体的cluster1、cluster2等目录下的文件。

修正方案 #

1. 优化Wildcard约束 #

保留正则约束，改用原始字符串避免转义问题：

wildcard_constraints:
    dirname=r"cluster\d+"

2. 从cluster_genomes.csv提取动态dirname列表 #

cluster_genomes.csv是define_ANI_species规则的输出，包含所有cluster信息，我们可以从这个文件中读取dirname，替代提前用glob_wildcards读取不存在的目录：

def get_dirnames():
    import csv
    dirnames = set()
    with open("Campylobacter/ANI_results/cluster_genomes.csv", "r") as f:
        reader = csv.reader(f)
        next(reader)  # 跳过表头，根据实际csv结构调整
        for row in reader:
            # 假设csv第一列是cluster名称（如cluster1），根据实际列位置修改
            cluster = row[0]
            dirnames.add(cluster)
    return sorted(dirnames)

3. 修改rule all的输入生成 #

用上面的函数动态获取dirname列表，替换原来的expand：

rule all:
    input:
        distmat_out = "Campylobacter/ANI_results/ani/ani.distmat",
        parse_distances_out = "Campylobacter/ANI_results/genome_pairs.csv",
        cluster_genomes_out = "Campylobacter/ANI_results/cluster_genomes.csv",
        liste_genomes = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt", dirname=get_dirnames()),
        core_genome_within_species = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/concat.fa", dirname=get_dirnames()),
        distances_between_genomes_r = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/distances.dist", dirname=get_dirnames())

4. 修正ReferenceDatabase规则 #

当前规则无法同时生成多个不同dirname的输出，需要修改脚本让它为每个cluster生成对应文件，同时更新规则输出为所有cluster对应的文件：

rule ReferenceDatabase:
    input:
        cluster_genomes = "Campylobacter/ANI_results/cluster_genomes.csv",
        root_dir = "Campylobacter"
    output:
        expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt", dirname=get_dirnames())
    shell:
        "python pipelines/ConSpecifix/create_Refdb.py {input.cluster_genomes} {input.root_dir}"

同时修改create_Refdb.py，接收cluster_genomes.csv路径和根目录作为参数，遍历每个cluster，创建对应的目录并生成path_to_genome_list.txt文件。

5. 可选：用Checkpoint处理完全动态的生成 #

如果cluster信息无法提前从csv获取，可使用Snakemake的checkpoint功能先生成动态目录，再收集输出路径：

checkpoint ReferenceDatabase:
    input:
        cluster_genomes = "Campylobacter/ANI_results/cluster_genomes.csv",
        root_dir = "Campylobacter"
    output:
        directory("Campylobacter/Gene_Flow/ReferenceDatabase")
    shell:
        "python pipelines/ConSpecifix/create_Refdb.py {input.root_dir}"

def gather_refdb_files():
    import os, glob
    cluster_dirs = glob.glob("Campylobacter/Gene_Flow/ReferenceDatabase/cluster*")
    dirnames = [os.path.basename(d) for d in cluster_dirs]
    return expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/path_to_genome_list.txt", dirname=dirnames)

rule all:
    input:
        distmat_out = "Campylobacter/ANI_results/ani/ani.distmat",
        parse_distances_out = "Campylobacter/ANI_results/genome_pairs.csv",
        cluster_genomes_out = "Campylobacter/ANI_results/cluster_genomes.csv",
        liste_genomes = gather_refdb_files,
        core_genome_within_species = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/concat.fa", dirname=lambda wc: [os.path.basename(d) for d in glob.glob("Campylobacter/Gene_Flow/ReferenceDatabase/cluster*")]),
        distances_between_genomes_r = expand("Campylobacter/Gene_Flow/ReferenceDatabase/{dirname}/core_genome/distances.dist", dirname=lambda wc: [os.path.basename(d) for d in glob.glob("Campylobacter/Gene_Flow/ReferenceDatabase/cluster*")])

Checkpoint会先运行生成动态目录，再收集实际生成的文件路径作为后续规则的输入。

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