首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >输入通配符约束Snakemake

输入通配符约束Snakemake
EN

Stack Overflow用户
提问于 2017-10-01 00:58:20
回答 1查看 2K关注 0票数 1

我对snakemake还不熟悉,我正在使用它将两个管道的步骤合并到一个更大的管道中。几个步骤创建同名文件的问题,而且我无法找到限制通配符的方法,因此我在一个无法解决的步骤上得到了一个Missing input files for rule错误。

完整的snakefile很长,可以在这里获得:https://mdd.li/snakefile

相关章节如下(档案中缺少部分):

代码语言:javascript
复制
wildcard_constraints:
    sdir="[^/]+",
    name="[^/]+",
    prefix="[^/]+"

# Several mapping rules here

rule find_intersecting_snps:
    input:
        bam="hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.bam"
    params:
        hornet=os.path.expanduser(config['hornet']),
        snps=config['snps']
    output:
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.remap.fq1.gz",
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.remap.fq2.gz",
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.keep.bam",
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.to.remap.bam",
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.to.remap.num.gz"
    shell:
        dedent(
            """\
            python2 {params.hornet}/find_intersecting_snps.py \
            -p {input.bam} {params.snps}
            """
        )

# Several remapping steps, similar to the first mapping steps, but in a different directory

rule wasp_merge:
    input:
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.keep.bam",
        "hic_mapping/wasp_results/{sdir}_{prefix}_filt_hg19.remap.kept.bam"
    output:
        "hic_mapping/wasp_results/{sdir}_{prefix}_filt_hg19.bwt2pairs.filt.bam"
    params:
        threads=config['job_cores']
    shell:
        dedent(
            """\
            {module}
            module load samtools
            samtools merge --threads {params.threads} {output} {input}
            """
        )

# All future steps use the name style wildcard, like below

rule move_results:
    input:
        "hic_mapping/wasp_results/{name}_filt_hg19.bwt2pairs.filt.bam"
    output:
        "hic_mapping/wasp_results/{name}_filt_hg19.bwt2pairs.bam"
    shell:
        dedent(
            """\
            mv {input} {output}
            """
        )

这个管道本质上是在一个类似于hic_mapping/bowtie_results/bwt2/<subdir>/<file>的目录结构中执行一些映射步骤(其中subdir是三个不同的目录),然后过滤结果,然后在hic_remap/bowtie_results/bwt2/<subdir>/<file>中执行另一个几乎相同的映射步骤,然后将结果合并到一个全新的目录中,然后将子目录折叠到文件名:hic_mapping/wasp_results/<subdir>_<file>中。

我遇到的问题是,如果wasp_merge步骤将子目录名称折叠到文件名中,则find_intersecting_snps步骤会中断。如果我只维护子目录结构,一切都正常。然而,这样做会破坏管道的未来步骤。

我得到的错误是:

代码语言:javascript
复制
MissingInputException in line 243 of /godot/quertermous/PooledHiChip/pipeline/Snakefile:
Missing input files for rule find_intersecting_snps:
hic_mapping/bowtie_results/bwt2/HCASMC5-8_HCASMC-8-CAGAGAGG-TACTCCTT_S8_L006/001_hg19.bwt2pairs.bam

正确的文件是:hic_mapping/bowtie_results/bwt2/HCASMC5-8/HCASMC-8-CAGAGAGG-TACTCCTT_S8_L006_001_hg19.bwt2pairs.bam

但它正在寻找:hic_mapping/bowtie_results/bwt2/HCASMC5-8_HCASMC-8-CAGAGAGG-TACTCCTT_S8_L006/001_hg19.bwt2pairs.bam

它不是在任何地方创建的,也不是由任何规则定义的。我认为wasp_merge步骤:hic_mapping/wasp_results/HCASMC5-8_HCASMC-8-CAGAGAGG-TACTCCTT_S8_L006_001_filt_hg19.bwt2pairs.filt.bam所创建的文件的存在使人感到困惑。

或者可能是下游文件(在创建此错误的目标之后):hic_mapping/wasp_results/HCASMC5-8_HCASMC-8-CAGAGAGG-TACTCCTT_S8_L006_001_filt_hg19.bwt2pairs.bam

但是,我不知道为什么这两个文件都会混淆find_intersecting_snps规则,因为目录结构完全不同。

我觉得我一定是遗漏了一些显而易见的东西,因为这个错误太荒谬了,但我不知道它是什么。

EN

回答 1

Stack Overflow用户

发布于 2017-10-02 01:11:18

问题是目录名和文件名都包含下划线,而在最后的文件名中,我用下划线分隔这两个组件。

通过更改分离字符,或者用从其他地方获取名称的python函数替换规则,我可以解决这个问题。

这样做是可行的:

代码语言:javascript
复制
rule wasp_merge:
    input:
        "hic_mapping/bowtie_results/bwt2/{sdir}/{prefix}_hg19.bwt2pairs.keep.bam",
        "hic_mapping/wasp_results/{sdir}{prefix}_filt_hg19.remap.kept.bam"
    output:
        "hic_mapping/wasp_results/{sdir}{prefix}_filt_hg19.bwt2pairs.filt.bam"
    params:
        threads=config['job_cores']
    shell:
        dedent(
            """\
            {module}
            module load samtools
            samtools merge --threads {params.threads} {output} {input}
            """
        )
票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/46508335

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档