python - 当两个不同的规则路径可以生成给定的输出时,snakemake 能否避免歧义?

标签 python bioinformatics snakemake

初始工作流程

我有一个蛇文件,可以从配对端数据生成一些输出。

在此蛇文件中,我有一条规则“安装”配置文件中存储的给定信息的数据 (get_raw_data)。

然后,我有一个规则,使用该数据生成其余工作流程所依赖的中间文件 (run_tophat)。

以下是这些规则的输入和输出(OPJ 代表 os.path.join):

rule get_raw_data:
    output:
        OPJ(raw_data_dir, "{lib}_1.fastq.gz"),
        OPJ(raw_data_dir, "{lib}_2.fastq.gz"),

(稍后详细介绍此规则的实现)

rule run_tophat:
    input:
        transcriptome = OPJ(annot_dir, "dmel-all-r5.9.gff"),
        fq1 = OPJ(raw_data_dir, "{lib}_1.fastq.gz"),
        fq2 = OPJ(raw_data_dir, "{lib}_2.fastq.gz"),
    output:
        junctions = OPJ(output_dir, "{lib}", "junctions.bed"),
        bam = OPJ(output_dir, "{lib}", "accepted_hits.bam"),

(简化)我的主要规则是这样的:

rule all:
    input:
        expand(OPJ(output_dir, "{lib}", "junctions.bed"), lib=LIBS),

将工作流程扩展到单端数据

我现在必须在单端数据上运行我的工作流程。

我想避免最终输出根据数据是单端还是成对端而具有不同的名称模式。

我可以轻松地制作适用于单端数据的上述两个规则的变体(get_raw_data_single_endrun_tophat_single_end),其输入和输出如下:

rule get_raw_data_single_end:
    output:
        OPJ(raw_data_dir, "{lib}.fastq.gz")

rule run_tophat_single_end:
    input:
        transcriptome = OPJ(annot_dir, "dmel-all-r5.9.gff"),
        fq = OPJ(raw_data_dir, "{lib}.fastq.gz"),
    output:
        junctions = OPJ(output_dir, "{lib}", "junctions.bed"),
        bam = OPJ(output_dir, "{lib}", "accepted_hits.bam"),

如何为snakemake提供足够的信息来选择正确的规则路径?

配置文件包含有关lib通配符是否与单端或成对端数据关联的信息,其方式如下:库名称是lib2raw<中的键lib2raw_single_end 字典(两个字典都是从配置文件中读取的)。

我不希望相同的库名称成为两个字典中的键。因此,从某种意义上来说,我想要执行工作流的单端分支还是双端分支是不明确的。

函数lib2data(使用这些字典)由get_raw_dataget_raw_data_single_end使用来确定运行哪个shell命令来“安装” “数据。

这是该函数的简化版本(实际版本包含一个额外的分支,用于为来自 SRR 标识符的数据生成命令):

def lib2data(wildcards):
    lib = wildcards.lib
    if lib in lib2raw:
        raw = lib2raw[lib]
        link_1 = "ln -s %s %s_1.fastq.gz" % (raw.format(mate="1"), lib)
        link_2 = "ln -s %s %s_2.fastq.gz" % (raw.format(mate="2"), lib)
        return "%s\n%s\n" % (link_1, link_2)
    elif lib in lib2raw_single_end:
        raw = lib2raw_single_end[lib]
        return "ln -s %s %s.fastq.gz\n" % (raw, lib)
    else:
        raise ValueError("Procedure to get raw data for %s unknown." % lib)

除了输出之外,这两个 get_raw_data* 规则是相同的,并按以下方式工作:

params:
    shell_command = lib2data,
shell:
    """
    (
    cd {raw_data_dir}
    {params.shell_command}
    )
    """

给定未在规则输入和输出中编码、而仅在配置文件和函数中编码的信息,snakemake 是否能够确定正确的规则路径?

看来事实并非如此。事实上,我正在尝试测试我的新 Snakefile(添加了 *_single_end 规则),但是在执行 get_raw_data 期间发生 KeyError 规则,而执行该规则的库与单端数据关联

如何实现所需的行为(两个分支工作流程能够使用配置中的信息来选择正确的分支)?

编辑:KeyError是由于lib2data中的错误造成的

使用正确的字典获取与库名称关联的数据后,我最终遇到以下错误:

AmbiguousRuleException:
Rules run_tophat and run_tophat_single_end are ambiguous for the file tophat_junction_discovery_revision_supplement/HWT3/junctions.bed.
Expected input files:
        run_tophat: ./HWT3_1.fastq.gz ./HWT3_2.fastq.gz Annotations/dmel-all-r5.9.gff
        run_tophat_single_end: ./HWT3.fastq.gz Annotations/dmel-all-r5.9.gff

编辑 2:向 get_raw_data* 规则添加输入

阅读后this post on the snakemake mailing list ,我尝试在我的规则中添加一些输入以避免歧义。

def lib2data_input(wildcards):
    lib = wildcards.lib
    if lib in lib2raw:
        raw = lib2raw[lib]
        return [raw.format(mate="1"), raw.format(mate="2")]
    elif lib in lib2raw_single_end:
        raw = lib2raw_single_end[lib]
        return [raw]
    else:
        raise ValueError("Procedure to get raw data for %s unknown." % lib)

rule get_raw_data:
    input:
        lib2data_input
# [same output, params and shell as before]
# [same modification for the single-end case]

这会导致MissingInputException。奇怪的是,据报道丢失的文件确实存在。这个技巧应该有效吗?(无法重现这个,现在结果是:)

AmbiguousRuleException:
Rules run_tophat_single_end and run_tophat are ambiguous for the file tophat_junction_discovery_revision_supplement/HTW2/junctions.bed.
Expected input files:
        run_tophat_single_end: ./HTW2.fastq.gz Annotations/dmel-all-r5.9.gff
        run_tophat: ./HTW2_1.fastq.gz ./HTW2_2.fastq.gz Annotations/dmel-all-r5.9.gff

我指定“数据安装”规则输入的方式显然不足以引导snakemake遵循正确的规则。

最佳答案

我不知道这是否有帮助,但您可以使用函数来定义规则的输入。这样,您就可以使用相同的规则来处理单端或成对端数据,前提是规则的输出相同......

def my_inputs(wildcards):
    data_type = config["data_type"]
    if (data_type == "pe"):
        input = ...
    elif (data_type == "se"):
        input = ...
    return input

rule my_rule:
    input: my_inputs
    ...

关于python - 当两个不同的规则路径可以生成给定的输出时,snakemake 能否避免歧义?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/40510347/

相关文章:

snakemake - 是否可以在蛇形试运行中打印命令而不是规则?

conda - Snakemake 包装器中(生物)conda 版本的最佳实践?

python - Google Docs python gdata 2.0.16 将文件上传到现有集合

python - MySQL VarChar 255 及以上

python - 沿轴连接 numpy 字符串数组?

python - BioPython,如何从 .fasta 转换为 .aln 以进行簇比对?

javascript - 将数学公式转换为 node.js

python - 根据行和列确定字母 - python

r - 使用 R 绘制分类数据

python - Snakefile 的各种迭代都会给出相同的错误