在snakemake中组合shell命令行

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中组合多个shell行，我将不胜感激。

Answer 1

您可能会得到一个无效的语法，因为您在shell中使用的"是：Description="Gene name">。这将关闭您的shell。您可以转义这些引号或使用"""语法：

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。打印提示符时，输出的格式很好。使用"""语法时，必须在每个换行符的末尾使用\对其进行转义，并在打印时保留行首的空格。

如果您有很多管道工作要做，请查看pipe flag。您将第一步编写为规则，然后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}"

这样做的好处是，您可以在整个管道中重用每个规则以进行相交或注释，同时避免临时文件。