Snakemake-使用检查点从输出目录创建通配符

Question

我正在将一个多 fasta 文件解析为单个 fasta 文件，我想为每个文件创建通配符，因为下一个规则需要为每个文件并行化。我的问题是我无法从生成的 fasta 文件创建通配符，因为输出会根据我拥有的多 fasta 文件动态变化。这是我的代码：


rule all:
    input:
        final = "kmers/{sample}/results.fasta",
        merge = "merge.fasta",
        arc = "All_fasta/"
       
checkpoint pyt:
    input:
        mm = "fasta.fasta"
    output:
        arc = directory("All_fasta/")
    script:
        "pyt.py"
        
def some_func(wildcards):
    checkpoint_in = checkpoints.pyt.get(**wildcards).output[0]
    # I don't know what to do here
    return "{sample}"


rule not_working:
    input:
        new_input = "{sample}"
    output:
        final = "kmers/{sample}/results.fasta"
    shell:
        "somecmd {input} > {output}"

rule merge:
    input:
        final = expand("kmers/{sample}/results.fasta",sample=sample)
    output:
        merge = "merge.fasta"
    shell:
        "cat {input} > {ouput}"

简而言之，我想为我的 multi-fasta 输入文件中的每个 fasta 条目创建一个通配符。

谢谢！

Answer 1

我想这就是你想要的...

输入文件fasta.fasta是：

cat fasta.fasta 
>a
AAAA
>b
CCCC
>c
TTTT

sequence_names里面生成的aggregate_input是来自原始fasta文件的序列名列表。

import glob

rule all:
    input:
        'merge.fasta',

checkpoint pyt:
    input:
        "fasta.fasta"
    output:
        all_fasta= directory("All_fasta")
    shell:
        # Split multifasta into single files. I assume this is what pyt.py
        # does. This is just a dummy script
        r"""
        mkdir {output}

        grep -A 1 '>a' {input} > {output}/a.fasta
        grep -A 1 '>b' {input} > {output}/b.fasta
        grep -A 1 '>c' {input} > {output}/c.fasta
        """

rule process_sequences:
    input:
        new_input= 'All_fasta/{sequence_name}.fasta',
    output:
        final= 'kmers/{sequence_name}/results.fasta',
    shell:
        "cp {input} {output}"

def aggregate_input(wildcards):
    checkpoint_output = checkpoints.pyt.get().output['all_fasta']
    fasta_files = sorted(glob.glob(os.path.join(checkpoint_output, '*.fasta')))
    sequence_names = [re.sub('\.fasta$', '', os.path.basename(x)) for x in fasta_files]
    return expand('kmers/{sequence_name}/results.fasta', sequence_name= sequence_names)

rule merge:
    input:
        final = aggregate_input,
    output:
        merge = "merge.fasta"
    shell:
        "cat {input} > {output}"

我避免使用 glob_wildcards 因为我觉得它太神秘了。我更喜欢使用 glob.glob 收集输出文件并使用 python 内置函数解析它们。

Snakemake-使用检查点从输出目录创建通配符

Snakemake-create wildcards from output directory using checkpoints

python

bioinformatics

snakemake