Змеиный мастер управляет парой семплов

Question

Я хочу связать процесс выравнивания с независимым повторным выравниванием. Это правила:

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.

Я полагаю, мне нужно использовать также {опухоль}_{нормальный}, но я не могу использовать. Snakemake:

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"

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

conditions.csv

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

Спасибо за любое предложение

Answer 1

Я не уверен, что вам нужно включить два входных файла в GATK IndelRealigner. Исходя из этого предположения, вы можете изменить правило, чтобы стать безразличным к "типу (опухоль против нормального)" файла, который он обрабатывает. Я читал спецификации здесь. Пожалуйста, если я ошибаюсь, перестаньте читать и поправьте меня.

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

Изменив правило на независимость от бэма (вымышленное слово), вы получаете два преимущества, и есть один главный недостаток.

Преимущества:

1. Теперь у нас есть только один шаблон

2. Мы можем выполнить выравнивание каждого файла.bam независимо, что, как мы надеемся, с выделенным процессором должно ускорить процесс.

Недостаток:

1. Теперь мы, вероятно, помещаем две копии генома в память где-то, поскольку потоки теперь выполняются как отдельные процессы, больше нет совместного использования памяти файлом генома. (На моей предыдущей должности доступность оборудования обычно не была проблемой, поэтому я сильно склонен к тому, чтобы все разделить)

Причина, по которой я думаю, что в документации GATK есть настройка для приема нескольких файлов 'bam', заключается в том, что, если вы просто используете это как однократный вызов, вы хотите перечислить все файлы одновременно. Нам это не нужно, поскольку мы автоматизируем процесс вызова. Нам безразличен 1 звонок или 100 звонков.