Форматирование входных файлов Snakemake в оболочке - PullRequest
1 голос
/ 03 марта 2020

Я использую конвейер snakemake для запуска команды GATK MarkDuplicate с несколькими входными файлами bam из разных групп чтения.

rule mark_duplicates:
    input:
        get_dedup_input
    output:
        bam=temp("bams/{patient}.{sample_type}.markdups.bam"),
        md5=temp("bams/{patient}.{sample_type}.markdups.bam.md5"),
        metrics="qc/gatk/{patient}_{sample_type}_dup_metrics.txt"
    conda:
        "../envs/gatk.yml"
    shell:
        """
        gatk MarkDuplicates -I {input} -O {output.bam} -M {output.metrics} \
            --CREATE_MD5_FILE true --ASSUME_SORT_ORDER "queryname"
        """

get_dedup_input возвращает список файлов ввода bam. MarkDuplicates требует, чтобы каждый входной файл был указан с флагом -I. Если бы у меня был только один файл bam, я мог бы просто написать -I {input}, но это не удалось, потому что это указывает -I file1.bam file2.bam, это должно быть -I file1.bam -I file2.bam. Каков наилучший способ отформатировать ввод, чтобы каждый входной файл указывался как -I [input file]?

Вот два сценария ios ниже, чтобы уточнить, как будут выглядеть команды ввода, вывода и оболочки, если я должны были запустить команду вручную. Для краткости я опустил некоторые несущественные MarkDuplicate флаги:

Одна группа чтения

Inputs: patient101.normal.rg1.bam
Output: patient101.normal.markdups.bam
Shell: 
gatk MarkDuplicates -I patient101.normal.rg1.bam \
-O patient101.normal.markdups.bam \
-M metrics.txt

Две группы чтения

Inputs: patient101.normal.rg1.bam, patient101.normal.rg2.bam
Output: patient101.normal.markdups.bam
Shell: 
gatk MarkDuplicates -I patient101.normal.rg1.bam \
-I patient101.normal.rg2.bam \
-O patient101.normal.markdups.bam \
-M metrics.txt

1 Ответ

1 голос
/ 03 марта 2020

Спасибо за обновление вашего ответа.

Так что, вероятно, лучше всего сделать, это сделать параметр, который станет нужной нам строкой, например:

rule mark_duplicates:
    input:
        get_dedup_input
    output:
        bam=temp("bams/{patient}.{sample_type}.markdups.bam"),
        md5=temp("bams/{patient}.{sample_type}.markdups.bam.md5"),
        metrics="qc/gatk/{patient}_{sample_type}_dup_metrics.txt"
    conda:
        "../envs/gatk.yml"
    params:
        input=lambda wildcards, input: " -I ".join(input)
    shell:
        """
        gatk MarkDuplicates -I {params.input} -O {output.bam} -M {output.metrics} \
            --CREATE_MD5_FILE true --ASSUME_SORT_ORDER "queryname"
        """

Если наш ввод это одиночный бам, params.input будет просто patient101.normal.rg1.bam и -I добавляется вперед как обычно.

Если у нас есть два входных бамса, наша лямбда-функция помещает -I между ними и команда оболочки добавляет -I впереди.

...