Как обрабатывать пространственные пары файлов с помощью Snakemake?
Я немного новичок в Snakemake и сейчас ломаю голову над тем, как обрабатывать определенные пары файлов с помощью конвейера Snakemake (извиняюсь, если вопрос неубедительный).
У меня есть набор файлов fastq из опухолей и соответствующих нормальных образцов. Каждый образец опухоли необходимо будет обрабатывать вместе с соответствующим ему нормальным образцом при вызове вариантов и т. д. Итак, я подготовил файл конфигурации со следующими образцами:
sample_list:
- sample: 1
tumor: AO1_04_RN_1_T_4_S4
control: AO2_07_C007558T1Wa_S37
- sample: 2
tumor: AO2_01_C007589T1FTa_S2
control: AO2_07_C007589T1Wa_S34
- sample: 3
tumor: AO9_09_FM_1_T_7_S13
control: AO2_07_C007558T1Wa_S37
Однако первая проблема, с которой я столкнулся, заключается в том, что эти образцы содержат по 4 файла fastq.gz каждый, которые мне нужно объединить в 2 - например, файлы для образца опухоли AO1_04_RN_1_T_4_S4:
- AO1_04_RN_1_T_4_S4_L001_R1_001.fastq.gz
- AO1_04_RN_1_T_4_S4_L002_R1_001.fastq.gz
- AO1_04_RN_1_T_4_S4_L001_R2_001.fastq.gz
- AO1_04_RN_1_T_4_S4_L002_R2_001.fastq.gz
и мне понадобится первое правило (скажем, rule concat_fastq), чтобы просто объединить для каждого образца файлы R1 и R2 в два отдельных fastq AO1_04_RN_1_T_4_S4_R1.fastq.gz и AO1_04_RN_1_T_4_S4_R2.fastq.gz, так как мне нужно будет предоставить их в качестве входных данных для выравнивания.
На данный момент я обнаружил, что указание входных данных выглядит следующим образом:
r1 = expand("{path}/{sample}_L{num}_R1_001.fastq.gz",
path = config["input_path"],
num = ["001","002"],
sample = [sample["tumor"] for sample in config["sample_list"]] + [sample["control"] for sample in config["sample_list"]]),
r2 = expand("{path}/{sample}_L{num}_R2_001.fastq.gz",
path = config["input_path"],
num = ["001","002"],
sample = [sample["tumor"] for sample in config["sample_list"]] + [sample["control"] for sample in config["sample_list"]])
очевидно, каждый раз объединяет слишком много файлов, а что-то вроде этого, без использования expand:
r1 = ["{path}/{sample}_L001_R1_001.fastq.gz", "{path}/{sample}_L002_R1_001.fastq.gz"],
r2 = ["{path}/{sample}_L001_R2_001.fastq.gz", "{path}/{sample}_L002_R2_001.fastq.gz"]
приводит к следующей ошибке:
Wildcards in input files cannot be determined from output files:
'path'
поскольку я не могу указать, что такое {path}.
Может кто-нибудь дать совет, как обращаться с парными файлами с помощью Snakemake? Я также надеюсь, что мой файл конфигурации структурирован правильно, чтобы затем обрабатывать пары опухоль-норма.
Большое спасибо за вашу помощь
Ответы 1
Использование правила all для настройки примеров подстановочных знаков немного упрощает эту задачу. Затем в concat_fastq мы снова используем расширение, чтобы получить две полосы для каждой пары чтения. Я также могу предложить использовать CSV для указания образцов и метаданных, поскольку CSV легче анализировать (по моему мнению).
configfile: "config.yaml"
TUMORS = [d["tumor"] for d in config["sample_list"]]
CONTROLS = [d["control"] for d in config["sample_list"]]
SAMPLES = TUMORS + CONTROLS
READS_PATH = "data/reads"
rule all:
input:
reads=expand(
"results/concat_fastq/{sample}_R{num}.fastq.gz", sample=SAMPLES, num=[1, 2]
),
rule concat_fastq:
input:
# Use {{ }} to escape sample wildcard in expand.
r1=expand(READS_PATH + "{{sample}}_L00{lane}_R1.fastq.gz", lane=[1, 2]),
r2=expand(READS_PATH + "{{sample}}_L00{lane}_R2.fastq.gz", lane=[1, 2]),
output:
r1 = "results/concat_fastq/{sample}_R1.fastq.gz",
r2 = "results/concat_fastq/{sample}_R2.fastq.gz",
shell:
"""
cat {input.r1} > {output.r1}
cat {input.r2} > {output.r2}
"""
Еще раз спасибо за очень полезный ответ. Я продвинулся в разработке последующих правил, но теперь снова застрял в получении правила для вызова вариантов с помощью Mutect2 для обработки n-го образца из TUMORS (скажем,
{input.tumorbam}), в частности, вместе с n-ным из CONTROLS ({input.controlbam}). Другими словами, мне нужно, чтобы Mutect2 выполнился с помощью-I nth_tumor -I nth_control, но я постоянно получаю (даже используяzip) то, что он пытается выполнить-I 1st_tumor 2nd_tumor (etc) -I 1st_control 2nd_control (etc)и терпит неудачу. Буду очень признателен за дополнительную помощь, если это возможно...!