Snakemake：将多个输入用于具有多个子组的一个输出的规则

Question

我有一个工作流水线用于下载、比对和对 public 测序数据执行变异检出。问题是它目前只能在每个样本的基础上工作（即样本作为每个单独的测序实验）。如果我想对一组实验（例如样本的生物学 and/or 技术复制）执行变体调用，则它不起作用。我试图解决它，但我无法让它工作。

这里是对齐规则的简化：

rule alignment:
    input:
        rules.download.output.fastq
    output:
        '{group}/alignment/{sample}.bam'
    shell:
        "bash scripts/02_alignment.sh {wildcards.group} {wildcards.samples}"

变异调用也一样：

rule variant_calling:
    input:
        rules.alignment.output
    output:
        '{group}/variants/{sample}.vcf.gz'
    shell:
        "bash scripts/03_variant_calling.sh {wildcards.sample} {wildcards.group}"

这很好用，因为每个对齐的 .bam 文件都会生成一个 .vcf 文件。但我想做的是从任意数量的 .bam 文件生成一个 .vcf 文件。我有一个 pandas 数据框，其中包含所有 sample 名称及其对应的 group。我基本上想将第二条规则的 output 更改为 '{group}/variants/{group}.vcf'，但我所做的一切都以某种方式失败了。

我的想法是为规则提供所有按组对齐的 .bam 文件作为输入，然后只给脚本运行它们所在的目录。问题是我找不到以这种每组方式提供输入的方法：要么我给出每个样本（作为工作管道），要么我为每个组变量调用提供所有 .bam 文件，无论他们实际上属于哪个组。我不能只使用通配符，因为最后的输出中不存在 {sample} 通配符。我也尝试过使用函数作为输入，但这会导致与上述相同的问题。

问题的关键似乎是分组的层次：如果我想对整个数据集中所有对齐的 .bam 文件执行变体调用，那可能会很好，给出我的以上提到的问题。问题来自整个数据集的子组：

  sample1      sample2             sample1      sample2      sample3
     |            |                   |            |            |
     |            |                   |            |            |
     --------------                   ---------------------------
            |                                      |
            |                                      |
          group1                                 group2

有什么解决办法吗？

Answer 1

您必须使用某种结构来将样本分组：

GROUPS = {
    "group1":["sample1","sample2"],
    "group2":["sample1","sample2","sample3"]
}

然后是这样的：

rule all:
    input:
         expand("{group}/variants/{group}.vcf.gz", group=list(GROUPS.keys()))

rule alignment:
    input:
        rules.download.output.fastq
    output:
        '{group}/alignment/{sample}.bam'
    shell:
        "bash scripts/02_alignment.sh {wildcards.group} {wildcards.samples}"

rule variant_calling:
    input:
        lambda wildcards: expand("{group}/alignment/{sample}.bam", group=wildcards.group, sample=GROUPS[wildcards.group]
    output:
        '{group}/variants/{group}.vcf.gz'
    shell:
        "bash scripts/03_variant_calling.sh {input} {output}"

当然还有一些你没有展示的遗漏规则，但我想你会明白的！

来自规则 variant_calling 的 shell 命令可能难以处理，但您始终可以将目录定义为参数，例如：

params: groupAlignDir = "{group}/alignment"

并在 shell 中使用它：

"bash scripts/03_variant_calling.sh {params.groupAlignDir} {output}"

并获取脚本目录中的所有 bam 文件 "variant_calling.sh"

Snakemake：将多个输入用于具有多个子组的一个输出的规则

Snakemake: rule for using many inputs for one output with multiple sub-groups

python

bioinformatics

snakemake

vcf-variant-call-format