3

我有一个工作管道，用于下载、对齐和执行公共测序数据的变体调用。问题是它目前只能在每个样本的基础上工作（即样本作为每个单独的测序实验）。如果我想对一组实验（例如样本的生物和/或技术复制）执行变体调用，它不起作用。我试图解决它，但我无法让它工作。

这是对齐规则的简化：

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}"

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

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

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

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

关于如何解决这个问题的任何想法？

4

1 回答 1

3

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

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"

并在外壳中使用它：

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

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

于 2018-02-27T13:56:07.880 回答