feat(rule): add GATK's HaplotypeCaller and filtering

Author: fcaretti

Diff for: .test/config/config.yml

@@ -7,4 +7,7 @@ reference:
 
 known_sites:
   folder: ".test/data"
-  filename: "placeholder.vcf"
+  filename: "placeholder.vcf"
+
+filtering:
+  params: ""

Diff for: config/config.yml

@@ -7,4 +7,7 @@ reference:
 
 known_sites:
   folder: "/home/federico/Desktop/RNA_SNPs_calling/data/reference"
-  filename: "resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf"
+  filename: "resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf"
+
+filtering:
+  params: "-i 'QUAL > 30 && INFO/DP > 10'"

Diff for: workflow/Snakefile

@@ -5,13 +5,14 @@ import os
 configfile: "config/config.yml"
 
 
-# Access the data folder from the config
+# Define data and reference files
 data_folder = config["data"]["folder"]
 reference_folder = config["reference"]["folder"]
 genome_name = config["reference"]["genome"]
 reference = os.path.join(reference_folder, genome_name)
 reference_idx = f"{reference}.fai"
 reference_dict = f"{reference}.dict"
+# Used only in FreeBayes
 known_sites_folder = config["known_sites"]["folder"]
 known_filename = config["known_sites"]["filename"]
 known_sites = os.path.join(known_sites_folder, known_filename)
@@ -20,19 +21,23 @@ known_sites_idx = f"{known_sites}.idx"
 sample_files = glob.glob(os.path.join(data_folder, "*.bam"))
 samples = [os.path.basename(f).replace(".bam", "") for f in sample_files]
 
-recalibrated_files = [f"results/recal/{sample}.bam" for sample in samples]
-
 first_summaries = [
     f"results/stats/{sample}_before_recal.summary.txt" for sample in samples
 ]
 second_summaries = [
     f"results/stats/{sample}_after_recal.summary.txt" for sample in samples
 ]
 
+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]
+haplo_calls = "results/calls/calls_gatk.vcf"
+
 
 rule all:
     input:
-        recalibrated_files,
+        "results/calls/filtered_calls.vcf",
         first_summaries,
         second_summaries,
 
@@ -43,3 +48,5 @@ include: "rules/index_genome.smk"
 include: "rules/split_n_cigar_reads.smk"
 include: "rules/recalibration.smk"
 include: "rules/alignment_summary.smk"
+include: "rules/gatk_haplocaller.smk"
+include: "rules/filter.smk"

Diff for: workflow/rules/filter.smk

@@ -0,0 +1,12 @@
+rule bcf_filter_o_vcf:
+    input:
+        haplo_calls,
+    output:
+        "results/calls/filtered_calls.vcf",
+    log:
+        "logs/filter/filter.vcf.log",
+    params:
+        filter=config["filtering"]["params"],
+        extra="",
+    wrapper:
+        "v3.12.1/bio/bcftools/filter"

Diff for: workflow/rules/gatk_haplocaller.smk

@@ -0,0 +1,32 @@
+rule haplotype_caller:
+    input:
+        # single or list of bam files
+        bam="results/recal/{sample}.bam",
+        ref=reference,
+    output:
+        vcf="results/calls_gatk/{sample}.vcf",
+    log:
+        "logs/gatk/haplotypecaller/{sample}.log",
+    params:
+        extra="",  # optional
+        java_opts="",  # optional
+    threads: 4
+    resources:
+        mem_mb=1024,
+    wrapper:
+        "v3.12.1/bio/gatk/haplotypecaller"
+
+
+rule merge_vcfs:
+    input:
+        vcfs=vcfs,
+    output:
+        "results/calls/calls_gatk.vcf",
+    log:
+        "logs/picard/mergevcfs.log",
+    params:
+        extra="",
+    resources:
+        mem_mb=1024,
+    wrapper:
+        "v3.12.1/bio/picard/mergevcfs"