documentation of MaxEntScan support

Author: sigven

pcgr/annoutils.py

@@ -236,6 +236,7 @@ def assign_cds_exon_intron_annotations(csq_record, grantham_scores, logger):
     csq_record['CDS_RELATIVE_POSITION'] = '.'
     csq_record['LOSS_OF_FUNCTION'] = False
     csq_record['LOF_FILTER'] = '.'
+    csq_record['MAXENTSCAN'] = '.'
 
     splice_variant = False
     #print(csq_record.keys())
@@ -262,25 +263,33 @@ def assign_cds_exon_intron_annotations(csq_record, grantham_scores, logger):
     if re.search(pcgr_vars.CSQ_NULL_PATTERN, str(csq_record['Consequence'])) is not None:
         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:
+        csq_record['MAXENTSCAN'] = 'MES|' + str(csq_record['MaxEntScan_diff']) + '|' + \
+            str(csq_record['MaxEntScan_ref']) + '|' + str(csq_record['MaxEntScan_alt'])
+    
     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:
             csq_record['SPLICE_DONOR_RELEVANT'] = True
         if not csq_record['MaxEntScan_diff'] is None and re.search('splice_donor_(5th|variant)',str(csq_record['Consequence'])) is not None:
             if abs(csq_record['MaxEntScan_diff']) >= pcgr_vars.DONOR_DISRUPTION_MES_CUTOFF:
                 csq_record['LOSS_OF_FUNCTION'] = True
+                csq_record['MAXENTSCAN'] = str(csq_record['MAXENTSCAN']) + '|DONOR_DISRUPTING'
             else:
                 if csq_record['LOSS_OF_FUNCTION'] is True:
                     csq_record['LOSS_OF_FUNCTION'] = False
                     csq_record['LOF_FILTER'] = "NON_DONOR_DISRUPTING"
+                csq_record['MAXENTSCAN'] = str(csq_record['MAXENTSCAN']) + '|NON_DONOR_DISRUPTING'
 
     if re.search(pcgr_vars.CSQ_SPLICE_ACCEPTOR_PATTERN, str(csq_record['Consequence'])) is not None:
         if not csq_record['MaxEntScan_diff'] is None and re.search('splice_acceptor', str(csq_record['Consequence'])) is not None:
             if abs(csq_record['MaxEntScan_diff']) >= pcgr_vars.ACCEPTOR_DISRUPTION_MES_CUTOFF:
                 csq_record['LOSS_OF_FUNCTION'] = True
+                csq_record['MAXENTSCAN'] = str(csq_record['MAXENTSCAN']) + '|ACCEPTOR_DISRUPTING'
             else:
                 if csq_record['LOSS_OF_FUNCTION'] is True:
                     csq_record['LOSS_OF_FUNCTION'] = False
                     csq_record['LOF_FILTER'] = "NON_ACCEPTOR_DISRUPTING"
+                csq_record['MAXENTSCAN'] = str(csq_record['MAXENTSCAN']) + '|NON_ACCEPTOR_DISRUPTING'
 
     if re.search(pcgr_vars.CSQ_SPLICE_REGION_PATTERN, str(csq_record['Consequence'])) is not None:
         match = re.search(

pcgrr/R/input_data.R

@@ -640,6 +640,21 @@ load_dna_variants <- function(
       )
   }
 
+  if ("MAXENTSCAN" %in% colnames(results[['variant']]) &
+      "SPLICE_EFFECT" %in% colnames(results[['variant']])) {
+    results[['variant']] <-
+      results[['variant']] |>
+      dplyr::mutate(
+        SPLICE_EFFECT = paste(
+          SPLICE_EFFECT, MAXENTSCAN, sep = ", "
+      )) |>
+      dplyr::mutate(
+        SPLICE_EFFECT = stringr::str_replace_all(
+          .data$SPLICE_EFFECT, "^(NA, )|^(NA, NA)$", ""
+        )
+      )
+  }
+
   ## Rename annotations for more clarity
   if ("TSG_SUPPORT" %in% colnames(results[['variant']])) {
     results[['variant']] <-

pcgrr/data-raw/data-raw.R

@@ -206,6 +206,8 @@ data_coltype_defs[['snv_indel_somatic_raw']] <- readr::cols_only(
   IMPACT = readr::col_character(),
   LOSS_OF_FUNCTION = readr::col_logical(),
   LOF_FILTER = readr::col_character(),
+  MAXENTSCAN = readr::col_character(),
+  MaxEntScan_diff = readr::col_number(),
   SPLICE_DONOR_RELEVANT = readr::col_logical(),
   SPLICE_EFFECT_MUTSPLICEDB = readr::col_character(),
   NULL_VARIANT = readr::col_logical(),
@@ -313,6 +315,8 @@ data_coltype_defs[['snv_indel_germline_raw']] <- readr::cols_only(
   IMPACT = readr::col_character(),
   LOSS_OF_FUNCTION = readr::col_logical(),
   LOF_FILTER = readr::col_character(),
+  MAXENTSCAN = readr::col_character(),
+  MaxEntScan_diff = readr::col_number(),
   SPLICE_DONOR_RELEVANT = readr::col_logical(),
   SPLICE_EFFECT_MUTSPLICEDB = readr::col_character(),
   NULL_VARIANT = readr::col_logical(),
@@ -787,6 +791,7 @@ dt_display[['snv_indel_gene_actionable']] <-
     'HGVSc',
     'HGVSc_RefSeq',
     'PREDICTED_EFFECT',
+    'SPLICE_EFFECT',
     'LOSS_OF_FUNCTION',
     'LOF_FILTER',
     'ONCOGENICITY',

pcgrr/data/data_coltype_defs.rda

@@ -107,15 +107,17 @@ For SNVs and InDels, PCGR considers the following consequence types as candidate
  * Frameshift variants - [SO:0001589](http://www.sequenceontology.org/miso/current_svn/term/SO:0001589)
  * Splice-site disruptions (2bp donor/acceptor site) - [SO:0001574](http://www.sequenceontology.org/miso/current_svn/term/SO:0001574),
  [SO:0001575](http://www.sequenceontology.org/miso/current_svn/term/SO:0001575)
- * Start losses - [SO:0001574](http://www.sequenceontology.org/miso/current_svn/term/SO:0002012)
+ * Splice donor 5th base disruptions - [SO:0001787](http://www.sequenceontology.org/miso/current_svn/term/SO:0001787)
+ * Start losses - [SO:0002012](http://www.sequenceontology.org/miso/current_svn/term/SO:0002012)
 
-If variants of other consequence types (e.g. synonymous variants, or splice site variants beyond the canonical 2bp site) are found to affect splicing, specifically through records in [MutSpliceDB](https://brb.nci.nih.gov/splicing/), these are also marked as loss-of-function candidates.
+If variants of other consequence types (e.g. synonymous variants, or splice site variants beyond the canonical 2bp site/5th donor base) are found to affect splicing, specifically through records in [MutSpliceDB](https://brb.nci.nih.gov/splicing/), these are also marked as loss-of-function candidates.
 
 A collection of filters is next applied, which can potentialy remove the loss-of-function property for candidates identified above:
 
  * Frameshifts/stop gains within the last 5% of the CDS
  * Splice site variants that are not predicted to affect a donor site (GC -> GT)
  * Variants where [MaxEntScan](https://pubmed.ncbi.nlm.nih.gov/15285897/) does not predict an effect on splicing
+    - A MaxEntScan score difference (between reference and alternative allele) of at least 6 is considered disruptive for donor sites, while 7 is required for acceptor sites (thresholds adopted from [LOFTEE](https://github.com/konradjk/loftee))
     - Annotations in MutSpliceDB will have precedence if any conflicting evidence with MaxEntScan output is found
 
 If a variant is filtered as non-LoF through any of these criteria, this will be evident from the `LOF_FILTER` variable (found in the interactive tables of the HTML report as well as the TSV/Excel output).

pcgrr/inst/templates/pcgr_quarto_report/documentation.qmd

@@ -239,11 +239,19 @@ noncoding_variant_set <-
          "PROTEIN_DOMAIN",
          "MUTATION_HOTSPOT",
          "MUTATION_HOTSPOT_CANCERTYPE",
-         "HGVSc_RefSeq",
-         "PREDICTED_EFFECT",
          "TARGETED_INHIBITORS",
          "TARGETED_INHIBITORS_ALL",
          "VEP_ALL_CSQ"))) |>
+  dplyr::select(
+    c("SYMBOL",
+      "ALTERATION",
+      "GENENAME",
+      "CONSEQUENCE",
+      "ONCOGENICITY",
+      "COSMIC_ID",
+      "LOSS_OF_FUNCTION"),
+    dplyr::everything()
+  ) |>
  dplyr::arrange(
     dplyr::desc(.data$ONCOGENICITY_SCORE),
     dplyr::desc(.data$TISSUE_ASSOC_RANK),

pcgrr/inst/templates/pcgr_quarto_report/snv_indel/oncogenicity.qmd

@@ -67,6 +67,7 @@ for oncogenicity classification
 - Added oncogenicity documentation in variant tables of HTML report (indicates which criteria that was matched for a given variant)
 - Removed `VEP_ALL_CSQ` from variant tables in HTML report to reduce file size - still available in TSV output
 - Added MutSpliceDB splice site effects (used e.g. for loss-of-function annotation)
+- Added MaxEntScan plugin in VEP for splice site disruption prediction (used e.g. for loss-of-function annotation)
 - Added more prediction algorithms from dbNSFP
 - Fixed bug that caused crash for missing values in DP/AF values of input VCF
 - Fixed erroneous re-formatting of MAF integer columns to floats

pcgrr/vignettes/CHANGELOG.Rmd

@@ -189,7 +189,8 @@ A VCF file containing annotated, somatic calls (single nucleotide variants and i
 | `PFAM_DOMAIN` | Pfam domain identifier (from VEP) |
 | `INTOGEN_DRIVER_MUT` | Indicates if existing variant is predicted as driver mutation from IntoGen Catalog of Driver Mutations |
 | `EFFECT_PREDICTIONS` | Insilico predictions variant effect on protein function and pre-mRNA splicing from [database of non-synonymous functional predictions - dbNSFP v5.0](https://www.dbnsfp.org/). Predicted effects are provided by different sources/algorithms (separated by `&`), `T` = Tolerated, `N` = Neutral, `D` = Damaging |
-| `SPLICE_EFFECT` | Effect of splicing, from MutSpliceDB. Format: <ENTREZGENE>|<SYMBOL>|<REFSEQ_TRANSCRIPT_ID>|<HGVSc>|<EFFECT><AFFECTED_EXONS>|<SOURCE> |
+| `SPLICE_EFFECT` | Effect of splicing, from MutSpliceDB and/or MaxEntScan. Format: <ENTREZGENE>|<SYMBOL>|<REFSEQ_TRANSCRIPT_ID>|<HGVSc>|<EFFECT><AFFECTED_EXONS>|<SOURCE> (MutSpliceDB), 
+MES|<SCORE_DIFF>|SCORE_REF>|<SCORE_ALT>|EFFECT (MaxEntScan). |
 | `DBNSFP_BAYESDEL_ADDAF` | predicted effect from BayesDel (dbNSFP) |
 | `DBNSFP_LIST_S2` | predicted effect from LIST-S2 (dbNSFP) |
 | `DBNSFP_SIFT` | predicted effect from SIFT (dbNSFP) |
@@ -322,7 +323,7 @@ The following variables are included in the TSV file (VCF tags issued by the use
 | 48. `TUMOR_SUPPRESSOR_SUPPORT` | Tumor suppressor annotation support (CGC/CancerMine/NCG) |
 | 49. `TARGETED_INHIBITORS2` | Targeted inhibitors |
 | 50. `EFFECT_PREDICTIONS` | Variant effect predictions - from dbNSFP |
-| 51. `SPLICE_EFFECT` | Splice effect annotation from MutSpliceDB |
+| 51. `SPLICE_EFFECT` | Splice effect annotations from MutSpliceDB and MaxEntScan (see details above) |
 | 52. `REGULATORY_ANNOTATION` | Regulatory annotation |
 | 53. `VEP_ALL_CSQ` | VEP consequence - all transcripts |
 | 54. `gnomADe_AF` | gnomAD exomes allele frequency - globally |