Объединение нескольких файлов vcf с помощью snakemake - PullRequest
Новая
       66

Объединение нескольких файлов vcf с помощью snakemake

0 голосов
/ 29 мая 2020

Я пытаюсь объединить несколько файлов vcf по хромосомам с помощью snakemake. Мои файлы такие, и, как видите, имеют разные координаты. Как лучше всего объединить все chr1A и все chr1B? 

chr1A:0-2096.filtered.vcf
chr1A:2096-7896.filtered.vcf
chr1B:0-3456.filtered.vcf
chr1B:3456-8796.filtered.vcf

Мой псевдокод: 

chromosomes=["chr1A","chr1B"]

rule all:
    input: 
        expand("{sample}.vcf", sample=chromosomes)

rule merge:
    input:
        I1="path/to/file/{sample}.xxx.filtered.vcf",
        I2="path/to/file/{sample}.xxx.filtered.vcf",
    output:
        outf ="{sample}.vcf"
    shell:
        """
        java -jar picard.jar GatherVcfs I={input.I1} I={input.I2} O={output.outf}
        """

EDIT: 

workdir: "/media/prova/Maxtor2/vcf2/merged/"

import subprocess

d = {"chr1A": ["chr1A:0-2096.flanking.view.filtered.vcf", "chr1A:2096-7896.flanking.view.filtered.vcf"],
     "chr1B": ["chr1B:0-3456.flanking.view.filtered.vcf", "chr1B:3456-8796.flanking.view.filtered.vcf"]}


rule all:
    input:
        expand("{sample}.vcf", sample=d)

def f(w):
    return d.get(w.chromosome, "")


rule merge:
    input:
        f
    output:
        outf ="{chromosome}.vcf"
    params:
        lambda w: "I=" + " I=".join(d[w.chromosome])
    shell:
        "java -jar /home/Documents/Tools/picard.jar GatherVcfs {params[0]} O={output.outf}"

1 Ответ

1 голос
/ 29 мая 2020

Мне удалось воспроизвести вашу ошибку. При ограничении подстановочных знаков он работает: 

d = {"chr1A": ["chr1A:0-2096.flanking.view.filtered.vcf", "chr1A:2096-7896.flanking.view.filtered.vcf"],
     "chr1B": ["chr1B:0-3456.flanking.view.filtered.vcf", "chr1B:3456-8796.flanking.view.filtered.vcf"]}

chromosomes = list(d)

rule all:
    input:
        expand("{sample}.vcf", sample=chromosomes)

# these tell Snakemake exactly what values the wildcards may take
# we use "|" to create the regex chr1A|chr1B
wildcard_constraints:
    chromosome = "|".join(chromosomes)

rule merge:
    input:
        # a lambda is an unnamed function
        # the first argument is the wildcards
        # we merely use it to look up the appropriate files in the dict d
        lambda w: d[w.chromosome]
    output:
        outf = "{chromosome}.vcf"
    params:
        # here we create the string 
        # "I=chr1A:0-2096.flanking.view.filtered.vcf I=chr1A:2096-7896.flanking.view.filtered.vcf"
        # for use in our command
        lambda w: "I=" + " I=".join(d[w.chromosome])
    shell:
        "java -jar /home/Documents/Tools/picard.jar GatherVcfs {params[0]} O={output.outf}"

Он должен был работать и без ограничений; это похоже на ошибку в Snakemake.

...