Snakemake workflow for the vg toolkit.

public 1yr ago 0 bookmarks

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output, operation, topic

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

This repository contains a Snakemake workflow for the vg toolkit . In the current stage it's simply a place to host a usable Snakefile and config.yaml files. With it, you can build a variation graph from a reference FASTA file and variants in a VCF file, map reads to the graph and genotype those variants.

For now, this workflow has been used locally or on a large AWS/Kubernetes instance. We will soon test it on a HPC system and update the examples below.

This workflow was tested on vg 1.21.0 .

Input files
Examples
Overview of the workflow
Overview of the workflow when using the Giraffe mapper
FAQ

Input files

To use the workflow you should first copy the Snakefile and config.yaml files in the working directory.

Graph construction

For construction, two files are necessary in the working directory:

a reference file named {GENOME}.fa
a VCF file called {VARIANTS}.vcf.gz (sorted and indexed)

The {GENOME} and {VARIANTS} labels should match the graph label in the configuration ( config.yaml file or in the command line). The corresponding graph label is {GENOME}-{VARIANTS} . These labels shouldn't contain the character - .

Mapping and genotyping samples

The reads for each sample should be in a folder named with the sample name. In the {sample} folder, FASTQ files should be named {sample}_1.fastq.gz and {sample}_2.fastq.gz . The sample names shouldn't contain the character - .

Examples

The two rules used most in practice are:

construct_all to build the graph and all the necessary indexes for mapping and variant calling.
genotype to map reads and genotype samples.

See an example on a small dataset in the testdata folder .

Construct a variation graph

To build a graph from a hg38.fa and a hgsvc.vcf.gz , you would use graph: "hg38-hgsvc" in config.yaml and:

snakemake --configfile config.yaml --resources mem_mb=200000 --cores 16 construct_all

Or specify the graph label in the command line:

snakemake --configfile config.yaml --config graph="hg38-hgsvc" --resources mem_mb=200000 --cores 16 construct_all

Note the --resources mem_mb=200000 --cores 16 which specifies the maximum amount memory or cores to use.

Map reads and genotype variants

Let's assume the graph indexes constructed (see above) and you want to analyze reads for 2 samples in SAMP1/SAMP1_*.fastq.gz and SAMP2/SAMP2_*.fastq.gz .

snakemake --configfile config.yaml --config samples="SAMP1 SAMP2" --resources mem_mb=100000 --cores 16 genotype

(Assuming that the graph label, e.g. hg38-hgsvc , is defined in the config.yaml file)

Overview of the workflow

We can easily visualize the workflow with Snakemake . Her we use the small test data with two chromosomes, one sample, and reads split in 2 chunks (more info in the testdata folder ).

Rule graph

DAG of jobs

File graph

Overview of the workflow when using the Giraffe mapper

The newer and faster mapper, giraffe, uses different graph indexes.

Rule graph

DAG of jobs

File graph

FAQ

How do I specify the resources?

When running locally or on an instance, we use --resources mem_mb=200000 --cores 16 to specify the maximum memory and cores to use.

To run on an HPC, you would need to add the HPC specific configuration (more details soon).

The current resources listed in config.yaml were tuned for a human genome.

Where can I run this?

You can run this on any resource supported by Snakemake . For example, locally and on HPC.

Can I use the workflow for non-human organisms?

Yes, juts remember to update the chromosome names in the config.yaml .

Why are there AWS/S3 commands in the workflow?

We sometimes run the workflow on an AWS instance and use these command to save the outputs in a S3 bucket. If you don't want to do this, use the default config value s3save: False .

Code Snippets

shell:
    "rm -rf {params.ff}"

SnakeMake From line 59 of master/Snakefile

run:
    shell('(vg construct -t {threads} -r {input.ref} -v {input.vcf} -a -f -S -R {wildcards.chr} -C | vg ids --sort - > {output}) 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 76 of master/Snakefile

run:
    shell('vg ids --join --mapping {output.mapping} {input} 2> {log}')
    if config['s3save']:
        shell("for ff in {input}; do aws s3 cp --quiet $ff {SROOT}/; done")
        shell("aws s3 cp --quiet {output.mapping} {SROOT}/")

SnakeMake From line 89 of master/Snakefile

run:
    shell('vg index -L -t {threads} -x {output} {input.vg} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 106 of master/Snakefile

run:
    shell('vg snarls -t {threads} {input} > {output} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 120 of master/Snakefile

run:
    shell('vg prune -t {threads} -M 32 --restore-paths {input.vg} > {output} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")        

SnakeMake From line 134 of master/Snakefile

run:
    shell('mkdir -p {params.tmp_dir}')
    shell('vg index --temp-dir {params.tmp_dir} -p -t {threads} -g {output.gcsa} {input} 2> {log}')
    shell('rm -r {params.tmp_dir}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.gcsa} {SROOT}/")
        shell("aws s3 cp --quiet {output.gcsalcp} {SROOT}/")

SnakeMake From line 152 of master/Snakefile

run:
    shell('vg gbwt -n {wildcards.n} -g {output.gg} -o {output.gbwt} -x {input} -P 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.gbwt} {SROOT}/")
        shell("aws s3 cp --quiet {output.gg} {SROOT}/")

SnakeMake From line 171 of master/Snakefile

run:
    shell('vg minimizer -k {wildcards.k} -w {wildcards.w} -t {threads} -i {output} -g {input.gbwt} {input.xg} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 187 of master/Snakefile

run:
    shell('vg snarls -t {threads} --include-trivial {input} > {output} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 200 of master/Snakefile

run:
    shell('vg index -t {threads} -j {output} -x {input.xg} -s {input.snarls} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 215 of master/Snakefile

run:
    NLINES=config['nb_split_reads'] * 4
    shell('mkdir -p {output.dir}')
    shell("gzip -cd {input} | split -d -l {NLINES} --filter='pigz -p {threads} > ${{FILE}}.fastq.gz' - \"{output.dir}/{wildcards.sample}_1.part\"")

SnakeMake From line 233 of master/Snakefile

run:
    NLINES=config['nb_split_reads'] * 4
    shell('mkdir -p {output.dir}')
    shell("gzip -cd {input} | split -d -l {NLINES} --filter='pigz -p {threads} > ${{FILE}}.fastq.gz' - \"{output.dir}/{wildcards.sample}_2.part\"")

SnakeMake From line 243 of master/Snakefile

run:
    shell("vg map -t {threads} -x {input.xg} -g {input.gcsa} -f {input.r1} -f {input.r2} > {output} 2> {log}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{output}")

SnakeMake From line 273 of master/Snakefile

run:
    shell("vg mpmap -S -t {threads} -x {input.xg} -g {input.gcsa} -f {input.r1} -f {input.r2} > {output} 2> {log}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{output}")

SnakeMake From line 292 of master/Snakefile

run:
    shell("vg gaffe -p -t {threads} -m {input.min} -d {input.dist} --gbwt-name {input.gbwt} -x {input.xg} -N {wildcards.sample} -f {input.r1} -f {input.r2} > {output} 2> {log}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{output}")

SnakeMake From line 312 of master/Snakefile

run:
    shell("cat {input} > {output}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{wildcards.sample}/")

SnakeMake From line 331 of master/Snakefile

run:
    shell("vg pack -x {input.xg} -g {input.gam} -Q {wildcards.minq} -t {threads} -o {output} 2> {log}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{wildcards.sample}/")

SnakeMake From line 348 of master/Snakefile

run:
    shell("vg call -k {input.pack} -t {threads} -s {wildcards.sample} --snarls {input.snarls} {input.xg} > {params.tmp_raw_vcf} 2> {log}")
    shell("bcftools sort {params.tmp_raw_vcf} | bgzip > {output.vcf}")
    shell("tabix -f -p vcf {output.vcf}")
    shell('rm {params.tmp_raw_vcf}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.vcf} {SROOT}/{wildcards.sample}/")
        shell("aws s3 cp --quiet {output.idx} {SROOT}/{wildcards.sample}/")

SnakeMake BCFtools tabix From line 369 of master/Snakefile

run:
    shell("vg call -k {input.pack} -t {threads} -s {wildcards.sample} --snarls {input.snarls} -v {input.vcf} {input.xg} > {params.tmp_raw_vcf} 2> {log}")
    shell("bcftools view -e 'GT=\"0/0\" || GT=\"./.\"' {params.tmp_raw_vcf} | bcftools sort | bgzip > {output.vcf}")
    shell("tabix -f -p vcf {output.vcf}")
    shell('rm {params.tmp_raw_vcf}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.vcf} {SROOT}/{wildcards.sample}/")
        shell("aws s3 cp --quiet {output.idx} {SROOT}/{wildcards.sample}/")

SnakeMake BCFtools tabix From line 396 of master/Snakefile

run:
    shell('vg convert {input} -p > {output} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/")

SnakeMake From line 418 of master/Snakefile

run:
    shell('vg augment {input.pg} {input.gam} -t {threads} -m 4 -q 5 -Q 5 -A {output.gam} > {output.pg} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.gam} {SROOT}/{wildcards.sample}/")
        shell("aws s3 cp --quiet {output.pg} {SROOT}/{wildcards.sample}/")

SnakeMake From line 436 of master/Snakefile

run:
    shell('vg snarls -t {threads} {input} > {output} 2> {log}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{wildcards.sample}/")

SnakeMake From line 451 of master/Snakefile

run:
    shell("vg pack -x {input.pg} -g {input.gam} -Q {wildcards.minq} -t {threads} -o {output} 2> {log}")
    if config['s3save']:
        shell("aws s3 cp --quiet {output} {SROOT}/{wildcards.sample}/")

SnakeMake From line 467 of master/Snakefile

run:
    shell("vg call -k {input.pack} -t {threads} -s {wildcards.sample} --snarls {input.snarls} {input.xg} > {params.tmp_raw_vcf} 2> {log}")
    shell("bcftools sort {params.tmp_raw_vcf} | bgzip > {output.vcf}")
    shell("tabix -f -p vcf {output.vcf}")
    shell('rm {params.tmp_raw_vcf}')
    if config['s3save']:
        shell("aws s3 cp --quiet {output.vcf} {SROOT}/{wildcards.sample}/")
        shell("aws s3 cp --quiet {output.idx} {SROOT}/{wildcards.sample}/")