Snakemake，使用通配符时如何更改输出文件名答案

【问题标题】：Snakemake, how to change output filename when using wildcardsSnakemake，使用通配符时如何更改输出文件名
【发布时间】：2020-03-18 23:11:24
【问题描述】：

我想我有一个简单的问题，但我不知道如何解决它。

我的输入文件夹包含如下文件：

AAAAA_S1_R1_001.fastq
AAAAA_S1_R2_001.fastq
BBBBB_S2_R1_001.fastq
BBBBB_S2_R2_001.fastq

我的蛇形代码：

import glob

samples = [os.path.basename(x) for x in sorted(glob.glob("input/*.fastq"))]
name = []
for x in samples:
    if "_R1_" in x:
        name.append(x.split("_R1_")[0])
NAME = name

rule all:
    input:
        expand("output/{sp}_mapped.bam", sp=NAME),

rule bwa:
    input:
        R1 = "input/{sample}_R1_001.fastq",
        R2 = "input/{sample}_R2_001.fastq"
    output:
        mapped = "output/{sample}_mapped.bam"
    params:
        ref = "refs/AF086833.fa"
    run:
        shell("bwa mem {params.ref} {input.R1} {input.R2} | samtools sort > {output.mapped}")

输出文件名是：

AAAAA_S1_mapped.bam
BBBBB_S2_mapped.bam

我希望输出文件是：

AAAAA_mapped.bam
BBBBB_mapped.bam

如何在 bwa 规则之前或之后更改输出名称或重命名文件。

【问题讨论】：

在 python 交互式解释器中尝试以下操作："AAAAA_S1_R1_001.fastq".split("_R1_")。您会看到“_S1”包含在结果列表的第一个单词中。由于这是您构建NAME 示例名称列表的方式，因此那些“_S1”部分将保留在您的sample 通配符中。

标签： bioinformatics snakemake

【解决方案1】：

试试这个：

import pathlib

indir = pathlib.Path("input")
paths = indir.glob("*_S?_R?_001.fastq")
samples = set([x.stem.split("_")[0] for x in paths])

rule all:
    input:
        expand("output/{sample}_mapped.bam", sample=samples)


def find_fastqs(wildcards):
    fastqs = [str(x) for x in indir.glob(f"{wildcards.sample}_*.fastq")]
    return sorted(fastqs)


rule bwa:
    input:
        fastqs = find_fastqs
    output:
        mapped = "output/{sample}_mapped.bam"
    params:
        ref = "refs/AF086833.fa"
    shell:
        "bwa mem {params.ref} {input.fastqs} | samtools sort > {output.mapped}"

使用输入函数查找rule bwa 的正确样本。可能有更优雅的解决方案，但我现在看不到。不过，我认为这应该可行。

（已编辑以反映 OP 的编辑。）

【讨论】：

非常感谢！我只需要改变我的问题。它并不总是“S1”。我更新了我的问题。
太棒了，现在我需要弄清楚所有代码的含义 =) 谢谢！
@BioFrank 该解决方案构建了一个“正确”的样本名称列表：通过拆分"_"，它排除了"_S*" 部分。然后，它使用一个函数作为输入。这使得以比仅在字符串中替换通配符更复杂的方式确定输入文件名成为可能。在这种情况下，它使用glob 来查找所有以所需通配符开头的fastq 文件。如果您的示例名称包括"_"，或者如果您在input 文件夹中有以示例名称和"_" 开头并以".fastq" 结尾的其他文件，则此解决方案可能会失败。

【解决方案2】：

不幸的是，我也遇到了具有以下逻辑的文件名问题：{batch}/{seq_run}_{index}_{flowcell}_{lane}_{read_orientation}.fastq.gz。

我认为核心问题是没有一个通配符是唯一的。此外，并非所有通配符的所有值都可以组合； seq_run1 在 lane1 上运行，而不是在 lane2 上。因此，expand() 不起作用。

在多次尝试 Snakemake（见下文）后，我的解决方案是使用 mv / sed / rename 标准化输入。删除{batch}、{flowcell} 和{lane} 后，可以使用{sample}，这是{seq_run} 和{index} 的独特组合。

没有有什么作用（但在相同情况下值得为其他人尝试）：

Adding the zip argument to expand()
使用以下语法重命名output：output: "_".join(re.split("[/_]", "{full_filename}")[1,2]+".fastq.gz"

【讨论】：