Snakemake workflow: DNA-seq variant calling with Varlociraptor

public public 1yr ago 0 bookmarks
View Workflow
Report
dna-seq-varlociraptor
dna-seq-varlociraptor
View Workflow

This is the template for a new Snakemake workflow. Replace this text with a comprehensive description covering the purpose and domain. Insert your code into the respective folders, i.e. scripts , rules , and envs . Define the entry point of the workflow in the Snakefile and the main configuration in the config.yaml file.

Authors

  • Felix Mölder (@FelixMoelder)

  • Johannes Köster (@johanneskoester)

Usage

In any case, if you use this workflow in a paper, don't forget to give credits to the authors by citing the URL of this (original) repository and, if available, its DOI (see above).

Step 1: Obtain a copy of this workflow

  1. Create a new github repository using this workflow as a template .

  2. Clone the newly created repository to your local system, into the place where you want to perform the data analysis.

Step 2: Configure workflow

Configure the workflow according to your needs via editing the file config.yaml . Add samples and sequencing units (lanes, replicates) to samples.tsv and units.tsv . For each sample, sample_name , alias , platform , and group has to be defined. Samples within the same group will be called jointly. Aliases represent the name of the sample within its group (they can be the same as the sample name, or something simpler, e.g. tumor or normal). For each group, a scenario is rendered via Jinja . Therefore, edit the template scenario ( scenario.yaml ) according to your needs. The sample sheet is available for jinja rendering as a pandas data frame in the variable samples . This allows to customize the scenario according to the contents of the sample sheet. You can therefore add additional columns to the sample sheet (e.g. purity) and access them in the scenario template, in order to pass the information to Varlociraptor.

Step 3: Execute workflow

Test your configuration by performing a dry-run via 

snakemake --use-conda -n

Execute the workflow locally via 

snakemake --use-conda --cores $N

using $N cores or run it in a cluster environment via 

snakemake --use-conda --cluster qsub --jobs 100

or 

snakemake --use-conda --drmaa --jobs 100

If you not only want to fix the software stack but also the underlying OS, use 

snakemake --use-conda --use-singularity

in combination with any of the modes above. See the Snakemake documentation for further details.

Step 4: Investigate results

After successful execution, you can create a self-contained interactive HTML report with all results via: 

snakemake --report report.html

This report can, e.g., be forwarded to your collaborators. An example (using some trivial test data) can be seen here .

Step 5: Commit changes

Whenever you change something, don't forget to commit the changes back to your github copy of the repository: 

git commit -a
git push

Step 6: Obtain updates from upstream

Whenever you want to synchronize your workflow copy with new developments from upstream, do the following.

  1. Once, register the upstream repository in your local copy: git remote add -f upstream [email protected]:snakemake-workflows/dna-seq-varlociraptor.git or git remote add -f upstream https://github.com/snakemake-workflows/dna-seq-varlociraptor.git if you do not have setup ssh keys.

  2. Update the upstream version: git fetch upstream .

  3. Create a diff with the current version: git diff HEAD upstream/master workflow > upstream-changes.diff .

  4. Investigate the changes: vim upstream-changes.diff .

  5. Apply the modified diff via: git apply upstream-changes.diff .

  6. Carefully check whether you need to update the config files: git diff HEAD upstream/master config . If so, do it manually, and only where necessary, since you would otherwise likely overwrite your settings and samples.

Step 7: Contribute back

In case you have also changed or added steps, please consider contributing them back to the original repository:

  1. Fork the original repo to a personal or lab account.

  2. Clone the fork to your local system, to a different place than where you ran your analysis.

  3. Copy the modified files from your analysis to the clone of your fork, e.g., cp -r envs rules scripts path/to/fork . Make sure to not accidentally copy config file contents or sample sheets.

  4. Commit and push your changes to your fork.

  5. Create a pull request against the original repository.

Testing

Test cases are in the subfolder .test . They are automtically executed via continuous integration with Github actions.

Code Snippets

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


shell(
    "bcftools concat {snakemake.params} -o {snakemake.output[0]} "
    "{snakemake.input.calls}")
Python Snakemake BCFtools From line 1 of concat/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


shell(
    "bcftools concat {snakemake.params} -o {snakemake.output[0]} "
    "{snakemake.input.calls}")
Python Snakemake BCFtools From line 1 of concat/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
__author__ = "Johannes Köster, Julian de Ruiter"
__copyright__ = "Copyright 2016, Johannes Köster and Julian de Ruiter"
__email__ = "[email protected], [email protected]"
__license__ = "MIT"


from os import path

from snakemake.shell import shell


# Extract arguments.
extra = snakemake.params.get("extra", "")

sort = snakemake.params.get("sort", "none")
sort_order = snakemake.params.get("sort_order", "coordinate")
sort_extra = snakemake.params.get("sort_extra", "")

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

# Check inputs/arguments.
if not isinstance(snakemake.input.reads, str) and len(snakemake.input.reads) not in {
    1,
    2,
}:
    raise ValueError("input must have 1 (single-end) or " "2 (paired-end) elements")

if sort_order not in {"coordinate", "queryname"}:
    raise ValueError("Unexpected value for sort_order ({})".format(sort_order))

# Determine which pipe command to use for converting to bam or sorting.
if sort == "none":

    # Simply convert to bam using samtools view.
    pipe_cmd = "samtools view -Sbh -o {snakemake.output[0]} -"

elif sort == "samtools":

    # Sort alignments using samtools sort.
    pipe_cmd = "samtools sort {sort_extra} -o {snakemake.output[0]} -"

    # Add name flag if needed.
    if sort_order == "queryname":
        sort_extra += " -n"

    prefix = path.splitext(snakemake.output[0])[0]
    sort_extra += " -T " + prefix + ".tmp"

elif sort == "picard":

    # Sort alignments using picard SortSam.
    pipe_cmd = (
        "picard SortSam {sort_extra} INPUT=/dev/stdin"
        " OUTPUT={snakemake.output[0]} SORT_ORDER={sort_order}"
    )

else:
    raise ValueError("Unexpected value for params.sort ({})".format(sort))

shell(
    "(bwa mem"
    " -t {snakemake.threads}"
    " {extra}"
    " {snakemake.params.index}"
    " {snakemake.input.reads}"
    " | " + pipe_cmd + ") {log}"
)
Python Snakemake SAMtools Picard From line 1 of mem/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


shell(
    "picard MarkDuplicates {snakemake.params} INPUT={snakemake.input} "
    "OUTPUT={snakemake.output.bam} METRICS_FILE={snakemake.output.metrics} "
    "&> {snakemake.log}"
)
Python Snakemake Picard From line 1 of markduplicates/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


shell("samtools index {snakemake.params} {snakemake.input[0]} {snakemake.output[0]}")
Python Snakemake SAMtools From line 1 of index/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


n = len(snakemake.input)
assert n == 2, "Input must contain 2 (paired-end) elements."

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

shell(
    "cutadapt"
    " {snakemake.params.adapters}"
    " {snakemake.params.others}"
    " -o {snakemake.output.fastq1}"
    " -p {snakemake.output.fastq2}"
    " -j {snakemake.threads}"
    " {snakemake.input}"
    " > {snakemake.output.qc} {log}"
)
Python Snakemake Cutadapt From line 3 of pe/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


log = snakemake.log_fmt_shell(stdout=False, stderr=True)

shell(
    "cutadapt"
    " {snakemake.params}"
    " -j {snakemake.threads}"
    " -o {snakemake.output.fastq}"
    " {snakemake.input[0]}"
    " > {snakemake.output.qc} {log}"
)
Python Snakemake Cutadapt From line 3 of se/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


try:
    exclude = "-x " + snakemake.input.exclude
except AttributeError:
    exclude = ""


extra = snakemake.params.get("extra", "")
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

shell(
    "OMP_NUM_THREADS={snakemake.threads} delly call {extra} "
    "{exclude} -g {snakemake.input.ref} "
    "-o {snakemake.output[0]} {snakemake.input.samples} {log}"
)
Python Snakemake delly From line 1 of delly/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2017, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell

shell.executable("bash")

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

params = snakemake.params.get("extra", "")

pipe = ""
if snakemake.output[0].endswith(".bcf"):
    pipe = "| bcftools view -Ob -"

if snakemake.threads == 1:
    freebayes = "freebayes"
else:
    chunksize = snakemake.params.get("chunksize", 100000)
    freebayes = (
        "freebayes-parallel <(fasta_generate_regions.py "
        "{snakemake.input.ref}.fai {chunksize}) "
        "{snakemake.threads}"
    ).format(snakemake=snakemake, chunksize=chunksize)

shell(
    "({freebayes} {params} -f {snakemake.input.ref}"
    " {snakemake.input.samples} {pipe} > {snakemake.output[0]}) {log}"
)
Python Snakemake BCFtools FreeBayes From line 1 of freebayes/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
__author__ = "Patrik Smeds"
__copyright__ = "Copyright 2016, Patrik Smeds"
__email__ = "[email protected]"
__license__ = "MIT"

from os import path

from snakemake.shell import shell

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

# Check inputs/arguments.
if len(snakemake.input) == 0:
    raise ValueError("A reference genome has to be provided!")
elif len(snakemake.input) > 1:
    raise ValueError("Only one reference genome can be inputed!")

# Prefix that should be used for the database
prefix = snakemake.params.get("prefix", "")

if len(prefix) > 0:
    prefix = "-p " + prefix

# Contrunction algorithm that will be used to build the database, default is bwtsw
construction_algorithm = snakemake.params.get("algorithm", "")

if len(construction_algorithm) != 0:
    construction_algorithm = "-a " + construction_algorithm

shell(
    "bwa index" " {prefix}" " {construction_algorithm}" " {snakemake.input[0]}" " {log}"
)
Python Snakemake BWA From line 1 of index/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2018, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


extra = snakemake.params.get("extra", "")
log = snakemake.log_fmt_shell(stdout=False, stderr=True)

shell(
    "picard "
    "CreateSequenceDictionary "
    "{extra} "
    "R={snakemake.input[0]} "
    "O={snakemake.output[0]} "
    "{log}"
)
Python Snakemake Picard From line 1 of createsequencedictionary/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2019, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"

import subprocess as sp
from snakemake.shell import shell

species = snakemake.params.species.lower()
release = snakemake.params.release
build = snakemake.params.build

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

suffixes = ""
datatype = snakemake.params.get("datatype", "")
if datatype == "dna":
    suffixes = ["dna.primary_assembly.fa.gz", "dna.toplevel.fa.gz"]
elif datatype == "cdna":
    suffixes = ["cdna.all.fa.gz"]
elif datatype == "cds":
    suffixes = ["cds.all.fa.gz"]
elif datatype == "ncrna":
    suffixes = ["ncrna.fa.gz"]
elif datatype == "pep":
    suffixes = ["pep.all.fa.gz"]
else:
    raise ValueError("invalid datatype, must be one of dna, cdna, cds, ncrna, pep")

success = False
for suffix in suffixes:
    url = "ftp://ftp.ensembl.org/pub/release-{release}/fasta/{species}/{datatype}/{species_cap}.{build}.{suffix}".format(
        release=release,
        species=species,
        datatype=datatype,
        build=build,
        suffix=suffix,
        species_cap=species.capitalize(),
    )

    try:
        shell("curl -sSf {url} > /dev/null 2> /dev/null")
    except sp.CalledProcessError:
        continue

    shell("(curl -L {url} | gzip -d > {snakemake.output[0]}) {log}")
    success = True
    break

if not success:
    raise ValueError(
        "Requested sequence does not seem to exist on ensembl FTP servers or servers are unavailable (url {})".format(
            url
        )
    )
Python Snakemake From line 1 of ensembl-sequence/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2019, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"

import tempfile
from snakemake.shell import shell

species = snakemake.params.species.lower()
release = snakemake.params.release
type = snakemake.params.type

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

if type == "all":
    if species == "homo_sapiens":
        suffixes = [
            "-chr{}".format(chrom) for chrom in list(range(1, 23)) + ["X", "Y", "MT"]
        ]
    else:
        suffixes = [""]
elif type == "somatic":
    suffixes = ["_somatic"]
elif type == "structural_variations":
    suffixes = ["_structural_variations"]
else:
    raise ValueError(
        "Unsupported type {} (only all, somatic, structural_variations are allowed)".format(
            type
        )
    )

urls = [
    "ftp://ftp.ensembl.org/pub/release-{release}/variation/vcf/{species}/{species}{suffix}.vcf.gz".format(
        release=release, species=species, suffix=suffix
    )
    for suffix in suffixes
]

if snakemake.input.get("fai"):
    # in case of a given .fai, reheader the VCF such that contig lengths are defined
    with tempfile.TemporaryDirectory() as tmpdir:
        shell(
            "(bcftools concat -Ob {urls} > {tmpdir}/out.bcf && "
            " bcftools reheader --fai {snakemake.input.fai} {tmpdir}/out.bcf | bcftools view -Oz -o {snakemake.output[0]}) {log}"
        )
else:
    # without .fai, just concatenate
    shell("bcftools concat -Oz {urls} > {snakemake.output[0]} {log}")
Python Snakemake BCFtools From line 1 of ensembl-variation/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
__author__ = "Michael Chambers"
__copyright__ = "Copyright 2019, Michael Chambers"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell


shell("samtools faidx {snakemake.params} {snakemake.input[0]} > {snakemake.output[0]}")
Python Snakemake SAMtools From line 1 of faidx/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

shell("tabix {snakemake.params} {snakemake.input[0]} {log}")
Python Snakemake tabix From line 1 of tabix/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2018, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


from tempfile import TemporaryDirectory
import os

from snakemake.shell import shell

extra = snakemake.params.get("extra", "")
java_opts = snakemake.params.get("java_opts", "")

with TemporaryDirectory() as tmpdir:
    recal_table = os.path.join(tmpdir, "recal_table.grp")
    log = snakemake.log_fmt_shell(stdout=True, stderr=True)
    known = snakemake.input.get("known", "")
    if known:
        known = "--known-sites {}".format(known)

    shell(
        "gatk --java-options '{java_opts}' BaseRecalibrator {extra} "
        "-R {snakemake.input.ref} -I {snakemake.input.bam} "
        "-O {recal_table} {known} {log}"
    )

    log = snakemake.log_fmt_shell(stdout=True, stderr=True, append=True)
    shell(
        "gatk --java-options '{java_opts}' ApplyBQSR -R {snakemake.input.ref} -I {snakemake.input.bam} "
        "--bqsr-recal-file {recal_table} "
        "-O {snakemake.output.bam} {log}"
    )
Python Snakemake gatk From line 1 of baserecalibrator/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
__author__ = "Bradford Powell"
__copyright__ = "Copyright 2018, Bradford Powell"
__email__ = "[email protected]"
__license__ = "BSD"


from snakemake.shell import shell
from os import path
import shutil
import tempfile

shell.executable("bash")

outcalls = snakemake.output.calls
if outcalls.endswith(".vcf.gz"):
    outprefix = "| bcftools view -Oz"
elif outcalls.endswith(".bcf"):
    outprefix = "| bcftools view -Ob"
else:
    outprefix = ""

incalls = snakemake.input[0]
if incalls.endswith(".bcf"):
    incalls = "<(bcftools view {})".format(incalls)

log = snakemake.log_fmt_shell(stdout=False, stderr=True)

extra = snakemake.params.get("extra", "")

data_dir = snakemake.params.get("data_dir", "")
if data_dir:
    data_dir = '-dataDir "%s"' % data_dir

stats = snakemake.output.get("stats", "")
csvstats = snakemake.output.get("csvstats", "")
csvstats_opt = "" if not csvstats else "-csvStats {}".format(csvstats)
stats_opt = "-noStats" if not stats else "-stats {}".format(stats)

shell(
    "(bcftools view {incalls} | "
    "snpEff {data_dir} {stats_opt} {csvstats_opt} {extra} "
    "{snakemake.params.reference} "
    "{outprefix} > {outcalls}) {log}"
)
Python Snakemake BCFtools snpEff From line 1 of snpeff/wrapper.py 
13
14
wrapper:
    "d90abc388ed573be3cf9c3d011e164b694ee509c/bio/snpeff"
SnakeMake From line 13 of annotation/pass1.smk 
32
33
shell:
    "bcftools view {input.bcf} {params.pipes} | bcftools view -Ob > {output}"
SnakeMake BCFtools From line 32 of annotation/pass2.smk 
46
47
shell:
    "rbt vcf-annotate-dgidb {input} > {output}"
SnakeMake rust-bio-tools From line 46 of annotation/pass2.smk 
57
58
shell:
    "wget {params.zip} -O {output}"
SnakeMake From line 57 of annotation/pass2.smk 
69
70
71
72
73
74
75
76
shell:
    """
    unzip {input} -d resources/dbnsfp/ &&
    (zcat resources/dbnsfp/*_variant.chr1.gz |
    head -n 1 ; zcat resources/dbnsfp/*_variant.chr* |
    grep -v '^#' ) | bgzip -@ {threads} > {output} &&
    rm -r resources/dbnsfp
    """
SnakeMake From line 69 of annotation/pass2.smk 
90
91
92
shell:
    "bcftools view {input.bcf} | SnpSift dbnsfp -db {input.db} -f {params.fields} {params.extra} /dev/stdin | "
    "sed 's/\\(^##INFO=<ID=dbNSFP_\\w*,Number=\\)A/\\1./g' | bcftools view -Ob > {output}"
SnakeMake BCFtools SnpSift From line 90 of annotation/pass2.smk 
102
103
shell:
    "tabix -s 1 -b 2 -e 2 {input}"
SnakeMake tabix From line 102 of annotation/pass2.smk 
10
11
script:
    "../scripts/render-scenario.py"
SnakeMake From line 10 of rules/calling.smk 
26
27
28
shell:
    "varlociraptor preprocess variants --candidates {input.candidates} "
    "{input.ref} --bam {input.bam} --output {output} 2> {log}"
SnakeMake Varlociraptor From line 26 of rules/calling.smk 
42
43
44
45
shell:
    "varlociraptor "
    "call variants generic --obs {params.obs} "
    "--scenario {input.scenario} > {output} 2> {log}"
SnakeMake Varlociraptor From line 42 of rules/calling.smk 
61
62
wrapper:
    "0.36.0/bio/bcftools/concat"
SnakeMake From line 61 of rules/calling.smk 
14
15
wrapper:
    "0.43.0/bio/freebayes"
SnakeMake FreeBayes From line 14 of rules/candidate_calling.smk 
29
30
wrapper:
    "0.43.0/bio/delly"
SnakeMake delly From line 29 of rules/candidate_calling.smk 
7
8
shell:
    "fasterq-dump {wildcards.accession}"
SnakeMake From line 7 of rules/download.smk 
10
11
shell:
    "bcftools view {input} | SnpSift filter \"{params.filter}\" | bcftools view -Ob > {output}"
SnakeMake BCFtools SnpSift From line 10 of rules/filtering.smk 
24
25
26
shell:
    "varlociraptor filter-calls control-fdr {input} --var {wildcards.vartype} "
    "--events {params.events} --fdr {params.threshold} > {output}"
SnakeMake Varlociraptor From line 24 of rules/filtering.smk 
46
47
wrapper:
    "0.37.1/bio/bcftools/concat"
SnakeMake From line 46 of rules/filtering.smk 
15
16
wrapper:
    "0.39.0/bio/bwa/mem"
SnakeMake From line 15 of rules/mapping.smk 
29
30
wrapper:
    "0.39.0/bio/picard/markduplicates"
SnakeMake From line 29 of rules/mapping.smk 
47
48
wrapper:
    "0.47.0/bio/gatk/baserecalibrator"
SnakeMake From line 47 of rules/mapping.smk 
8
9
script:
    "../scripts/build_oncoprint_matrix.py"
SnakeMake From line 8 of rules/oncoprint.smk 
18
19
script:
    "../scripts/oncoprint.R"
SnakeMake From line 18 of rules/oncoprint.smk 
 9
10
wrapper:
    "0.45.1/bio/reference/ensembl-sequence"
SnakeMake From line 9 of rules/ref.smk 
18
19
wrapper:
    "0.45.1/bio/samtools/faidx"
SnakeMake From line 18 of rules/ref.smk 
29
30
wrapper:
    "0.45.1/bio/picard/createsequencedictionary"
SnakeMake From line 29 of rules/ref.smk 
43
44
wrapper:
    "0.45.1/bio/reference/ensembl-variation"
SnakeMake From line 43 of rules/ref.smk 
54
55
shell:
    "rbt vcf-fix-iupac-alleles < {input} | bcftools view -Oz > {output}"
SnakeMake BCFtools rust-bio-tools From line 54 of rules/ref.smk 
67
68
wrapper:
    "0.45.1/bio/tabix"
SnakeMake tabix From line 67 of rules/ref.smk 
80
81
wrapper:
    "0.45.1/bio/bwa/index"
SnakeMake From line 80 of rules/ref.smk 
13
14
15
16
17
18
shell:
    "bcftools view {input.bcf} > {output}.vcf;"
    "create_report {output}.vcf {input.ref} --flanking 100 "
    "--info-columns ANN dgiDB_drugs cosmic_LEGACY_ID --info-columns-prefixes PROB_ dbNSFP_ --sample-columns DP AF OBS"
    " --template {params} --tracks {input.bams} --output {output} --standalone; "
    "rm {output}.vcf"
SnakeMake BCFtools From line 13 of rules/report.smk 
11
12
13
14
15
16
shell:
    "varlociraptor estimate tmb "
    "--coding-genome-size {params.coding_genome_size} "
    "--somatic-tumor-events {params.somatic_events} "
    "--tumor-sample {params.tumor_sample} "
    "< {input} > {output}"
SnakeMake Varlociraptor From line 11 of rules/tmb.smk 
15
16
wrapper:
    "0.42.0/bio/cutadapt/pe"
SnakeMake From line 15 of rules/trimming.smk 
30
31
wrapper:
    "0.42.0/bio/cutadapt/se"
SnakeMake From line 30 of rules/trimming.smk 
38
39
40
41
run:
    if input[0].endswith(".gz"):
        shell("cat {input} > {output}")
    else:
SnakeMake From line 38 of rules/trimming.smk 
8
9
shell:
    "bcftools index {input}"
SnakeMake BCFtools From line 8 of rules/utils.smk 
17
18
wrapper:
    "0.39.0/bio/samtools/index"
SnakeMake From line 17 of rules/utils.smk 
8
9
shell:
    "vl2svg {input} {output}"
SnakeMake From line 8 of rules/vega.smk 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import pysam
import pandas as pd
import os


input_files = snakemake.input

df = pd.DataFrame(columns=["Sample"])
for sample_file in input_files:
    variant_file = pysam.VariantFile(sample_file)
    sample_name = os.path.basename(sample_file).split(".")[0]
    gene_variant_dict = {"Sample": [sample_name]}
    for rec in variant_file.fetch():
        for sample in rec.samples:
            allele_frequencies = rec.samples[sample]["AF"] #Can be multiple entries
            for allele_frequency in allele_frequencies:
                variant = rec.info["SVLEN"]
                if variant[0]:
                    variant_type = "INDEL"
                else:
                    variant_type = "SNV" 
                transcripts = rec.info["ANN"]
                for transcript in transcripts:
                    gene = transcript.split("|")[3]
                    if gene not in gene_variant_dict:
                        gene_variant_dict[gene] = set()
                    gene_variant_dict[gene].add(variant_type)
                break
    for key, value in gene_variant_dict.items():
        gene_variant_dict[key] = ','.join(value)
    sample_df = pd.DataFrame(gene_variant_dict, index=[0])
    df = pd.concat([df, sample_df], join="outer", ignore_index=False, sort=False)
df.set_index("Sample", inplace=True)
with open(snakemake.output[0], 'w') as output_f:
    print(df.to_csv(sep="\t", index=True), file=output_f)
Python Pandas pysam From line 1 of scripts/build_oncoprint_matrix.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
library(ComplexHeatmap)
library(ggplot2)

table = read.table(snakemake@input[[1]], sep="\t", header=TRUE, row.names=1)
mat = as.matrix(table)
mat = t(mat)


col = c(SNV = "blue", INDEL = "red")

alter_fun = list(
        SNV = function(x, y, w, h) grid.rect(x, y, w*0.9, h*0.9, 
            gp = gpar(fill = col["SNV"], col = NA)),
        INDEL = function(x, y, w, h) grid.rect(x, y, w*0.9, h*0.4, 
            gp = gpar(fill = col["INDEL"], col = NA))
    )


heatmap_legend_param = list(title = "Alterations", at = c("SNV", "INDEL"), 
        labels = c("SNV", "INDEL"))

mat <- mat[order(apply(mat, 1, function(row) sum(row != "")), decreasing = T), ]

if (nrow(mat) > 2000) {
    mat <- mat[1:2000,]
}
rows_matrix <- nrow(mat)
height_plot <- (rows_matrix/5)
if (height_plot < 4) {
    height_plot <- 4
}
pdf(file = snakemake@output[[1]], height=height_plot)
if (rows_matrix > 0) {
    oncoprint <- oncoPrint(mat,
        alter_fun = alter_fun, col = col, 
        remove_empty_columns = FALSE, remove_empty_rows = TRUE,
        pct_side = "right", row_names_side = "left",
        show_column_names=T,
        column_title = "OncoPrint", heatmap_legend_param = heatmap_legend_param)
    draw(oncoprint, newpage=F)
}
dev.off()
R ggplot2 ComplexHeatmap From line 1 of scripts/oncoprint.R 
1
2
3
4
5
6
7
8
from jinja2 import Template
import pandas as pd

with open(snakemake.input[0]) as template, open(snakemake.output[0], "w") as out:
    samples = snakemake.params.samples
    out.write(Template(template.read()).render(
        samples=samples[samples["group"] == snakemake.wildcards.group]
    ))
Python Pandas Jinja2 From line 1 of scripts/render-scenario.py
ShowHide 42 more snippets with no or duplicated tags.

Login to post a comment if you would like to share your experience with this workflow.

Do you know this workflow well? If so, you can request seller status , and start supporting this workflow.

Free

Created: 1yr ago
Updated: 1yr ago
Maitainers: public
URL: https://github.com/vsoch/dna-seq-varlociraptor
Name: dna-seq-varlociraptor
Version: 1
Badge:
workflow icon

Insert copied code into your website to add a link to this workflow.

Other Versions:
Version 2 Version v5.0.1
Downloaded: 0
Copyright: Public Domain
License: MIT License
  • Future updates

Related Workflows

psychip_snakemake
Show Details View Workflow
ENCODE pipeline for histone marks developed for the psychENCODE project
public

public

psychip pipeline is an improved version of the ENCODE pipeline for histone marks developed for the psychENCODE project. The o...
raw sequence reads Alignment Sequence alignment report macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
ncov_2
Show Details View Workflow
Near-real time tracking of SARS-CoV-2 in Connecticut
public

public

Repository containing scripts to perform near-real time tracking of SARS-CoV-2 in Connecticut using genomic data. This pipeli...
JSON nextclade Augur Biopython FOCUS Pandas Snakemake bs4 epiweeks geopy matplotlib numpy pycountry pycountry-convert uszipcode
Free
28
cellranger-snakemake-gke
Show Details View Workflow
snakemake workflow to run cellranger on a given bucket using gke.
public

public

A Snakemake workflow for running cellranger on a given bucket using Google Kubernetes Engine. The usage of this workflow ...
macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
atlas
Show Details View Workflow
ATLAS - Three commands to start analyzing your metagenome data
public

public

Metagenome-atlas is a easy-to-use metagenomic pipeline based on snakemake. It handles all steps from QC, Assembly, Binning, t...
raw sequence reads Genome assembly Annotation track checkm2 gunc prodigal snakemake-wrapper-utils MEGAHIT Atlas BBMap Biopython BioRuby Bwa-mem2 cd-hit CheckM DAS Diamond eggNOG-mapper v2 MetaBAT 2 Minimap2 MMseqs MultiQC Pandas Picard pyfastx SAMtools SemiBin Snakemake SPAdes SqueezeMeta TADpole VAMB CONCOCT ete3 gtdbtk h5py networkx numpy plotly psutil utils metagenomics
Free
175
rna-seq-star-deseq2
Show Details View Workflow
RNA-seq workflow using STAR and DESeq2
public

public

This workflow performs a differential gene expression analysis with STAR and Deseq2. The usage of this workflow is described ...
raw sequence reads Gene expression profiling ID list snakemake-wrapper-utils biomaRt PCAtools BiocParallel BWA Cutadapt DESeq2 GFFutils MultiQC Pandas SAMtools Snakemake STAR cli stringr tidyverse RNA-Seq
Free
41
dna-seq-gatk-variant-calling
Show Details View Workflow
This Snakemake pipeline implements the GATK best-practices workflow
public

public

This Snakemake pipeline implements the GATK best-practices workflow for calling small germline variants. The usage of thi...
VCF raw sequence reads Variant calling genetic variants gatk rust-bio-tools snakemake-wrapper-utils tabix BCFtools BWA FastQC MultiQC Pandas Picard SAMtools Snakemake Trimmomatic Variant Effect Predictor (VEP) common matplotlib numpy seaborn DNA
Free
57