Комбинируйте командные строки оболочки в snakemake

Question

Я хотел бы объединить две командные строки в одну, чтобы избежать промежуточных файлов.

workdir: "/path/to/workdir/"

rule all:
    input: 
        "my.filtered.vcf.gz"

rule bedtools:
    input:
        invcf="/path/to/my.vcf.gz",
        bedgz="/path/to/my.bed.gz"
    output:
        outvcf="my.filtered.vcf.gz"
    shell:
        "/Tools/bedtools2/bin/bedtools intersect -a {input.invcf} -b {input.bedgz} -header -wa |"
        "/Tools/bcftools/bcftools annotate -c CHROM,FROM,TO,GENE -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description="Gene name">') > {output.outvcf}"

Я получаю неверную синтаксическую ошибку. Я был бы признателен, если бы вы могли объяснить, как объединить несколько линий оболочки в snakemake.

Answer 1

Вероятно, вы получили неверный синтаксис из-за " вы используете здесь в своей оболочке: Description="Gene name">. Это закрывает вашу оболочку. Вы можете либо избежать этих кавычек, либо использовать""" синтаксис:

rule bedtools:
    input:
        invcf="/path/to/my.vcf.gz",
        bedgz="/path/to/my.bed.gz"
    output:
        outvcf="my.filtered.vcf.gz"
    shell:
        "/Tools/bedtools2/bin/bedtools intersect -a {input.invcf} -b {input.bedgz} -header -wa |"
        "/Tools/bcftools/bcftools annotate -c CHROM,FROM,TO,GENE -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description=\"Gene name\">') > {output.outvcf}"

или

rule bedtools:
    input:
        invcf="/path/to/my.vcf.gz",
        bedgz="/path/to/my.bed.gz"
    output:
        outvcf="my.filtered.vcf.gz"
    shell:
        """
        /Tools/bedtools2/bin/bedtools intersect -a {input.invcf} -b {input.bedgz} -header -wa | /Tools/bcftools/bcftools annotate -c CHROM,FROM,TO,GENE -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description="Gene name">') > {output.outvcf}
        """

Обратите внимание, что вы можете использовать многострочность с """. Пример без труб:

shell:
    """
    bedtools .... {input} > tempFile 
    bcftools .... tempFile > tempFile2
    whatever .... tempFile2 > {output}
    """

Answer 2

Проблема заключается в том, чтобы избежать двойных кавычек, но нужно добавить немного больше о форматировании и каналах.

Я предпочитаю синтаксис, заключающийся в переносе каждой строки в " так что строки могут быть лучше расположены:

rule bedtools:
    input:
        invcf="/path/to/my.vcf.gz",
        bedgz="/path/to/my.bed.gz"
    output:
        outvcf="my.filtered.vcf.gz"
    shell:
        "/Tools/bedtools2/bin/bedtools "
           "intersect "
           "-a {input.invcf} "
           "-b {input.bedgz} "
           "-header -wa "
        "| /Tools/bcftools/bcftools "
           "annotate "
           "-c CHROM,FROM,TO,GENE "
           "-h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description=\"Gene name\">') "
        "> {output.outvcf}"

Я считаю, что этот аргумент легче увидеть, и его легче изменить, перемещая строки. Но обратите внимание, что конечный пробел каждой строки требуется, и вы должны использовать явную новую строку,\n, если вам нужна отдельная команда. Когда приглашение напечатано, вывод будет красиво оформлен. С""" синтаксис, вы должны экранировать каждую новую строку с помощью \ в конце и пробелы в начале строки сохраняются при печати.

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

rule bedtools_intersect:
    input:
        invcf="/path/to/my.vcf.gz",
        bedgz="/path/to/my.bed.gz"
    output:
        outvcf=pipe("my.intersected.vcf.gz")
    shell:
        "/Tools/bedtools2/bin/bedtools "
           "intersect "
           "-a {input.invcf} "
           "-b {input.bedgz} "
           "-header -wa "
        "> {output.outvcf}"

rule bcftools_annotate:
    input:
        invcf="my.intersected.vcf.gz"
    output:
        outvcf="my.filtered.vcf.gz"
    shell:
        "/Tools/bcftools/bcftools "
           "annotate "
           "-c CHROM,FROM,TO,GENE "
           "-h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description=\"Gene name\">') "
           "{input.invcf} "
        "> {output.outvcf}"

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