Merge branch 'main' of github.com:sigven/pcgr

Author: sigven

README.md

@@ -28,6 +28,9 @@ PCGR originates from the [Norwegian Cancer Genomics Consortium (NCGC)](https://c
 
 ### Top News
 
+- *July 15th 2025:* **2.2.2 release**
+  - [CHANGELOG](https://sigven.github.io/pcgr/articles/CHANGELOG.html)
+
 - *March 23rd 2025:* **2.2.1 release**
   - fix bug in CPSR for ClinVar variants with non-standard significance levels
 

pcgr/annoutils.py

@@ -230,7 +230,7 @@ def assign_cds_exon_intron_annotations(csq_record, grantham_scores, logger):
     csq_record['CDS_CHANGE'] = '.'
     csq_record['HGVSp_short'] = '.'
     csq_record['PROTEIN_CHANGE'] = '.'
-    csq_record['GRANTHAM_DISTANCE'] = '.'
+    csq_record['GRANTHAM_DISTANCE'] = -1
     csq_record['ALTERATION'] = '.'
     csq_record['EXON_AFFECTED'] = '.'
     csq_record['CDS_RELATIVE_POSITION'] = '.'
@@ -264,8 +264,14 @@ def assign_cds_exon_intron_annotations(csq_record, grantham_scores, logger):
         csq_record['NULL_VARIANT'] = True
 
     if not csq_record['MaxEntScan_diff'] is None and not csq_record['MaxEntScan_ref'] is None and not csq_record['MaxEntScan_alt'] is None:
+        fraction_drop = 0.0
+        #if float(csq_record['MaxEntScan_ref']) > 0:
+        #    fraction_drop = float(csq_record['MaxEntScan_diff']) / float(csq_record['MaxEntScan_ref']).round(4)
+        #else:
+        #    fraction_drop = 0.0
         csq_record['MAXENTSCAN'] = 'MES|' + str(csq_record['MaxEntScan_diff']) + '|' + \
-            str(csq_record['MaxEntScan_ref']) + '|' + str(csq_record['MaxEntScan_alt'])
+            str(csq_record['MaxEntScan_ref']) + '|' + str(csq_record['MaxEntScan_alt']) #+ \
+            #'|' + str(fraction_drop)
 
     if re.search(pcgr_vars.CSQ_SPLICE_DONOR_PATTERN, str(csq_record['Consequence'])) is not None:
         if re.search(r'(\+3(A|G)>|\+4A>|\+5G>)', str(csq_record['HGVSc'])) is not None:
@@ -435,10 +441,14 @@ def assign_cds_exon_intron_annotations(csq_record, grantham_scores, logger):
                         if 'Amino_acids' in csq_record.keys():
                             if not csq_record['Amino_acids'] is None:
                                 if '/' in str(csq_record['Amino_acids']):
-                                    key = str(csq_record['Amino_acids']).split('/')[0] + '_' + str(csq_record['Amino_acids']).split('/')[1]
-                                    if key in grantham_scores.keys():
-                                        #print('GRANTHAM\t' + str(key) + ' ' + str(grantham_scores[key]))
-                                        csq_record['GRANTHAM_DISTANCE'] = grantham_scores[key]
+                                    aaref = str(csq_record['Amino_acids']).split('/')[0]
+                                    aalt = str(csq_record['Amino_acids']).split('/')[1]
+                                    key = str(aaref) + '_' + str(aalt)
+                                    ## check cases for double aminio acid changes in 'Amino_acids', e.g. GQ/GY - only consider the second amino acid change
+                                    if len(aaref) == 2 and len(aalt) == 2 and aaref[0] == aalt[0]:                                        
+                                        key = aaref[1] + '_' + aalt[1]
+                                    if key in grantham_scores.keys():                                        
+                                        csq_record['GRANTHAM_DISTANCE'] = grantham_scores[key]                                   
 
 
                     csq_record['PROTEIN_CHANGE'] = protein_change_VEP

pcgr/cna.py

@@ -145,7 +145,7 @@ def annotate_cna_segments(output_segment_gene_fname: str,
     ## load copy-number biomarker evidence
 
     cna_query_segment_df['segment_length_mb'] = \
-        ((cna_query_segment_df['segment_end'] - cna_query_segment_df['segment_start']) / 1e6).astype(float).round(5)
+        ((cna_query_segment_df['segment_end'] - cna_query_segment_df['segment_start']) / 1e6).astype(float).round(4)
 
     biomarkers = {}
     cna_actionable_dict = {}
@@ -170,29 +170,32 @@ def annotate_cna_segments(output_segment_gene_fname: str,
 
     ## Mark copy number amplifications (threshold defined by user) in input
     cna_query_segment_df['aberration_key'] = 'nan'
-    cna_query_segment_df['amp_cond'] = True
-    cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] < n_copy_amplifications,"amp_cond"] = False
+    cna_query_segment_df['amp_cond'] = False
     cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] >= n_copy_amplifications,"amp_cond"] = True
 
     cna_query_segment_df.loc[cna_query_segment_df.amp_cond, 'aberration_key'] =  \
         cna_query_segment_df.loc[cna_query_segment_df.amp_cond, 'entrezgene'].astype(str) + '_amplification'
 
     ## Mark homozygous deletions in input
-    cna_query_segment_df['loss_cond'] = True
-    cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] > 0,"loss_cond"] = False
-    cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] == 0,"loss_cond"] = True
+    cna_query_segment_df['homloss_cond'] = False
+    cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] == 0,"homloss_cond"] = True
+    
+    ## Mark heterozygous deletions in input
+    cna_query_segment_df['hetloss_cond'] = False
+    cna_query_segment_df.loc[cna_query_segment_df['n_major'] + cna_query_segment_df['n_minor'] == 1,"hetloss_cond"] = True
 
     cna_query_segment_df['variant_class'] = 'undefined'
     cna_query_segment_df.loc[cna_query_segment_df.amp_cond, 'variant_class'] = 'gain'
-    cna_query_segment_df.loc[cna_query_segment_df.loss_cond, 'variant_class'] = 'homdel'
+    cna_query_segment_df.loc[cna_query_segment_df.homloss_cond, 'variant_class'] = 'homdel'
+    cna_query_segment_df.loc[cna_query_segment_df.hetloss_cond, 'variant_class'] = 'hetdel'    
 
-    cna_query_segment_df.loc[cna_query_segment_df.loss_cond, 'aberration_key'] =  \
-        cna_query_segment_df.loc[cna_query_segment_df.loss_cond, 'entrezgene'].astype(str) + '_ablation'
+    cna_query_segment_df.loc[cna_query_segment_df.homloss_cond, 'aberration_key'] =  \
+        cna_query_segment_df.loc[cna_query_segment_df.homloss_cond, 'entrezgene'].astype(str) + '_ablation'
 
     ## Append actionability evidence to input amplifications (column 'biomarker_match')
     cna_query_segment_df = cna_query_segment_df.merge(
         cna_actionable_df, left_on=["aberration_key"], right_on=["aberration_key"], how="left")
-    cna_query_segment_df.drop(['amp_cond', 'loss_cond', 'aberration_key'], axis=1, inplace=True)    
+    cna_query_segment_df.drop(['amp_cond', 'hetloss_cond', 'homloss_cond','aberration_key'], axis=1, inplace=True)    
     cna_query_segment_df.loc[cna_query_segment_df['biomarker_match'].isnull(),"biomarker_match"] = '.'
 
     ## remove all temporary files

pcgr/oncogenicity.py

@@ -460,7 +460,12 @@ def load_oncogenic_variants(oncogenic_variants_fname: str, logger: Logger):
          gene = str(row['entrezgene'])
          if not 'oncogenic' in str(row['oncogenicity']).lower():
             continue         
-         oncogenic_variants[str(gene) + '-' + str(row['var_id'])] = row
+         oncogenic_variants[str(gene) + '-' + str(row['var_id'])] = row         
+         if 'grantham_distance' in row.keys():
+            if row['grantham_distance'] == '':
+               row['grantham_distance'] = -1
+            else:
+               row['grantham_distance'] = float(row['grantham_distance'])
          if not len(row['hgvsp']) == 0:
             oncogenic_variants[str(gene) + '-' + str(row['hgvsp'])] = row
          if not len(row['hgvs_c']) == 0:
@@ -530,10 +535,11 @@ def match_oncogenic_variants(transcript_csq_elements, oncogenic_variants, rec, p
             if oncogenic_varkeys[oncogenic_varkey].startswith('by_codon'):
                grantham_distance = rec.INFO.get('GRANTHAM_DISTANCE')
                if not grantham_distance is None:
-                  if float(grantham_distance / float(oncogenic_variants[oncogenic_varkey]['grantham_distance'])) > 0.8:                  
-                     if not oncogenic_info in known_oncogenic_sites:
-                        known_oncogenic_sites[oncogenic_info] = []
-                     known_oncogenic_sites[oncogenic_info].append(oncogenic_varkeys[oncogenic_varkey])
+                  if grantham_distance > 0 and oncogenic_variants[oncogenic_varkey]['grantham_distance'] > 0:                     
+                     if float(grantham_distance / float(oncogenic_variants[oncogenic_varkey]['grantham_distance'])) >= 0.8:                  
+                        if not oncogenic_info in known_oncogenic_sites:
+                           known_oncogenic_sites[oncogenic_info] = []
+                        known_oncogenic_sites[oncogenic_info].append(oncogenic_varkeys[oncogenic_varkey])
             else:
                if not oncogenic_info in known_oncogenic_matches:
                   known_oncogenic_matches[oncogenic_info] = []

pcgrr/R/cna.R

@@ -86,7 +86,7 @@ plot_cna_segments <- function(chrom_coordinates = NULL,
 
 
   ## Identify segments that involve oncogene gain or
-  ## tumor suppressor loss
+  ## tumor suppressor loss (homozygous or heterozygous)
   onc_gain_tsg_loss <- cna_gene |>
     dplyr::select(
       c("CHROM", "SEGMENT_START", "SEGMENT_END",
@@ -97,9 +97,11 @@ plot_cna_segments <- function(chrom_coordinates = NULL,
       (.data$ONCOGENE == TRUE &
          .data$VARIANT_CLASS == "gain") |
         (.data$TUMOR_SUPPRESSOR == TRUE &
-           .data$VARIANT_CLASS == "homdel"))
+           (.data$VARIANT_CLASS == "homdel" |
+              .data$VARIANT_CLASS == "hetdel")))
 
   tsg_loss <- data.frame()
+  tsg_het_loss <- data.frame()
   onc_gain <- data.frame()
 
   ## If there are oncogene gains or tumor suppressor losses,
@@ -146,6 +148,27 @@ plot_cna_segments <- function(chrom_coordinates = NULL,
           .groups = "drop")
     }
 
+    tsg_het_loss <- onc_gain_tsg_loss |>
+      dplyr::filter(
+        .data$TUMOR_SUPPRESSOR == TRUE &
+          .data$VARIANT_CLASS == "hetdel")
+
+    ## For now, if multiple TSGs are involved in a lost segment, we will only
+    ## show the top three in the plot (hover)
+    if(NROW(tsg_het_loss) > 0){
+      tsg_het_loss <- tsg_het_loss |>
+        dplyr::arrange(
+          dplyr::desc(.data$TUMOR_SUPPRESSOR_RANK)) |>
+        dplyr::group_by(
+          .data$CHROM,
+          .data$SEGMENT_START,
+          .data$SEGMENT_END) |>
+        dplyr::summarise(
+          TSG_HET_LOSS = paste(
+            utils::head(.data$SYMBOL, 3), collapse = ", "),
+          .groups = "drop")
+    }
+
   }
 
   ## Prepare data for plotting
@@ -199,13 +222,29 @@ plot_cna_segments <- function(chrom_coordinates = NULL,
         !is.na(.data$TSG_LOSS),
         paste0(
           .data$SegmentInfo,
-          "<br> - Tumor suppressor loss: ",
+          "<br> - Tumor suppressor loss (homozygous del): ",
           .data$TSG_LOSS),
         .data$SegmentInfo))
   }else{
     cna_segments_global$TSG_LOSS <- as.character(NA)
   }
 
+  if(NROW(tsg_het_loss) > 0){
+    cna_segments_global <- cna_segments_global |>
+      dplyr::left_join(
+        tsg_het_loss,
+        by = c("CHROM", "SEGMENT_START", "SEGMENT_END")
+      ) |>
+      dplyr::mutate(SegmentInfo = dplyr::if_else(
+        !is.na(.data$TSG_HET_LOSS),
+        paste0(.data$SegmentInfo,
+               "<br> - Tumor suppressor loss (heterozygous del): ",
+               .data$TSG_HET_LOSS),
+        .data$SegmentInfo))
+  }else{
+    cna_segments_global$TSG_HET_LOSS <- as.character(NA)
+  }
+
   if(NROW(onc_gain) > 0){
     cna_segments_global <- cna_segments_global |>
       dplyr::left_join(
@@ -379,10 +418,22 @@ get_oncogenic_cna_events <- function(cna_df_display = NULL){
       )
     )
 
+  tsgene_hetloss_variants <-
+    dplyr::filter(
+      cna_df_display,
+        .data$TUMOR_SUPPRESSOR == TRUE &
+        .data$VARIANT_CLASS == "hetdel") |>
+    dplyr::select(
+      dplyr::any_of(
+        pcgrr::dt_display$cna_other_oncogenic
+      )
+    )
+
   cna_oncogenic_events <-
     dplyr::bind_rows(
       oncogene_gain_variants,
-      tsgene_loss_variants
+      tsgene_loss_variants,
+      tsgene_hetloss_variants
     ) |>
     dplyr::select(
       dplyr::any_of(
@@ -394,6 +445,12 @@ get_oncogenic_cna_events <- function(cna_df_display = NULL){
       dplyr::desc(.data$GLOBAL_ASSOC_RANK),
     )
 
+  if("SEGMENT_LENGTH_MB" %in% colnames(cna_oncogenic_events)){
+    cna_oncogenic_events <- cna_oncogenic_events |>
+      dplyr::mutate(SEGMENT_LENGTH_MB = round(
+        .data$SEGMENT_LENGTH_MB, digits = 2))
+  }
+
   return(cna_oncogenic_events)
 
 

pcgrr/R/input_data.R

@@ -65,7 +65,9 @@ load_somatic_cna <- function(
     )) |>
     tidyr::separate(
       col = "SEGMENT_NAME",
-      into = c("SEGMENT_ID", "N_MAJOR","N_MINOR","ARM","CYTOBAND","EVENT_TYPE"),
+      into = c("SEGMENT_ID", "N_MAJOR",
+               "N_MINOR","ARM",
+               "CYTOBAND","EVENT_TYPE"),
       sep = "\\|",
       remove = T
     ) |>
@@ -90,7 +92,6 @@ load_somatic_cna <- function(
     variant_origin = "Somatic")
 
   callset_cna[['segment']] <- segments
-
   if (NROW(callset_cna$variant) > 0) {
     callset_cna[['variant']] <- callset_cna[['variant']] |>
       dplyr::mutate(CN_TOTAL =

pcgrr/R/report.R

@@ -220,9 +220,11 @@ init_cna_vstats <- function() {
 
   vstats <- list()
   for (t in c("n_tsg_loss",
+              "n_tsg_hetloss",
               "n_oncogene_gain",
               "n_other_drugtarget_gain",
               "n_segments_loss",
+              "n_segments_hetloss",
               "n_segments_gain",
               "n_actionable_tier1",
               "n_actionable_tier2")) {

pcgrr/R/utils.R

@@ -661,6 +661,7 @@ variant_stats_report <- function(
      "VAR_ID" %in% colnames(callset$variant) &
      "VARIANT_CLASS" %in% colnames(callset$variant)){
     for (n in c("n_tsg_loss",
+                "n_tsg_hetloss",
                 "n_oncogene_gain",
                 "n_other_drugtarget_gain")) {
       call_stats[[name]][[n]] <- 0
@@ -672,6 +673,13 @@ variant_stats_report <- function(
           .data$ACTIONABILITY_TIER == 3 &
           .data$VARIANT_CLASS == "homdel") |>
       nrow()
+    call_stats[[name]][["n_tsg_hetloss"]] <-
+      callset$variant |>
+      dplyr::filter(
+        !is.na(.data$ACTIONABILITY_TIER) &
+          .data$ACTIONABILITY_TIER == 3 &
+          .data$VARIANT_CLASS == "hetdel") |>
+      nrow()
     call_stats[[name]][["n_oncogene_gain"]] <-
       callset$variant |>
       dplyr::filter(
@@ -711,6 +719,13 @@ variant_stats_report <- function(
       dplyr::select(.data$VAR_ID) |>
       dplyr::distinct() |>
       NROW()
+
+    call_stats[[name]][["n_segments_hetloss"]] <-
+      callset$variant |>
+      dplyr::filter(.data$VARIANT_CLASS == "hetdel") |>
+      dplyr::select(.data$VAR_ID) |>
+      dplyr::distinct() |>
+      NROW()
   }
 
   if(vartype == 'snv_indel' &

pcgrr/R/variant_classification.R

@@ -173,7 +173,8 @@ assign_amp_asco_tiers <- function(
           dplyr::mutate(AMP_ASCO_TIER_OTHER_VARS = dplyr::if_else(
             (!is.na(.data$TUMOR_SUPPRESSOR) &
                .data$TUMOR_SUPPRESSOR == TRUE &
-               .data$VARIANT_CLASS == "homdel") |
+               (.data$VARIANT_CLASS == "homdel" |
+                  .data$VARIANT_CLASS == "hetdel")) |
               (!is.na(.data$ONCOGENE) &
                  .data$ONCOGENE == TRUE &
                  .data$VARIANT_CLASS == "gain"),
@@ -226,7 +227,8 @@ assign_amp_asco_tiers <- function(
           dplyr::mutate(AMP_ASCO_TIER = dplyr::if_else(
             (!is.na(.data$TUMOR_SUPPRESSOR) &
                .data$TUMOR_SUPPRESSOR == TRUE &
-               .data$VARIANT_CLASS == "homdel") |
+               (.data$VARIANT_CLASS == "homdel" |
+                  .data$VARIANT_CLASS == "hetdel")) |
               (!is.na(.data$ONCOGENE) &
                  .data$ONCOGENE == TRUE &
                  .data$VARIANT_CLASS == "gain"),

pcgrr/data-raw/data-raw.R

@@ -101,8 +101,8 @@ for (c in c("pathogenicity",
     color_palette[[c]][["values"]] <- c("#9B3297", "#0073C2")
   }
   if (c == "cna_variant_class") {
-    color_palette[[c]][["levels"]] <- c("gain", "homdel")
-    color_palette[[c]][["values"]] <- c("#00a65a", "#CD534C")
+    color_palette[[c]][["levels"]] <- c("gain", "hetdel", "homdel")
+    color_palette[[c]][["values"]] <- c("#00a65a", "#E58B85", "#CD534C")
   }
   if (c == "warning") {
     color_palette[[c]] <- "#ff7518"