fix a bug for using csto as the cds length error in stoploss snv muta…

…tion annotation, which leads the error extent length
ltnetcase · Jul 8, 2021 · eb6f0d1 · eb6f0d1
1 parent 0f5cc57
commit eb6f0d1
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 12 deletions.
diff --git a/Changes b/Changes
@@ -309,3 +309,6 @@
 
 1.23  Tue Apr  6 17:52:33 2021 
     - add standard_cHGVS for gene flanking region cHGVS annotation.
+
+1.24  Thu Jul  8 13:29:19 2021
+    - fix a bug for using csto as the cds length error in stoploss snv mutation annotation, which leads the error extent length.
diff --git a/VERSION b/VERSION
@@ -1,2 +1,2 @@
-1.2.3
+1.2.4
 
diff --git a/lib/BedAnno.pm b/lib/BedAnno.pm
@@ -9,13 +9,13 @@ use Time::HiRes qw(gettimeofday tv_interval);
 
 use Tabix;
 
-our $VERSION = '1.23';
+our $VERSION = '1.24';
 
 =head1 NAME
 
 BedAnno - Perl module for annotating variation depend on the BED format database.
 
-=head2 VERSION v1.23
+=head2 VERSION v1.24
 
 From version 0.32 BedAnno will change to support CG's variant shell list
 and use ncbi annotation release 104 as the annotation database
@@ -3270,9 +3270,11 @@ sub getTrChange {
         my $chgvs_5 =
           ($cdsOpt) ? $trannoEnt->{cdsBegin} : $trannoEnt->{rnaBegin};
         my $chgvs_3 = ($cdsOpt) ? $trannoEnt->{cdsEnd} : $trannoEnt->{rnaEnd};
+        my $cds_len = 0;
         if ($cdsOpt) {
             ( $trannoEnt->{protBegin}, $trannoEnt->{protEnd} ) =
               _getCoveredProd( $trdbEnt, $chgvs_5, $chgvs_3 );
+            $cds_len = $trdbEnt->{csto} - $trdbEnt->{csta};
         }
 
         # [ aaPos, codon, aa, polar, frame, [framealt] ]
@@ -3381,15 +3383,14 @@ sub getTrChange {
             my $lesser_len = ($real_rl < $real_al) ? $real_rl : $real_al;
             if (
                 ( $chgvs_5 =~ /^\-(\d+)$/ and $1 <= $lesser_len )
-                or ( $chgvs_5 =~ /^(\d+)$/
-                    and ( $trdbEnt->{csto} - $1 ) <= $lesser_len )
+                or ( $chgvs_5 =~ /^(\d+)$/ and $cdsOpt
+                    and ( $cds_len - $1 ) <= $lesser_len )
               )
             {
-                # 5utr and first cds
                 my $offset = $1;
                 my $opt_53 = ( $chgvs_5 =~ /^\-/ ) ? 1 : 0;
                 if ( !$opt_53 ) {
-                    $offset = $trdbEnt->{csto} - $offset;
+                    $offset = $cds_len - $offset - 1;
                 }
                 my $change_cn = 0;
                 $change_cn++
@@ -3668,8 +3669,7 @@ sub getTrChange {
                     if (
                         (
                                 $cdsOpt
-                            and $chgvs_3 + 1 <=
-                            ( $trdbEnt->{csto} - $trdbEnt->{csta} )
+                            and $chgvs_3 + 1 <= $cds_len
                         )
                         or ( !$cdsOpt and $chgvs_3 + 1 <= $trdbEnt->{len} )
                       )
@@ -4018,7 +4018,7 @@ sub getTrChange {
                     next;
                 }
                 elsif ( $chgvs_5 =~ /^(\d+)\+1/
-                    and $1 == ( $trdbEnt->{csto} - $trdbEnt->{csta} ) )
+                    and $1 == $cds_len )
                 {
                     $trannoEnt->{func} = 'utr-3';
                     next;
@@ -4222,12 +4222,12 @@ sub getTrChange {
                     # probably frame shift flag
                     my $frameshift_flag = ( $diff_ra % 3 > 0 ) ? 1 : 0;
 
-                    my $start_in_cds_flag = ($chgvs_5 > ( $trdbEnt->{csto} - 3 )) ? 0 : 1;
+                    my $start_in_cds_flag = ($chgvs_5 > ( $cds_len - 3 )) ? 0 : 1;
 
                     # end in cds or not.
                     my $end_in_cds_flag = (
                              $chgvs_3 =~ /^\*/
-                          or $chgvs_3 > ( $trdbEnt->{csto} - 3 )
+                          or $chgvs_3 > ( $cds_len - 3 )
                     ) ? 0 : 1;
 
                     my $ready_to_add_3;

diff --git a/t/BedAnno.t b/t/BedAnno.t
@@ -2054,6 +2054,7 @@ my $rep_span_cds_utr3 = {
     },
     'r_End'    => 'C6E',
     'c'	       => 'c.833_*12delAAAAAAAAAAAAAAinsAAAAAAAAAAAAA',
+    'uncurated_cHGVS' => 'c.833_*12delAAAAAAAAAAAAAAinsAAAAAAAAAAAAA',
     'standard_cHGVS'=> 'c.833_*12delinsAAAAAAAAAAAAA',
     'rnaBegin' => 908,
     'ei_Begin' => 'EX6E',
@@ -3924,10 +3925,15 @@ my $stoploss_var_anno = $bare_beda->anno("chr17", 71334725, 71334726, "T", "C");
 my $stoploss_varname = $stoploss_var_anno->{var}->{varName};
 my $stoploss_delvar_anno = $bare_beda->anno("chr17", 71334725, 71334726, "T", "");
 my $stoploss_del_varname = $stoploss_delvar_anno->{var}->{varName};
+my $stoploss_snv_anno = $bare_beda->anno("chr2", 152236045, 152236046, "A", "T");
+my $stoploss_snv_varname = $stoploss_snv_anno->{var}->{varName};
 
 ok( $stoploss_varname eq "NM_001144952.1(SDK2): c.6519A>G (p.*2173Wext*61)",
     "for [ stop-loss snv ]" )
   or explain "The anno info: ", $stoploss_var_anno;
+ok( $stoploss_snv_varname eq "NM_007115.3(TNFAIP6): c.833A>T (p.*278Lext*6)",
+    "for [ stop-loss snv2 ]" )
+  or explain "The anno info: ", $stoploss_snv_anno;
 ok(
     $stoploss_del_varname eq
       "NM_001144952.1(SDK2): c.6519delA (p.*2173Cext*39)",