feat(rule): add BQSR

Author: fcaretti

Diff for: .test/config/config.yml

@@ -3,4 +3,8 @@ data:
 
 reference:
   folder: ".test/data"
-  genome: "placeholder.fa"
+  genome: "placeholder.fa"
+
+known_sites:
+  folder: ".test/data"
+  filename: "placeholder.vcf"

Diff for: .test/data/placeholder.vcf

@@ -1,6 +1,10 @@
 data:
-  folder: "data_folder"
+  folder: "/home/federico/Desktop/RNA_SNPs_calling/data"
 
 reference:
-  folder: "genome_folder"
-  genome: "genome.fa"
+  folder: "/home/federico/Desktop/RNA_SNPs_calling/data/reference"
+  genome: "GRCh38.primary_assembly.genome.fa"
+
+known_sites:
+  folder: "/home/federico/Desktop/RNA_SNPs_calling/data/reference"
+  filename: "resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf"

Diff for: .test/data/placeholder.vcf.idx

@@ -12,21 +12,27 @@ genome_name = config["reference"]["genome"]
 reference = os.path.join(reference_folder, genome_name)
 reference_idx = f"{reference}.fai"
 reference_dict = f"{reference}.dict"
+known_sites_folder = config["known_sites"]["folder"]
+known_filename = config["known_sites"]["filename"]
+known_sites = os.path.join(known_sites_folder, known_filename)
+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]
 
 read_groups = [f"results/grouped/{sample}.bam" for sample in samples]
 deduped_files = [f"results/dedup/{sample}.bam" for sample in samples]
 splitted_files = [f"results/split/{sample}.bam" for sample in samples]
+recalibrated_files = [f"results/recal/{sample}.bam" for sample in samples]
 
 
 rule all:
     input:
-        splitted_files,
+        recalibrated_files,
 
 
 include: "rules/add_or_replace_rg.smk"
 include: "rules/mark_duplicates.smk"
 include: "rules/index_genome.smk"
 include: "rules/split_n_cigar_reads.smk"
+include: "rules/recalibration.smk"

Diff for: config/config.yml

@@ -0,0 +1,37 @@
+rule gatk_baserecalibrator:
+    input:
+        bam="results/split/{sample}.bam",
+        ref=reference,
+        dict=reference_dict,
+        known=known_sites,
+    output:
+        recal_table="results/recal_tables/{sample}.grp",
+    log:
+        "logs/baserecalibrator/{sample}.log",
+    params:
+        extra="",  # optional
+        java_opts="",  # optional
+    resources:
+        mem_mb=1024,
+    wrapper:
+        "v3.12.1/bio/gatk/baserecalibrator"
+
+
+rule gatk_applybqsr:
+    input:
+        bam="results/split/{sample}.bam",
+        ref=reference,
+        dict=reference_dict,
+        recal_table="results/recal_tables/{sample}.grp",
+    output:
+        bam="results/recal/{sample}.bam",
+    log:
+        "logs/gatk_applybqsr/{sample}.log",
+    params:
+        extra="",  # optional
+        java_opts="",  # optional
+        embed_ref=True,  # embed the reference in cram output
+    resources:
+        mem_mb=1024,
+    wrapper:
+        "v3.12.1/bio/gatk/applybqsr"