Любопытный случай создания змеи

У меня такая же цель, как и в Snakemake: неизвестные выходные / входные файлы однако после расщепления по хромосомам, как уже указывалось, я заранее знаю, что у меня есть например, 5 хромосом в моем sample.bam файле. На примере игрушки:

$ 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 касается их пустыми версиями, стирая то, что было помещено туда. Думаю, я удивлен таким поведением и задаюсь вопросом, есть ли для этого веская причина. Обратите внимание, что это поведение не ограничивается touch() snakemake, поскольку, если я заменю touch('sample.REF_{chromosome}.bam') на просто 'sample.REF_{chromosome}.bam', а затем добавлю shell :touch {output} `, я получу тот же результат. Теперь, конечно, я нашел вполне приемлемый обходной путь:

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
person Murray Patterson    schedule 01.03.2019    source источник
comment
Пожалуйста, найдите время, чтобы отформатировать свой вопрос.   -  person matt    schedule 01.03.2019

Ответы (1)


arrow_upward
0
arrow_downward

Я думаю, что ваш пример несколько надуманный по нескольким причинам. Правило 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')
person dariober    schedule 01.03.2019