时间戳

首页 > 解决方案 > Snakemake 使用自定义脚本按 contig 拆分床位图

python - Snakemake 使用自定义脚本按 contig 拆分床位图

问题描述

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

我编写了一个脚本，它读取 .bedgraph 文件并为每个 contig 输出单独的文件，因为我想分别分析每个 contig。我已经编写了代码来输出具有输入文件的基本名称 + 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 中。我也很不清楚如何将该输出合并到我的脚本中。

编辑：我之前说过我认为我应该使用动态，但在环顾四周后，它似乎已被弃用。

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")

如果有人能指出我正确的方向，将不胜感激！该教程非常适合入门，并且在此之前我有很多可以很好地工作的规则，但是在这一点上我有点迷茫。

编辑2：

我已经确认，只需将以下内容与扩展 split_beds 文件夹的“全部”规则一起使用，就可以制作一个带有单个 contig 的单个文件，因此脚本工作正常......只需要能够进行多个输出到不同的文件夹...

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

标签： pythonbioinformaticssnakemake

解决方案

在从知道自己在做什么的人那里得到一些帮助之后，问题变得更容易解决了。

我最终确实使用了检查点（感谢@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"])

推荐阅读