snakemake 管理对样本和 indelrealigner

Question

我想将 realigner 进程与 indel realignment 连接起来。 这是规则：

rule gatk_IndelRealigner:
    input:
        tumor="mapped_reads/merged_samples/{tumor}.sorted.dup.reca.bam",
        normal="mapped_reads/merged_samples/{normal}.sorted.dup.reca.bam",
        id="mapped_reads/merged_samples/operation/{tumor}_{normal}.realign.intervals"

    output:
        "mapped_reads/merged_sample/CoClean/{tumor}.sorted.dup.reca.cleaned.bam",
        "mapped_reads/merged_sample/CoClean/{normal}.sorted.dup.reca.cleaned.bam",
    params:
        genome=config['reference']['genome_fasta'],
        mills= config['mills'],
        ph1_indels= config['know_phy'],
    log:
        "mapped_reads/merged_samples/logs/{tumor}.indel_realign_2.log"
    threads: 8
    shell:
        "gatk -T IndelRealigner -R {params.genome}  "
        "-nt {threads} "
        "-I {input.tumor} -I {input.normal}  -known {params.ph1_indels} -known {params.mills} -nWayOut  .cleaned.bam --maxReadsInMemory 500000 --noOriginalAligmentTags --targetIntervals {input.id} >& {log}  "

这是错误：

Not all output files of rule gatk_IndelRealigner contain the same wildcards.

我想我也需要使用 {tumor}_{normal} 但我不能使用。 蛇人：

rule all:
      input:expand("mapped_reads/merged_samples/CoClean/{sample}.sorted.dup.reca.cleaned.bam",sample=config['samples']),
           expand("mapped_reads/merged_samples/operation/{sample[1][tumor]}_{sample[1][normal]}.realign.intervals", sample=read_table(config["conditions"], ",").iterrows())

config.yml

conditions: "conditions.csv"

条件.csv

tumor,normal
411,412

在这里您可以看到代码示例（用于测试目的）给出了相同的错误：

目录

$ tree  prova/
prova/
├── condition.csv
├── config.yaml
├── output
│   ├── ABC.bam
│   ├── pippa.bam
│   ├── Pippo.bam
│   ├── TimBorn.bam
│   ├── TimNorm.bam
│   ├── TimTum.bam
│   └── XYZ.bam
└── Snakefile

这是蛇神

$ cat prova/Snakefile 
from pandas import read_table

configfile: "config.yaml"

rule all:
    input:
         expand("{pathDIR}/{sample[1][tumor]}_{sample[1][normal]}.bam", pathDIR=config["pathDIR"], sample=read_table(config["sampleFILE"], " ").iterrows()),
         expand("CoClean/{sample[1][tumor]}.bam",  sample=read_table(config["sampleFILE"], " ").iterrows()),
         expand("CoClean/{sample[1][normal]}.bam", sample=read_table(config["sampleFILE"], " ").iterrows())

rule gatk_RealignerTargetCreator:
     input:
         "{pathGRTC}/{normal}.bam",
         "{pathGRTC}/{tumor}.bam",
     output:
         "{pathGRTC}/{tumor}_{normal}.bam"
 #    wildcard_constraints:
 #        tumor = '[^_|-|\/][0-9a-zA-Z]*',
 #        normal = '[^_|-|\/][0-9a-zA-Z]*'
     run:
         call('touch ' + str(wildcard.tumor) + '_' + str(wildcard.normal) + '.bam', shell=True)



rule gatk_IndelRealigner:
    input:
        t1="output/{tumor}.bam",
        n1="output/{normal}.bam",

    output:
        "CoClean/{tumor}.sorted.dup.reca.cleaned.bam",
        "CoClean/{normal}.sorted.dup.reca.cleaned.bam",
    log:
        "mapped_reads/merged_samples/logs/{tumor}.indel_realign_2.log"
    threads: 8
    shell:
        "gatk -T IndelRealigner -R {params.genome}  "
        "-nt {threads} -I {input.t1} -I {input.n1}  & {log}  "

条件.csv

$ more condition.csv 
tumor normal
TimTum TimBorn
XYZ ABC
Pippo pippa

感谢您的建议

Answer 1

我不认为您必须在 GATK IndelRealigner 中包含两个输入文件。基于该假设，您可以更改规则以使其对处理的文件的“类型（肿瘤与正常）”无动于衷。 I read the specs here。如果我错了，请停止阅读并纠正我。

rule gatk_IndelRealigner:
    input:
        inputBAM="output/{sampleGATKIR}.bam",

    output:
        "CoClean/{sampleGATKIR}.sorted.dup.reca.cleaned.bam",
    log:
        "mapped_reads/merged_samples/logs/{sampleGATKIR}.indel_realign_2.log"
    params:
        genome="**DONT FORGET TO ADD THIS""
    threads: 8
    shell:
        "gatk -T IndelRealigner -R {params.genome}  "
        "-nt {threads} -I {input.inputBAM}   & {log}  "

通过将规则更改为 bam 类型不可知论（组成词），您可以获得两个优势，并且有一个主要劣势。

优点：

1. 现在我们只有一个通配符

2. 我们可以独立运行每个 .bam 文件的对齐方式，使用专用 CPU 有望使事情变得更快。

缺点：

1. 我们现在可能会将基因组的两个副本放在内存中的某个位置，因为线程现在作为单独的进程运行，不再共享基因组文件的内存。 （在我之前的职位上，硬件可用性通常不是问题，所以我非常倾向于将所有内容拆分）

我认为 GATK 文档将其设置为接受多个“bam”文件的原因是，如果您只是将其用作一次性调用，您希望同时列出所有文件。我们不需要那个，因为我们正在自动化呼叫过程。我们对 1 个电话或 100 个电话无动于衷。