Snakemake带有自定义脚本，可通过contig分割床图

Question

我对生物信息学和Snakemake都很陌生，但我正在尝试为Tn-seq数据分析建立一个自动化管道

我编写了一个脚本，读取.bedgraph文件并为每个挫伤输出单独的文件，因为我想分别分析每个挫伤。我已经编写了输出一个文件的代码，该文件的基本名称为输入文件+contig名称：

input_handle = FILE
path = PATH

import csv
import re

contigs = {}

with open(input_handle) as data:
    data_reader = csv.reader(data, delimiter='\t')
    contigs = {row[0] for row in data_reader}

for c in contigs:
    with open(input_handle) as data:
        data_reader = csv.reader(data, delimiter='\t')
        out_file = path + re.search(r".+\/(.+)(?=.bedgraph)", input_handle).group(1) + "-" + c + ".bedgraph"
        f_out = open(out_file, 'w')
        for row in data_reader:
            if row[0] == c:
                f_out.write("\t".join(row)+"\n")

我正在努力找出如何将其适当地纳入Snakemake中。我还不清楚如何在脚本中合并该输出

编辑：我之前说过我认为我应该使用dynamic，但在仔细查看之后，似乎不赞成使用dynamic

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
---->   "split_beds/{sample}/?????"
    script:
        "scripts/split_bed.py"

电话：

input_bedgraph = snakemake.input[0]

import csv
import re

contigs = {}

with open(input_bedgraph) as data:
    data_reader = csv.reader(data, delimiter='\t')
    contigs = {row[0] for row in data_reader}

for c in contigs:
    with open(input_bedgraph) as data:
        data_reader = csv.reader(data, delimiter='\t')
---->   out_file = snakemake.output[0]
        f_out = open(out_file, 'w')
        for row in data_reader:
            if row[0] == c:
                f_out.write("\t".join(row)+"\n")

如果有人能为我指出正确的方向，我将不胜感激！本教程非常适合入门，在此之前，我有很多规则可以很好地使用，但现在我有点迷失了方向

EDIT2:

我已经证实，简单地使用下面的“全部”规则来扩展split_beds文件夹，就可以生成一个文件，每个文件都有一个contig，所以脚本可以很好地工作。。。只需要能够做多个输出到不同的文件夹

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        "split_beds/{sample}.bedgraph"
    script:
        "scripts/split_bed.py"

Answer 1

在得到一些知道他们在做什么的人的帮助后，问题变得更容易解决

我最终使用了检查点（谢谢@Luigi）。我相信使用awk one liner而不是脚本可能会简化问题，但使用下面的代码以及修改的all规则最终会得到我需要的所有文件

checkpoint split_bed:
    input:
        "bedgraphs/{sample}.bedgraph"
    output:
        directory("split_beds/{sample}/")
    shell:
        "mkdir split_beds/{wildcards.sample}; awk \'{{print $0 > \"split_beds/{wildcards.sample}/{wildcards.sample}_\"$1\".bedgraph\"}}\' {input}"

def aggregate_input(wildcards):
    checkpoint_output = checkpoints.split_bed.get(**wildcards).output[0]
    return expand("split_beds/{sample}/{contig}.bedgraph",
           sample=wildcards.sample,
           contig=glob_wildcards(os.zpath.join(checkpoint_output, "{contig}.bedgraph")).contig)

rule aggregate:
    input:
        aggregate_input
    output:
        "split_beds/{sample}/{contig}"

“所有”规则是：

rule all:
    input:
        expand("split_beds/{sample}/", sample=config["samples"])

Answer 2

您是否需要在蛇形管道下游使用contigs？如果是这样的话，dynamic()无疑是一条路要走This related question是dynamic()的一个非常容易访问的介绍。它可能看起来像这样：

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        dynamic("split_beds/{sample}/{contig}.fasta")
    script:
        "scripts/split_bed.py"

rule process_contigs:
    input:
        dynamic("split_beds/{sample}/{contig}.fasta")
    ...

如果我的理解是正确的，snakemake将其理解为“process_contigs需要一些名称与模式'split_beds/{sample}/{contig}.fasta'匹配的文件，这些文件将由split_bed生成”

或者，如果您只想散开分割的contigs，而不需要在snakemake中单独处理这些文件，那么一种不太优雅的方法是在处理完contigs后只需触摸一个文件，然后将其用作snakemake规则的输出（使用pathlib）

您的样本处理代码：

input_handle = FILE
path = PATH

...

from pathlib import Path

...

for c in contigs:

    ...

Path('split_beds/status/{sample}_completed.txt').touch()

你的规则：

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        "split_beds/status/{sample}_completed.txt"
    script:
        "scripts/split_bed.py"

如果分割的重叠是管道的终点，并且您不想麻烦制定另一个规则以将动态重叠文件作为输入，那么我只会真正执行第二个操作

编辑：

我刚刚想起checkpoints的存在。我认为这些被替换的dynamic()检查点基本上允许snakemake在工作流的某些点重新评估DAG，我认为这是您所需要的，因为您将为每个样本生成数量可变的单个重叠文件。因此，DAG不能在运行开始时完成（我认为这就是设计检查点要解决的问题）。我自己没有使用过它们，所以我不会试图通过伪代码编写一个示例来误导您，但希望这能为您指明正确的方向