Snakemake：关于如何正确访问配置文件的困惑

Question

这个问题是我之前提出的一个问题，它涉及了解如何使用Snakemake正确访问配置文件。我有一个特定的问题，我需要首先解决，而在理解索引工作原理方面则是一个普遍的问题，其次我要问。

我正在使用snakemake进行运行，并从Alignment / QC到动机分析运行ATAC-seq管道。

答：具体问题

我正在尝试添加一条规则，trim_galore_pe以在对齐之前从我的fastq文件中修剪适配器，并从snakemake引发错误语句，因为生成的输出文件的名称与snakemake期望的名称trim galore不匹配。这是因为我无法解决如何在snakemake文件中正确写入输出文件语句以使名称匹配。

由TRIM GALORE包含SRA号生成的名称的示例，例如：

trimmed_fastq_files/SRR2920475_1_val_1.fq.gz

snakemake期望的文件包含示例引用，并且应显示为：

trimmed_fastq_files/Corces2016_4983.7B_Mono_1_val_1.fq.gz

这也会影响规则之后的后续规则trim_galore_pe。我需要找到一种方法来使用配置文件中的信息来生成所需的输出文件。

对于Snakefile中显示的规则之后的所有规则，我需要使用示例名称（即）来命名文件Corces2016_4983.7A_Mono。对于下面的Snakefile中显示的所有FAST_QC和MULTIQC规则，在输出文件名结构中具有示例名称也将很有用，它们在当前Snakefile中都已经执行了。

但是，Bowtie2，FASTQC规则的输入以及规则的输入和输出trim_galore_pe需要包含SRA编号。问题始于trim_galore并影响所有下游规则。

尽管我在以前的规则中已经提取了SRA编号，但是当不使用fastq_files配置文件中明确说明的文件夹时，我不确定如何执行此操作。通过引入trim_galore_pe规则，我已将一组新的SRA文件有效地移到了新文件trimmed_fastq_files夹中。如何从包含旧文件夹名称的SRA文件配置文件列表中仅提取SRA编号，同时引用trimmed_fastq_filesSnakefile中的新文件夹是我问题的症结所在。

我希望这很清楚。

这是我的配置文件：

samples:
    Corces2016_4983.7A_Mono: fastq_files/SRR2920475
    Corces2016_4983.7B_Mono: fastq_files/SRR2920476
cell_types:
    Mono:
    - Corces2016_4983.7A
index: /home/genomes_and_index_files/hg19

这是我的Snakefile：

# read config info into this namespace
configfile: "config.yaml"
print (config['samples'])

rule all:
    input:
        expand("FastQC/PRETRIM/{sample}_{num}_fastqc.zip", sample=config["samples"], num=['1', '2']),
        expand("bam_files/{sample}.bam", sample=config["samples"]),
        "FastQC/PRETRIM/fastq_multiqc.html",
        "FastQC/POSTTRIM/fastq_multiqc.html"

rule fastqc_pretrim:
    input:
        sample=lambda wildcards: f"{config['samples'][wildcards.sample]}_{wildcards.num}.fastq.gz"
    output:
        # Output needs to end in '_fastqc.html' for multiqc to work
        html="FastQC/PRETRIM/{sample}_{num}_fastqc.html",
        zip="FastQC/PRETRIM/{sample}_{num}_fastqc.zip"
    wrapper:
        "0.23.1/bio/fastqc"

rule multiqc_fastq_pretrim:
    input:
        expand("FastQC/PRETRIM/{sample}_{num}_fastqc.html", sample=config["samples"], num=['1', '2'])
    output:
        "FastQC/PRETRIM/fastq_multiqc.html"
    wrapper:
        "0.23.1/bio/multiqc"

rule trim_galore_pe:
    input:
        sample=lambda wildcards: expand(f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])
    output:
        "trimmed_fastq_files/{sample}_1_val_1.fq.gz",
        "trimmed_fastq_files/{sample}_1.fastq.gz_trimming_report.txt",
        "trimmed_fastq_files/{sample}_2_val_2.fq.gz",
        "trimmed_fastq_files/{sample}_2.fastq.gz_trimming_report.txt"
    params:
        extra="--illumina -q 20"
    log:
        "logs/trim_galore/{sample}.log"
    wrapper:
        "0.23.1/bio/trim_galore/pe"

rule fastqc_posttrim:
    input:
        "trimmed_fastq_files/{sample}_1_val_1.fq.gz", "trimmed_fastq_files/{sample}_2_val_2.fq.gz"
    output:
        # Output needs to end in '_fastqc.html' for multiqc to work
        html="FastQC/POSTTRIM/{sample}_{num}_fastqc.html",
        zip="FastQC/POSTTRIM/{sample}_{num}_fastqc.zip"
    wrapper:
        "0.23.1/bio/fastqc"

rule multiqc_fastq_posttrim:
    input:
        expand("FastQC/POSTTRIM/{sample}_{num}.trim_fastqc.html", sample=config["samples"], num=['1', '2'])
    output:
        "FastQC/POSTTRIM/fastq_multiqc.html"
    wrapper:
        "0.23.1/bio/multiqc"

rule bowtie2:
    input:
        "trimmed_fastq_files/{sample}_1_val_1.fq.gz", "trimmed_fastq_files/{sample}_2_val_2.fq.gz"
    output:
        "bam_files/{sample}.bam"
    log:
        "logs/bowtie2/{sample}.txt"
    params:
       index=config["index"],  # prefix of reference genome index (built with bowtie2-build),
       extra=""
    threads: 8
    wrapper:
        "0.23.1/bio/bowtie2/align"

这当前正在运行，并使用提供完整的作业列表snakemake -np，但会引发上述错误。

B：一般问题

是否有在线资源简要说明了如何使用python来引用配置文件，尤其是关于蛇形的引用？在线文档还远远不够，并假设您具有python的先验知识。

我的编程经验主要是在bash和R中，但是我喜欢Snakemake，并且通常了解字典和列表在python中的工作方式以及如何引用其中存储的项目。但是，我发现上述某些Snakemake规则中括号，通配符和反逗号的复杂用法令人困惑，因此当尝试在配置文件中引用文件名的不同部分时，往往会遇到困难。我想充分了解如何利用这些元素。

例如，在上面发布的Snakefile中的规则中：

sample=lambda wildcards: expand(f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2]) 

该命令实际发生了什么？我的理解是，我们正在使用来访问配置文件，config['samples']并且我们正在使用该[wildcards.sample]部分来显式访问fastq_files/SRR2920475配置文件的一部分。扩展允许我们遍历配置文件中适合命令参数的每个项目，即所有SRA文件，并且需要lambda通配符才能使用wildcards命令中的调用。我不确定的是：

1. f扩展后该怎么办？为什么需要它？
2. 为什么config['samples']在方括号内包含反向逗号，但周围不需要反向逗号[wildcards.sample]？
3. 为什么使用单双括号？
4. 在上面的Snakefile中，一些规则包含将数字序列分配给的部分num，但是这些数字有时又被包围在反引号中，有时却不是...为什么？

任何建议，技巧，指针将不胜感激。

C：澄清@bli在下面提出的建议

我已按照您在注释中的建议编辑了配置文件，并省略了文件夹名称，仅保留了SRA编号。这对我来说很有意义，但是我还有其他一些问题使我无法运行此Snakefile。

新的配置文件：

samples:
    Corces2016_4983.7A_Mono: SRR2920475
    Corces2016_4983.7B_Mono: SRR2920476
cell_types:
    Mono:
    - Corces2016_4983.7A
index: /home/c1477909/genomes_and_index_files/hg19

新的Snakefile：

# read config info into this namespace
configfile: "config.yaml"
print (config['samples'])

rule all:
    input:
        expand("FastQC/PRETRIM/{sample}_{num}_fastqc.zip", sample=config["samples"], num=['1', '2']),
        expand("bam_files/{sample}.bam", sample=config["samples"]),
        "FastQC/PRETRIM/fastq_multiqc.html",
        "FastQC/POSTTRIM/fastq_multiqc.html",

rule fastqc_pretrim:
    input:
      lambda wildcards: f"fastq_files/{config['samples'][wildcards.sample]}_{wildcards.num}.fastq.gz"
    output:
        # Output needs to end in '_fastqc.html' for multiqc to work
        html="FastQC/PRETRIM/{sample}_{num}_fastqc.html",
        zip="FastQC/PRETRIM/{sample}_{num}_fastqc.zip"
    wrapper:
        "0.23.1/bio/fastqc"

rule multiqc_fastq_pretrim:
    input:
        expand("FastQC/PRETRIM/{sample}_{num}_fastqc.html", sample=config["samples"], num=['1', '2'])
    output:
        "FastQC/PRETRIM/fastq_multiqc.html"
    wrapper:
        "0.23.1/bio/multiqc"

rule trim_galore_pe:
    input:
        lambda wildcards: expand(f"fastq_files/{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])
    output:
        "trimmed_fastq_files/{wildcards.sample}_1_val_1.fq.gz",
        "trimmed_fastq_files/{wildcards.sample}_1.fastq.gz_trimming_report.txt",
        "trimmed_fastq_files/{wildcards.sample}_2_val_2.fq.gz",
        "trimmed_fastq_files/{wildcards.sample}_2.fastq.gz_trimming_report.txt"
    params:
        extra="--illumina -q 20"
    log:
        "logs/trim_galore/{sample}.log"
    wrapper:
        "0.23.1/bio/trim_galore/pe"

rule fastqc_posttrim:
    input:
        lambda wildcards: expand(f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}_val_{{num}}.fq.gz", num=[1,2])
    output:
        # Output needs to end in '_fastqc.html' for multiqc to work
        html="FastQC/POSTTRIM/{sample}_{num}_fastqc.html",
        zip="FastQC/POSTTRIM/{sample}_{num}_fastqc.zip"
    wrapper:
        "0.23.1/bio/fastqc"

rule multiqc_fastq_posttrim:
    input:
        expand("FastQC/POSTTRIM/{sample}_{num}.trim_fastqc.html", sample=config["samples"], num=['1', '2'])
    output:
        "FastQC/POSTTRIM/fastq_multiqc.html"
    wrapper:
        "0.23.1/bio/multiqc"

rule bowtie2:
    input:
        lambda wildcards: expand(f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}_val_{{num}}.fq.gz", num=[1,2])
    output:
        "bam_files/{sample}.bam"
    log:
        "logs/bowtie2/{sample}.txt"
    params:
        index=config["index"],  # prefix of reference genome index (built with bowtie2-build),
        extra=""
    threads: 8
    wrapper:
        "0.23.1/bio/bowtie2/align"

使用这些新文件，最初一切正常，由创建了部分作业列表snakemake -np。但是，这是因为完整作业列表的一半已经运行；即trimmed_fastq_files生成了文件夹，并且其中已正确命名了修剪过的fastq文件。当我删除所有先前创建的文件以查看Snakefile的整个新版本是否正常运行时，snakemake -np失败，并指出该规则下游的trim_galore_pe规则缺少输入文件。

如您所见，我正在尝试在输出部分{wildcard.sample}的trim_galore_pe规则的输入部分中调用变量集，但是snakemake不喜欢这样。有可能这样做吗？

我还使用以下答案中的技巧尝试了此操作，但这也不起作用：

rule trim_galore_pe:
    input:
        sample=lambda wildcards: expand(f"fastq_files/{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])
    output:
        expand(f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}_val_{{num}}.fq.gz", num=[1,2]),
        expand(f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}.fastq.gz_trimming_report.txt", num=[1,2])
    params:
        extra="--illumina -q 20"
    log:
        "logs/trim_galore/{sample}.log"
    wrapper:
        "0.23.1/bio/trim_galore/pe"

错误然后指出wildcards not defined。因此，从逻辑上讲，我尝试lambda wildcards:在输出部分的两个expand语句之前尝试定义通配符，但这会引发语法错误Only input files can be specified as functions。我还尝试使用下面的一些索引建议，但无法获得正确的组合。

这可能是由我不确定有关Snakefiles的另一件事引起的，这就是作用域确定的方式。

• 如果在中定义变量，rule all所有其他规则都可以访问它吗？
• 如果我在规则的输入部分中定义了一个变量，该变量是否可用于该规则的所有其他部分（即输出，shell命令等），但仅适用于该规则？
• 如果是，{wildcard.sample}如果在输入部分定义了变量，为什么不能访问该变量？是因为该变量包含在“封闭”范围的lambda函数中？

任何（进一步的）建议将不胜感激。

Answer 1

我会尽量回答你的问题 B，并提供额外的细节，希望对你和其他人有用。

编辑：我在最后添加了一些回答问题 C 的尝试。

关于报价

首先，关于你所说的“倒置”逗号，它们通常被称为“单引号”，它们在python中用于构建字符串。双引号也可以用于相同的目的。主要区别在于当您尝试创建包含引号的字符串时。使用双引号可以创建包含单引号的字符串，反之亦然。否则，您需要使用反斜杠（“\”）“转义”引号：

s1 = 'Contains "double quotes"'
s1_bis = "Contains \"double quotes\""
s2 = "Contains 'single quotes'"
s2_bis = 'Contains \'single quotes\''

（我倾向于更喜欢双引号，这只是个人品味。）

分解示例

rule trim_galore_pe:
    input:
        sample=lambda wildcards: expand(f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])

您正在将函数 ( lambda wildcards: ...)分配给变量 ( sample)，该变量恰好属于规则的输入部分。

这将导致snakemake在确定规则的特定实例的输入时使用此函数，基于通配符的当前值（从它想要生成的输出的当前值推断）。

为清楚起见，很可能通过将函数定义与规则声明分开来重写这一点，而不使用lambda构造，并且它的工作方式相同：

def determine_sample(wildcards):
    return expand(
        f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz",
        num=[1,2])

rule trim_galore_pe:
    input:
        sample = determine_sample

expand是一个snakemake特定的函数（但你可以在任何python程序或交互式解释器中使用 导入它from snakemake.io import expand），这使得生成字符串列表更容易。在下面的交互式 python3.6 会话中，我们将尝试使用不同的本机 python 构造重现使用它时发生的情况。

访问配置

# We'll try to see how `expand` works, we can import it from snakemake
from snakemake.io import expand
    ?
# We want to see how it works using the following example
# expand(f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])

# To make the example work, we will first simulate the reading
# of a configuration file
import yaml

config_text = """
samples:
    Corces2016_4983.7A_Mono: fastq_files/SRR2920475
    Corces2016_4983.7B_Mono: fastq_files/SRR2920476
cell_types:
    Mono:
    - Corces2016_4983.7A
index: /home/genomes_and_index_files/hg19
"""
# Here we used triple quotes, to have a readable multi-line string.
?
# The following is equivalent to what snakemake does with the configuration file:
config = yaml.load(config_text)
config

输出：

{'cell_types': {'Mono': ['Corces2016_4983.7A']},
 'index': '/home/genomes_and_index_files/hg19',
 'samples': {'Corces2016_4983.7A_Mono': 'fastq_files/SRR2920475',
  'Corces2016_4983.7B_Mono': 'fastq_files/SRR2920476'}}

我们获得了一个字典，其中键“样本”与嵌套字典相关联。

# We can access the nested dictionary as follows
config["samples"]
# Note that single quotes could be used instead of double quotes
# Python interactive interpreter uses single quotes when it displays strings

输出：

{'Corces2016_4983.7A_Mono': 'fastq_files/SRR2920475',
 'Corces2016_4983.7B_Mono': 'fastq_files/SRR2920476'}

# We can access the value corresponding to one of the keys
# again using square brackets
config["samples"]["Corces2016_4983.7A_Mono"]

输出：

'fastq_files/SRR2920475'

# Now, we will simulate a `wildcards` object that has a `sample` attribute
# We'll use a namedtuple for that
# https://docs.python.org/3/library/collections.html#collections.namedtuple
from collections import namedtuple
Wildcards = namedtuple("Wildcards", ["sample"])
wildcards = Wildcards(sample="Corces2016_4983.7A_Mono")
wildcards.sample

输出：

'Corces2016_4983.7A_Mono'

编辑（15/11/2018）：我发现了一种更好的创建通配符的方法：

from snakemake.io import Wildcards
wildcards = Wildcards(fromdict={"sample": "Corces2016_4983.7A_Mono"})

# We can use this attribute as a key in the nested dictionary
# instead of using directly the string
config["samples"][wildcards.sample]
# No quotes here: `wildcards.sample` is a string variable

输出：

'fastq_files/SRR2920475'

解构 expand

# Now, the expand of the example works, and it results in a list with two strings
expand(f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])    
# Note: here, single quotes are used for the string "sample",
# in order not to close the opening double quote of the whole string

输出：

['fastq_files/SRR2920475_1.fastq.gz', 'fastq_files/SRR2920475_2.fastq.gz']

# Internally, I think what happens is something similar to the following:
filename_template = f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz"

# This template is then used for each element of this "list comprehension"    
[filename_template.format(num=num) for num in [1, 2]]

输出：

['fastq_files/SRR2920475_1.fastq.gz', 'fastq_files/SRR2920475_2.fastq.gz']

# This is equivalent to building the list using a for loop:
filenames = []
for num in [1, 2]:
    filename = filename_template.format(num=num)
    filenames.append(filename)
filenames

输出：

['fastq_files/SRR2920475_1.fastq.gz', 'fastq_files/SRR2920475_2.fastq.gz']

字符串模板和格式

# It is interesting to have a look at `filename_template`    
filename_template

输出：

'fastq_files/SRR2920475_{num}.fastq.gz'

# The part between curly braces can be substituted
# during a string formatting operation:
"fastq_files/SRR2920475_{num}.fastq.gz".format(num=1)

输出：

'fastq_files/SRR2920475_1.fastq.gz'

现在让我们进一步展示如何使用字符串格式。

# In python 3.6 and above, one can create formatted strings    
# in which the values of variables are interpreted inside the string    
# if the string is prefixed with `f`.
# That's what happens when we create `filename_template`:
filename_template = f"{config['samples'][wildcards.sample]}_{{num}}.fastq.gz"    
filename_template    

输出：

'fastq_files/SRR2920475_{num}.fastq.gz'

在格式化字符串期间发生了两次替换：

1. 的值config['samples'][wildcards.sample]用于制作字符串的第一部分。（使用单引号是sample因为这个 python 表达式位于用双引号构建的字符串中。）

2. num作为格式化操作的一部分，周围的双括号被减少为单个括号。这就是为什么我们可以在涉及num.

# Equivalently, without using 3.6 syntax:    
filename_template = "{filename_prefix}_{{num}}.fastq.gz".format(
    filename_prefix = config["samples"][wildcards.sample])
filename_template

输出：

'fastq_files/SRR2920475_{num}.fastq.gz'

# We could achieve the same by first extracting the value
# from the `config` dictionary    
filename_prefix = config["samples"][wildcards.sample]
filename_template = f"{filename_prefix}_{{num}}.fastq.gz"
filename_template

输出：

'fastq_files/SRR2920475_{num}.fastq.gz'

# Or, equivalently:
filename_prefix = config["samples"][wildcards.sample]
filename_template = "{filename_prefix}_{{num}}.fastq.gz".format(
    filename_prefix=filename_prefix)
filename_template

输出：

'fastq_files/SRR2920475_{num}.fastq.gz'

# We can actually perform string formatting on several variables
# at the same time:
filename_prefix = config["samples"][wildcards.sample]
num = 1
"{filename_prefix}_{num}.fastq.gz".format(
    filename_prefix=filename_prefix, num=num)

输出：

'fastq_files/SRR2920475_1.fastq.gz'

# Or, using 3.6 formatted strings
filename_prefix = config["samples"][wildcards.sample]
num = 1
f"{filename_prefix}_{num}.fastq.gz"

输出：

'fastq_files/SRR2920475_1.fastq.gz'

# We could therefore build the result of the expand in a single step:
[f"{config['samples'][wildcards.sample]}_{num}.fastq.gz" for num in [1, 2]]

输出：

['fastq_files/SRR2920475_1.fastq.gz', 'fastq_files/SRR2920475_2.fastq.gz']

对问题 C 的评论

就 Python 将如何构建字符串而言，以下内容有点复杂：

input:
    lambda wildcards: expand(f"fastq_files/{config['samples'][wildcards.sample]}_{{num}}.fastq.gz", num=[1,2])

但它应该可以工作，正如我们在以下模拟中看到的那样：

from collections import namedtuple
from snakemake.io import expand

Wildcards = namedtuple("Wildcards", ["sample"])
wildcards = Wildcards(sample="Corces2016_4983.7A_Mono")
config = {"samples": {
    "Corces2016_4983.7A_Mono": "SRR2920475",
    "Corces2016_4983.7B_Mono": "SRR2920476"}}
expand(
    f"fastq_files/{config['samples'][wildcards.sample]}_{{num}}.fastq.gz",
    num=[1,2])

输出：

['fastq_files/SRR2920475_1.fastq.gz', 'fastq_files/SRR2920475_2.fastq.gz']

trim_galore_pe规则中的问题实际上出在它的output部分：你不应该{wildcards.sample}在那里使用，而只是{sample}.

output规则的部分是您通过将要获取的文件与给定的模式进行匹配来通知snakemake 通配符属性将用于规则的给定实例的地方。与花括号匹配的部分将用于设置相应属性名称的值。

例如，如果 snakemake 想要一个名为 的文件"trimmed_fastq_files/Corces2016_4983.7A_Mono_1_val_1.fq.gz"，它会尝试将它与所有规则输出部分中存在的所有模式进行匹配，并最终找到这个："trimmed_fastq_files/{sample}_1_val_1.fq.gz"

幸运的是，它可以通过在Corces2016_4983.7A_Mono和{sample}部分之间建立对应关系来将文件名与模式匹配。然后它会sample在本地通配符实例中放置一个属性，有点像我手动执行以下操作：

Wildcards = namedtuple("Wildcards", ["sample"])
wildcards = Wildcards(sample="Corces2016_4983.7A_Mono")

我不知道如果你使用{wildcards.sample}而不是在snakemake中会发生什么{wildcards}，但让我们尝试使用我的模拟框架：

Wildcards = namedtuple("Wildcards", ["sample"])
wildcards = Wildcards(wildcards.sample="Corces2016_4983.7A_Mono")
  File "<ipython-input-12-c02ce12bff85>", line 1
    wildcards = Wildcards(wildcards.sample="Corces2016_4983.7A_Mono")
                         ^
SyntaxError: keyword can't be an expression

你接下来的尝试呢？

output:
    expand(f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}_val_{{num}}.fq.gz", num=[1,2]),

在这里，我的理解是 Python 首先尝试将f字符串格式应用于f"trimmed_fastq_files/{config['samples'][wildcards.sample]}_{{num}}_val_{{num}}.fq.gz". 为此，它需要能够评估config['samples'][wildcards.sample]，但该wildcards对象尚不存在。因此wildcards not defined。 wildcards只有在将“下游”规则所需的文件名与包含{attribute_name}模式的字符串匹配后才会生成。但这是snakemake目前正在尝试构建的字符串。

以下是一些需要记住的要点：

• wildcards 在将其输出与另一个“下游”规则实例所需的文件匹配后，实际上仅存在于本地，在规则实例中。
• 您不在输入部分中定义变量。您使用变量来构建规则实例将需要的文件的具体名称（或者，更准确地说，您说您希望在运行规则实例之前存在：规则不需要实际使用这些文件）。这些变量是那些在规则范围之外定义的变量，直接在蛇文件的底层、纯 Python 模式和本地wildcards对象。默认情况下，{attribute_name}占位符将被本地wildcards对象的属性替换（"{sample}"变成"Corces2016_4983.7A_Mono"），但是如果你想做更复杂的事情来构建文件名，你需要通过一个必须显式处理这个wildcards对象的函数来完成（lambda wildcards: f"{wildcards.sample}"变成"Corces2016_4983.7A_Mono"）。