Fixes for increased stability

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: SingleMoleculeGenomicsIO
 Title: Tools for reading single molecule genomics data 
-Version: 0.1.5
+Version: 0.1.6
 Authors@R: c(
     person("Michael", "Stadler", , "michael.stadler@fmi.ch",
            role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2269-4934")),

NEWS.md

@@ -1,3 +1,8 @@
+# SingleMoleculeGenomicsIO 0.1.6
+
+* Fix reporting of total number of processed modifications
+* Avoid processing of prior modification data upon removing of soft-masked bases (may improve performance of readModBam several-fold)
+
 # SingleMoleculeGenomicsIO 0.1.5
 
 * Breaking change: Require specification of column names for read regrouping by to include the name of the nested `colData` column where they should be found

src/filter_modbam_cpp.cpp

@@ -136,11 +136,11 @@ Rcpp::NumericVector filter_modbam_cpp(std::string infile,
     Rcpp::RObject bar;
 
     // determine outmode based on extension of outfile
-    if (outfile.compare(outfile.size() - 4, 4, ".bam") == 0 ||
-        outfile.compare(outfile.size() - 4, 4, ".BAM") == 0) {
+    if (outfile.size() >= 4 && (outfile.compare(outfile.size() - 4, 4, ".bam") == 0 ||
+        outfile.compare(outfile.size() - 4, 4, ".BAM") == 0)) {
         outmode = "wb";
-    } else if (outfile.compare(outfile.size() - 4, 4, ".sam") == 0 ||
-        outfile.compare(outfile.size() - 4, 4, ".SAM") == 0) {
+    } else if (outfile.size() >= 4 && (outfile.compare(outfile.size() - 4, 4, ".sam") == 0 ||
+        outfile.compare(outfile.size() - 4, 4, ".SAM") == 0)) {
         outmode = "w";
     } else {
         had_error = true;

src/read_modbam_cpp.cpp

@@ -155,9 +155,10 @@ int process_bam_record(bam1_t *bamdata,        // bam record
 
     // ... remove unaligned (e.g. soft-masked) read-bases
     //     (iterate backwards to avoid messing up indices
-    //      when removing elements)
-    n_total += ref_position.size();
-    for (size_t e = ref_position.size(); e-- > 0;) {
+    //     when removing elements). Only iterate over the
+    //     elements added for this read
+    n_total += aligned_read_position_converted.size();
+    for (size_t e = ref_position.size(); e-- > size_before_this_read;) {
         if (ref_position[e] == -1) {
             n_unaligned++;
             read_id.erase(read_id.begin() + e);

src/utils.cpp

@@ -92,11 +92,13 @@ std::string concatenate_hts_files(std::vector<std::string> input_files,
     }
 
     // determine outmode based on extension of output_file
-    if (output_file.compare(output_file.size() - 4, 4, ".bam") == 0 ||
-        output_file.compare(output_file.size() - 4, 4, ".BAM") == 0) {
+    if (output_file.size() >= 4 &&
+        (output_file.compare(output_file.size() - 4, 4, ".bam") == 0 ||
+        output_file.compare(output_file.size() - 4, 4, ".BAM") == 0)) {
         outmode = "wb";
-    } else if (output_file.compare(output_file.size() - 4, 4, ".sam") == 0 ||
-        output_file.compare(output_file.size() - 4, 4, ".SAM") == 0) {
+    } else if (output_file.size() >= 4 &&
+        (output_file.compare(output_file.size() - 4, 4, ".sam") == 0 ||
+        output_file.compare(output_file.size() - 4, 4, ".SAM") == 0)) {
         outmode = "w";
     } else {
         had_error = true;