Подпроцесс: не удается неявно преобразовать объект _io.BufferedReader в str

Я работаю над написанием скрипта, который представляет собой комбинацию кода snakemake и python для автоматизации большого количества файлов, которые поставляются в паре. Точнее, я работаю над согласованием чтения с BWA MEM с парными конечными чтениями (http://bio-bwa.sourceforge.net/bwa.shtml). В первой части сценария я перебирал список имен в моем файле (которые являются файлами fastq bunzip), а затем отсортировал их соответствующим образом в списке. Вот несколько файлов:

['NG-8653_ 1A _lib95899_4332_7_ 1', 'NG-8653_ 1A _lib95899_4332_7_ 2', 'NG- 8653_ 1B _lib95900_4332_7_ 1 ',' NG-8653_ 1B _lib95900_4332_7_ 2 ',' NG-8653_ 1N _lib95898_4332_7_ 1 ',' NG-8653_ 1N _lib95898_4332_7_ 2 ']

Как видите, чтения отсортированы по два (1A _... 1 и 1A ..._ 2 и т. Д.). Теперь, используя подпроцесс, я хочу выровнять их, распаковав их с помощью bunzip2, а затем передав их через стандартный ввод в bwa mem. Команда bwa mem преобразует файлы формата fastq в файлы .sam, а затем я должен использовать samtools для их преобразования в формат .bam. Вот сценарий:

import re, os, subprocess, bz2

WDIR = "/home/alaa/Documents/snakemake"
workdir: WDIR
SAMPLESDIR = "/home/alaa/Documents/snakemake/fastq/"
REF = "/home/alaa/Documents/inputs/reference/hg19_ref_genome.fa"

FILE_FASTQ = glob_wildcards("fastq/{samples}.fastq.bz2")
LIST_FILE_SAMPLES = []

for x in FILE_FASTQ[0]:
    LIST_FILE_SAMPLES.append(x)

LIST_FILE_SAMPLES = sorted(LIST_FILE_SAMPLES)
print(LIST_FILE_SAMPLES)

rule fastq_to_bam:
    run:
        for x in range(0, len(LIST_FILE_SAMPLES), 2):
            # get the name of the sample (1A, 1B ...)
            samp = ""
            samp += LIST_FILE_SAMPLES[x].split("_")[1]

            # get the corresponding read (1 or 2)
            r1 = SAMPLESDIR + LIST_FILE_SAMPLES[x] + ".fastq.bz2"
            r2 = SAMPLESDIR + LIST_FILE_SAMPLES[x+1] + ".fastq.bz2"

            # gunzipping the files and pipping them
            p1 = subprocess.Popen(['bunzip2', '-kc', r1], stdout=subprocess.PIPE)
            p2 = subprocess.Popen(['bunzip2', '-kc', r2], stdout=subprocess.PIPE)           


            # now write the output file to .bam format after aligning them
            with open("sam/" + samp + ".bam", "w") as stdout:
                fastq2sam = subprocess.Popen(["bwa", "mem", "-T 1", REF, p1.stdout, p2.stdout], stdout=subprocess.PIPE)
                fastq2samOutput = subprocess.Popen(["samtools", "view", "-Sb", "-"], shell = True, stdin=fastq2sam.stdout, stdout=stdout)

Я пытался отладить сценарий, пытаясь построчно. При записи bunzip2 в выходной файл он работал нормально. Теперь, если я попытаюсь передать его по конвейеру, я получаю сообщение об ошибке:

Error in job fastq_to_bam while creating output file .
RuleException:
TypeError in line 39 of /home/alaa/Documents/snakemake/Snakefile:
Can't convert '_io.BufferedReader' object to str implicitly
  File "/home/alaa/Documents/snakemake/Snakefile", line 39, in __rule_fastq_to_bam
  File "/usr/lib/python3.5/subprocess.py", line 947, in __init__
  File "/usr/lib/python3.5/subprocess.py", line 1490, in _execute_child
  File "/usr/lib/python3.5/concurrent/futures/thread.py", line 55, in run
 Exiting because a job execution failed. Look above for error message
 Will exit after finishing currently running jobs.
 Exiting because a job execution failed. Look above for error message

Подскажите, пожалуйста, в чем проблема со скриптом? Я пытаюсь найти проблему с сегодняшнего утра и, кажется, не могу ее понять. Любая помощь горячо приветствуется. Заранее спасибо.

ИЗМЕНИТЬ 1:

Прочитав больше об отзывах от @bli и @Johannes, я зашел так далеко:

import re, os, subprocess, bz2, multiprocessing
from os.path import join
from contextlib import closing

WDIR = "/home/alaa/Documents/snakemake"
workdir: WDIR
SAMPLESDIR = "fastq/"
REF = "/home/alaa/Documents/inputs/reference/hg19_ref_genome.fa"


FILE_FASTQ = glob_wildcards("fastq/{samples, NG-8653_\d+[a-zA-Z]+_.+}")
LIST_FILE_SAMPLES = []

for x in FILE_FASTQ[0]:
    LIST_FILE_SAMPLES.append("_".join(x.split("_")[0:5]))

LIST_FILE_SAMPLES = sorted(LIST_FILE_SAMPLES)
print(LIST_FILE_SAMPLES)


rule final:
    input:
        expand('bam/' + '{sample}.bam', sample = LIST_FILE_SAMPLES)

rule bunzip_fastq:
    input:
        r1 = SAMPLESDIR + '{sample}_1.fastq.bz2',
        r2 = SAMPLESDIR + '{sample}_2.fastq.bz2'
    output:
        o1 = SAMPLESDIR + '{sample}_r1.fastq.gz',
        o2 = SAMPLESDIR + '{sample}_r2.fastq.gz'
    shell:
        """
        bunzip2 -kc < {input.r1} | gzip -c > {output.o1}
        bunzip2 -kc < {input.r2} | gzip -c > {output.o2}
        """

rule fastq_to_bam:
    input:
        r1 = SAMPLESDIR + '{sample}_r1.fastq.gz',
        r2 = SAMPLESDIR + '{sample}_r2.fastq.gz',
        ref = REF
    output:
        'bam/' + '{sample}.bam'
    shell:
        """
        bwa mem {input.ref} {input.r1} {input.r2} | samtools -b > {output}
        """

Большое спасибо за вашу помощь! Думаю, с этого момента я смогу справиться.

С уважением, Алаа


person Zen    schedule 10.04.2017    source источник
comment
Я повторяю замечание Йоханнеса Кестера в его ответе. Вы можете подумать о том, чтобы создать отдельное правило для пакетирования, в котором вы могли бы использовать раздел оболочки вместо того, чтобы запускать что-то вручную с subprocess. Затем вы делаете выходные данные этого правила входными данными вашего правила сопоставления (и удаляете цикл и вместо этого используете подстановочные знаки).   -  person bli    schedule 11.04.2017


Ответы (1)


Ваша проблема здесь:

["bwa", "mem", "-T 1", REF, p1.stdout, p2.stdout]

p1.stdout и p2.stdout относятся к типу BufferedReader, но subprocess.Popen ожидает список строк. Возможно, вы захотите использовать, например, p1.stdout.read().

Однако имейте в виду, что ваш подход не является идиоматическим способом использования Snakemake, на самом деле в настоящее время в скрипте нет ничего, что действительно использовало бы функции Snakemake.

В Snakemake вы бы предпочли иметь правило, которое обрабатывает одиночный сэмпл с помощью bwa mem, принимая fastq в качестве входных данных и сохраняя bam в качестве выходных. См. этот пример в официальном руководстве по Snakemake. Он делает именно то, что вы пытаетесь достичь здесь, но с гораздо менее необходимым шаблоном. Просто позвольте Snakemake сделать эту работу, не пытайтесь реализовать это самостоятельно.

person Johannes Köster    schedule 11.04.2017
comment
Спасибо за ответ Йоханнес. Ваш комментарий полезен. И то, что я опубликовал, было только началом сценария. Я полностью понимаю, как работает Snakemake, потому что у меня есть дополнительные правила ниже по течению! Еще раз спасибо за вашу помощь. Я все еще работаю над сценарием, и я опубликую свой обходной путь, когда закончу с ним. - person Zen; 11.04.2017