Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
155 changes: 85 additions & 70 deletions workflow/rules/vaf.smk
Original file line number Diff line number Diff line change
@@ -1,23 +1,25 @@
rule snps_to_ancestor:
threads: 2
retries: 3
shadow: "minimal"
conda: "../envs/var_calling.yaml"
shadow:
"minimal"
conda:
"../envs/var_calling.yaml"
params:
mpileup_depth = config["VC"]["MAX_DEPTH"],
mpileup_quality = 0,
ivar_quality = config["VC"]["MIN_QUALITY"],
ivar_freq = config["VC"]["MIN_FREQ"],
ivar_depth = config["VC"]["MIN_DEPTH"],
mpileup_depth=config["VC"]["MAX_DEPTH"],
mpileup_quality=0,
ivar_quality=config["VC"]["MIN_QUALITY"],
ivar_freq=config["VC"]["MIN_FREQ"],
ivar_depth=config["VC"]["MIN_DEPTH"],
input:
reference_fasta = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
bam = get_input_bam,
gff = OUTDIR/"reference.gff3"
reference_fasta=OUTDIR / f"{OUTPUT_NAME}.ancestor.fasta",
bam=get_input_bam,
gff=OUTDIR / "reference.gff3",
output:
tsv = temp(OUTDIR/"vaf"/"{sample}.tsv"),
reference_fasta_renamed = temp(OUTDIR/"vaf"/"{sample}.reference.fasta"),
tsv=temp(OUTDIR / "vaf" / "vc" / "{sample}.tsv"),
reference_fasta_renamed=temp(OUTDIR / "vaf" / "{sample}.reference.fasta"),
log:
LOGDIR / "snps_to_ancestor" / "{sample}.log.txt"
LOGDIR / "snps_to_ancestor" / "{sample}.log.txt",
shell:
"""
set -e
Expand Down Expand Up @@ -55,65 +57,70 @@ rule snps_to_ancestor:

rule mask_tsv:
threads: 1
conda: "../envs/biopython.yaml"
conda:
"../envs/biopython.yaml"
params:
mask_class = ["mask"]
input:
tsv = OUTDIR/"vaf"/"{sample}.tsv",
vcf = lambda wildcards: select_problematic_vcf()
mask_class=["mask"],
input:
tsv=OUTDIR / "vaf" / "vc" / "{sample}.tsv",
vcf=lambda wildcards: select_problematic_vcf(),
output:
masked_tsv = temp(OUTDIR/"vaf"/"{sample}.masked.tsv")
masked_tsv=temp(OUTDIR / "vaf" / "masked" / "{sample}.tsv"),
log:
LOGDIR / "mask_tsv" / "{sample}.log.txt"
LOGDIR / "mask_tsv" / "{sample}.log.txt",
script:
"../scripts/mask_tsv.py"


rule filter_tsv:
threads: 1
conda: "../envs/renv.yaml"
conda:
"../envs/renv.yaml"
params:
min_depth = 20,
min_alt_rv = 2,
min_alt_dp = 2,
input:
tsv = OUTDIR/"vaf"/"{sample}.masked.tsv"
min_depth=20,
min_alt_rv=2,
min_alt_dp=2,
input:
tsv=OUTDIR / "vaf" / "masked" / "{sample}.tsv",
output:
filtered_tsv = temp(OUTDIR/"vaf"/"{sample}.masked.prefiltered.tsv")
filtered_tsv=temp(OUTDIR / "vaf" / "filtered" / "{sample}.tsv"),
log:
LOGDIR / "filter_tsv" / "{sample}.log.txt"
LOGDIR / "filter_tsv" / "{sample}.log.txt",
script:
"../scripts/filter_tsv.R"


rule tsv_to_vcf:
threads: 1
conda: "../envs/biopython.yaml"
conda:
"../envs/biopython.yaml"
params:
ref_name = config["ALIGNMENT_REFERENCE"],
input:
tsv = OUTDIR/"vaf"/"{sample}.masked.prefiltered.tsv",
ref_name=config["ALIGNMENT_REFERENCE"],
input:
tsv=OUTDIR / "vaf" / "filtered" / "{sample}.tsv",
output:
vcf = temp(OUTDIR/"vaf"/"{sample}.vcf")
vcf=temp(OUTDIR / "vaf" / "vcf" / "{sample}.vcf"),
log:
LOGDIR / "tsv_to_vcf" / "{sample}.log.txt"
LOGDIR / "tsv_to_vcf" / "{sample}.log.txt",
script:
"../scripts/tsv_to_vcf.py"


rule variants_effect:
threads: 1
shadow: "minimal"
conda: "../envs/snpeff.yaml"
shadow:
"minimal"
conda:
"../envs/snpeff.yaml"
params:
ref_name = config["ALIGNMENT_REFERENCE"],
snpeff_data_dir = (BASE_PATH / "config" / "snpeff").resolve()
ref_name=config["ALIGNMENT_REFERENCE"],
snpeff_data_dir=(BASE_PATH / "config" / "snpeff").resolve(),
input:
vcf = OUTDIR/"vaf"/"{sample}.vcf"
vcf=OUTDIR / "vaf" / "vcf" / "{sample}.vcf",
output:
ann_vcf = OUTDIR/"vaf"/"{sample}.annotated.vcf"
ann_vcf=OUTDIR / "vaf" / "annotated" / "{sample}.vcf",
log:
LOGDIR / "variants_effect" / "{sample}.log.txt"
LOGDIR / "variants_effect" / "{sample}.log.txt",
retries: 2
shell:
"""
Expand All @@ -133,74 +140,82 @@ rule variants_effect:

rule extract_vcf_fields:
threads: 1
conda: "../envs/snpeff.yaml"
conda:
"../envs/snpeff.yaml"
params:
extract_columns = [f"'{col}'" for col in config["ANNOTATION"]["SNPEFF_COLS"].values()],
sep = ",",
extract_columns=[
f"'{col}'" for col in config["ANNOTATION"]["SNPEFF_COLS"].values()
],
sep=",",
input:
vcf = OUTDIR/"vaf"/"{sample}.annotated.vcf"
vcf=OUTDIR / "vaf" / "annotated" / "{sample}.vcf",
output:
tsv = OUTDIR/"vaf"/"{sample}.vcf_fields.tsv"
tsv=OUTDIR / "vaf" / "fields" / "{sample}.tsv",
log:
LOGDIR / "tsv_to_vcf" / "{sample}.log.txt"
LOGDIR / "tsv_to_vcf" / "{sample}.log.txt",
shell:
"SnpSift extractFields -e 'NA' -s {params.sep:q} {input.vcf:q} {params.extract_columns} >{output.tsv:q} 2>{log:q}"


rule format_vcf_fields_longer:
conda: "../envs/renv.yaml"
conda:
"../envs/renv.yaml"
params:
sample = "{sample}",
colnames_mapping = config["ANNOTATION"]["SNPEFF_COLS"],
filter_include = config["ANNOTATION"]["FILTER_INCLUDE"],
filter_exclude = config["ANNOTATION"]["FILTER_EXCLUDE"],
variant_name_pattern = lambda wildcards: config["ANNOTATION"]["VARIANT_NAME_PATTERN"], # lambda to deactivate automatic wildcard expansion in pattern
sep = ",",
sample="{sample}",
colnames_mapping=config["ANNOTATION"]["SNPEFF_COLS"],
filter_include=config["ANNOTATION"]["FILTER_INCLUDE"],
filter_exclude=config["ANNOTATION"]["FILTER_EXCLUDE"],
variant_name_pattern=lambda wildcards: config["ANNOTATION"][
"VARIANT_NAME_PATTERN"
], # lambda to deactivate automatic wildcard expansion in pattern
sep=",",
input:
tsv = OUTDIR/"vaf"/"{sample}.vcf_fields.tsv",
tsv=OUTDIR / "vaf" / "fields" / "{sample}.tsv",
output:
tsv = OUTDIR/"vaf"/"{sample}.vcf_fields.longer.tsv",
tsv=OUTDIR / "vaf" / "fields_longer" / "{sample}.tsv",
log:
LOGDIR / "format_vcf_fields_longer" / "{sample}.log.txt"
LOGDIR / "format_vcf_fields_longer" / "{sample}.log.txt",
script:
"../scripts/format_vcf_fields_longer.R"


rule concat_vcf_fields:
params:
sep = "\t",
sep="\t",
input:
expand(OUTDIR/"vaf"/"{sample}.vcf_fields.longer.tsv", sample=iter_samples()),
expand(OUTDIR / "vaf" / "fields_longer" / "{sample}.tsv", sample=iter_samples()),
output:
OUTDIR/f"{OUTPUT_NAME}.vcf_fields.longer.tsv",
OUTDIR / f"{OUTPUT_NAME}.vcf_fields.longer.tsv",
run:
import pandas as pd
from functools import reduce

reduce(
lambda a, b: pd.concat((a, b), axis="rows", ignore_index=True),
(pd.read_csv(path, sep=params.sep) for path in input)
(pd.read_csv(path, sep=params.sep) for path in input),
).to_csv(output[0], sep=params.sep, index=False)


rule merge_annotation:
threads: 1
conda: "../envs/renv.yaml"
conda:
"../envs/renv.yaml"
params:
sample = "{sample}",
ref_name = config["ALIGNMENT_REFERENCE"],
sample="{sample}",
ref_name=config["ALIGNMENT_REFERENCE"],
input:
tsv = OUTDIR/"vaf"/"{sample}.masked.prefiltered.tsv",
annot = OUTDIR/"vaf"/"{sample}.vcf_fields.longer.tsv",
tsv=OUTDIR / "vaf" / "filtered" / "{sample}.tsv",
annot=OUTDIR / "vaf" / "fields_longer" / "{sample}.tsv",
output:
tsv = OUTDIR/"vaf"/"{sample}.variants.tsv"
tsv=OUTDIR / "vaf" / "variants" / "{sample}.tsv",
log:
LOGDIR / "merge_annotation" / "{sample}.log.txt"
LOGDIR / "merge_annotation" / "{sample}.log.txt",
script:
"../scripts/merge_annotation.R"


use rule concat_vcf_fields as concat_variants with:
input:
expand(OUTDIR/"vaf"/"{sample}.variants.tsv", sample=iter_samples()),
expand(OUTDIR / "vaf" / "variants" / "{sample}.tsv", sample=iter_samples()),
output:
OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
OUTDIR / f"{OUTPUT_NAME}.variants.tsv",
Loading