Сценарий snakemake имеет доступ к stdin / stdout для обработки потока - PullRequest
0 голосов
/ 15 января 2019

Для рабочего процесса Snakemake мне нужно манипулировать тегами во многих файлах BAM, и я хотел бы обработать их, передав их через сценарий ( с использованием скрипта Snakemake: директива ).Конкретным образом я делаю это с обработкой потока pysam .

infile = pysam.AlignmentFile("-", "rb")
outfile = pysam.AlignmentFile("-", "wb", template=infile)

for s in infile:
    (flowcell, lane) = s.query_name.split(':')[0:2]
    rg_id = ".".join([flowcell, lane])
    s.set_tag('RG',rg_id,'Z')
    outfile.write(s)

Этот скрипт работает хорошо автономно, но я не смог понять, как его интегрировать с помощьюдиректива о змейке script.Я предпочитаю этот способ минимизировать использование ввода-вывода и оперативной памяти.

Редактировать: прибегать к прямой загрузке для исправления тега RG.

# parameters passed from snakemake
bam_file = snakemake.input[0]
fixed_bam_file = snakemake.output[0]

bamfile  = pysam.AlignmentFile(bam_file, "rb")
fixed_bamfile = pysam.AlignmentFile(fixed_bam_file, "wb", template = bamfile)

for i, read in enumerate(bamfile.fetch()):
    (flowcell, lane) = read.query_name.split(':')[0:2]
    rg_id = ".".join([flowcell, lane])
    read.set_tag('RG',rg_id,'Z')
    fixed_bamfile.write(read)
    if not (i % 100000):
        print("Updated the read group for {} reads to {}".format(i, rg_id))

bamfile.close()
fixed_bamfile.close()

РЕДАКТИРОВАТЬ: установлены директивы Snakemakes run: и shell:каталог workdir:, в то время как директива script: работает относительно каталога, в котором был запущен Snakefile (сохраняя все аккуратно и аккуратно).Отсюда проблема размещения потокового процессора под script:.

Ответы [ 2 ]

0 голосов
/ 17 января 2019

@ Кришан, кажется, вы уже нашли решение, и если это так, может быть, хорошо, чтобы опубликовать его в качестве ответа.

В качестве альтернативы, вы можете использовать объект {workflow}, чтобы получить каталог Snakefile и оттуда построить путь к вашему скрипту Python. Если ваша структура каталогов:

./
├── Snakefile
├── data
│   └── sample.bam
└── scripts
    └── edit_bam.py

Змеиный файл может выглядеть так:

rule all:
    input:
        'test.tmp',

rule one:
    input:
        'sample.bam',
    output:
        'test.tmp',
    shell:
        r"""
        cat {input} \
        | {workflow.basedir}/scripts/edit_bam.py >  {output}
        """

Выполнено с snakemake -d data ...

Кажется, объект workflow не задокументирован, но проверьте эту ветку Есть ли способ получить полный путь к Snake-файлу внутри Snake-файла?

0 голосов
/ 16 января 2019

Использование shell вместо script директивы:

rule all:
    input:
        expand('{sample}_edited.bam'), sample=['a', 'b', 'c']


rule somename:
    input:
        '{sample}.bam'
    output:
        '{sample}_edited.bam'
    shell:
        '''
        cat {input} > python edit_bam.py > {output}
        '''
...