Блог

Я обрабатываю большое количество файлов fastq. У меня есть несколько папок, по одной для каждого образца, каждая из которых содержит несколько копий (repl), например, папка может содержать:

 Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R{r}_{repl}.fastq.gz
{r} = 1 or 2
{repl} = [1...n 1] is a three digit integer with leading zeros
 

в каждой папке разное количество файлов, например, n может быть 13, например, в этом случае будут следующие файлы:

 Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_001.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_002.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_003.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_004.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_005.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_006.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_007.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_008.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_009.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_010.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_011.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_012.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_013.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_014.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_001.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_002.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_003.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_004.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_005.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_006.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_007.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_008.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_009.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_010.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_011.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_012.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_013.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_014.fastq.gz
 

…

Мне нужно передать salmon в моей командной оболочке r1 и r2, где r1 — это список файлов для образца, содержащего ‘R1’ в имени файла, аналогично для R2.

Поскольку я не знаю, какие могут быть реплики в каждой папке, я подумал проверить, существует ли файл, и вернуть список из моей собственной функции repl_expand . Однако он создает список всех файлов во всех образцах, а не ограничивается для каждого образца, поскольку ему передается список expand(["{sample}"], sample=DATASETS) .

Казалось бы, я мог бы просто передать «{sample»} в свою функцию, но это всего лишь строка… замена {sample} фактическим именем образца, по-видимому, выполняется snakemake после вызова моей функции, но перед вызовом, если я использую функцию расширения snakemake.

Чтобы было ясно, я хочу, чтобы salmon запускался три раза, каждый раз, когда r1 является списком файлов ‘R1’ для образца, аналогично для r2.

Поэтому я хотел бы иметь возможность передавать только одно повторяющееся имя образца в мою функцию в каждом из трех вызовов, чтобы создать r1 и r2 для моей командной строки salmon для этого образца.

Вот мой Snakefile, который производит три прогона лосося, но все они одинаковы и содержат данные из всех выборок.

 import os

MAX_REPLICATES = 30
DATASETS = ["Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001",
 "Sample_BG-1058-ME/BG-1058-ME_ACAGTG_L002",
 "Sample_BO-0712-ME/BO-0712-ME_GTGAAA_L002"]
# and many more...

def repl_expand(templates, direction):
    ss = []
    for templ in templates:
        templ = templ   "_R{dir}_{repl:03d}.fastq.gz"
        s = []
        for i in range (1, MAX_REPLICATES):
            name = templ.format(dir=direction, repl=i)
            if os.path.exists(templ.format(dir=direction, repl=i)):
                s.append(name)
        ss.append(s)
    return ss

SALMON = "/path/to/salmon/salmon-1.4.0_linux_x86_64/salmon-latest_linux_x86_64/bin/salmon"

rule all:
    input: expand(["quants/{dataset}.sf"], dataset=DATASETS)


rule salmon_quant:
    input:
        r1 = repl_expand(expand(["{sample}"], sample=DATASETS), 1),
        r2 = repl_expand(expand(["{sample}"], sample=DATASETS), 2),
        index = "path/to/salmon/salmon_sa_index"
    output:
        "quants/{sample}.sf"
    params:
        dir = "quants",
        libtype = 'A',
        threads = 12
    shell:
        "{SALMON} quant 
         -i {input.index} 
         --libtype {params.libtype} 
         --threads {params.threads} 
         --validateMappings 
         --gcBias 
         -o {params.dir} 
         -1 {input.r1} -2 {input.r2}"
 

Я новичок в snakemake, поэтому, возможно, подходил к этому совершенно неправильно и был бы рад, если бы меня поставили на правильный путь…

Спасибо

Первое замечание: вам не нужно предоставлять список expand функции. Эта функция принимает строку шаблона и возвращает список строк как результат подстановки переменных. Итак, правило all должно быть:

 rule all:
    input: expand("quants/{dataset}.sf", dataset=DATASETS)
 

Следующее замечание касается именования подстановочных знаков. Однако Snakemake не волнует, вызываете ли вы один и тот же объект с разными именами в разных правилах, я бы рекомендовал вам быть последовательным и явным, чтобы избежать ошибок, поэтому вывод rule salmon_quant может звучать как:

 rule salmon_quant:
    output:
        "quants/{dataset}.sf"
 

Интересно, что вы используете имя «{sample}» в своем rule salmon_quant в двух разных контекстах, поэтому они означают совершенно разные вещи. {sample} В разделе вывода (тот, который я переименовал) является шаблоном, в то время {sample} как на входе является переменной для expand функции. Вот почему вы наблюдаете эту «аномалию» подстановки.

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

 rule salmon_quant:
    input:
        r1 = lambda wildcards: repl_expand([wildcards.dataset], 1),
        r2 = lambda wildcards: repl_expand([wildcards.dataset], 2),
        index = "path/to/salmon/salmon_sa_index"
    output:
        "quants/{dataset}.sf"
 

На самом деле вам не нужен этот уродливый список из одного элемента, если вы упростите реализацию repl_expand удаления ненужного цикла for templ in templates: .