Promote secondary when primary alignment fails to surject

.gitmodules

@@ -66,7 +66,7 @@
 	url = https://github.com/vgteam/libhandlegraph.git
 [submodule "deps/libvgio"]
 	path = deps/libvgio
-	url = https://github.com/vgteam/libvgio.git
+	url = https://github.com/gaoj66-roche/libvgio.git
 [submodule "deps/jemalloc"]
 	path = deps/jemalloc
 	url = https://github.com/jemalloc/jemalloc.git

src/alignment.cpp

@@ -689,6 +689,11 @@ string alignment_to_sam_internal(const Alignment& alignment,
     if (has_annotation(alignment, "nearest_ref_pos")) {
         sam << "\tNR:Z:" << get_annotation<string>(alignment, "nearest_ref_pos");
     }
+    if (has_annotation(alignment, "rescued_secondary")) {
+        if (get_annotation<bool>(alignment, "rescued_secondary")) {
+            sam << "\tYF:i:1";
+        }
+    }
 
     sam << "\n";
     return sam.str();
@@ -1041,14 +1046,20 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header,
         string pos = get_annotation<string>(alignment, "nearest_ref_pos");
         bam_aux_append(bam, "NR", 'Z', pos.size() + 1, (uint8_t*) pos.c_str());
     }
-
+    if (has_annotation(alignment, "rescued_secondary")) {
+        if (get_annotation<bool>(alignment, "rescued_secondary")) {
+            int32_t val = 1;
+            bam_aux_append(bam, "YF", 'i', sizeof(int32_t), (uint8_t*) &val);
+        }
+    }
+
     // TODO: it would be nice wrap htslib and set the other tags this way as well
     if (has_annotation(alignment, "tags")) {
         // encode the alignments SAM tags
         auto parsed_tags = parse_sam_tags(get_annotation<string>(alignment, "tags"));
         for (const auto& tag : parsed_tags) {
 
-            if (get<0>(tag) == "AS" || get<0>(tag) == "RG" || get<0>(tag) == "SS" || get<0>(tag) == "GR" || get<0>(tag) == "NR") {
+            if (get<0>(tag) == "AS" || get<0>(tag) == "RG" || get<0>(tag) == "SS" || get<0>(tag) == "GR" || get<0>(tag) == "NR" || get<0>(tag) == "YF") {
                 // we handle these tags separately
                 continue;
             }

src/hts_alignment_emitter.cpp

@@ -64,7 +64,8 @@ unique_ptr<AlignmentEmitter> get_alignment_emitter(const string& filename, const
                 flags & ALIGNMENT_EMITTER_FLAG_HTS_ADD_GRAPH_ALIGNMENT_TAG,
                 flags & ALIGNMENT_EMITTER_FLAG_HTS_SUPPLEMENTARY,
                 flags & ALIGNMENT_EMITTER_FLAG_HTS_OFF_REF_POSITION,
-                flags & ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN);
+                flags & ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN,
+                flags & ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY);
         }
 
     } else {

src/hts_alignment_emitter.hpp

@@ -53,7 +53,11 @@ enum alignment_emitter_flags_t {
     /// When surjecting, annote off-reference reads with their nearest reference position
     ALIGNMENT_EMITTER_FLAG_HTS_OFF_REF_POSITION = 64,
     /// When surjecting, attempt to left align
-    ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN = 128
+    ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN = 128,
+    /// When surjecting, if a read's primary alignment fails to surject, tag
+    /// the best-scoring secondary that does surject with YF:i:1 so downstream
+    /// tools may treat it as primary. The secondary flag is not changed.
+    ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY = 256
 };
 
 /// Represents a path or subpath's sequence dictionary information. Holds

src/multipath_mapper.cpp

@@ -535,7 +535,7 @@ namespace vg {
         cerr << "computing mapping quality and sorting mappings" << endl;
 #endif
         sort_and_compute_mapping_quality(multipath_alns_out, cluster_idxs, &multiplicities_out);
-        
+
         if (!multipath_alns_out.empty() && likely_mismapping(multipath_alns_out.front())) {
             multipath_alns_out.front().set_mapping_quality(0);
         }
@@ -2261,7 +2261,7 @@ namespace vg {
         // Now compute the MAPQ for the best alignment
         auto placement_mapqs = compute_raw_mapping_qualities_from_scores(scores,
                                                                          !multipath_aln.quality().empty());
-        // And min it in with what;s there already.
+        // And min it in with what's there already.
         alns_out[0].set_mapping_quality(min(alns_out[0].mapping_quality(), placement_mapqs.front()));
         for (size_t i = 1; i < alns_out.size(); i++) {
             // And zero all the others

src/subcommand/giraffe_main.cpp

@@ -732,6 +732,11 @@ void help_giraffe(char** argv, const BaseOptionGroup& parser, const std::map<std
              << "                                in the GR tag as a cs-style difference string" << endl
              << "      --off-ref-position        annotate off-reference mappings in HTSLib output" << endl
              << "                                with the nearest reference position in NR tag" << endl
+             << "      --rescue-secondary        in HTSlib output, tag the best secondary alignment" << endl
+             << "                                (tag with YF:i:1) when primary alignment is not" << endl
+             << "                                surjectable; downstream tools may treat secondary" << endl
+             << "                                alignments tagged with YF:i:1 as if they were" << endl
+             << "                                primary (needs --max-multimaps > 1)" << endl
              << "  -n, --discard                 discard all output alignments (for profiling)" << endl
              << "      --output-basename NAME    write output to a GAM file with the given prefix" << endl
              << "                                for each setting combination. Setting values for" << endl
@@ -797,6 +802,7 @@ int main_giraffe(int argc, char** argv) {
     constexpr int OPT_OFF_REF_POSITION = 1014;
     constexpr int OPT_LEFT_ALIGN = 1015;
     constexpr int OPT_NO_REC_MODE = 1016;
+    constexpr int OPT_RESCUE_SECONDARY = 1017;
     constexpr int OPT_HAPLOTYPE_NAME = 1100;
 
     constexpr int OPT_KFF_NAME = 1101;
@@ -908,6 +914,9 @@ int main_giraffe(int argc, char** argv) {
     // When surjecting, should we annotate the off-reference reads with the nearest reference position?
     bool annotate_off_ref_position = false;
 
+    // When surjecting, if a read's primary fails to surject, should we rescue its best surjectable secondary?
+    bool rescue_secondary = false;
+
     // For GAM format, should we report in named-segment space instead of node ID space?
     bool named_coordinates = false;
 
@@ -1166,6 +1175,7 @@ int main_giraffe(int argc, char** argv) {
         {"ref-name", required_argument, 0, OPT_REF_NAME},
         {"add-graph-aln", no_argument, 0, OPT_ADD_GRAPH_ALIGNMENT},
         {"off-ref-position", no_argument, 0, OPT_OFF_REF_POSITION},
+        {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY},
         {"left-align", no_argument, 0, OPT_LEFT_ALIGN},
         {"named-coordinates", no_argument, 0, OPT_NAMED_COORDINATES},
         {"discard", no_argument, 0, 'n'},
@@ -1357,6 +1367,10 @@ int main_giraffe(int argc, char** argv) {
             case OPT_OFF_REF_POSITION:
                 annotate_off_ref_position = true;
                 break;
+
+            case OPT_RESCUE_SECONDARY:
+                rescue_secondary = true;
+                break;
 
             case OPT_LEFT_ALIGN:
                 left_align = true;
@@ -2057,6 +2071,7 @@ int main_giraffe(int argc, char** argv) {
         report_flag("interleaved", interleaved);
         report_flag("add-graph-aln", add_graph_alignment);
         report_flag("off-ref-position", annotate_off_ref_position);
+        report_flag("rescue-secondary", rescue_secondary);
         report_flag("left-align", left_align);
         report_flag("set-refpos", set_refpos);
         minimizer_mapper.set_refpos = set_refpos;
@@ -2202,6 +2217,23 @@ int main_giraffe(int argc, char** argv) {
                     // When surjecting, attempt to left align
                     flags |= ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN;
                 }
+                if (rescue_secondary && minimizer_mapper.max_multimaps < 2) {
+                    logger.warn() << "--rescue-secondary requires --max-multimaps > 1; "
+                                  << "with the current setting only one alignment is produced per read, "
+                                  << "so there are no secondary alignments to rescue in case of an "
+                                  << "unsurjectable primary alignment. Ignoring." << endl;
+                    rescue_secondary = false;
+                }
+                if (rescue_secondary) {
+                    if (paired && !interleaved) {
+                        logger.warn() << "--rescue-secondary with two-file paired input uses paired-end "
+                                      << "promotion semantics: promotion fires only when both mates of the "
+                                      << "primary pair fail to surject, and the demoted pair is kept in the "
+                                      << "output as secondary records. Use -i if your input is interleaved." << endl;
+                    }
+                    // When surjecting, rescue a secondary if the primary fails to surject
+                    flags |= ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY;
+                }
 
                 // We send along the positional graph when we have it, and otherwise we send the GBWTGraph which is sufficient for GAF output.
                 // TODO: What if we need both a positional graph and a NamedNodeBackTranslation???

src/subcommand/map_main.cpp

@@ -132,6 +132,11 @@ void help_map(char** argv) {
          << "      --ref-paths FILE             ordered list of paths in graph, one per line" << endl
          << "                                   or HTSlib .dict, for HTSLib @SQ headers" << endl
          << "      --ref-name NAME              reference assembly in graph for HTSlib output" << endl
+         << "      --rescue-secondary          in HTSlib output, tag the best secondary alignment" << endl
+         << "                                   (tag with YF:i:1) when primary alignment is not" << endl
+         << "                                   surjectable; downstream tools may treat secondary" << endl
+         << "                                   alignments tagged with YF:i:1 as if they were" << endl
+         << "                                   primary (needs --max-multimaps > 1)" << endl
          << "  -X, --compare                    realign -G GAM input, writing alignment with" << endl
          << "                                   \"correct\" field set to overlap with input" << endl
          << "  -v, --refpos-table               for efficient testing output a table of" << endl
@@ -165,6 +170,7 @@ int main_map(int argc, char** argv) {
     constexpr int OPT_COMMENTS_AS_TAGS = 1005;
     constexpr int OPT_MAX_GAP_LENGTH = 1006;
     constexpr int OPT_XDROP_ALIGNMENT = 1007;
+    constexpr int OPT_RESCUE_SECONDARY = 1008;
     string matrix_file_name;
     string seq;
     string qual;
@@ -243,6 +249,7 @@ int main_map(int argc, char** argv) {
     uint32_t max_gap_length = 40;
     bool log_time = false;
     bool comments_as_tags = false;
+    bool rescue_secondary = false;
 
     int c;
     optind = 2; // force optind past command positional argument
@@ -322,6 +329,7 @@ int main_map(int argc, char** argv) {
                 {"gaf", no_argument, 0, '%'},
                 {"log-time", no_argument, 0, '^'},
                 {"comments-as-tags", no_argument, 0, OPT_COMMENTS_AS_TAGS},
+                {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY},
                 {"help", no_argument, 0, 'h'},
                 {0, 0, 0, 0}
             };
@@ -650,6 +658,10 @@ int main_map(int argc, char** argv) {
             comments_as_tags = true;
             break;
 
+        case OPT_RESCUE_SECONDARY:
+            rescue_secondary = true;
+            break;
+
         case 'h':
         case '?':
             /* getopt_long already printed an error message. */
@@ -800,9 +812,35 @@ int main_map(int argc, char** argv) {
         paths = get_sequence_dictionary(ref_paths_name, {}, reference_assembly_names, *xgidx);
     }
 
+    if (rescue_secondary && !hts_output) {
+        logger.warn() << "--rescue-secondary has no effect unless surjecting to SAM, BAM, or CRAM "
+                      << "(--surject-to); ignoring." << endl;
+        rescue_secondary = false;
+    }
+
+    if (rescue_secondary && max_multimaps < 2) {
+        logger.warn() << "--rescue-secondary requires --max-multimaps > 1; "
+                      << "with the current setting only one alignment is produced per read, "
+                      << "so there are no secondary alignments to rescue in case of an "
+                      << "unsurjectable primary alignment. Ignoring." << endl;
+        rescue_secondary = false;
+    }
+
+    if (rescue_secondary && !interleaved_input && !fastq2.empty()) {
+        logger.warn() << "--rescue-secondary with two-file paired input uses paired-end "
+                      << "promotion semantics: promotion fires only when both mates of the "
+                      << "primary pair fail to surject, and the demoted pair is kept in the "
+                      << "output as secondary records. Use -i if your input is interleaved." << endl;
+    }
+
     // Set up output to an emitter that will handle serialization and surjection
+    int emitter_flags = ALIGNMENT_EMITTER_FLAG_NONE;
+    if (rescue_secondary) {
+        // When surjecting, rescue a secondary if the primary fails to surject.
+        emitter_flags |= ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY;
+    }
     unique_ptr<vg::io::AlignmentEmitter> alignment_emitter = get_alignment_emitter("-", output_format, paths, 
-                                                                                   thread_count, xgidx);
+                                                                                   thread_count, xgidx, emitter_flags);
 
     // We have one function to dump alignments into
     auto output_alignments = [&](vector<Alignment>& alns1, vector<Alignment>& alns2) {