【问题标题】：Is it possible to use wildcards in config files for a Snakemake pipeline?是否可以在 Snakemake 管道的配置文件中使用通配符？
【发布时间】：2020-10-24 01:49:15
【问题描述】：

我是构建 Snakefiles 的新手，对于我的生物信息学研究，我正在尝试在多个样本上循环我的规则。我寻找了类似的问题和答案，但我似乎无法解决这个问题。这可能是因为我仍然不太了解 Snakemake 是如何工作的。如果你们能帮助我，那就太好了。

目前我有多个规则，目前适用于一个示例：

# variables for every species
SAMPLE = "SRR8528338"
SAMPLES = "SRR8528338 SRR8528339 SRR8528340".split()

configfile: "./envs/contigs/" + SAMPLE + ".yaml"

var_variables = expand("results/4_mapped_contigs/" + SAMPLE + "/var/Contig{nr}_AT_sort.var", nr = config["contig_nrs"])
#make_contig_consensus = expand("results/5_consensus_contigs/{sample}", sample = SAMPLES)

rule all:
    input:
         var_variables
#        make_contig_consensus

rule convert_to_fBAM:
    input:
        "results/4_mapped_contigs/" + SAMPLE + "/sam/Contig{nr}_AT.sam"
    output:
        "results/4_mapped_contigs/" + SAMPLE + "/bam/Contig{nr}_AT.bam"
    shell:
        "samtools view -bS {input} > {output}"

rule sort_fBAM:
    input:
        "results/4_mapped_contigs/" + SAMPLE + "/bam/Contig{nr}_AT.bam"
    output:
        "results/4_mapped_contigs/" + SAMPLE + "/sorted_bam/Contig{nr}_AT_sort.bam"
    shell:
        "samtools sort -m5G {input} -o {output}"

rule convert_to_fpileup:
    input:
        "results/4_mapped_contigs/" + SAMPLE + "/sorted_bam/Contig{nr}_AT_sort.bam"
    output:
        "results/4_mapped_contigs/" + SAMPLE + "/pileup/Contig{nr}_AT_sort.pileup"
    shell:
        "samtools mpileup -B {input} > {output}"

rule SNP_calling:
    input:
        "results/4_mapped_contigs/" + SAMPLE + "/pileup/Contig{nr}_AT_sort.pileup"
    output:
        "results/4_mapped_contigs/" + SAMPLE + "/var/Contig{nr}_AT_sort.var"
    shell:
        "varscan pileup2cns {input} "
        "--min-freq-for-hom 0.6 "
        "--min-coverage 5 "
        "--min-var-freq 0.6 "
        "--p-value 0.1 "
        "--min-reads2 5 "
        "> {output}"

rule make_contig_consensus:
    input:
        "src/read_var.py"
    output:
        "results/5_consensus_contigs/{sample}"
    params:
        "{sample}"
    shell:
        "python3 {input} {params}"

每个样本的配置文件都不同（重叠群的数量）。对于 SRR8528338，它看起来像这样：

contig_nrs:
    1: ./results/4_mapped_contigs/SRR8528338/var/Contig1_AT_sort.var
    2: ./results/4_mapped_contigs/SRR8528338/var/Contig2_AT_sort.var
    3: ./results/4_mapped_contigs/SRR8528338/var/Contig3_AT_sort.var
    ...
    2146: ./results/4_mapped_contigs/SRR8528338/var/Contig2146_AT_sort.var

但是，我想在“SAMPLES”变量中引用的多个样本上循环所有这些规则。现在我之前尝试使用双括号，它适用于多个样本。（将所有 'SAMPLES' 更改为 {{sample}} 并添加：，sample = SAMPLES）。那么我的代码应该是这样的：

# variables for every species
SAMPLES = "SRR8528338 SRR8528339 SRR8528340".split()

for sample in SAMPLES:
    configfile: "./envs/contigs/" + sample + ".yaml"

var_variables = expand("results/4_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var", sample = SAMPLES, nr = config["contig_nrs"])
make_contig_consensus = expand("results/5_consensus_contigs/{sample}", sample = SAMPLES)

rule all:
    input:
         var_variables
#        make_contig_consensus

rule convert_to_fBAM:
    input:
        "results/4_mapped_contigs/{{sample}}/sam/Contig{nr}_AT.sam"
    output:
        "results/4_mapped_contigs/{{sample}}/bam/Contig{nr}_AT.bam"
    shell:
        "samtools view -bS {input} > {output}"

rule sort_fBAM:
    input:
        "results/4_mapped_contigs/{{sample}}/bam/Contig{nr}_AT.bam"
    output:
        "results/4_mapped_contigs/{{sample}}/sorted_bam/Contig{nr}_AT_sort.bam"
    shell:
        "samtools sort -m5G {input} -o {output}"

rule convert_to_fpileup:
    input:
        "results/4_mapped_contigs/{{sample}}/sorted_bam/Contig{nr}_AT_sort.bam"
    output:
        "results/4_mapped_contigs/{{sample}}/pileup/Contig{nr}_AT_sort.pileup"
    shell:
        "samtools mpileup -B {input} > {output}"

rule SNP_calling:
    input:
        "results/4_mapped_contigs/{{sample}}/pileup/Contig{nr}_AT_sort.pileup"
    output:
        "results/4_mapped_contigs/{{sample}}/var/Contig{nr}_AT_sort.var"
    shell:
        "varscan pileup2cns {input} "
        "--min-freq-for-hom 0.6 "
        "--min-coverage 5 "
        "--min-var-freq 0.6 "
        "--p-value 0.1 "
        "--min-reads2 5 "
        "> {output}"

rule make_contig_consensus:
    input:
        "src/read_var.py"
    output:
        "results/5_consensus_contigs/{sample}"
    params:
        "{sample}"
    shell:
        "python3 {input} {params}"

但是，当我运行它时，我得到一个错误。我不太确定，但我认为这是因为 for 循环（SAMPLES 中的示例）：

Missing input files for rule all:
results/4_mapped_contigs/SRR8528338/var/Contig1266_AT_sort.var
results/4_mapped_contigs/SRR8528338/var/Contig1299_AT_sort.var
...

现在我想知道：有没有办法通过使用通配符来扩展配置文件？比如：

configfile: expand("./envs/contigs/{sample}.yaml", sample = SAMPLES)

这样做会给我错误：

TypeError in line 4
expected str, bytes or os.PathLike object, not list

或者你有其他解决这个问题的方法吗？

谢谢！

更新：

我一直在尝试一些事情，我认为将配置文件更改为嵌套字典而不是单独的字典会很有用。它应该看起来像这样：

    contigs:
         SRR8528336: - 1
                     - 2
                     - ...
                     - 2113
         SRR8528337: - 1
                      ...
          ...
    exons:
         SRR8528336: - 1
                      ...
                     - 1827
         SRR8528337: - 1
                       ...
                     - 1826
          ...

例如，如果我想运行样本：SRR8528338 直到 SRR8528340，我将其作为输入：

SAMPLES = "SRR8528338 SRR8528339 SRR8528340".split()

并按样本名称调用重叠群：

var_variables = expand("results/4_mapped_contigs/{{sample}}/var/Contig{nr}_AT_sort.var", nr = config["contigs"][wildcards.sample])

或外显子：

expand("results/7_exons/{{sample}}/var/exon{nr}_AT_sort.var", nr = config["exons"][wildcards.sample])

如果我只想获取值，“wildcards.sample”是如何工作的？

解决方案（以及下一个问题）2020 年 7 月 31 日

我根据 bli 进行了更改，现在正在运行：

# variables for every species
SAMPLES = "SRR8528347 SRR8528355 SRR8528356".split()

configfile: "./envs/config_contigs.yaml"

# bam = []
# sort_bam = []
# fpileup = []
var_variables = []
make_contig_consensus = []
blat_variables = []
extract_hits_psl = []
for sample in SAMPLES:
    contig_nrs = config[sample]
    for nr in contig_nrs:
        # bam.append("results/A04_mapped_contigs/{sample}/bam/Contig{nr}_AT.bam".format(sample=sample, nr=nr))
        # sort_bam.append("results/A04_mapped_contigs/{sample}/sorted_bam/Contig{nr}_AT_sort.bam".format(sample=sample, nr=nr))
        # fpileup.append("results/A04_mapped_contigs/{sample}/pileup/Contig{nr}_AT_sort.pileup".format(sample=sample, nr=nr))
        var_variables.append("results/A04_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var".format(sample=sample, nr=nr))
        make_contig_consensus.append("results/A05_consensus_contigs/{sample}/Contig{nr}.fasta".format(sample=sample, nr=nr))
        blat_variables.append("results/A06_identified_contigs_blat/{sample}/contig{nr}_AT.psl".format(sample=sample, nr=nr))
        extract_hits_psl.append("results/A07_mapped_exons/{sample}/".format(sample=sample, nr=nr))

rule all:
    input:
        # bam,
        # sort_bam,
        # fpileup,
        var_variables,
        make_contig_consensus,
        blat_variables,
        extract_hits_psl

rule convert_to_fBAM:
    input:
        "results/A04_mapped_contigs/{sample}/sam/Contig{nr}_AT.sam"
    output:
        "results/A04_mapped_contigs/{sample}/bam/Contig{nr}_AT.bam"
    shell:
        "samtools view -bS {input} > {output}"

rule sort_fBAM:
    input:
        "results/A04_mapped_contigs/{sample}/bam/Contig{nr}_AT.bam"
    output:
        "results/A04_mapped_contigs/{sample}/sorted_bam/Contig{nr}_AT_sort.bam"
    shell:
        "samtools sort -m5G {input} -o {output}"

rule convert_to_fpileup:
    input:
        "results/A04_mapped_contigs/{sample}/sorted_bam/Contig{nr}_AT_sort.bam"
    output:
        "results/A04_mapped_contigs/{sample}/pileup/Contig{nr}_AT_sort.pileup"
    shell:
        "samtools mpileup -B {input} > {output}"

rule SNP_calling:
    input:
        "results/A04_mapped_contigs/{sample}/pileup/Contig{nr}_AT_sort.pileup"
    output:
        "results/A04_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var"
    shell:
        "varscan pileup2cns {input} "
        "--min-freq-for-hom 0.6 "
        "--min-coverage 5 "
        "--min-var-freq 0.6 "
        "--p-value 0.1 "
        "--min-reads2 5 "
        "> {output}"

rule make_contig_consensus:
    input:
        script = "src/read_var.py",
        file = "results/A04_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var"
    output:
        "results/A05_consensus_contigs/{sample}/Contig{nr}.fasta"
    params:
        "{sample}"
    shell:
        "python3 {input.script} {params}"

rule BLAT_assembled:
    input:
        "data/exons/exons_AT.fasta",
        "results/A05_consensus_contigs/{sample}/Contig{nr}.fasta"
    output:
        "results/A06_identified_contigs_blat/{sample}/contig{nr}_AT.psl"
    shell:
        "blat "
        "-t=dnax "
        "-q=dnax "
        "-stepSize=5 "
        "-repMatch=2253 "
        "-minScore=0 "
        "-minIdentity=0 "
        "{input} {output}"

rule extract_hits_psl:
    input:
        script = "src/extract_hits_psl.py"
        # file = "results/A06_identified_contigs_blat/{sample}/contig{nr}_AT.psl"
    output:
        "results/A07_mapped_exons/{sample}/"
    params:
        "{sample}"
    shell:
        "python {input.script} {params}"

config_contigs.yaml：

SRR8528347:
    - 1
    - ...
    - 5
SRR8528348:
    - 1
    - ...
    - 5
...

现在从 .yaml 调用它们是有效的，但规则应该按照与编写的顺序相同的顺序运行（从上到下）。运行此命令时，规则以不同的顺序运行，因此由于文件尚不存在而出现错误。我读到之前订单的输出应该和之后的输入一样，但是它不起作用。

【问题讨论】：

关于您的最后一个错误：expand 创建一个字符串列表，但 configfile 应该是文件的路径（可以表示为字符串），而不是字符串列表.如果您实际上每个样本都有一个配置文件，这意味着您的工作流程一次只能处理一个样本。我写了这个答案，我试图解释各种有用的东西来了解蛇的工作原理：stackoverflow.com/questions/50198628/…也许它对你有用。
谢谢！我看了看，确实帮助了我更多的理解，但我仍然卡住了。我看到你使用类似'config [sample] [wildcards.sample]'的东西。我想如果我可以像这样从配置文件中获取我的值，那会有所帮助。但我不明白如何获取这些值（请参阅上面的编辑帖子）。
我不确定您关于执行顺序的问题。你可以再详细一点吗？在 Snakemake 中，执行顺序将根据规则之间的显式依赖关系（即输入和输出部分）确定。如果没有这样的依赖，就无法保证执行顺序。
所以当我运行这个蛇文件时，它会运行规则 convert_to_fBAM 直到 convert_to_fpileup。但它不是在之后运行 SNP_calling，而是跳过此规则并转到下一个规则“make_contig_consensus”。然后会发生错误，因为“results/A04_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var”不存在。我将当前代码添加到主要问题 (31/07/2020)
这很奇怪。据我所知，make_contig_consensus 的输入依赖于SNP_calling 的输出，所以 Snakemake 不应该直接跳到make_contig_consensus。您可以在帖子中添加确切的错误消息吗？

标签： python config wildcard bioinformatics snakemake

【解决方案1】：

这里尝试回答您在上次编辑（2020 年 7 月 23 日）中提出的问题。

如果您的配置文件如您所见，您应该有一个 config Python 字典，其中包含一个 "contigs" 条目，该条目又是一个可以将数字列表与示例标识符相关联的字典。

您还生成一个样本标识符列表，并将其存储在列表SAMPLES中。

为了“驱动”管道，我们决定哪些文件将成为最终目标，这将作为all 规则的input 给出。在这里，您要生成遵循此模式的输出文件列表：

"results/4_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var"

{sample} 占位符将采用 SAMPLE 列表中的值，并且其中{nr} 占位符的取值取决于{sample} 占位符的值。

（注意：我使用“占位符”而不是“通配符”来帮助明确这里的蛇形 wildcards 对象不涉及“魔法”。wildcards 对象在每个规则实例中定义，但在这里我们只是在“普通” Python 中准备文件列表，超出规则。）

让我们用普通的 Python 来做，而不是使用 expand。

从一个空列表开始：

var_variables = []

现在为样本创建一个主循环：

for sample in SAMPLES:
    # Note: No `wildcards` here. We are not in a rule.
    # `sample` is just a (loop) variable in Python,
    # referring to one of the character strings in `SAMPLES`
    contig_nrs = config["contigs"][sample]
    for nr in contig_nrs:
        # I'm using python 3.6+ "formatted strings" to have the placeholders
        # replaced with the corresponding variable's values
        var_variables.append(f"results/4_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var")

现在，在这个双循环完成后，我们可以将这个文件列表提供给驱动规则：

rule all:
    input:
         var_variables

可能有一些方法可以在某处使用嵌套的expands 和双花括号语法 ({{}}) 来做同样的事情，但是如果您了解 Python 的一些基础知识（列表、字典、字符串...），你可以不用它，而且它可能更容易理解。（expand 是 Snakemake 提供的实用程序，标准 Python 中不存在。）

对于wildcards 对象，它们将仅存在于规则内部，现有属性将通过将output 部分的模式与规则实例必须生成的实际文件的名称相匹配来确定.

关于这方面的更多信息，您可以查看我为改进 Snakemake 的文档而提出的这个建议：https://bitbucket.org/snakemake/snakemake/pull-requests/307/documentation-to-help-understand-expand/diff

【讨论】：

非常感谢您帮助我。它与主循环一起工作。虽然，这适用于我的笔记本电脑。但是，当我尝试在 ssh 服务器上运行此脚本时（通过 git clone 完全相同的脚本），我在“expand”行上收到“无效语法”错误，当我从 expand( f"...") 行。虽然，如果我删除它，该行将不再用作通配符。你知道如何解决这个问题吗？
@Elboo 语法无效可能是因为服务器上的 Python 版本低于 3.6。不要在字符串前加上“f”，而是使用.format 方法："results/4_mapped_contigs/{sample}/var/Contig{nr}_AT_sort.var".format(sample=sample, nr=nr)
其实我还有一个问题。现在可以通过 de config 文件调用它们，但规则应该以正确的顺序运行。但是，它不再这样做了。我的脚本发布在主要问题中。我知道之前的规则的输出应该和下一条规则的输入一样才能让它们以正确的顺序运行，但是它不起作用。