Snakefile

import pandas

# import sample information
sinfo = pandas.read_csv(config['sample_info'], sep='\t', dtype={"sample": str}).set_index('sample', drop=False).sort_index()

if len(config['sample']) == 0:
    config['sample'] = list(sinfo['sample'])
elif not isinstance(config['sample'], list):
    config['sample'] = config['sample'].split(',')

# optional: TSV with (sample, hpo) columns. hpo is a comma-sep list of HPO terms
hpo = {}
hpo_samples = []
if 'hpo_tsv' in config:
    hpo = pandas.read_csv(config['hpo_tsv'], sep='\t', dtype={"sample": str}).set_index('sample', drop=False).sort_index()
    for samp in config['sample']:
        if samp in hpo['sample']:
            hpo_samples.append(samp)

# enumerate probands (samples with parents)
trio_probands = []
if 'bam_p1' in sinfo:
    for samp in sinfo[(sinfo.bam_p1.notnull()) & (sinfo.bam_p2.notnull())]['sample']:
        if samp in config['sample']:
            trio_probands.append(samp)

# AnnotSV output
annotsv_tsvs = []
for samp in config['sample']:
    if samp in hpo_samples:
        annotsv_tsvs.append("results/{sample}/snpeff.{sample}.svs_napu.sveval.annotsv_pheno.tsv.gz".format(sample=samp))
    else:
        annotsv_tsvs.append("results/{sample}/snpeff.{sample}.svs_napu.sveval.annotsv.tsv.gz".format(sample=samp))

# main rules
rule all:
    input:
        vcf=expand("results/{sample}/{sample}.annot.vcf.gz", sample=config['sample']),
        annotsv=annotsv_tsvs,
        tr=expand("results/{sample}/{sample}.tr.tsv.gz", sample=config['sample'])

rule annot_tr:
    input: expand("results/{sample}/{sample}.tr.tsv.gz", sample=config['sample'])

rule stats:
    input: expand("results/{sample}/{sample}.vcf_stats.tsv", sample=config['sample'])

rule denovo_svs:
    input: expand("results/{sample}/{sample}.denovo_svs.vcf.gz", sample=trio_probands)

rule reports:
    input: expand("results/{sample}/{sample}.variant_report_napu.html", sample=config['sample'])

rule reports_sniffles:
    input: expand("results/{sample}/{sample}.variant_report_napu_sniffles.html", sample=config['sample'])

# RULES

def gethvcf(wildcards):
    return(sinfo.hvcf[wildcards.sample])

def getsniffles(wildcards):
    return(sinfo.sniffles[wildcards.sample])

def gethvcf_p1(wildcards):
    return(sinfo.hvcf_p1[wildcards.sample])

def gethvcf_p2(wildcards):
    return(sinfo.hvcf_p2[wildcards.sample])

def getbam(wildcards):
    return(sinfo.bam[wildcards.sample])

def getbamp1(wildcards):
    return(sinfo.bam_p1[wildcards.sample])

def getbamp2(wildcards):
    return(sinfo.bam_p2[wildcards.sample])

rule split_multiallelic_vcf:
    input: gethvcf
    output: temp("results/{sample}/{sample}.norm.vcf.gz")
    threads: 4
    resources:
        mem='8G',
        runtime='1h'
    benchmark: 'benchmark/split_multiallelic_vcf.{sample}.tsv'
    shell:
        """
        bcftools norm -m -both --threads {threads} -Oz -o {output} {input}
        """

rule split_multiallelic_vcf_sniffles:
    input: getsniffles
    output: temp("results/{sample}/sniffles.{sample}.norm.vcf.gz")
    threads: 4
    resources:
        mem='8G',
        runtime='1h'
    benchmark: 'benchmark/split_multiallelic_vcf.{sample}.tsv'
    shell:
        """
        zcat {input} | sed "s/FORMAT\tSAMPLE/FORMAT\t{wildcards.sample}/" | bcftools view -i "INFO/SVLEN>-20000000 & INFO/SVLEN<20000000" | bcftools norm -m -both --threads {threads} -Oz -o {output}
        """

rule keep_nonref_vcf:
    input: "results/{sample}/{sample}.norm.vcf.gz"
    output: temp("results/{sample}/{sample}.norm.nonref.vcf.gz")
    resources:
        mem='8G',
        runtime='1h'
    benchmark: 'benchmark/keep_nonref_vcf.{sample}.tsv'
    shell:
        """
        bcftools view -c 1 -Oz -o {output} {input}
        """

rule sort_vcf:
    input: "results/{file}.vcf.gz"
    resources:
        mem='8G',
        runtime='3h'
    output:
        vcf="results/{file}.sorted.vcf.gz",
        tbi="results/{file}.sorted.vcf.gz.tbi"
    params:
        tmp_dir=config['tmp_dir']
    shell:
        """
        bcftools sort -T {params.tmp_dir} -Oz -o {output.vcf} {input}
        bcftools index -t -o {output.tbi} {output.vcf}
        """

rule annotate_with_snpeff:
    input:
        vcf="results/{sample}/{sample}.norm.nonref.sorted.vcf.gz",
        db=config['snpeff'] + "/GRCh38.105/cytoBand.txt.gz"
    output: "results/{sample}/snpeff.{sample}.vcf.gz"
    params:
        db=config['snpeff']
    container: 'docker://quay.io/biocontainers/snpeff@sha256:7ac091da707f5d63f307eef4ee57c3f0e94eed49f86bbdace3d4be3a514ed410'
    benchmark: 'benchmark/annotate_with_snpeff.{sample}.tsv'
    resources:
        mem='16G',
        runtime='3h'
    shell:
        """
        zcat {input.vcf} | snpEff -Xmx16g -nodownload -no-intergenic -dataDir "${{PWD}}/{params.db}" GRCh38.105 | gzip > {output}
        """

rule annotate_sniffles_with_snpeff:
    input:
        vcf="results/{sample}/sniffles.{sample}.norm.sorted.vcf.gz",
        db=config['snpeff'] + "/GRCh38.105/cytoBand.txt.gz"
    output: "results/{sample}/snpeff.{sample}.svs_sniffles.vcf.gz"
    params:
        db=config['snpeff']
    container: 'docker://quay.io/biocontainers/snpeff@sha256:7ac091da707f5d63f307eef4ee57c3f0e94eed49f86bbdace3d4be3a514ed410'
    benchmark: 'benchmark/annotate_sniffles_with_snpeff.{sample}.tsv'
    resources:
        mem='16G',
        runtime='3h'
    shell:
        """
        zcat {input.vcf} | snpEff -Xmx16g -nodownload -no-intergenic -dataDir "${{PWD}}/{params.db}" GRCh38.105 | gzip > {output}
        """

rule split_small_large_variants:
    input: "results/{sample}/snpeff.{sample}.vcf.gz"
    output:
        small=temp("results/{sample}/snpeff.{sample}.smallvars.vcf.gz"),
        sv="results/{sample}/snpeff.{sample}.svs_napu.vcf.gz"
    resources:
        mem='8G',
        runtime='6h'
    params:
        tmp_dir=config['tmp_dir']
    benchmark: 'benchmark/split_small_large_variants.{sample}.tsv'
    shell:
        """
        bcftools view -i "STRLEN(REF)>=30 | MAX(STRLEN(ALT))>=30" -Oz -o {output.sv} {input}
        bcftools view -i "STRLEN(REF)<30 & MAX(STRLEN(ALT))<30" {input} | bcftools sort -T {params.tmp_dir} -Oz -o {output.small}
        """

rule annotate_svs_sveval:
    input:
        script="scripts/annotate_svs.R",
        vcf="results/{sample}/snpeff.{sample}.svs{svmeth}.vcf.gz",
        db=config['svdb']
    output: "results/{sample}/snpeff.{sample}.svs{svmeth}.sveval.vcf.gz"
    container: "docker://quay.io/jmonlong/svannotate_sveval@sha256:c694bd3db95a07a49a231856c4662326976cc7f1d9875f746adb122ffc094ad6"
    params:
        tmp_dir=config['tmp_dir'],
        tmp_vcf=config['tmp_dir'].rstrip('/') + "/temp_{sample}_svs{svmeth}.annotated.vcf"
    resources:
        mem='8G',
        runtime='3h'
    benchmark: 'benchmark/annotate_svs_sveval.svs{svmeth}.{sample}.tsv'
    shell:
        """
        # Rscript {input.script} {input.vcf} {input.db} {params.tmp_vcf}
        gunzip -c {input.vcf} | bgzip > {params.tmp_vcf}
        R CMD BATCH '--args {params.tmp_vcf} {input.db} {params.tmp_vcf}' {input.script} {params.tmp_vcf}.out
        bcftools sort -T {params.tmp_dir} -Oz -o {output} {params.tmp_vcf}
        rm -f {params.tmp_vcf}
        """

rule annotate_svs_annotsv:
    input:
        vcf="results/{sample}/snpeff.{sample}.svs{svmeth}.sveval.vgval.vcf.gz",
        db=config['annotsv'] + '/Annotations_Human/Genes/GRCh38/genes.RefSeq.sorted.bed'
    output: "results/{sample}/snpeff.{sample}.svs{svmeth}.sveval.annotsv.tsv.gz"
    log: 'logs/annotate_svs{svmeth}_annotsv.{sample}.log'
    container: "docker://quay.io/jmonlong/annotsv:3.4"
    params:
        tmp_dir=config['tmp_dir'].rstrip('/') + "/temp_annotsv.svs{svmeth}.{sample}",
        annotsvdb=config['annotsv']
    resources:
        mem='32G',
        runtime='3h'
    benchmark: 'benchmark/annotate_svs{svmeth}_annotsv.{sample}.tsv'
    shell:
        """
        AnnotSV -annotationsDir {params.annotsvdb} -SvinputFile {input.vcf} -outputDir {params.tmp_dir} > {log}
        gzip -c {params.tmp_dir}/*.annotated.tsv > {output}
        rm  -rf {params.tmp_dir}
        """

rule annotate_svs_annotsv_pheno:
    input:
        vcf="results/{sample}/snpeff.{sample}.svs{svmeth}.sveval.vgval.vcf.gz",
        db=config['annotsv'] + '/Annotations_Human/Genes/GRCh38/genes.RefSeq.sorted.bed'
    output: "results/{sample}/snpeff.{sample}.svs{svmeth}.sveval.annotsv_pheno.tsv.gz"
    log: 'logs/annotate_svs{svmeth}_annotsv_pheno.{sample}.log'
    container: "docker://quay.io/jmonlong/annotsv:3.4"
    params:
        tmp_dir=config['tmp_dir'].rstrip('/') + "/temp_annotsv.svs{svmeth}.{sample}",
        annotsvdb=config['annotsv'],
        hpo=lambda wildcards: hpo.hpo[wildcards.sample]
    resources:
        mem='32G',
        runtime='3h'
    benchmark: 'benchmark/annotate_svs{svmeth}_annotsv.{sample}.tsv'
    shell:
        """
        AnnotSV -annotationsDir {params.annotsvdb} -SvinputFile {input.vcf} -hpo {params.hpo} -outputDir {params.tmp_dir} > {log}
        gzip -c {params.tmp_dir}/*.annotated.tsv > {output}
        rm  -rf {params.tmp_dir}
        """

rule merge_annotsv_tsvs:
    input:
        a_napu="results/{sample}/snpeff.{sample}.svs_napu.sveval.annotsv{pheno}.tsv.gz",
        a_sniffles="results/{sample}/snpeff.{sample}.svs_sniffles.sveval.annotsv{pheno}.tsv.gz"
    output: "results/{sample}/snpeff.{sample}.svs_napu_sniffles.sveval.annotsv{pheno}.tsv.gz"
    params:
        tmp_tsv=config['tmp_dir'].rstrip('/') + "/temp_merge_annotsv.{sample}.{pheno}.tsv"
    localrule: True
    shell:
        """
        gunzip -c {input.a_napu} > {params.tmp_tsv}
        gunzip -c {input.a_sniffles} | sed 1d >> {params.tmp_tsv}
        gzip -c {params.tmp_tsv} > {output}
        rm -f {params.tmp_tsv}
        """

rule annotate_smallvars_snpsift:
    input:
        vcf="results/{sample}/snpeff.{sample}.smallvars.vcf.gz",
        gnomad=config['gnomad'],
        clinvar=config['clinvar'],
        dbnsfp=config['dbnsfp']
    output: "results/{sample}/snpeff.{sample}.smallvars.annot.vcf.gz"
    container: "docker://quay.io/biocontainers/snpsift@sha256:049babfac841d15a92d8febfc10a25f5aa109c9fe6670af35ea79583a1c78402"
    params:
        tmp_vcf=config['tmp_dir'].rstrip('/') + '/temp_snpsift.{sample}.vcf'
    threads: 4
    resources:
        mem="8G",
        runtime='3h'
    benchmark: 'benchmark/annotate_smallvars_snpsift.{sample}.tsv'
    shell:
        """
        zcat {input.vcf} | \
        SnpSift -Xmx4g annotate -info CLNSIG {input.clinvar} | \
        SnpSift -Xmx1g filter "(ANN[*].IMPACT has 'HIGH') | (ANN[*].IMPACT has 'MODERATE') | ((exists LOF[*].PERC) & (LOF[*].PERC > 0.9)) | (exists CLNSIG)" | \
        SnpSift -Xmx8g annotate -noId {input.gnomad} > {params.tmp_vcf}
        SnpSift -Xmx8g dbnsfp -db {input.dbnsfp} -f GERP++_RS,CADD_raw,CADD_phred,MetaRNN_score,MetaRNN_pred,ALFA_Total_AF {params.tmp_vcf} | gzip > {output}
        rm -f {params.tmp_vcf}
        """

rule validate_svs_with_vg:
    input:
        vcf="results/{sample}/snpeff.{sample}.svs_napu.sveval.vcf.gz",
        bam=getbam,
        ref=config['ref'],
        script='scripts/validate-svs.py'
    output: "results/{sample}/snpeff.{sample}.svs_napu.sveval.vgval.vcf.gz"
    container: 'docker://quay.io/jmonlong/svvalidate_vgcall@sha256:8020ea9a9fd784c5dab74fd366cd089df5bb7f17154b79e4dcc780596c30b0ed'
    threads: 8
    params:
        tmp_dir=config['tmp_dir'].rstrip('/') + '/temp.{sample}.svs.sveval'
    resources:
        mem="8G",
        runtime='6h',
        tmpdir=config['tmp_dir']
    benchmark: 'benchmark/validate_svs_with_vg.{sample}.tsv'
    shell:
        """
        python3 {input.script} -b {input.bam} -f {input.ref} -v {input.vcf} -d {params.tmp_dir} -t {threads} -o {output}
        rm -rf {params.tmp_dir}
        """

rule dont_validate_sniffles_svs_with_vg:
    input: "results/{sample}/snpeff.{sample}.svs_sniffles.sveval.vcf.gz",
    output: "results/{sample}/snpeff.{sample}.svs_sniffles.sveval.vgval.vcf.gz"
    localrule: True
    shell: "cp {input} {output}"

rule merge_small_svs:
    input:
        svs="results/{sample}/snpeff.{sample}.svs_napu.sveval.vgval.vcf.gz",
        small="results/{sample}/snpeff.{sample}.smallvars.annot.vcf.gz"
    output: "results/{sample}/{sample}.annot.vcf.gz"
    params:
        tmp_dir=config['tmp_dir'],
        tmp_small=config['tmp_dir'].rstrip('/') + '/temp.merge_small_svs.{sample}.small.vcf.gz',
        tmp_svs=config['tmp_dir'].rstrip('/') + '/temp.merge_small_svs.{sample}.svs.vcf.gz'
    benchmark: 'benchmark/merge_small_svs.{sample}.tsv'
    resources:
        mem='8G',
        runtime='1h'
    shell:
        """
        bcftools sort -T {params.tmp_dir} -Oz -o {params.tmp_small} {input.small}
        bcftools index -t {params.tmp_small}
        bcftools sort -T {params.tmp_dir} -Oz -o {params.tmp_svs} {input.svs}
        bcftools index -t {params.tmp_svs}
        bcftools concat -a -Oz -o {output} {params.tmp_small} {params.tmp_svs}
        rm -f {params.tmp_small} {params.tmp_svs} {params.tmp_small}.tbi {params.tmp_svs}.tbi
        """

rule merge_small_svs_with_sniffles:
    input:
        svs_napu="results/{sample}/snpeff.{sample}.svs_napu.sveval.vgval.vcf.gz",
        svs_sniffles="results/{sample}/snpeff.{sample}.svs_sniffles.sveval.vcf.gz",
        small="results/{sample}/snpeff.{sample}.smallvars.annot.vcf.gz"
    output: "results/{sample}/{sample}.annot.with_sniffles.vcf.gz"
    params:
        tmp_dir=config['tmp_dir'],
        tmp_small=config['tmp_dir'].rstrip('/') + '/temp.merge_small_svs_with_sniffles.{sample}.small.vcf.gz',
        tmp_svs_napu=config['tmp_dir'].rstrip('/') + '/temp.merge_small_svs_with_sniffles.{sample}.svs_napu.vcf.gz',
        tmp_svs_sniffles=config['tmp_dir'].rstrip('/') + '/temp.merge_small_svs_with_sniffles.{sample}.svs_sniffles.vcf.gz'
    benchmark: 'benchmark/merge_small_svs.{sample}.tsv'
    resources:
        mem='8G',
        runtime='1h'
    shell:
        """
        bcftools sort -T {params.tmp_dir} -Oz -o {params.tmp_small} {input.small}
        bcftools index -t {params.tmp_small}
        bcftools sort -T {params.tmp_dir} -Oz -o {params.tmp_svs_napu} {input.svs_napu}
        bcftools index -t {params.tmp_svs_napu}
        bcftools sort -T {params.tmp_dir} -Oz -o {params.tmp_svs_sniffles} {input.svs_sniffles}
        bcftools index -t {params.tmp_svs_sniffles}
        bcftools concat -a -Oz -o {output} {params.tmp_small} {params.tmp_svs_napu} {params.tmp_svs_sniffles}
        rm -f {params.tmp_small} {params.tmp_svs_napu} {params.tmp_svs_sniffles} {params.tmp_small}.tbi {params.tmp_svs_napu}.tbi {params.tmp_svs_sniffles}.tbi
        """

rule prepare_tr:
    input:
        vcf=gethvcf,
        script='scripts/prepare_tr_from_vcf.py'
    output: "results/{sample}/{sample}.tr.tsv.gz"
    benchmark: 'benchmark/prepare_tr.{sample}.tsv'
    shell: "zcat {input.vcf} | python3 {input.script} | gzip > {output}"

rule get_vcf_stats:
    input:
        vcf=gethvcf,
        script='scripts/get_vcf_stats.py'
    output: "results/{sample}/{sample}.vcf_stats.tsv"
    benchmark: 'benchmark/get_vcf_stats.{sample}.tsv'
    resources:
        runtime=600
    shell: "python3 {input.script} -v {input.vcf} -o {output}"

rule annotate_denovo_svs_sveval:
    input:
        script="scripts/annotate_denovo_svs.R",
        vcf="results/{sample}/snpeff.{sample}.svs.vcf.gz",
        vcf_p1=gethvcf_p1,
        vcf_p2=gethvcf_p2,
        db=config['svdb']
    output: "results/{sample}/snpeff.{sample}.svs.sveval.denovo_cand.vcf.gz"
    container: "docker://quay.io/jmonlong/svannotate_sveval@sha256:c694bd3db95a07a49a231856c4662326976cc7f1d9875f746adb122ffc094ad6"
    params:
        tmp_dir=config['tmp_dir'],
        tmp_vcf=config['tmp_dir'].rstrip('/') + "/temp_{sample}_svs.denovo.annotated.vcf"
    resources:
        mem='8G',
        runtime='3h'
    benchmark: 'benchmark/annotate_svs_sveval.{sample}.tsv'
    shell:
        """
        Rscript {input.script} {input.vcf} {input.vcf_p1} {input.vcf_p2} {input.db} {params.tmp_vcf}
        bcftools sort -T {params.tmp_dir} -Oz -o {output} {params.tmp_vcf}
        rm -f {params.tmp_vcf}
        """

rule validate_denovo_svs_with_vg:
    input:
        vcf="results/{sample}/snpeff.{sample}.svs.sveval.denovo_cand.vcf.gz",
        bam=getbam,
        bam_p1=getbamp1,
        bam_p2=getbamp2,
        ref=config['ref'],
        script='scripts/validate-svs-trio.py'
    output: "results/{sample}/{sample}.denovo_svs.vcf.gz"
    container: 'docker://quay.io/jmonlong/svvalidate_vgcall@sha256:8020ea9a9fd784c5dab74fd366cd089df5bb7f17154b79e4dcc780596c30b0ed'
    threads: 8
    params:
        tmp_dir=config['tmp_dir'].rstrip('/') + '/temp.{sample}.validate.denovo'
    resources:
        mem="8G",
        runtime='6h',
        tmpdir=config['tmp_dir']
    benchmark: 'benchmark/validate_denovo_svs_with_vg.{sample}.tsv'
    shell:
        """
        python3 {input.script} -b {input.bam},{input.bam_p1},{input.bam_p2} -f {input.ref} -v {input.vcf} -d {params.tmp_dir} -t {threads} -o {output}
        rm -rf {params.tmp_dir}
        """

def getannotsv(wildcards):
    if wildcards.sample in hpo_samples:
        return("results/{sample}/snpeff.{sample}.svs_{svmeth}.sveval.annotsv_pheno.tsv.gz")
    else:
        return("results/{sample}/snpeff.{sample}.svs_{svmeth}.sveval.annotsv.tsv.gz")

def getVcfForReport(wildcards):
    if wildcards.svmeth == 'napu':
        return("results/{sample}/{sample}.annot.vcf.gz")
    elif wildcards.svmeth == 'napu_sniffles':
        return("results/{sample}/{sample}.annot.with_sniffles.vcf.gz")

if 'gene_list' in config:
    rule make_report:
        input:
            rmd="scripts/variant-report-sample.Rmd",
            vcf=getVcfForReport,
            annotsv=getannotsv,
            tr="results/{sample}/{sample}.tr.tsv.gz",
            tr_control=config['tr_catalog'],
            db=config['svdb'],
            gl=config['gene_list']
        output:
            html="results/{sample}/{sample}.variant_report_{svmeth}.html",
            rds="results/{sample}/{sample}.variant_report_{svmeth}.rds"
        params:
            tmp_rmd="{sample}.variant_report_{svmeth}.Rmd",
            tmp_html="{sample}.variant_report_{svmeth}.html"
        resources:
            mem='16G',
            runtime='3h'
        container: 'docker://quay.io/jmonlong/svannotate_sveval:0.5'
        benchmark: 'benchmark/make_report_{svmeth}.{sample}.tsv'
        shell:
            """
            cp {input.rmd} {params.tmp_rmd}
            Rscript -e 'rmarkdown::render("{params.tmp_rmd}")' -g {input.gl} -v {input.vcf} -r {input.tr} -rc {input.tr_control} -d {input.db} -a {input.annotsv} -o {output.rds}
            mv {params.tmp_html} {output.html}
            rm -f {params.tmp_rmd}
            """
else:
    rule make_report:
        input:
            rmd="scripts/variant-report-sample.Rmd",
            vcf=getVcfForReport,
            annotsv=getannotsv,
            tr="results/{sample}/{sample}.tr.tsv.gz",
            tr_control=config['tr_catalog'],
            db=config['svdb']
        output:
            html="results/{sample}/{sample}.variant_report_{svmeth}.html",
            rds="results/{sample}/{sample}.variant_report_{svmeth}.rds"
        params:
            tmp_rmd="{sample}.variant_report_{svmeth}.Rmd",
            tmp_html="{sample}.variant_report_{svmeth}.html"
        resources:
            mem='16G',
            runtime='3h'
        benchmark: 'benchmark/make_report_{svmeth}.{sample}.tsv'
        container: 'docker://quay.io/jmonlong/svannotate_sveval:0.5'
        shell:
            """
            cp {input.rmd} {params.tmp_rmd}
            Rscript -e 'rmarkdown::render("{params.tmp_rmd}")' -v {input.vcf} -r {input.tr} -rc {input.tr_control} -d {input.db} -a {input.annotsv} -o {output.rds}
            mv {params.tmp_html} {output.html}
            rm -f {params.tmp_rmd}
            """