一个奇怪的蛇形案例答案

【问题标题】：A curious case of snakemake一个奇怪的蛇形案例
【发布时间】：2019-03-01 14:05:19
【问题描述】：

我的目标与Snakemake: unknown output/input files after splitting by chromosome 中的目标相似，但是，正如我所指出的，我事先知道我的sample.bam 文件中有例如，5 条染色体。用作玩具示例：

$ cat sample.bam 
chromosome 1
chromosome 2
chromosome 3
chromosome 4
chromosome 5

我希望“拆分”这个 bam 文件，然后在生成的染色体上执行一堆 每个染色体 下游作业。我能想到的最简单的解决方案是：

chromosomes = '1 2 3 4 5'.split()

rule master :
    input :
        expand('sample.REF_{chromosome}.bam',
            chromosome = chromosomes)


rule chromosome :
    output :
        touch('sample.REF_{chromosome}.bam')

    input : 'split.done'


rule split_bam :
    output :
        touch('split.done')

    input : 'sample.bam'

    run :
        print('splitting bam..')
        chromosome = 1
        for line in open(input[0]) :
            outfile = 'sample.REF_{}.bam'.format(chromosome)
            print(line, end = '', file = open(outfile, 'w'))
            chromosome += 1

导致空的sample_REF_{chromosome}.bam 文件。我理解为什么会发生这种情况，而且确实 snakemake 甚至会发出警告，例如，

Warning: the following output files of rule chromosome were not present when the DAG was created:
{'sample.REF_3.bam'}
Touching output file sample.REF_3.bam.

也就是说，这些文件一开始就没有在 DAG 中，snakemake 会用空版本来处理它们，从而删除放在那里的内容。我想我对这种行为感到惊讶，并想知道这是否有充分的理由。请注意，此行为不仅限于snakemake 的touch()，因为我是否应该将touch('sample.REF_{chromosome}.bam') 简单地替换为'sample.REF_{chromosome}.bam'，然后使用shell :touch {输出}`，我会得到相同的结果。现在，当然，我找到了一个完全可以接受的解决方法：

chromosomes = '1 2 3 4 5'.split()

rule master :
    input :
        expand('sample.REF_{chromosome}.bam',
            chromosome = chromosomes)


rule chromosome :
    output : 'sample.REF_{chromosome}.bam'

    input : 'split_dir'

    shell : 'mv {input}/{output} {output}'


rule split_bam :
    output :
        temp(directory('split_dir'))

    input : 'sample.bam'

    run :
        print('splitting bam..')
        shell('mkdir {output}')
        chromosome = 1
        for line in open(input[0]) :
            outfile = '{}/sample.REF_{}.bam'.format(output[0], chromosome)
            print(line, end = '', file = open(outfile, 'w'))
            chromosome += 1

但我很惊讶我必须通过这些体操来完成一项看似简单的任务。出于这个原因，我想知道是否有更好的设计，或者我没有问正确的问题。欢迎任何建议/想法。

【问题讨论】：

请花一些时间来格式化您的问题。

标签： snakemake

【解决方案1】：

我认为你的例子有点做作，有几个原因。规则split_bam 已经产生了最终输出sample.REF_{chromosome}.bam。此外，规则master 使用从变量chromosomes 中获取的染色体，而规则split_bam 遍历bam 文件以获取染色体。

我的印象是你想要的可能是这样的：

chromosomes= '1 2 3 4 5'.split()

rule master:
    input:
        expand('sample.REF_{chromosome}.bam',
            chromosome = chromosomes)

rule split_bam :
    input:
        'sample.bam'
    output:
        expand('sample.split.{chromosome}.bam', chromosome= chromosomes)
    run:
        print('splitting bam..')
        for chromosome in chromosomes:
            outfile = 'sample.split.{}.bam'.format(chromosome)
            print(chromosome, end = '', file = open(outfile, 'w'))

rule chromosome:
    input:
        'sample.split.{chromosome}.bam'
    output:
        touch('sample.REF_{chromosome}.bam')

【讨论】：