Merge pull request #281 from sigven/read_support_patch

Read support patch

Author: sigven

README.md

@@ -28,6 +28,12 @@ PCGR originates from the [Norwegian Cancer Genomics Consortium (NCGC)](https://c
 
 ### Top News
 
+- *September 17th 2025:* **2.2.5 release**
+  - fixing missing support for dp/af filtering, adding ad filtering
+  - skip processing when no PASS variants detected in input VCF
+  - more streamlined plotting functions in quarto report templates
+  - [CHANGELOG](https://sigven.github.io/pcgr/articles/CHANGELOG.html)
+
 - *September 8th 2025:* **2.2.4 release**
   - various minor bug fixes, addition of `--sex` option for sex-adjusted CNA annotation
   - [CHANGELOG](https://sigven.github.io/pcgr/articles/CHANGELOG.html)
@@ -87,7 +93,7 @@ PCGR originates from the [Norwegian Cancer Genomics Consortium (NCGC)](https://c
 
 ### Example reports
 
-[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.15068347.svg)](https://doi.org/10.5281/zenodo.15068347)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17140659.svg)](https://doi.org/10.5281/zenodo.17140659)
 
 ### Why use PCGR?
 

pcgr/annoutils.py

@@ -100,6 +100,7 @@ def write_pass_vcf(annotated_vcf, logger):
     """
     #out_vcf = re.sub(r'\.annotated\.vcf\.gz$','.annotated.pass.vcf',annotated_vcf)
     out_vcf = re.sub(r'\.vcf\.gz$', '.pass.vcf', annotated_vcf)
+
     vcf = VCF(annotated_vcf)
     w = Writer(out_vcf, vcf)
 
@@ -118,16 +119,21 @@ def write_pass_vcf(annotated_vcf, logger):
     logger.info('Number of non-PASS/REJECTED variant calls: ' +
                 str(num_rejected))
     logger.info('Number of PASSed variant calls: ' + str(num_pass))
+
+    
+    vcf_no_pass_variants = False
     if num_pass == 0:
+        vcf_no_pass_variants = True
         logger.warning(
             'There are zero variants with a \'PASS\' filter in the VCF file')
-        os.system('bgzip -dc ' + str(annotated_vcf) +
-                  ' egrep \'^#\' > ' + str(out_vcf))
+        os.system(f"bgzip -dc {annotated_vcf} | egrep '^#' >  {out_vcf}")
     # else:
-    os.system('bgzip -f ' + str(out_vcf))
-    os.system('tabix -f -p vcf ' + str(out_vcf) + '.gz')
+    os.system(f'bgzip -f {out_vcf}')
+    os.system(f'tabix -f -p vcf {out_vcf}.gz')
+
+    #exit(-1)
 
-    return
+    return(vcf_no_pass_variants)
 
 
 def map_regulatory_variant_annotations(vep_csq_records):

pcgr/arg_checker.py

@@ -40,51 +40,118 @@ def verify_args(arg_dict):
 
     # check that tumor purity and tumor ploidy is set correctly
     if arg_dict['tumor_purity'] is not None:
-        if not (arg_dict['tumor_purity'] > 0 and arg_dict['tumor_purity'] <= 1):
+        if not (float(arg_dict['tumor_purity']) > 0 and float(arg_dict['tumor_purity']) <= 1):
             err_msg = f"Tumor purity value ('--tumor_purity' = {arg_dict['tumor_purity']}) must be within (0, 1]"
             error_message(err_msg, logger)
 
     if arg_dict['tumor_ploidy'] is not None:
-        if not arg_dict['tumor_ploidy'] > 0:
+        if not float(arg_dict['tumor_ploidy']) > 0:
             err_msg = f"Tumor ploidy value ('--tumor_ploidy' = {arg_dict['tumor_ploidy']}) must be > 0"
             error_message(err_msg, logger)
 
-    # check that minimum/maximum depth/allelic fractions are set correctly
-    if int(arg_dict['tumor_dp_min']) < 0:
-        err_msg = f"Minimum sequencing depth tumor ('tumor_dp_min' = {arg_dict['tumor_dp_min']}) must be >= 0"
+    ## check that allelic support tags are set correctly when minimum/maximum depth/allelic fractions are set
+    if arg_dict['tumor_dp_tag'] == "_NA_" and arg_dict['tumor_dp_min'] is not None:
+        err_msg = f"Minimum sequencing depth tumor ('tumor_dp_min' = {arg_dict['tumor_dp_min']}) requires '--tumor_dp_tag' to be set"
         error_message(err_msg, logger)
-
-    if float(arg_dict['tumor_af_min']) < 0 or float(arg_dict['tumor_af_min']) > 1:
-        err_msg = f"Minimum AF tumor ('tumor_af_min' = {arg_dict['tumor_af_min']}) must be within [0, 1]"
+    
+    if arg_dict['tumor_af_tag'] == "_NA_" and (arg_dict['tumor_af_min'] is not None or arg_dict['tumor_ad_min'] is not None):
+        err_msg = f"Minimum AF/AD tumor ('tumor_af_min' = {arg_dict['tumor_af_min']}, 'tumor_ad_min' = {arg_dict['tumor_ad_min']}) requires '--tumor_af_tag' to be set"
         error_message(err_msg, logger)
-
-    if int(arg_dict['control_dp_min']) < 0:
-        err_msg = f"Minimum sequencing depth control ('control_dp_min' = {arg_dict['control_dp_min']}) must be >= 0"
+    
+    if arg_dict['control_dp_tag'] == "_NA_" and arg_dict['control_dp_min'] is not None:
+        err_msg = f"Minimum sequencing depth control ('control_dp_min' = {arg_dict['control_dp_min']}) requires '--control_dp_tag' to be set"
         error_message(err_msg, logger)
-
-    if float(arg_dict['control_af_max']) < 0 or float(arg_dict['control_af_max']) > 1:
-        err_msg = f"Maximum AF control ('control_af_max' = {arg_dict['control_af_max']}) must be within [0, 1]"
+    
+    if arg_dict['control_af_tag'] == "_NA_" and (arg_dict['control_af_max'] is not None or arg_dict['control_ad_max'] is not None):
+        err_msg = f"Maximum AF/AD control ('control_af_max' = {arg_dict['control_af_max']}, 'control_ad_max' = {arg_dict['control_ad_max']}) requires '--control_af_tag' to be set"
         error_message(err_msg, logger)
+
+    # check that minimum/maximum depth/allelic fractions are set correctly
+    dp_tumor_set = False
+    if arg_dict['tumor_dp_min'] is not None:
+        dp_tumor_set = True
+        if int(arg_dict['tumor_dp_min']) <= 0:
+            err_msg = f"Minimum sequencing depth tumor ('tumor_dp_min' = {arg_dict['tumor_dp_min']}) must be > 0"
+            error_message(err_msg, logger)
+
+    if arg_dict['tumor_af_min'] is not None:
+        if (float(arg_dict['tumor_af_min']) > 1 or float(arg_dict['tumor_af_min']) < 0):
+            err_msg = f"Minimum AF tumor ('tumor_af_min' = {arg_dict['tumor_af_min']}) must be within [0, 1]"
+            error_message(err_msg, logger)
+
+    dp_control_set = False
+    if arg_dict['control_dp_min'] is not None:
+        dp_control_set = True
+        if int(arg_dict['control_dp_min']) <= 0:            
+            err_msg = f"Minimum sequencing depth control ('control_dp_min' = {arg_dict['control_dp_min']}) must be > 0"
+            error_message(err_msg, logger)
+
+    if arg_dict['tumor_ad_min'] is not None:
+        err = 0
+        if int(arg_dict['tumor_ad_min']) <= 0:
+            err = 1
+        if dp_tumor_set is True and int(arg_dict['tumor_ad_min']) > int(arg_dict['tumor_dp_min']):
+            err = 1
+        if err == 1:            
+            err_msg = (
+                f"Minimum allelic depth tumor - ('tumor_ad_min' = {arg_dict['tumor_ad_min']}) must be > 0 "
+                f"and less than or equal to minimum sequencing depth tumor ('tumor_dp_min' = {arg_dict['tumor_dp_min']})"
+            )
+            error_message(err_msg, logger)
+
+    if arg_dict['control_af_max'] is not None:
+        if float(arg_dict['control_af_max']) < 0 or float(arg_dict['control_af_max']) > 1:
+            err_msg = f"Maximum AF control ('control_af_max' = {arg_dict['control_af_max']}) must be within [0, 1]"
+            error_message(err_msg, logger)
 
+    if arg_dict['control_ad_max'] is not None:
+        err = 0
+        if int(arg_dict['control_ad_max']) < 0:
+            err = 1
+        if dp_control_set is True and int(arg_dict['control_ad_max']) > int(arg_dict['control_dp_min']):
+            err = 1
+        if err == 1:
+            err_msg = (
+                f"Maximum allelic depth control - ('control_ad_max' = {arg_dict['control_ad_max']}) must be >= 0 "
+                f"and less than or equal to minimum sequencing depth control ('control_dp_min' = {arg_dict['control_dp_min']})"
+            )
+            error_message(err_msg, logger)
 
-     # TMB: check that minimum/maximum depth/allelic fractions are set correctly
-    if int(arg_dict['tmb_dp_min']) < 0:
-        err_msg = f"Minimum sequencing depth tumor - TMB calculation ('tmb_dp_min' = {arg_dict['tmb_dp_min']}) must be >= 0"
-        error_message(err_msg, logger)
-        
-    if int(arg_dict['tmb_dp_min']) > 0 and (int(arg_dict['tmb_dp_min']) < int(arg_dict['tumor_dp_min'])):
-        err_msg = f"Minimum sequencing depth (tumor) for TMB calculation ('tmb_dp_min' = {str(arg_dict['tmb_dp_min'])}) must be "
-        err_msg += f"greater or equal to minimum sequencing depth tumor {str(arg_dict['tumor_dp_min'])} (i.e. filter for variant inclusion in report)"
-        error_message(err_msg, logger)
+    # TMB: check that minimum/maximum depth/allelic fractions are set correctly
+    if arg_dict['tmb_dp_min'] is not None:
+        if int(arg_dict['tmb_dp_min']) <= 0:
+            err_msg = f"Minimum sequencing depth tumor - TMB calculation ('tmb_dp_min' = {arg_dict['tmb_dp_min']}) must be > 0"
+            error_message(err_msg, logger)
+        if dp_tumor_set is True and (int(arg_dict['tmb_dp_min']) < int(arg_dict['tumor_dp_min'])):
+            err_msg = f"Minimum sequencing depth (tumor) for TMB calculation ('tmb_dp_min' = {str(arg_dict['tmb_dp_min'])}) must be "
+            err_msg += f"greater or equal to minimum sequencing depth tumor {str(arg_dict['tumor_dp_min'])} (i.e. global filter for variant inclusion in report)"
+            error_message(err_msg, logger)
+    
+    if arg_dict['tmb_ad_min'] is not None:
+        if int(arg_dict['tmb_ad_min']) <= 0:
+            err_msg = f"Minimum allelic depth tumor - TMB calculation ('tmb_ad_min' = {arg_dict['tmb_ad_min']}) must be > 0"
+            error_message(err_msg, logger)
+        if arg_dict['tumor_ad_min'] is not None:
+            if int(arg_dict['tmb_ad_min']) > 0 and int(arg_dict['tumor_ad_min']) > 0:
+                if int(arg_dict['tmb_ad_min']) < int(arg_dict['tumor_ad_min']):
+                    err_msg = f"Minimum allelic depth (tumor) for TMB calculation ('tmb_ad_min' = {str(arg_dict['tmb_ad_min'])}) must be "
+                    err_msg += f"greater or equal to minimum allelic depth tumor {str(arg_dict['tumor_ad_min'])} (i.e. global filter for variant inclusion in report)"
+                    error_message(err_msg, logger)
+                if dp_tumor_set is True and (int(arg_dict['tmb_ad_min']) > int(arg_dict['tmb_dp_min']) and int(arg_dict['tmb_dp_min']) > 0):
+                    err_msg = f"Minimum allelic depth (tumor) for TMB calculation ('tmb_ad_min' = {str(arg_dict['tmb_ad_min'])}) must be "
+                    err_msg += f"less than or equal to minimum sequencing depth (tumor) for TMB calculation ('tmb_dp_min' = {str(arg_dict['tmb_dp_min'])})"
+                    error_message(err_msg, logger)
 
-    if float(arg_dict['tmb_af_min']) < 0 or float(arg_dict['tmb_af_min']) > 1:
-        err_msg = f"Minimum AF (tumor) for TMB calculation ('tmb_af_min' = {arg_dict['tmb_af_min']}) must be within [0, 1]"
-        error_message(err_msg, logger)
+    if arg_dict['tmb_af_min'] is not None:
+        if float(arg_dict['tmb_af_min']) < 0 or float(arg_dict['tmb_af_min']) > 1:
+            err_msg = f"Minimum AF (tumor) for TMB calculation ('tmb_af_min' = {arg_dict['tmb_af_min']}) must be within [0, 1]"
+            error_message(err_msg, logger)
 
-    if float(arg_dict['tmb_af_min']) > 0 and (float(arg_dict['tmb_af_min']) < float(arg_dict['tumor_af_min'])):
-        err_msg = f"Minimum AF (tumor) for TMB calculation ('tmb_af_min' = {str(arg_dict['tmb_af_min'])}) must be "
-        err_msg += f"greater or equal to minimum AF tumor {str(arg_dict['tumor_dp_min'])} (i.e. filter for variant inclusion in report)"
-        error_message(err_msg, logger)
+        if arg_dict['tumor_af_min'] is not None:
+            if float(arg_dict['tmb_af_min']) > 0 and (float(arg_dict['tmb_af_min']) < float(arg_dict['tumor_af_min'])):
+                err_msg = f"Minimum AF (tumor) for TMB calculation ('tmb_af_min' = {str(arg_dict['tmb_af_min'])}) must be "
+                err_msg += f"greater or equal to minimum AF tumor ({str(arg_dict['tumor_af_min'])}, i.e. global filter for variant inclusion in report)"
+                error_message(err_msg, logger)
 
     # Check that coding target size region of sequencing assay is set correctly
     if float(arg_dict['effective_target_size_mb']) < 0 or float(arg_dict['effective_target_size_mb']) > float(pcgr_vars.CODING_EXOME_SIZE_MB):

pcgr/config.py

@@ -1,13 +1,13 @@
 #!/usr/bin/env python
 
 from pcgr import pcgr_vars
-from pcgr.utils import check_file_exists, error_message
+from pcgr.utils import check_file_exists, error_message, warn_message
 
 import pandas as pd
 import os
 import csv
 
-def create_config(arg_dict, workflow = "PCGR"):
+def create_config(arg_dict, workflow = "PCGR", logger=None):
 
     conf_options = {}
     if workflow == "PCGR" or workflow == "CPSR":
@@ -65,7 +65,21 @@ def create_config(arg_dict, workflow = "PCGR"):
         if arg_dict['tumor_ploidy'] is not None:
             conf_options['sample_properties']['tumor_ploidy'] = float(arg_dict['tumor_ploidy'])
         if arg_dict['sex'] is not None:
+            sex = str(arg_dict['sex'])
+            ## Breast, Ovary/Fallopian Tube, Uterus, Vulva/Vagina
+            if arg_dict['tsite'] == 6 or arg_dict['tsite'] == 18 or arg_dict['tsite'] == 29 or arg_dict['tsite'] == 30:
+                if sex == 'MALE':
+                    if arg_dict['tsite'] == 6:
+                        warn_message(f"Tumor site '{conf_options['sample_properties']['site']}' is typically observed in females - ensure '--sex' option is correctly set", logger)
+                    else:
+                        error_message(f"Tumor site '{conf_options['sample_properties']['site']}' is not observed in males - please check the '--sex' option", logger)
+            ## Prostate, Testis
+            if arg_dict['tsite'] == 23 or arg_dict['tsite'] == 26:
+                if sex == 'FEMALE':
+                    error_message(f"Tumor site '{conf_options['sample_properties']['site']}' is not observed in females - please check the '--sex' option", logger)
+
             conf_options['sample_properties']['sex'] = str(arg_dict['sex'])
+
 
         #conf_options['clinicaltrials'] = {
         #    'run': int(arg_dict['include_trials'])
@@ -97,10 +111,12 @@ def create_config(arg_dict, workflow = "PCGR"):
 
         conf_options['somatic_snv'] = {}
         conf_options['somatic_snv']['allelic_support'] = {
-            'tumor_dp_min': int(arg_dict['tumor_dp_min']),
-            'control_dp_min': int(arg_dict['control_dp_min']),
-            'tumor_af_min': float(arg_dict['tumor_af_min']),
-            'control_af_max': float(arg_dict['control_af_max']),
+            'tumor_dp_min': arg_dict['tumor_dp_min'],
+            'control_dp_min': arg_dict['control_dp_min'],
+            'tumor_af_min': arg_dict['tumor_af_min'],
+            'control_af_max': arg_dict['control_af_max'],
+            'tumor_ad_min': arg_dict['tumor_ad_min'],
+            'control_ad_max': arg_dict['control_ad_max'],
             'control_dp_tag': str(arg_dict['control_dp_tag']),
             'control_af_tag': str(arg_dict['control_af_tag']),
             'tumor_dp_tag': str(arg_dict['tumor_dp_tag']),
@@ -119,11 +135,16 @@ def create_config(arg_dict, workflow = "PCGR"):
         conf_options['somatic_snv']['msi'] = {
             'run': int(arg_dict['estimate_msi'])
         }
+
+        tmb_dp_min = 0 if arg_dict['tmb_dp_min'] is None else int(arg_dict['tmb_dp_min'])
+        tmb_af_min = 0.0 if arg_dict['tmb_af_min'] is None else float(arg_dict['tmb_af_min'])
+        tmb_ad_min = 0 if arg_dict['tmb_ad_min'] is None else int(arg_dict['tmb_ad_min'])
         conf_options['somatic_snv']['tmb'] = {
             'run': int(arg_dict['estimate_tmb']), 
             'tmb_display': arg_dict['tmb_display'],           
-            'tmb_dp_min': arg_dict['tmb_dp_min'],
-            'tmb_af_min': arg_dict['tmb_af_min']
+            'tmb_dp_min': tmb_dp_min,
+            'tmb_af_min': tmb_af_min,
+            'tmb_ad_min': tmb_ad_min
         }
 
         for pop in ['nfe', 'fin', 'amr', 'eas', 'sas', 'asj', 'oth', 'afr', 'global']: