feat: Post qual gnomad filtration

BALSAMIC/commands/run/analysis.py

@@ -3,7 +3,6 @@
 import logging
 import os
 import subprocess
-import shlex
 import sys
 from pathlib import Path
 from typing import List

BALSAMIC/constants/variant_filters.py

@@ -159,7 +159,7 @@ class BaseSNVFilters:
     @classmethod
     def filter_criteria(
         cls,
-        category: Literal["clinical", "research", "quality"],
+        category: Literal["clinical", "research", "post_quality", "quality"],
         analysis_type: Optional[Enum] = None,
         variant_caller: Optional[Enum] = None,
         exclude_variantcaller_filters: Optional[bool] = False,
@@ -169,7 +169,7 @@ def filter_criteria(
         Shared filtering logic to get filters based on criteria.
 
         Args:
-            category (Literal["clinical", "research", "quality"]): The filter category to use.
+            category (Literal["clinical", "research", "post_quality", "quality"]): The filter category to use.
             analysis_type (Optional[Enum]): Filter based on analysis type (default: None).
             variant_caller (Optional[Enum]): Filter based on variant caller (default: None).
             exclude_variantcaller_filters (Optional[bool]): If True, excludes the variantcaller filters.
@@ -208,7 +208,7 @@ def filter_matches(f: VCFFilter) -> bool:
     @classmethod
     def get_bcftools_filter_string(
         cls,
-        category: Literal["clinical", "research", "quality"],
+        category: Literal["clinical", "research", "post_quality", "quality"],
         analysis_type: Optional[Enum] = None,
         variant_caller: Optional[Enum] = None,
         soft_filter_normals: Optional[bool] = None,
@@ -218,7 +218,7 @@ def get_bcftools_filter_string(
         Get a set of filter names based on various attributes.
 
         Args:
-            category (Literal["clinical", "research", "quality"]): The filter category to use.
+            category (Literal["clinical", "research", "post_quality", "quality"]): The filter category to use.
             analysis_type (Optional[Enum]): Filter based on analysis type (default: None).
             variant_caller (Optional[Enum]): Filter based on variant caller (default: None).
             soft_filter_normals (Optional[bool]): If True, removes filters in MATCHED_NORMAL_FILTER_NAMES.
@@ -256,7 +256,7 @@ def get_bcftools_filter_string(
     @classmethod
     def get_filters(
         cls,
-        category: Literal["clinical", "research", "quality"],
+        category: Literal["clinical", "research", "post_quality", "quality"],
         analysis_type: Optional[Enum] = None,
         variant_caller: Optional[Enum] = None,
         exclude_variantcaller_filters: Optional[bool] = True,
@@ -266,7 +266,7 @@ def get_filters(
         Get a list of filters matching the specified attributes.
 
         Args:
-            category (Literal["clinical", "research", "quality"]): The filter category to use.
+            category (Literal["clinical", "research",, "post_quality", "quality"]): The filter category to use.
             analysis_type (Optional[Enum]): Filter based on analysis type (default: None).
             variant_caller (Optional[Enum]): Filter based on variant caller (default: None).
             exclude_variantcaller_filters (Optional[bool]): If True, excludes the variantcaller filters.
@@ -285,12 +285,15 @@ def get_filters(
 
 
 class WgsSNVFilters(BaseSNVFilters):
+    post_quality = [
+        VCFFilter(tag_value=0.1, filter_name="gnomad_hard_threshold", field="INFO"),
+    ]
     research = [
         VCFFilter(tag_value=0.01, filter_name="SWEGENAF", field="INFO"),
         VCFFilter(tag_value=0.001, filter_name="balsamic_high_pop_freq", field="INFO"),
     ]
     clinical = research + [
-        VCFFilter(tag_value=0.007, filter_name="Frq", field="INFO"),
+        VCFFilter(tag_value=0.005, filter_name="Frq", field="INFO"),
         VCFFilter(tag_value=0.1, filter_name="ArtefactFrq", field="INFO"),
     ]
     quality = [
@@ -331,15 +334,18 @@ class WgsSNVFilters(BaseSNVFilters):
 
 
 class TgaSNVFilters(BaseSNVFilters):
+    post_quality = [
+        VCFFilter(tag_value=0.1, filter_name="gnomad_hard_threshold", field="INFO"),
+    ]
     research = [
         VCFFilter(
             filter_name="MERGED", Description="SNV Merged with neighboring variants"
         ),
         VCFFilter(tag_value=0.01, filter_name="SWEGENAF", field="INFO"),
-        VCFFilter(tag_value=0.005, filter_name="balsamic_high_pop_freq", field="INFO"),
+        VCFFilter(tag_value=0.001, filter_name="balsamic_high_pop_freq", field="INFO"),
     ]
     clinical = research + [
-        VCFFilter(tag_value=0.01, filter_name="Frq", field="INFO"),
+        VCFFilter(tag_value=0.005, filter_name="Frq", field="INFO"),
         VCFFilter(tag_value=0.1, filter_name="ArtefactFrq", field="INFO"),
         VCFFilter(tag_value=0.3, filter_name="HighOccurrenceFrq", field="INFO"),
     ]
@@ -403,10 +409,10 @@ class TgaSNVFilters(BaseSNVFilters):
 class TgaUmiSNVFilters(BaseSNVFilters):
     research = [
         VCFFilter(tag_value=0.01, filter_name="SWEGENAF", field="INFO"),
-        VCFFilter(tag_value=0.02, filter_name="balsamic_high_pop_freq", field="INFO"),
+        VCFFilter(tag_value=0.001, filter_name="balsamic_high_pop_freq", field="INFO"),
     ]
     clinical = research + [
-        VCFFilter(tag_value=0.01, filter_name="Frq", field="INFO"),
+        VCFFilter(tag_value=0.005, filter_name="Frq", field="INFO"),
         VCFFilter(tag_value=0.1, filter_name="ArtefactFrq", field="INFO"),
     ]
     quality = [

BALSAMIC/snakemake_rules/annotation/somatic_snv_annotation.rule

@@ -2,9 +2,37 @@
 # coding: utf-8
 # Somatic SNV annotation module.
 
+rule snv_gnomad_annotation_and_filtration:
+  input:
+    vcf_snv_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research.vcf.gz",
+  output:
+    research_gnomadfiltered = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research_gnomadfiltered.vcf.gz",
+    vcfanno_gnomad_annotation_toml = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.vcfanno_gnomadsnv.toml"
+  benchmark:
+    Path(benchmark_dir, "snv_gnomad_annotation_and_filtration.{case_name}.{var_caller}.tsv").as_posix()
+  singularity:
+    Path(singularity_image, config["bioinfo_tools"].get("ensembl-vep") + ".sif").as_posix()
+  params:
+    pop_freq = get_tag_and_filtername(snv_post_quality_filters, "gnomad_hard_threshold"),
+    gnomad_annotation = config_model.retrieve_annotations_toml(category=AnnotationCategory.POST_QUALITY),
+    tmpdir = tempfile.mkdtemp(prefix=tmp_dir)
+  message:
+    "Annotating and filtering based on gnomad frequencies on SNVs and INDELs"
+  shell:
+        """
+echo \'{params.gnomad_annotation}\' > {output.vcfanno_gnomad_annotation_toml};
+
+vcfanno -p {threads} {output.vcfanno_gnomad_annotation_toml} {input.vcf_snv_research} \
+| bcftools view --threads {threads} -O z -o {params.tmpdir}/SNV.somatic.research.gnomad_annotated.vcf.gz ;
+
+bcftools filter --threads {threads} --include 'INFO/GNOMADAF_popmax <= {params.pop_freq[0]} || INFO/GNOMADAF_popmax == \".\"' -O z -o {output.research_gnomadfiltered} {params.tmpdir}/SNV.somatic.research.gnomad_annotated.vcf.gz
+
+tabix -p vcf -f {output.research_gnomadfiltered};
+        """
+
 rule bcftools_get_somaticINDEL_research:
   input:
-    vcf_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research.vcf.gz",
+    vcf_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research_gnomadfiltered.vcf.gz",
   output:
     vcf_indel_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.indel.research.vcf.gz",
   benchmark:
@@ -43,7 +71,7 @@ CADD.sh -g GRCh37 -o {output.cadd_indel_research} {input.vcf_indel_research}
 
 rule bcftools_annotate_somaticINDEL_research:
   input:
-    vcf_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research.vcf.gz",
+    vcf_research = vcf_dir + "SNV.somatic.{case_name}.{var_caller}.research_gnomadfiltered.vcf.gz",
     cadd_indel_research = vep_dir + "SNV.somatic.{case_name}.{var_caller}.cadd_indel.research.tsv.gz",
   output:
     vcf_indel_research = vep_dir + "SNV.somatic.{case_name}.{var_caller}.cadd_indel.research.vcf.gz",

BALSAMIC/snakemake_rules/annotation/somatic_sv_annotation.rule

@@ -56,7 +56,7 @@ rule svdb_annotate_somatic_research_sv:
     shell:
         """
 if [[ -f "{params.swegen_sv_frequency}" ]]; then
-    svdb --query --bnd_distance 5000 --overlap 0.80 \
+    svdb --query --bnd_distance 1000 --overlap 0.80 \
     --in_occ OCC --out_occ swegen_obs --in_frq FRQ --out_frq SWEGENAF \
     --db {params.swegen_sv_frequency} --query_vcf {input.vcf_sv_research_vep} > {params.tmp_vcf};
     bgzip -l 9 -c {params.tmp_vcf} > {output.vcf_research};
@@ -86,7 +86,7 @@ rule svdb_annotate_clinical_obs_somatic_clinical_sv:
     shell:
         """
 if [[ -f "{params.clinical_sv_observations}" ]]; then
-    svdb --query --bnd_distance 10000 --overlap 0.80 \
+    svdb --query --bnd_distance 1000 --overlap 0.80 \
     --in_occ Obs --out_occ clin_obs --in_frq Frq --out_frq Frq \
     --db {params.clinical_sv_observations} --query_vcf {input.vcf_sv_research} > {params.vcf_clinical_obs}
     bgzip -l 9 -c {params.vcf_clinical_obs} > {output.vcf_sv_clinical_obs};
@@ -118,7 +118,7 @@ rule svdb_annotate_somatic_obs_somatic_clinical_sv:
     shell:
         """
 if [[ -f "{params.somatic_sv_observations}" ]]; then
-    svdb --query --bnd_distance 10000 --overlap 0.80 \
+    svdb --query --bnd_distance 1000 --overlap 0.80 \
     --in_occ Obs --out_occ Cancer_Somatic_Obs --in_frq Frq --out_frq Cancer_Somatic_Frq \
     --db {params.somatic_sv_observations} --query_vcf {input.vcf_sv_clinical_obs} > {params.vcf_somatic_obs}
     bgzip -l 9 -c {params.vcf_somatic_obs} > {output.vcf_sv_clinical};

BALSAMIC/snakemake_rules/annotation/varcaller_filter_tumor_normal_wgs.rule

@@ -69,7 +69,5 @@ python {params.rescue_script} --vcf {params.tmpdir}/clinical.filtered.pass_befor
 
 bcftools view --threads {threads} --exclude '{params.hard_filters}' -O z -o {output.vcf_pass_tnscope} {params.tmpdir}/clinical.filtered.pass_after_rescue.vcf;
 
-tabix -p vcf -f {output.vcf_pass_tnscope};
-
 bcftools +counts {output.vcf_pass_tnscope} > {output.bcftools_counts_clinical};
     """

BALSAMIC/workflows/balsamic.smk

@@ -128,10 +128,12 @@ else:
 if sequencing_type == "targeted":
     exome = config_model.panel.exome
     snv_quality_filters = SNV_FILTERS.get_filters(category="quality", analysis_type=analysis_type, exome=exome)
+    snv_post_quality_filters = SNV_FILTERS.get_filters(category="post_quality", analysis_type=analysis_type, exome=exome)
     snv_research_filters = SNV_FILTERS.get_filters(category="research", analysis_type=analysis_type, exome=exome)
     snv_clinical_filters = SNV_FILTERS.get_filters(category="clinical", analysis_type=analysis_type, exome=exome)
 else:
     snv_quality_filters = SNV_FILTERS.get_filters(category="quality", analysis_type=analysis_type)
+    snv_post_quality_filters = SNV_FILTERS.get_filters(category="post_quality",analysis_type=analysis_type)
     snv_research_filters = SNV_FILTERS.get_filters(category="research", analysis_type=analysis_type)
     snv_clinical_filters = SNV_FILTERS.get_filters(category="clinical", analysis_type=analysis_type)
 

CHANGELOG.rst

@@ -65,6 +65,7 @@ Fixed:
 ^^^^^^
 
 
+
 [17.0.1]
 --------