snakemake - 用匿名函数替换输入指令中的通配符

Question

我正在编写一个snakemake，它将为多个输入样本运行生物信息学管道。这些输入文件（每个分析两个，一个与部分字符串匹配R1，第二个与部分字符串匹配R2）以一个模式开始，以扩展名.fastq.gz 结束。不过，最终我想执行多个操作，对于这个示例，我只想使用 bwa mem 将 fastq 读取与参考基因组对齐。所以对于这个例子，我的输入文件是NIPT-N2002394-LL_S19_R1_001.fastq.gz，我想生成NIPT-N2002394-LL.bam（参见下面的代码，指定输入和输出所在的目录）。

我的config.yaml 文件如下所示：

# Run_ID
run: "200311_A00154_0454_AHHHKMDRXX"

# Base directory: the analysis directory from which I will fetch the samples
bd: "/nexusb/nipt/"


# Define the prefix
# will be used to subset the folders in bd
prefix: "NIPT"

# Reference:
ref: "/nexus/bhinckel/19/ONT_projects/PGD_breakpoint/ref_hg19_local/hg19_chr1-y.fasta"

下面是我的蛇文件

import os
import re
#############
# config file
#############
configfile: "config.yaml"


#######################################
# Parsing variables from config.yaml
#######################################
RUN = config['run']

BD = config['bd']

PREFIX = config['prefix']

FQDIR = f'/nexusb/Novaseq/{RUN}/Unaligned/'

BASEDIR = BD + RUN
SAMPLES = [sample for sample in os.listdir(BASEDIR) if sample.startswith(PREFIX)]
# explanation: in BASEDIR I have multiple subdirectories. The names of the subdirectories starting with PREFIX will be the name of the elements I want to have in the list SAMPLES, which eventually shall be my {sample} wildcard

#############
# RULES
#############
rule all:
    input:
        expand("aligned/{sample}.bam", sample = SAMPLES)


rule bwa_map:
    input:
        REF = config['ref'],
        R1 = FQDIR + "{sample}_S{s}_R1_001.fastq.gz",
        R2 = FQDIR + "{sample}_S{s}_R2_001.fastq.gz"
    output:
        "aligned/{sample}.bam"
    shell:
        "bwa mem {input.REF} {input.R1} {input.R2}| samtools view -Sb - > {output}"

但我得到了：

Building DAG of jobs...
WildcardError in line 55 of /nexusb/nipt/200311_A00154_0454_AHHHKMDRXX/testMetrics/snakemake/Snakefile:
Wildcards in input files cannot be determined from output files:
's'

拨打snakemake -np时

我相信我的错误在于输入指令中R1 和R2 的定义。我觉得它令人费解，因为根据official documentationsnakemake 应该将任何通配符解释为正则表达式.+。但是对于示例NIPT-PearlPPlasma-05-PPx 并没有这样做，其R1 和R2 应该分别为NIPT-PearlPPlasma-05-PPx_S5_R1_001.fastq.gz 和NIPT-PearlPPlasma-05-PPx_S5_R2_001.fastq.gz。

Answer 1

再看一下snakemake教程，了解如何从输出中推断输入，无论如何我认为问题在于这段代码：

output:
    expand("aligned/{sample}.bam", sample = SAMPLES)

需要改成

output:
    "aligned/{sample}.bam"

你没有工作，因为之前expand("aligned/{sample}.bam", sample = SAMPLES) 基本上变成了这样的列表["aligned/sample0.bam","aligned/sample1.bam"]。当您删除扩展时，您只需给出输出应该是什么样子的“描述”，因此snakemake 可以推断通配符和输入。

---

编辑：

由于我没有实际文件，因此很难对其进行测试，但是您应该这样做。如果存在多个 S-thingies，将无法正常工作。

def get_reads(wildcards):
    R1 = FQDIR + f"{wildcards.sample}_S{{s}}_R1_001.fastq.gz"
    R2 = FQDIR + f"{wildcards.sample}_S{{s}}_R2_001.fastq.gz"
    globbed = glob_wildcards(R1)
    R1, R2 = expand([R1, R2], s=globbed.s)
    return {"R1": R1, "R2": R2}


rule bwa_map:
    input:
        unpack(get_reads),
        REF = config['ref']
    output:
        "aligned/{sample}.bam"
    shell:
        "bwa mem {input.REF} {input.R1} {input.R2}| samtools view -Sb - > {output}"

Answer 2

问题出在这里：

rule bwa_map:
    input:
        REF = config['ref'],
        R1 = FQDIR + "{sample}_S{s}_R1_001.fastq.gz",
        R2 = FQDIR + "{sample}_S{s}_R2_001.fastq.gz"
    output:
        "aligned/{sample}.bam"

您的output 明确定义了{sample} 是通配符的模式。当 Snakemake 构建 DAG 并发现任何其他规则需要与此模式匹配的文件时，它会将具体值设置为 wildcard.sample。此时应定义所有输入，但您要引入一层间接性：未定义的通配符{s}。

{s} 的值应从output 中明确推断出来。如果你能在设计时做到这一点，用具体值替换它，否则你可以使用 Snakemake 的checkpoint 功能。