fix: check for length being a multiple of 3 before codon aware stripping

nextstrain · Feb 11, 2021 · a33dd6b · a33dd6b
1 parent 9a457c3
commit a33dd6b
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 6 deletions.
diff --git a/packages/nextalign/src/align/getGapOpenCloseScores.cpp b/packages/nextalign/src/align/getGapOpenCloseScores.cpp
@@ -33,7 +33,7 @@ std::vector<int> getGapOpenCloseScoresCodonAware(//
   for (const auto& [geneName, gene] : geneMap) {
 
     // TODO: might use std::fill()
-    for (int i = gene.start; i <= gene.end-2; i += 3) {
+    for (int i = gene.start; i < gene.end-2; i += 3) {
       gapOpenClose[i] = options.alignment.scoreGapOpenInFrame;
       gapOpenClose[i+1] = options.alignment.scoreGapOpenOutOfFrame;
       gapOpenClose[i+2] = options.alignment.scoreGapOpenOutOfFrame;

diff --git a/packages/nextalign/src/translate/extractGene.cpp b/packages/nextalign/src/translate/extractGene.cpp
@@ -9,7 +9,6 @@
 #include "../utils/safe_cast.h"
 #include "./removeGaps.h"
 
-
 NucleotideSequenceView extractGeneRef(const NucleotideSequenceView& ref, const Gene& gene) {
   precondition_less(gene.length, ref.size());
   precondition_less_equal(gene.length, ref.size());
@@ -86,6 +85,10 @@ NucleotideSequence extractGeneQuery(
 
 
   auto result = NucleotideSequence(query.substr(start, length));
+  const auto resultLengthPreStrip = safe_cast<int>(result.size());
+  if (resultLengthPreStrip % 3 != 0) {
+    throw ErrorExtractGeneLengthInvalid(gene.geneName, resultLengthPreStrip);
+  }
   stripGeneInPlace(result);
 
   const auto resultLength = safe_cast<int>(result.size());

diff --git a/packages/nextalign/src/translate/translateGenes.cpp b/packages/nextalign/src/translate/translateGenes.cpp
@@ -44,10 +44,6 @@ PeptidesInternal translateGenes(         //
 
   const auto coordMap = mapCoordinates(ref);
 
-  // Each position in the raw ref sequence should have a corresponding mapped position in aligned ref sequence
-  invariant_equal(coordMap.size(), ref.size());
-
-
   std::vector<PeptideInternal> queryPeptides;
   queryPeptides.reserve(geneMap.size());