Merge pull request #9 from fcaretti/Add-rules

feat(rule): GATK's HaplotypeCaller

config/config.yml

@@ -1,13 +1,20 @@
 data:
-  folder: "/home/federico/Desktop/RNA_SNPs_calling/data"
+  folder: "data_folder"
 
 reference:
-  folder: "/home/federico/Desktop/RNA_SNPs_calling/data/reference"
-  genome: "GRCh38.primary_assembly.genome.fa"
+  folder: "reference_folder"
+  genome: "genome.fa"
 
 known_sites:
-  folder: "/home/federico/Desktop/RNA_SNPs_calling/data/reference"
-  filename: "resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf"
+  folder: "known_sites_folder"
+  filename: "ex: resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf"
 
 filtering:
-  params: "-i 'QUAL > 30 && INFO/DP > 10'"
+  params: "ex: -i 'QUAL > 30 && INFO/DP > 10'"
+
+vep:
+  cache_dir: "cache_folder"
+  plugins_dir: "plugins_folder"
+  species: ""
+  build: ""
+  release: ""

workflow/Snakefile

@@ -31,6 +31,7 @@ second_summaries = [
 calls = [f"results/calls/{sample}.vcf" for sample in samples]
 alns = [f"results/recal/{sample}.bam" for sample in samples]
 idxs = [f"results/recal/{sample}.bai" for sample in samples]
+vcfs = [f"results/calls_gatk/{sample}.vcf" for sample in samples]
 vcf_zips = [f"results/calls_gatk/{sample}.vcf.gz" for sample in samples]
 vcf_idxs = [f"results/calls_gatk/{sample}.vcf.csi" for sample in samples]
 haplo_calls = "results/calls/calls_gatk.vcf"

workflow/rules/gatk_haplocaller.smk

@@ -17,14 +17,13 @@ rule haplotype_caller:
         "v3.12.1/bio/gatk/haplotypecaller"
 
 
-
 rule bgzip:
     input:
         "results/calls_gatk/{sample}.vcf",
     output:
         temp("results/calls_gatk/{sample}.vcf.gz"),
     params:
-        extra="", # optional
+        extra="",  # optional
     threads: 1
     log:
         "logs/bgzip/{sample}.log",

workflow/rules/vep.smk

@@ -0,0 +1,43 @@
+rule annotate_variants:
+    input:
+        calls="results/calls/calls_gatk.vcf",  # .vcf, .vcf.gz or .bcf
+        cache=config["vep"]["cache_dir"],  # can be omitted if fasta and gff are specified
+        plugins=config["vep"]["plugins_dir"],
+        fasta=reference,
+        fai=reference_idx,
+    output:
+        calls="results/calls/annotated_calls.vcf",  # .vcf, .vcf.gz or .bcf
+        stats="results/calls/variants.html",
+    params:
+        # Pass a list of plugins to use, see https://www.ensembl.org/info/docs/tools/vep/script/vep_plugins.html
+        # Plugin args can be added as well, e.g. via an entry "MyPlugin,1,FOO", see docs.
+        plugins=["LoFtool"],
+        extra="--everything",  # optional: extra arguments
+    log:
+        "logs/vep/annotate.log",
+    threads: 4
+    wrapper:
+        "v3.12.1/bio/vep/annotate"
+
+
+rule get_vep_cache:
+    output:
+        directory(config["vep"]["cache_dir"]),
+    params:
+        species=config["vep"]["species"],
+        build=config["vep"]["build"],
+        release=config["vep"]["release"],
+    log:
+        "logs/vep/cache.log",
+    cache: "omit-software"  # save space and time with between workflow caching (see docs)
+    wrapper:
+        "v3.12.1/bio/vep/cache"
+
+
+rule download_vep_plugins:
+    output:
+        temp(directory(config["vep"]["plugins_dir"])),
+    params:
+        release=config["vep"]["release"],
+    wrapper:
+        "v3.12.1/bio/vep/plugins"