snakemake - замена подстановочных знаков во входной директиве анонимной функцией - PullRequest
Новая
       4

snakemake - замена подстановочных знаков во входной директиве анонимной функцией

0 голосов
/ 30 марта 2020

Я пишу змейку, которая будет запускать конвейер биоинформатики для нескольких входных выборок. Эти входные файлы (два для каждого анализа, один с частичным совпадением строк R1 и второй с частичным совпадением строк R2) начинаются с шаблона и заканчиваются расширением .fastq.gz. В конце концов я хочу выполнить несколько операций, хотя для этого примера я просто хочу выровнять чтения fastq по эталонному геному, используя bwa mem . Таким образом, для этого примера мой входной файл NIPT-N2002394-LL_S19_R1_001.fastq.gz, и я хочу сгенерировать NIPT-N2002394-LL.bam (см. Код ниже, определяющий каталоги, в которых находятся вход и выход).

Мой config.yaml файл выглядит так: 

# Run_ID
run: "200311_A00154_0454_AHHHKMDRXX"

# Base directory: the analysis directory from which I will fetch the samples
bd: "/nexusb/nipt/"


# Define the prefix
# will be used to subset the folders in bd
prefix: "NIPT"

# Reference:
ref: "/nexus/bhinckel/19/ONT_projects/PGD_breakpoint/ref_hg19_local/hg19_chr1-y.fasta"

А ниже мой файл змеи 

import os
import re
#############
# config file
#############
configfile: "config.yaml"


#######################################
# Parsing variables from config.yaml
#######################################
RUN = config['run']

BD = config['bd']

PREFIX = config['prefix']

FQDIR = f'/nexusb/Novaseq/{RUN}/Unaligned/'

BASEDIR = BD + RUN
SAMPLES = [sample for sample in os.listdir(BASEDIR) if sample.startswith(PREFIX)]
# explanation: in BASEDIR I have multiple subdirectories. The names of the subdirectories starting with PREFIX will be the name of the elements I want to have in the list SAMPLES, which eventually shall be my {sample} wildcard

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


rule bwa_map:
    input:
        REF = config['ref'],
        R1 = FQDIR + "{sample}_S{s}_R1_001.fastq.gz",
        R2 = FQDIR + "{sample}_S{s}_R2_001.fastq.gz"
    output:
        "aligned/{sample}.bam"
    shell:
        "bwa mem {input.REF} {input.R1} {input.R2}| samtools view -Sb - > {output}"

Но я получаю: 

Building DAG of jobs...
WildcardError in line 55 of /nexusb/nipt/200311_A00154_0454_AHHHKMDRXX/testMetrics/snakemake/Snakefile:
Wildcards in input files cannot be determined from output files:
's'

При звонке snakemake -np

Я верю своей ошибке лежит в определениях R1 и R2 во входной директиве. Я нахожу это загадочным, потому что согласно официальной документации змеиный мастер должен интерпретировать любой подстановочный знак как регулярное выражение .+. Но это не так для образца NIPT-PearlPPlasma-05-PPx, у которого R1 и R2 должны быть NIPT-PearlPPlasma-05-PPx_S5_R1_001.fastq.gz и NIPT-PearlPPlasma-05-PPx_S5_R2_001.fastq.gz соответственно.

Ответы [ 2 ]

0 голосов
/ 30 марта 2020

Проблема здесь: 

rule bwa_map:
    input:
        REF = config['ref'],
        R1 = FQDIR + "{sample}_S{s}_R1_001.fastq.gz",
        R2 = FQDIR + "{sample}_S{s}_R2_001.fastq.gz"
    output:
        "aligned/{sample}.bam"

Ваш output четко определяет шаблон, где {sample} является подстановочным знаком. Когда Snakemake создает DAG и находит, что для любого другого правила требуется файл, соответствующий этому шаблону, он устанавливает конкретное значение wildcard.sample. В этот момент все входные данные должны быть определены, но вы вводите еще один уровень косвенности: подстановочный знак {s}, который не определен.

Значение {s} должно быть ясно выведено из output. Если вы можете сделать это во время разработки, замените его конкретными значениями, в противном случае вы можете использовать checkpoint функцию Snakemake.

0 голосов
/ 30 марта 2020

Посмотрите еще раз на учебник по змеиному мастерству о том, как вывод выводится из вывода, в любом случае, я думаю, что проблема заключается в следующем фрагменте кода: 

output:
    expand("aligned/{sample}.bam", sample = SAMPLES)

И его необходимо изменить на 

output:
    "aligned/{sample}.bam"

То, что у вас не работало, потому что раньше expand("aligned/{sample}.bam", sample = SAMPLES) в основном становилось списком, подобным этому ["aligned/sample0.bam","aligned/sample1.bam"]. Когда вы удаляете раскрытие, вы даете только «описание» того, как должен выглядеть вывод, и, следовательно, snakemake может вывести подстановочные знаки и ввод.

edit:

Это сложно чтобы проверить это, так как у меня нет реальных файлов, но вы должны сделать что-то вроде этого. Не будет работать, если существует несколько S -го вещей. 

def get_reads(wildcards):
    R1 = FQDIR + f"{wildcards.sample}_S{{s}}_R1_001.fastq.gz"
    R2 = FQDIR + f"{wildcards.sample}_S{{s}}_R2_001.fastq.gz"
    globbed = glob_wildcards(R1)
    R1, R2 = expand([R1, R2], s=globbed.s)
    return {"R1": R1, "R2": R2}


rule bwa_map:
    input:
        unpack(get_reads),
        REF = config['ref']
    output:
        "aligned/{sample}.bam"
    shell:
        "bwa mem {input.REF} {input.R1} {input.R2}| samtools view -Sb - > {output}"
...