Guard mappers from truncated quality strings

src/alignment.cpp

@@ -2847,6 +2847,30 @@ void alignment_set_distance_to_correct(Alignment& aln, const map<string ,vector<
     }
 }
 
+void check_quality_length(const Alignment& aln) {
+    size_t quality_length = aln.quality().length();
+    if (quality_length == 0 || quality_length == aln.sequence().length()) {
+        // This length is acceptable.
+        return;
+    }
+
+    bool too_short = quality_length < aln.sequence().length();
+
+    std::stringstream ss;
+    ss << "Read " << aln.name() << " has " << aln.sequence().length() << " bases of sequence but ";
+    if (too_short) {
+        ss << "only ";
+    }
+    ss << quality_length << " base quality values.";
+    if (too_short) {
+        ss << " Was the quality information truncated?";
+    }
+
+    #pragma omp critical (cerr)
+    std::cerr << "error [vg::alignment.cpp]: " << ss.str() << std::endl;
+    exit(1);
+}
+
 AlignmentValidity alignment_is_valid(const Alignment& aln, const HandleGraph* hgraph, bool check_sequence) {
     size_t read_idx = 0;
     for (size_t i = 0; i < aln.path().mapping_size(); ++i) {

src/alignment.hpp

@@ -328,6 +328,13 @@ map<string ,vector<pair<size_t, bool> > > alignment_refpos_to_path_offsets(const
 void alignment_set_distance_to_correct(Alignment& aln, const Alignment& base, const unordered_map<string, string>* translation = nullptr);
 void alignment_set_distance_to_correct(Alignment& aln, const map<string, vector<pair<size_t, bool>>>& base_offsets, const unordered_map<string, string>* translation = nullptr);
 
+/**
+ * Stop the program and print a useful error message if the alignment has
+ * quality values, but not the right number of them for the number of sequence
+ * bases. 
+ */
+void check_quality_length(const Alignment& aln);
+
 /**
  * Represents a report on whether an alignment makes sense in the context of a graph.
  */

src/subcommand/giraffe_main.cpp

@@ -2042,6 +2042,9 @@ int main_giraffe(int argc, char** argv) {
                             std::cerr << "Thread " << thread_num << " now mapping " << aln1.name() << ", " << aln2.name() << std::endl;
                         }
 
+
+                        check_quality_length(aln1);
+                        check_quality_length(aln2);
                         toUppercaseInPlace(*aln1.mutable_sequence());
                         toUppercaseInPlace(*aln2.mutable_sequence());
 
@@ -2150,6 +2153,7 @@ int main_giraffe(int argc, char** argv) {
                             std::cerr << "Thread " << thread_num << " now mapping " << aln.name() << std::endl;
                         }
 
+                        check_quality_length(aln);
                         toUppercaseInPlace(*aln.mutable_sequence());
 
                         // Map the read with the MinimizerMapper.

src/subcommand/map_main.cpp

@@ -978,6 +978,8 @@ int main_map(int argc, char** argv) {
                 return;
             }
 
+            check_quality_length(alignment);
+
             int tid = omp_get_thread_num();
             vector<Alignment> alignments = mapper[tid]->align_multi(alignment,
                                                                     kmer_size,
@@ -1016,6 +1018,8 @@ int main_map(int argc, char** argv) {
             };
 
             function<void(Alignment&,Alignment&)> lambda = [&](Alignment& aln1, Alignment& aln2) {
+                check_quality_length(aln1);
+                check_quality_length(aln2);
                 auto our_mapper = mapper[omp_get_thread_num()];
                 bool queued_resolve_later = false;
                 auto alnp = our_mapper->align_paired_multi(aln1,
@@ -1067,6 +1071,7 @@ int main_map(int argc, char** argv) {
         } else if (fastq2.empty()) {
             // single
             function<void(Alignment&)> lambda = [&](Alignment& alignment) {
+                        check_quality_length(alignment);
                         int tid = omp_get_thread_num();
                         vector<Alignment> alignments = mapper[tid]->align_multi(alignment,
                                                                                 kmer_size,
@@ -1092,6 +1097,8 @@ int main_map(int argc, char** argv) {
                 }
             };
             function<void(Alignment&,Alignment&)> lambda = [&](Alignment& aln1, Alignment& aln2) {
+                check_quality_length(aln1);
+                check_quality_length(aln2);
                 auto our_mapper = mapper[omp_get_thread_num()];
                 bool queued_resolve_later = false;
                 auto alnp = our_mapper->align_paired_multi(aln1,
@@ -1165,6 +1172,8 @@ int main_map(int argc, char** argv) {
                 }
             };
             function<void(Alignment&,Alignment&)> lambda = [&](Alignment& aln1, Alignment& aln2) {
+                check_quality_length(aln1);
+                check_quality_length(aln2);
                 auto our_mapper = mapper[omp_get_thread_num()];
                 bool queued_resolve_later = false;
                 auto alnp = our_mapper->align_paired_multi(aln1,
@@ -1214,6 +1223,7 @@ int main_map(int argc, char** argv) {
         } else {
             // Processing single-end GAM input
             function<void(Alignment&)> lambda = [&](Alignment& alignment) {
+                check_quality_length(alignment);
                 int tid = omp_get_thread_num();
                 vector<Alignment> alignments = mapper[tid]->align_multi(alignment,
                                                                         kmer_size,

src/subcommand/mpmap_main.cpp

@@ -2115,7 +2115,8 @@ int main_mpmap(int argc, char** argv) {
         if (watchdog) {
             watchdog->check_in(thread_num, alignment.name());
         }
-
+
+        check_quality_length(alignment);
         toUppercaseInPlace(*alignment.mutable_sequence());
 
         bool is_rna = uses_Us(alignment);
@@ -2181,6 +2182,8 @@ int main_mpmap(int argc, char** argv) {
             watchdog->check_in(thread_num, alignment_1.name());
         }
 
+        check_quality_length(alignment_1);
+        check_quality_length(alignment_2);
         toUppercaseInPlace(*alignment_1.mutable_sequence());
         toUppercaseInPlace(*alignment_2.mutable_sequence());