fix: problem with variants in multi-stop codon txs (#95)

src/mapper/altseq.rs

@@ -206,6 +206,10 @@ pub struct AltSeqBuilder {
     pub reference_data: RefTranscriptData,
     /// Whether the transcript reference amino acid sequence has multiple stop codons.
     pub ref_has_multiple_stops: bool,
+    /// Position of the first stop codon in the reference amino acid.  It is difficult
+    /// to handle transcripts with multiple stop codons, but we can make a reasonable
+    /// prediction up to the first stop codon.  This rescues changes to the `p.Met1`.
+    pub first_stop_pos: Option<usize>,
 }
 
 impl AltSeqBuilder {
@@ -214,11 +218,13 @@ impl AltSeqBuilder {
             panic!("Must initialize with HgvsVariant::CdsVariant");
         }
         let ref_has_multiple_stops = reference_data.aa_sequence.matches('*').count() > 1;
+        let first_stop_pos = reference_data.aa_sequence.find('*');
 
         Self {
             var_c,
             reference_data,
             ref_has_multiple_stops,
+            first_stop_pos,
         }
     }
 
@@ -433,7 +439,7 @@ impl AltSeqBuilder {
             &self.reference_data.protein_accession,
             &self.reference_data.aa_sequence,
             is_substitution,
-            self.ref_has_multiple_stops,
+            self.ref_has_multiple_stops && self.first_stop_pos.map(|p| p <= start).unwrap_or(false),
         )
     }
 
@@ -465,7 +471,7 @@ impl AltSeqBuilder {
             &self.reference_data.protein_accession,
             &self.reference_data.aa_sequence,
             false,
-            self.ref_has_multiple_stops,
+            self.ref_has_multiple_stops && self.first_stop_pos.map(|p| p <= start).unwrap_or(false),
         )
     }
 
@@ -496,7 +502,7 @@ impl AltSeqBuilder {
             &self.reference_data.protein_accession,
             &self.reference_data.aa_sequence,
             false,
-            self.ref_has_multiple_stops,
+            self.ref_has_multiple_stops && self.first_stop_pos.map(|p| p <= start).unwrap_or(false),
         )
     }
 
@@ -975,7 +981,7 @@ impl AltSeqToHgvsp {
             } else {
                 // is non-dup insertion
                 let start = *start as usize - 1;
-                let end = start as usize + 1;
+                let end = start + 1;
 
                 aa_start = Some(ProtPos {
                     aa: self.ref_seq()[(start - 1)..start].to_owned(),

src/mapper/assembly.rs

@@ -437,12 +437,18 @@ mod test {
         }
 
         #[test]
-        fn test_fail_c_to_p_brca2() -> Result<(), anyhow::Error> {
+        fn test_c_to_p_brca2() -> Result<(), anyhow::Error> {
             let mapper = build_mapper_38(true)?;
             let hgvs_c = "NM_000059.3:c.7790delAAG";
             let var_c = HgvsVariant::from_str(hgvs_c)?;
 
-            assert!(mapper.c_to_p(&var_c).is_err());
+            // NB: this used to fail in Python hgvs but works here now as we do not
+            // perform comprehensive validation yet (1 bp interval/position, but 3bp
+            // deleted).
+            assert_eq!(
+                "NP_000050.2:p.Glu2598LysfsTer50",
+                format!("{}", mapper.c_to_p(&var_c)?)
+            );
 
             Ok(())
         }

src/mapper/variant.rs

@@ -1626,6 +1626,9 @@ mod test {
         Ok(())
     }
 
+    /// Check the case with multiple stop codons.  We introduced a change in hgvs-rs
+    /// that does not handle multiple stop codons in the transcript sequence as
+    /// conservatively as the Python version.
     #[test]
     fn hgvs_c_to_p_multiple_stop_codons() -> Result<(), anyhow::Error> {
         let hgvsc = "NM_999992.1:c.4G>A";
@@ -1904,6 +1907,12 @@ mod test {
         run_gxp_test("tests/data/mapper/gcp/real.tsv", false)
     }
 
+    /// Check for issues with variants affecting `Met1` leading to `p.Met1?`.
+    #[test]
+    fn real_met1() -> Result<(), anyhow::Error> {
+        run_gxp_test("tests/data/mapper/gcp/real-met1.tsv", false)
+    }
+
     #[test]
     fn noncoding() -> Result<(), anyhow::Error> {
         run_gxp_test("tests/data/mapper/gcp/noncoding.tsv", false)

src/normalizer.rs

@@ -1384,9 +1384,15 @@ mod test {
 
         {
             // gene BRCA2
+            //
+            // NB: this used to fail in Python hgvs but works here now as we do not
+            // perform comprehensive validation yet (1 bp interval/position, but 3bp
+            // deleted).
             let raw = HgvsVariant::from_str("NM_000059.3:c.7790delAAG")?;
-            // with self.assertRaises(HGVSInvalidVariantError):
-            assert!(norm.normalize(&raw).is_err());
+            assert_eq!(
+                "NM_000059.3:c.7791delA",
+                format!("{}", norm.normalize(&raw)?)
+            );
         }
 
         Ok(())