использовать змеиную пару

Question

Я новичок в использовании snakemake, у меня есть простая проблема при выполнении картирования в snakemake. У меня есть пары _1.fastq.gz & _2.fastq.gz, и я хотел бы сделать отображение конца пары для примерно 10 пар fastq.gz. Поэтому я написал файл змеиного мейкера для этого.

import os
import snakemake.io
import glob

(SAMPLES,READS,) = glob_wildcards("raw/{sample}_{read}.fastq.gz")
READS=["1","2"]
REF="/data/data/reference/refs/ucsc.hg19.fasta.gz"

rule all:
    input: expand("raw/{sample}.bam",sample=SAMPLES)

rule bwa_map:
    input:
        ref=REF,
        r1=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS),
        r2=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS)

    output: "raw/{sample}.bam"

    shell: "bwa mem -M -t 8 {input.ref} {input.r1} {input.r2} | samtools view -Sbh - > {output}"

Ошибка:

RuleException:
CalledProcessError in line 17 of /data/data/Samples/snakemake-example/WGS-test/step2.smk:
Command ' set -euo pipefail;  bwa mem -M -t 8 /data/data/reference/refs/ucsc.hg19.fasta.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz | samtools view -Sbh - > raw/sub2.bam ' returned non-zero exit status 1.
  File "/data/data/Samples/snakemake-example/WGS-test/step2.smk", line 17, in __rule_bwa_map
  File "/root/miniconda3/envs/bioinfo/lib/python3.6/concurrent/futures/thread.py", line 56, in run

Мой желаемый результат будет похож на генерацию всех 10 файлов BAM, таких как

sub1.bam sub2.bam sub3.bam ...

Похоже, все команды fastq помещены в команду. Как я могу разделить их и запустить пару за парой автоматически без использования метода жесткого кода. Пожалуйста, совет.

Answer 1

Первое правило (здесь rule all) указывает файлы, которые вы хотели бы создать во время рабочего процесса snakemake.

Для данного файла, f, в rule all::input, snakemake просмотрит все правила и попытается найти то, которое может создать f (на основе сопоставления с образцом в сегменте output каждого правила) .

Предположим, f = raw/my_sample.bam

Как только snakemake найдет правило, которое может создать f, оно определит все входные файлы, необходимые для создания этого файла.

Итак, snakemake обнаруживает, что f = raw/my_sample.bam может быть создан с помощью rule bwa_map (поскольку f соответствует шаблону raw/<anything>.bam), а затем определяет, какие файлы необходимы для создания f на основе сегмента input .

Snakemake думает: я могу сделать raw/my_sample.bam, если у меня есть файл ref="/data/data/reference/refs/ucsc.hg19.fasta.gz" файлы r1=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS) и файлы r2=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS)

В expand, r1 расширяет sample до каждого значения в ОБРАЗЦАХ и read до каждого значения в ЧИТАНИЯХ. Но у вас есть 10 значений в SAMPLES и 2 значения в READS, поэтому r1 расширяется до 20 различных путей к файлам для каждого выходного файла, который он пытается создать. Он игнорирует подстановочный знак sample, присутствующий в предложении output (поскольку вы перезаписали его в вызове expand).

Необходимо разрешить использование подстановочных знаков, определенных в предложении вывода, до предложения ввода :

import os
import snakemake.io
import glob

(SAMPLES,READS,) = glob_wildcards("raw/{sample}_{read}.fastq.gz")
READS=["1","2"]
REF="/data/data/reference/refs/ucsc.hg19.fasta.gz"

rule all:
    input: expand("raw/{sample}.bam",sample=SAMPLES)

rule bwa_map:
    input:
        ref=REF,
        # determine `r1` based on the {sample} wildcard defined in `output`
        # and the fixed value `1` to indicate the read direction
        r1="raw/{sample}_1.fastq.gz",
        # determine `r2` based on the {sample} wildcard similarly
        r2="raw/{sample}_2.fastq.gz"

    output: "raw/{sample}.bam"

    # better to pass in the threads than to hardcode them in the shell command
    threads: 8

    shell: "bwa mem -M -t {threads} {input.ref} {input.r1} {input.r2} | samtools view -Sbh - > {output}"

Я бы посоветовал вам взглянуть на то, как правило для выравнивания bwa записано в ресурсе обертки snakemake (вы также можете подумать об использовании обертки): https://snakemake -wrappers.readthedocs. -й / о / стабильный / упаковщики / BWA / mem.html

Не по теме: С точки зрения обзора кода я спрашиваю, почему вы выводите выровненные данные в каталог raw? Имеет ли смысл выводить выровненные данные в align или aligned? Вы также, кажется, импортируете из пакетов, которые вы не используете.