From 874dd1ec9f0560dd6e9570cede1a7730fd060c0a Mon Sep 17 00:00:00 2001 From: gaoj66-roche Date: Wed, 15 Jul 2026 21:47:24 +0000 Subject: [PATCH 1/5] Promote best surjectable secondary when primary fails to surject Add --promote-secondary to vg surject, giraffe, map, and mpmap. When a read's primary alignment fails to surject (becomes unmapped) but a secondary surjects, the best-scoring mapped, non-supplementary secondary is promoted to primary instead of emitting an unmapped primary alongside mapped secondaries. Only alignments already present are considered with --max-multimaps > 1. For standalone vg surject, input must be collated by read name. A grouped parallel reader keeps each read's primary and secondaries in one batch, never splitting a group across worker threads (may exceed the batch size, enabled only when promotion is on). A one-time warning explains when no surjectable secondary is found: input not collated, or reads generated without secondaries (--max-multimaps). Point the libvgio submodule at the fork branch carrying the group-aware parallel iterator. Promote best surjectable secondary pair when paired primary fails to surject Extend --promote-secondary to paired-end HTSlib output. When both mates of a read pair's primary alignment fail to surject but a secondary pair surjects, the secondary pair with the highest summed mapped-mate score is promoted into the primary slot (ties broken toward fully-mapped pairs), instead of emitting an unmapped primary pair alongside mapped secondaries. A candidate pair is eligible if at least one mate surjected; half-mapped pairs are used only when no fully-mapped secondary pair exists. Mates are kept index-aligned and is_secondary flags updated on both promoted and demoted pairs. Only alignments already present are considered; no realignment is performed. Supports subcommands that can produce SAM/BAM/CRAM output. vg surject with GAF input is not supported because GAF does not handle secondary alignments. --- .gitmodules | 2 +- deps/libvgio | 2 +- src/alignment.cpp | 15 +- src/hts_alignment_emitter.cpp | 3 +- src/hts_alignment_emitter.hpp | 6 +- src/mapping_quality_calculator.cpp | 16 ++ src/minimizer_mapper.cpp | 10 +- src/minimizer_mapper_from_chains.cpp | 6 +- src/multipath_mapper.cpp | 4 +- src/subcommand/giraffe_main.cpp | 30 +++ src/subcommand/map_main.cpp | 38 ++- src/subcommand/mpmap_main.cpp | 219 ++++++++++++++++++ src/subcommand/surject_main.cpp | 330 +++++++++++++++++++++------ src/surjecting_alignment_emitter.cpp | 12 +- src/surjecting_alignment_emitter.hpp | 2 +- src/surjector.cpp | 197 ++++++++++++++++ src/surjector.hpp | 37 ++- src/unittest/surject.cpp | 278 ++++++++++++++++++++++ test/t/07_vg_map.t | 4 +- test/t/15_vg_surject.t | 32 ++- test/t/33_vg_mpmap.t | 4 +- test/t/50_vg_giraffe.t | 4 +- 22 files changed, 1160 insertions(+), 91 deletions(-) diff --git a/.gitmodules b/.gitmodules index fcf4e6d764a..2cfa9c1def6 100644 --- a/.gitmodules +++ b/.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 diff --git a/deps/libvgio b/deps/libvgio index 2dad5c163b1..f1d64e34115 160000 --- a/deps/libvgio +++ b/deps/libvgio @@ -1 +1 @@ -Subproject commit 2dad5c163b1137f5cf0c43d84bfc2bdf67e82e7d +Subproject commit f1d64e3411579ef9f2218ec49e57e87fb4e5a34c diff --git a/src/alignment.cpp b/src/alignment.cpp index efc6d029f30..f5e0e6e16da 100644 --- a/src/alignment.cpp +++ b/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(alignment, "nearest_ref_pos"); } + if (has_annotation(alignment, "promoted_from_secondary")) { + if (get_annotation(alignment, "promoted_from_secondary")) { + sam << "\tps: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(alignment, "nearest_ref_pos"); bam_aux_append(bam, "NR", 'Z', pos.size() + 1, (uint8_t*) pos.c_str()); } - + if (has_annotation(alignment, "promoted_from_secondary")) { + if (get_annotation(alignment, "promoted_from_secondary")) { + int32_t val = 1; + bam_aux_append(bam, "ps", '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(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) == "ps") { // we handle these tags separately continue; } diff --git a/src/hts_alignment_emitter.cpp b/src/hts_alignment_emitter.cpp index c0b8f576ca8..233d86cf209 100644 --- a/src/hts_alignment_emitter.cpp +++ b/src/hts_alignment_emitter.cpp @@ -64,7 +64,8 @@ unique_ptr 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_PROMOTE_SECONDARY); } } else { diff --git a/src/hts_alignment_emitter.hpp b/src/hts_alignment_emitter.hpp index d047dc262fc..a8973cd634d 100644 --- a/src/hts_alignment_emitter.hpp +++ b/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, promote + /// the best-scoring secondary that does surject to be the new primary + /// instead of emitting an unmapped primary alongside mapped secondaries. + ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_SECONDARY = 256 }; /// Represents a path or subpath's sequence dictionary information. Holds diff --git a/src/mapping_quality_calculator.cpp b/src/mapping_quality_calculator.cpp index 29f12819d56..b21f058c918 100644 --- a/src/mapping_quality_calculator.cpp +++ b/src/mapping_quality_calculator.cpp @@ -271,6 +271,22 @@ void MappingQualityCalculator::compute_mapping_quality(vector& alignm for (size_t i = 1; i < alignments.size(); ++i) { alignments[0].add_secondary_score(alignments[i].score()); } + + // Compute meaningful MAPQs for all non-primary alignments using the same + // score vector. The primary-specific adjustments (identity scaling, cluster + // blending, mq_estimate cap) are not applied to secondaries. + if (alignments.size() > 1) { + vector raw_scores(alignments.size()); + for (size_t i = 0; i < alignments.size(); ++i) { + raw_scores[i] = alignments[i].score(); + } + vector all_mapqs = compute_all_mapping_qualities(raw_scores); + for (size_t i = 0; i < alignments.size(); ++i) { + if (i == max_idx) continue; + int32_t mq = (i < all_mapqs.size()) ? all_mapqs[i] : 0; + alignments[i].set_mapping_quality(min(mq, max_mapping_quality)); + } + } } void MappingQualityCalculator::compute_paired_mapping_quality(pair, vector>& alignment_pairs, diff --git a/src/minimizer_mapper.cpp b/src/minimizer_mapper.cpp index 5f0248d348c..7024c43d07c 100644 --- a/src/minimizer_mapper.cpp +++ b/src/minimizer_mapper.cpp @@ -1136,9 +1136,9 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { crash_unless(!mappings.empty()); // Compute MAPQ if not unmapped. Otherwise use 0 instead of the 50% this would give us. - // Use exact mapping quality - double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : - get_regular_aligner()->mapq_calc->compute_max_mapping_quality(scores, false) ; + // Use exact mapping quality + double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : + get_regular_aligner()->mapq_calc->compute_max_mapping_quality(scores, false); #ifdef print_minimizer_table double uncapped_mapq = mapq; @@ -1180,7 +1180,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { // Make sure to clamp 0-60. mappings.front().set_mapping_quality(max(min(mapq, 60.0), 0.0)); - + if (!supplementaries.empty()) { // Estimate a mapping quality for the supplementaries // TODO: only count the score of the overlapping portion of other alignments @@ -1204,7 +1204,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { for (size_t i = 0; i < mappings.size(); i++) { // For each output alignment in score order auto& out = mappings[i]; - + // Assign primary and secondary status out.set_is_secondary(i > 0); } diff --git a/src/minimizer_mapper_from_chains.cpp b/src/minimizer_mapper_from_chains.cpp index 27325ab379b..a4190713604 100644 --- a/src/minimizer_mapper_from_chains.cpp +++ b/src/minimizer_mapper_from_chains.cpp @@ -987,8 +987,8 @@ vector MinimizerMapper::map_from_chains(Alignment& aln) { // Because the winning alignment won't necessarily *always* have the // maximum score, we need to use compute_first_mapping_quality and not // compute_max_mapping_quality. - double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : - get_regular_aligner()->mapq_calc->compute_first_mapping_quality(scaled_scores, false, &multiplicity_by_alignment) ; + double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : + get_regular_aligner()->mapq_calc->compute_first_mapping_quality(scaled_scores, false, &multiplicity_by_alignment); #ifdef debug_write_minimizers #pragma omp critical @@ -1106,7 +1106,7 @@ vector MinimizerMapper::map_from_chains(Alignment& aln) { for (size_t i = 0; i < mappings.size(); i++) { // For each output alignment in score order auto& out = mappings[i]; - + // Assign primary and secondary status out.set_is_secondary(i > 0); } diff --git a/src/multipath_mapper.cpp b/src/multipath_mapper.cpp index f90214123ca..575ac77d17f 100644 --- a/src/multipath_mapper.cpp +++ b/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 diff --git a/src/subcommand/giraffe_main.cpp b/src/subcommand/giraffe_main.cpp index 3dea6124aa9..7198f23758f 100644 --- a/src/subcommand/giraffe_main.cpp +++ b/src/subcommand/giraffe_main.cpp @@ -732,6 +732,9 @@ void help_giraffe(char** argv, const BaseOptionGroup& parser, const std::map 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 +800,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_PROMOTE_SECONDARY = 1017; constexpr int OPT_HAPLOTYPE_NAME = 1100; constexpr int OPT_KFF_NAME = 1101; @@ -908,6 +912,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 promote its best surjectable secondary? + bool promote_secondary = false; + // For GAM format, should we report in named-segment space instead of node ID space? bool named_coordinates = false; @@ -1166,6 +1173,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}, + {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, {"left-align", no_argument, 0, OPT_LEFT_ALIGN}, {"named-coordinates", no_argument, 0, OPT_NAMED_COORDINATES}, {"discard", no_argument, 0, 'n'}, @@ -1357,6 +1365,10 @@ int main_giraffe(int argc, char** argv) { case OPT_OFF_REF_POSITION: annotate_off_ref_position = true; break; + + case OPT_PROMOTE_SECONDARY: + promote_secondary = true; + break; case OPT_LEFT_ALIGN: left_align = true; @@ -2057,6 +2069,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("promote-secondary", promote_secondary); report_flag("left-align", left_align); report_flag("set-refpos", set_refpos); minimizer_mapper.set_refpos = set_refpos; @@ -2202,6 +2215,23 @@ int main_giraffe(int argc, char** argv) { // When surjecting, attempt to left align flags |= ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN; } + if (promote_secondary && minimizer_mapper.max_multimaps < 2) { + logger.warn() << "--promote-secondary requires --max-multimaps > 1; " + << "with the current setting only one alignment is produced per read, " + << "so there are no secondary alignments to promote in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + promote_secondary = false; + } + if (promote_secondary) { + if (paired && !interleaved) { + logger.warn() << "--promote-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, promote a mapped secondary if the primary fails to surject + flags |= ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_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??? diff --git a/src/subcommand/map_main.cpp b/src/subcommand/map_main.cpp index cc1114fc2fb..c36b4123d86 100644 --- a/src/subcommand/map_main.cpp +++ b/src/subcommand/map_main.cpp @@ -132,6 +132,9 @@ 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 + << " --promote-secondary in HTSlib output, if a read's primary fails to" << endl + << " surject, promote its best surjectable secondary" << endl + << " to 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 +168,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_PROMOTE_SECONDARY = 1008; string matrix_file_name; string seq; string qual; @@ -243,6 +247,7 @@ int main_map(int argc, char** argv) { uint32_t max_gap_length = 40; bool log_time = false; bool comments_as_tags = false; + bool promote_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -322,6 +327,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}, + {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, {"help", no_argument, 0, 'h'}, {0, 0, 0, 0} }; @@ -650,6 +656,10 @@ int main_map(int argc, char** argv) { comments_as_tags = true; break; + case OPT_PROMOTE_SECONDARY: + promote_secondary = true; + break; + case 'h': case '?': /* getopt_long already printed an error message. */ @@ -800,9 +810,35 @@ int main_map(int argc, char** argv) { paths = get_sequence_dictionary(ref_paths_name, {}, reference_assembly_names, *xgidx); } + if (promote_secondary && !hts_output) { + logger.warn() << "--promote-secondary has no effect unless surjecting to SAM, BAM, or CRAM " + << "(--surject-to); ignoring." << endl; + promote_secondary = false; + } + + if (promote_secondary && max_multimaps < 2) { + logger.warn() << "--promote-secondary requires --max-multimaps > 1; " + << "with the current setting only one alignment is produced per read, " + << "so there are no secondary alignments to promote in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + promote_secondary = false; + } + + if (promote_secondary && !interleaved_input && !fastq2.empty()) { + logger.warn() << "--promote-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 (promote_secondary) { + // When surjecting, promote a mapped secondary if the primary fails to surject. + emitter_flags |= ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_SECONDARY; + } unique_ptr 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& alns1, vector& alns2) { diff --git a/src/subcommand/mpmap_main.cpp b/src/subcommand/mpmap_main.cpp index b9e2361f5f1..41c4da1f7ef 100644 --- a/src/subcommand/mpmap_main.cpp +++ b/src/subcommand/mpmap_main.cpp @@ -98,6 +98,182 @@ static void error_if_negative(const Logger& logger, double value, const string& } } +/// Returns true if the given surjected multipath alignment is mapped (has at +/// least one aligned base). +static bool mp_aln_is_mapped(const multipath_alignment_t& mp_aln) { + for (size_t i = 0; i < mp_aln.subpath_size(); ++i) { + if (mp_aln.subpath(i).path().mapping_size() > 0) { + return true; + } + } + return false; +} + +/// Apply "promote secondary on failed surjection" to a single read's surjected +/// multipath alignments, keeping the parallel path_positions array in sync. +/// +/// mp_alns holds the read's surjected primary followed by its secondaries (and +/// possibly supplementaries appended at the end); path_positions holds the +/// matching (path name, is_reverse, offset) for each. If the primary (the first +/// non-supplementary entry) is unmapped, the best-scoring mapped, +/// non-supplementary secondary is promoted into its place (both the alignment +/// and its position). Only alignments already present are considered; no +/// realignment is performed. Supplementaries are never promoted. +/// +/// warned is a shared flag used to emit the "no promotable secondary" warning +/// only once. Returns true iff a promotion occurred. +static bool promote_secondary_mp(vector& mp_alns, + vector>& path_positions, + const string& read_name, + std::atomic_flag& warned) { + if (mp_alns.size() < 2) { + return false; + } + + // Find the primary: first non-supplementary entry. + int64_t primary_idx = -1; + for (size_t i = 0; i < mp_alns.size(); ++i) { + if (!is_supplementary(mp_alns[i])) { + primary_idx = (int64_t) i; + break; + } + } + if (primary_idx < 0 || mp_aln_is_mapped(mp_alns[primary_idx])) { + // No primary, or primary surjected fine. + return false; + } + + // Find the best-scoring mapped, non-supplementary secondary to promote. + int64_t best_idx = -1; + int32_t best_score = numeric_limits::min(); + for (size_t i = 0; i < mp_alns.size(); ++i) { + if ((int64_t) i == primary_idx || is_supplementary(mp_alns[i])) { + continue; + } + if (!mp_aln_is_mapped(mp_alns[i])) { + continue; + } + int32_t score = optimal_alignment_score(mp_alns[i], true); + if (score > best_score || + (score == best_score && best_idx >= 0 && + mp_alns[i].mapping_quality() > mp_alns[best_idx].mapping_quality())) { + best_score = score; + best_idx = (int64_t) i; + } + } + + if (best_idx < 0) { + // No mapped secondary to promote. Warn once. + if (!warned.test_and_set()) { + #pragma omp critical (cerr) + { + cerr << "warning:[vg mpmap] --promote-secondary was requested, but a read whose " + << "primary alignment failed to surject had no surjectable secondary alignment " + << "to promote (first seen for read \"" << read_name << "\"). This can happen " + << "because the reads were not generated with secondaries (run mpmap with " + << "--max-multimaps > 1). Affected reads are left with an unmapped primary. " + << "This warning is shown only once." << endl; + } + } + return false; + } + + // Promote: swap the chosen secondary into the primary slot, mark it primary, + // and mark the demoted (unmapped) alignment secondary. Keep positions in sync. + mp_alns[best_idx].set_annotation("secondary", false); + mp_alns[primary_idx].set_annotation("secondary", true); + mp_alns[best_idx].set_annotation("promoted_from_secondary", true); + std::swap(mp_alns[primary_idx], mp_alns[best_idx]); + std::swap(path_positions[primary_idx], path_positions[best_idx]); + return true; +} + +/// Paired counterpart of promote_secondary_mp for mpmap's surjected output. +/// +/// output_mp_aln_pairs holds one surjected pair per multimapping (index 0 is the +/// primary pair, later indices are secondary pairs); path_positions is the +/// matching parallel array of (mate1 pos, mate2 pos). If BOTH mates of the +/// primary pair are unmapped, the secondary pair with the highest summed +/// mapped-mate score is promoted into index 0 (ties broken toward fully-mapped +/// pairs), keeping path_positions in sync and updating the "secondary" +/// annotations. A candidate is eligible if at least one of its mates is mapped. +/// Only alignments already present are considered; no realignment is performed. +/// +/// warned is shared so the "no promotable secondary" warning is emitted once. +/// Returns true iff a promotion occurred. +static bool promote_secondary_pair_mp(vector>& output_mp_aln_pairs, + vector, tuple>>& path_positions, + const string& read_name, + std::atomic_flag& warned) { + if (output_mp_aln_pairs.size() < 2 || output_mp_aln_pairs.size() != path_positions.size()) { + return false; + } + + // Only act when BOTH mates of the primary pair failed to surject. + if (mp_aln_is_mapped(output_mp_aln_pairs[0].first) || + mp_aln_is_mapped(output_mp_aln_pairs[0].second)) { + return false; + } + + int64_t best_idx = -1; + int32_t best_score = numeric_limits::min(); + bool best_fully_mapped = false; + int32_t best_mapq = numeric_limits::min(); + for (size_t k = 1; k < output_mp_aln_pairs.size(); ++k) { + bool m1 = mp_aln_is_mapped(output_mp_aln_pairs[k].first); + bool m2 = mp_aln_is_mapped(output_mp_aln_pairs[k].second); + if (!m1 && !m2) { + continue; + } + int32_t score = (m1 ? optimal_alignment_score(output_mp_aln_pairs[k].first, true) : 0) + + (m2 ? optimal_alignment_score(output_mp_aln_pairs[k].second, true) : 0); + bool fully_mapped = m1 && m2; + int32_t mapq = (m1 ? output_mp_aln_pairs[k].first.mapping_quality() : 0) + + (m2 ? output_mp_aln_pairs[k].second.mapping_quality() : 0); + bool better = false; + if (score != best_score) { + better = score > best_score; + } + else if (fully_mapped != best_fully_mapped) { + better = fully_mapped; + } + else { + better = mapq > best_mapq; + } + if (best_idx < 0 || better) { + best_idx = (int64_t) k; + best_score = score; + best_fully_mapped = fully_mapped; + best_mapq = mapq; + } + } + + if (best_idx < 0) { + if (!warned.test_and_set()) { + #pragma omp critical (cerr) + { + cerr << "warning:[vg mpmap] --promote-secondary was requested, but a read pair whose " + << "primary alignment failed to surject had no surjectable secondary alignment " + << "to promote (first seen for read \"" << read_name << "\"). This can happen " + << "because the reads were not generated with secondaries (run mpmap with " + << "--max-multimaps > 1). Affected read pairs are left with an unmapped primary. " + << "This warning is shown only once." << endl; + } + } + return false; + } + + output_mp_aln_pairs[best_idx].first.set_annotation("secondary", false); + output_mp_aln_pairs[best_idx].second.set_annotation("secondary", false); + output_mp_aln_pairs[0].first.set_annotation("secondary", true); + output_mp_aln_pairs[0].second.set_annotation("secondary", true); + output_mp_aln_pairs[best_idx].first.set_annotation("promoted_from_secondary", true); + output_mp_aln_pairs[best_idx].second.set_annotation("promoted_from_secondary", true); + std::swap(output_mp_aln_pairs[0], output_mp_aln_pairs[best_idx]); + std::swap(path_positions[0], path_positions[best_idx]); + return true; +} + void help_mpmap(char** argv) { cerr << "usage: " << argv[0] << " mpmap [options] -x graph.xg -g index.gcsa " << "[-f reads1.fq [-f reads2.fq] | -G reads.gam] > aln.gamp" << endl @@ -142,6 +318,9 @@ void help_mpmap(char** argv) { << " [all reference paths, all generic paths]" << endl << " --ref-name NAME reference assembly in graph to use for" << endl << " HTSlib formats (see -F) [all references]" << endl + << " --promote-secondary in HTSlib output, if a read's primary fails to surject," << endl + << " promote its best surjectable secondary to primary" << endl + << " (needs --max-multimaps > 1)" << endl << " -N, --sample NAME add this sample name to output" << endl << " -R, --read-group NAME add this read group to output" << endl << " -p, --suppress-progress do not report progress to stderr" << endl @@ -271,6 +450,7 @@ int main_mpmap(int argc, char** argv) { constexpr int OPT_REF_NAME = 1039; constexpr int OPT_LINEAR_PATH = 1040; constexpr int OPT_LINEAR_INDEX = 1041; + constexpr int OPT_PROMOTE_SECONDARY = 1042; string matrix_file_name; string graph_name; string gcsa_name; @@ -357,6 +537,7 @@ int main_mpmap(int argc, char** argv) { int default_num_alt_alns = 16; int num_alt_alns = default_num_alt_alns; bool agglomerate_multipath_alns = false; + bool promote_secondary = false; double suboptimal_path_exponent = 1.25; double likelihood_approx_exp = 10.0; double likelihood_approx_exp_arg = numeric_limits::lowest(); @@ -456,6 +637,7 @@ int main_mpmap(int argc, char** argv) { {"same-strand", no_argument, 0, 'T'}, {"ref-paths", required_argument, 0, 'S'}, {"ref-name", required_argument, 0, OPT_REF_NAME}, + {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, {"output-fmt", required_argument, 0, 'F'}, {"snarls", required_argument, 0, 's'}, {"synth-tail-anchors", no_argument, 0, OPT_SUPPRESS_TAIL_ANCHORS}, @@ -629,6 +811,10 @@ int main_mpmap(int argc, char** argv) { case OPT_REF_NAME: reference_assembly_names.insert(optarg); break; + + case OPT_PROMOTE_SECONDARY: + promote_secondary = true; + break; case 's': snarls_name = require_exists(logger, optarg); @@ -1190,6 +1376,20 @@ int main_mpmap(int argc, char** argv) { << "when output format (-F) is SAM, BAM, or CRAM." << endl; ref_paths_name = ""; } + + if (promote_secondary && !hts_output) { + logger.warn() << "--promote-secondary has no effect unless output format (-F) is " + << "SAM, BAM, or CRAM; ignoring." << endl; + promote_secondary = false; + } + + if (promote_secondary && num_alt_alns < 2) { + logger.warn() << "--promote-secondary requires --alt-paths > 1; " + << "with the current setting only one alignment is produced per read, " + << "so there are no secondary alignments to promote in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + promote_secondary = false; + } if (!reference_assembly_names.empty() && !hts_output) { logger.warn() << "Reference assembly names (--ref-name) are only used " @@ -1957,6 +2157,9 @@ int main_mpmap(int argc, char** argv) { // during distribution estimation vector> ambiguous_pair_buffer; + // Shared flag so the "no promotable secondary" warning is emitted only once. + std::atomic_flag warned_no_promotable_secondary = ATOMIC_FLAG_INIT; + // do unpaired multipath alignment and write to buffer function do_unpaired_alignments = [&](Alignment& alignment) { #ifdef record_read_run_times @@ -2001,6 +2204,14 @@ int main_mpmap(int argc, char** argv) { get<2>(suppl_positions[j]), get<1>(suppl_positions[j])); } } + + if (promote_secondary) { + // If the primary failed to surject, promote the best mapped + // secondary in its place (positions kept in sync). Warns once if + // no surjectable secondary is available. + promote_secondary_mp(mp_alns, path_positions, alignment.name(), + warned_no_promotable_secondary); + } } if (is_rna) { @@ -2147,6 +2358,14 @@ int main_mpmap(int argc, char** argv) { } } } + + if (promote_secondary) { + // If both mates of the primary pair failed to surject, promote + // the best surjectable secondary pair into its place (positions + // kept in sync). Warns once if no surjectable secondary exists. + promote_secondary_pair_mp(output_mp_aln_pairs, path_positions, + alignment_1.name(), warned_no_promotable_secondary); + } } else { output_mp_aln_pairs = std::move(mp_aln_pairs); diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 54adfb0c176..4b7ddc1ac95 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -21,6 +21,8 @@ #include "../xg.hpp" #include #include +#include +#include #include "../utility.hpp" #include "../surjector.hpp" #include "../hts_alignment_emitter.hpp" @@ -86,6 +88,9 @@ void help_surject(char** argv) { << " of the pre-surjected graph alignment in GR tag" << endl << " --off-ref-position annotate SAM records that become unmapped during" << endl << " surject with the nearest ref. position in the NR tag" << endl + << " --promote-secondary if a read's primary fails to surject, promote its best" << endl + << " surjectable secondary to primary instead of emitting an" << endl + << " unmapped primary (input must be collated by read name)" << endl << " -C, --compression N level for compression [0-9]" << endl << " -V, --no-validate skip checking whether alignments plausibly are" << endl << " against the provided graph" << endl @@ -149,6 +154,7 @@ int main_surject(int argc, char** argv) { constexpr int OPT_NO_PRUNE_LOW_CPLX = 1000; constexpr int OPT_OFF_REF_POS = 1001; + constexpr int OPT_PROMOTE_SECONDARY = 1002; if (argc == 2) { help_surject(argv); @@ -187,6 +193,7 @@ int main_surject(int argc, char** argv) { bool validate = true; bool show_progress = false; bool left_align = false; + bool promote_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -213,6 +220,7 @@ int main_surject(int argc, char** argv) { {"sam-output", no_argument, 0, 's'}, {"supplementary", no_argument, 0, 'u'}, {"off-ref-position", no_argument, 0, OPT_OFF_REF_POS}, + {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, {"left-align", no_argument, 0, 'B'}, {"read-length", required_argument, 0, 'D'}, {"spliced", no_argument, 0, 'S'}, @@ -386,6 +394,10 @@ int main_surject(int argc, char** argv) { annotate_off_reference_pos = true; break; + case OPT_PROMOTE_SECONDARY: + promote_secondary = true; + break; + case 'h': case '?': help_surject(argv); @@ -504,6 +516,23 @@ int main_surject(int argc, char** argv) { surjector.report_supplementary = report_supplementary; surjector.left_align = left_align; surjector.multimap_to_all_paths = multimap; + surjector.promote_secondary_on_failed_surjection = promote_secondary; + surjector.warn_about_input_collation = true; + + bool hts_output = (output_format == "SAM" || output_format == "BAM" || output_format == "CRAM"); + if (promote_secondary && !hts_output) { + logger.warn() << "--promote-secondary has no effect unless output is SAM, BAM, or CRAM; " + << "ignoring." << endl; + promote_secondary = false; + surjector.promote_secondary_on_failed_surjection = false; + } + if (promote_secondary && input_format == "GAF") { + logger.warn() << "--promote-secondary is not supported with GAF input because GAF does not " + << "distinguish primary and secondary alignments; ignoring." << endl; + promote_secondary = false; + surjector.promote_secondary_on_failed_surjection = false; + } + // Count our threads int thread_count = vg::get_thread_count(); @@ -542,7 +571,178 @@ int main_surject(int argc, char** argv) { output_format, sequence_dictionary, thread_count, xgidx, ALIGNMENT_EMITTER_FLAG_HTS_RAW | (spliced * ALIGNMENT_EMITTER_FLAG_HTS_SPLICED)); - if (interleaved) { + // Emit one already-surjected read pair's multimappings, pairing mates by + // strand and handling supplementaries/unpaired mates. surjected1[k] and + // surjected2[k] need NOT be index-aligned here: mates are matched by + // reference strand, exactly as in the default interleaved path. + auto emit_surjected_pair = [&](vector& surjected1, vector& surjected2) { + // pair up non-supplementary alignments + unordered_map, size_t> strand_idx1, strand_idx2; + for (size_t i = 0; i < surjected1.size(); ++i) { + if (!is_supplementary(surjected1[i])) { + const auto& pos = surjected1[i].refpos(0); + strand_idx1[make_pair(pos.name(), pos.is_reverse())] = i; + } + } + for (size_t i = 0; i < surjected2.size(); ++i) { + if (!is_supplementary(surjected2[i])) { + const auto& pos = surjected2[i].refpos(0); + strand_idx2[make_pair(pos.name(), pos.is_reverse())] = i; + } + } + for (size_t i = 0; i < surjected1.size(); ++i) { + const auto& pos = surjected1[i].refpos(0); + auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse())); + if (!is_supplementary(surjected1[i]) && it != strand_idx2.end()) { + alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len); + } + else { + if (is_supplementary(surjected1[i]) && !has_annotation(surjected1[i], "mate_info")) { + string annotation; + if (!strand_idx2.empty()) { + const auto& mate = it != strand_idx2.end() ? surjected2[it->second] : surjected2[strand_idx2.begin()->second]; + annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), false)); + } + else { + annotation = std::move(mate_info("", -1, false, false)); + } + set_annotation(surjected1[i], "mate_info", annotation); + } + alignment_emitter->emit_single(std::move(surjected1[i])); + } + } + for (size_t i = 0; i < surjected2.size(); ++i) { + const auto& pos = surjected2[i].refpos(0); + auto it = strand_idx1.find(make_pair(pos.name(), !pos.is_reverse())); + if (is_supplementary(surjected2[i]) || it == strand_idx1.end()) { + if (is_supplementary(surjected2[i]) && !has_annotation(surjected2[i], "mate_info")) { + string annotation; + if (!strand_idx1.empty()) { + const auto& mate = it != strand_idx1.end() ? surjected1[it->second] : surjected1[strand_idx1.begin()->second]; + annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), true)); + } + else { + annotation = std::move(mate_info("", -1, false, true)); + } + set_annotation(surjected2[i], "mate_info", annotation); + } + alignment_emitter->emit_single(std::move(surjected2[i])); + } + } + }; + + if (interleaved && promote_secondary) { + // Paired secondary promotion needs a read pair's whole multimapping + // (its primary pair followed by its secondary pairs) together. We + // read grouped pairs (consecutive pair-records sharing a read name) + // so a primary pair and its secondary pairs are never split across + // worker threads. This requires the input to be collated by read + // name. If both mates of the primary pair fail to surject, the best + // surjectable secondary pair is promoted into its place. + using AlnPair = pair; + function&)> process_pair_group = [&](vector& group) { + if (group.empty()) { + return; + } + try { + set_crash_context(group.front().first.name()); + size_t thread_num = omp_get_thread_num(); + if (watchdog) { + watchdog->check_in(thread_num, group.front().first.name()); + } + // Surject every pair-record, building index-aligned mate + // vectors: promoted1[k] and promoted2[k] are the two mates of + // the same alignment pair, index 0 being the primary pair. + // Supplementaries are collected separately and emitted after. + vector promoted1, promoted2; + vector extra1, extra2; + for (auto& rec : group) { + Alignment& src1 = rec.first; + Alignment& src2 = rec.second; + if (validate) { + ensure_alignment_is_for_graph(logger, src1, *xgidx); + ensure_alignment_is_for_graph(logger, src2, *xgidx); + } + set_metadata(src1); + set_metadata(src2); + auto s1 = surjector.surject(src1, paths, subpath_global, spliced); + auto s2 = surjector.surject(src2, paths, subpath_global, spliced); + // Keep the non-supplementary (primary of this record) mate + // for pair-level promotion; route supplementaries to extra. + int64_t keep1 = -1, keep2 = -1; + for (size_t j = 0; j < s1.size(); ++j) { + if (keep1 < 0 && !is_supplementary(s1[j])) { + keep1 = (int64_t) j; + } else { + extra1.emplace_back(std::move(s1[j])); + } + } + for (size_t j = 0; j < s2.size(); ++j) { + if (keep2 < 0 && !is_supplementary(s2[j])) { + keep2 = (int64_t) j; + } else { + extra2.emplace_back(std::move(s2[j])); + } + } + // Fall back to a null (unmapped) alignment if somehow + // all results were supplementary, so downstream refpos + // access stays valid. + auto null_with_refpos = [&](const Alignment& src) { + Alignment a; + a.set_name(src.name()); + a.set_sequence(src.sequence()); + a.set_quality(src.quality()); + a.add_refpos(); + return a; + }; + promoted1.emplace_back(keep1 >= 0 ? std::move(s1[keep1]) : null_with_refpos(src1)); + promoted2.emplace_back(keep2 >= 0 ? std::move(s2[keep2]) : null_with_refpos(src2)); + } + // Promote the best surjectable secondary pair if both mates of + // the primary pair failed to surject. Warns once if none. + surjector.promote_secondary_pair_if_primary_unmapped(promoted1, promoted2); + // Emit the (possibly reordered) pairs plus any supplementaries. + for (auto& a : extra1) { + promoted1.emplace_back(std::move(a)); + } + for (auto& a : extra2) { + promoted2.emplace_back(std::move(a)); + } + emit_surjected_pair(promoted1, promoted2); + total_reads_surjected += 2 * group.size(); + if (watchdog) { + watchdog->check_out(thread_num); + } + clear_crash_context(); + } catch (const std::exception& ex) { + report_exception(ex); + } + }; + // Two pair-records belong to the same group if they share a read name. + function pairs_in_same_group = + [](const AlnPair& a, const AlnPair& b) { + return a.first.name() == b.first.name(); + }; + if (input_format == "GAM") { + get_input_file(file_name, [&](istream& in) { + vg::io::ProtobufIterator cursor(in); + function get_pair = [&](AlnPair& dest) { + if (!cursor.has_current()) { + return false; + } + dest.first = std::move(cursor.take()); + if (!cursor.has_current()) { + // Odd number of records in an interleaved GAM. + adjacent_but_not_paired_error(logger, dest.first.name(), ""); + return false; + } + dest.second = std::move(cursor.take()); + return true; + }; + vg::io::grouped_unpaired_for_each_parallel(get_pair, process_pair_group, pairs_in_same_group); + }); + } + } else if (interleaved) { // GAM input is paired, and for HTS output reads need to know their pair partners' mapping locations. // TODO: We don't preserve order relationships (like primary/secondary) beyond the interleaving. function lambda = [&](Alignment& src1, Alignment& src2) { @@ -594,72 +794,11 @@ int main_surject(int argc, char** argv) { // Surject auto surjected1 = surjector.surject(src1, paths, subpath_global, spliced); auto surjected2 = surjector.surject(src2, paths, subpath_global, spliced); - - // pair up non-supplementary alignments - unordered_map, size_t> strand_idx1, strand_idx2; - for (size_t i = 0; i < surjected1.size(); ++i) { - if (!is_supplementary(surjected1[i])) { - const auto& pos = surjected1[i].refpos(0); - strand_idx1[make_pair(pos.name(), pos.is_reverse())] = i; - } - } - for (size_t i = 0; i < surjected2.size(); ++i) { - if (!is_supplementary(surjected2[i])) { - const auto& pos = surjected2[i].refpos(0); - strand_idx2[make_pair(pos.name(), pos.is_reverse())] = i; - } - } - - for (size_t i = 0; i < surjected1.size(); ++i) { - const auto& pos = surjected1[i].refpos(0); - auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse())); - if (!is_supplementary(surjected1[i]) && it != strand_idx2.end()) { - // the alignments are paired on this strand - alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len); - } - else { - // supplementary or unpaired - if (is_supplementary(surjected1[i]) && !has_annotation(surjected1[i], "mate_info")) { - // we need to annotate this supplementary with mate info for SAM/BAM conversion - string annotation; - if (!strand_idx2.empty()) { - // there is a non-supplementary alignment available (prefer the one consistent with this path strand) - const auto& mate = it != strand_idx2.end() ? surjected2[it->second] : surjected2[strand_idx2.begin()->second]; - annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), false)); - } - else { - // we don't have access to the primary, but we can still record the read 1/2 identity - annotation = std::move(mate_info("", -1, false, false)); - } - set_annotation(surjected1[i], "mate_info", annotation); - } - alignment_emitter->emit_single(std::move(surjected1[i])); - } - } - for (size_t i = 0; i < surjected2.size(); ++i) { - const auto& pos = surjected2[i].refpos(0); - auto it = strand_idx1.find(make_pair(pos.name(), !pos.is_reverse())); - if (is_supplementary(surjected2[i]) || it == strand_idx1.end()) { - // this strand's surjection is unpaired or supplementary - if (is_supplementary(surjected2[i]) && !has_annotation(surjected2[i], "mate_info")) { - // we need to annotate this supplementary with mate info for SAM/BAM conversion - string annotation; - if (!strand_idx1.empty()) { - // there is a non-supplementary alignment available (prefer the one consistent with this path strand) - const auto& mate = it != strand_idx1.end() ? surjected1[it->second] : surjected1[strand_idx1.begin()->second]; - annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), true)); - } - else { - // we don't have access to the primary, but we can still record the read 1/2 identity - annotation = std::move(mate_info("", -1, false, true)); - } - set_annotation(surjected2[i], "mate_info", annotation); - } - alignment_emitter->emit_single(std::move(surjected2[i])); - } - } - + // Pair up mates by strand and emit (shared with the + // promotion path). + emit_surjected_pair(surjected1, surjected2); + total_reads_surjected += 2; if (watchdog) { watchdog->check_out(thread_num); @@ -681,6 +820,65 @@ int main_surject(int argc, char** argv) { }; vg::io::gaf_paired_interleaved_for_each_parallel(*xgidx, file_name, gaf_checking_lambda); } + } else if (promote_secondary) { + // Secondary promotion needs a read's whole multimapping (primary + // plus its secondaries) together so it can promote a mapped + // secondary if the primary fails to surject. We use a grouped + // parallel reader that keeps consecutive same-named reads in one + // group and never splits a group across worker threads. This + // requires the input to be collated by read name (as produced by + // vg giraffe / map / mpmap with --max-multimaps > 1). + function&)> process_group = [&](vector& group) { + if (group.empty()) { + return; + } + try { + set_crash_context(group.front().name()); + size_t thread_num = omp_get_thread_num(); + if (watchdog) { + watchdog->check_in(thread_num, group.front().name()); + } + // Surject every member of the group, collecting all results. + vector surjected_group; + for (auto& src : group) { + if (validate) { + ensure_alignment_is_for_graph(logger, src, *xgidx); + } + set_metadata(src); + auto surjected = surjector.surject(src, paths, subpath_global, spliced); + for (auto& s : surjected) { + surjected_group.emplace_back(std::move(s)); + } + } + // Promote the best mapped secondary if the primary failed to + // surject. Emits a one-time warning if none is available. + surjector.promote_secondary_if_primary_unmapped(surjected_group); + alignment_emitter->emit_singles(std::move(surjected_group)); + total_reads_surjected += group.size(); + if (watchdog) { + watchdog->check_out(thread_num); + } + clear_crash_context(); + } catch (const std::exception& ex) { + report_exception(ex); + } + }; + get_input_file(file_name, [&](istream& in) { + // Single-threaded reader; group boundaries are on read name. + vg::io::ProtobufIterator cursor(in); + function get_read = [&](Alignment& dest) { + if (!cursor.has_current()) { + return false; + } + dest = std::move(cursor.take()); + return true; + }; + function in_same_group = + [](const Alignment& a, const Alignment& b) { + return a.name() == b.name(); + }; + vg::io::grouped_unpaired_for_each_parallel(get_read, process_group, in_same_group); + }); } else { // We can just surject each Alignment by itself. // TODO: We don't preserve order relationships (like primary/secondary). diff --git a/src/surjecting_alignment_emitter.cpp b/src/surjecting_alignment_emitter.cpp index af5a789e4a3..5f564d118a9 100644 --- a/src/surjecting_alignment_emitter.cpp +++ b/src/surjecting_alignment_emitter.cpp @@ -15,7 +15,7 @@ using namespace std; SurjectingAlignmentEmitter::SurjectingAlignmentEmitter(const PathPositionHandleGraph* graph, unordered_set paths, unique_ptr&& backing, bool prune_suspicious_anchors, bool add_graph_alignment_tag, bool report_supplementary, - bool add_off_ref_position_tag, bool left_align) : surjector(graph), paths(paths), backing(std::move(backing)) { + bool add_off_ref_position_tag, bool left_align, bool promote_secondary) : surjector(graph), paths(paths), backing(std::move(backing)) { // Configure the surjector surjector.prune_suspicious_anchors = prune_suspicious_anchors; @@ -23,6 +23,7 @@ SurjectingAlignmentEmitter::SurjectingAlignmentEmitter(const PathPositionHandleG surjector.report_supplementary = report_supplementary; surjector.annotate_off_reference_pos = add_off_ref_position_tag; surjector.left_align = left_align; + surjector.promote_secondary_on_failed_surjection = promote_secondary; } void SurjectingAlignmentEmitter::surject_alignments_in_place(vector& alns) const { @@ -104,6 +105,10 @@ void SurjectingAlignmentEmitter::emit_mapped_singles(vector>&& for (auto& mappings : alns_batch_caught) { // Surject all mappings in place surject_alignments_in_place(mappings); + // Each inner vector holds one read's whole multimapping (its primary + // followed by its secondaries), so this is the right granularity to + // promote a mapped secondary if the primary failed to surject. + surjector.promote_secondary_if_primary_unmapped(mappings); } // Forward it along backing->emit_mapped_singles(std::move(alns_batch_caught)); @@ -129,6 +134,11 @@ void SurjectingAlignmentEmitter::emit_mapped_pairs(vector>&& a for (size_t i = 0; i < alns1_batch_caught.size(); ++i) { supplementary_batch.emplace_back(); surject_paired_alignments_in_place(alns1_batch_caught[i], alns2_batch_caught[i], supplementary_batch.back()); + // After surjection, alns1_batch_caught[i][k] and alns2_batch_caught[i][k] + // are the two mates of the same alignment pair (index 0 is the primary + // pair). If both mates of the primary pair failed to surject, promote + // the best surjectable secondary pair into its place. + surjector.promote_secondary_pair_if_primary_unmapped(alns1_batch_caught[i], alns2_batch_caught[i]); } // Forward it along backing->emit_mapped_pairs(std::move(alns1_batch_caught), std::move(alns2_batch_caught), std::move(tlen_limit_batch)); diff --git a/src/surjecting_alignment_emitter.hpp b/src/surjecting_alignment_emitter.hpp index f876aaa62ed..c65287ea2fe 100644 --- a/src/surjecting_alignment_emitter.hpp +++ b/src/surjecting_alignment_emitter.hpp @@ -37,7 +37,7 @@ class SurjectingAlignmentEmitter : public vg::io::AlignmentEmitter { unordered_set paths, unique_ptr&& backing, bool prune_suspicious_anchors = false, bool add_graph_alignment_tag = false, bool report_supplementary = false, bool add_off_ref_position_tag = false, - bool left_align = false); + bool left_align = false, bool promote_secondary = false); /// Force full length alignment in surjection resolution bool surject_subpath_global = true; diff --git a/src/surjector.cpp b/src/surjector.cpp index dc19d9f76fc..0380276915b 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -588,6 +588,10 @@ using namespace std; if (source_aln->is_secondary() || (i != 0 && !is_supplementary(alns_out->back()))) { alns_out->back().set_is_secondary(true); + // Non-promoted secondaries get mapq 0 in surjected output; + // promote_secondary_if_primary_unmapped restores mapq on the + // winner if promotion occurs. + alns_out->back().set_mapping_quality(0); } if (annotate_with_all_path_scores) { @@ -5311,6 +5315,199 @@ using namespace std; } } + bool Surjector::promote_secondary_if_primary_unmapped(vector& surjected_group) const { + + if (!promote_secondary_on_failed_surjection) { + return false; + } + + if (surjected_group.empty()) { + return false; + } + + // Find the primary: the unique non-secondary, non-supplementary member. + int64_t primary_idx = -1; + for (size_t i = 0; i < surjected_group.size(); ++i) { + if (!surjected_group[i].is_secondary() && !is_supplementary(surjected_group[i])) { + primary_idx = (int64_t) i; + break; + } + } + + if (primary_idx < 0) { + // No primary to replace (shouldn't normally happen). + return false; + } + + if (surjected_group[primary_idx].path().mapping_size() != 0) { + // The primary surjected successfully; nothing to do. + return false; + } + + // The primary is unmapped: look for the best-scoring mapped secondary to + // promote. Supplementary alignments are never eligible to become primary. + int64_t best_idx = -1; + int32_t best_score = numeric_limits::min(); + for (size_t i = 0; i < surjected_group.size(); ++i) { + if ((int64_t) i == primary_idx) { + continue; + } + const Alignment& candidate = surjected_group[i]; + if (candidate.path().mapping_size() == 0) { + // Still unmapped after surjection; not a usable candidate. + continue; + } + if (is_supplementary(candidate)) { + // Supplementaries can't stand in for a primary. + continue; + } + // Break ties toward higher mapping quality for determinism. + if (candidate.score() > best_score || + (candidate.score() == best_score && best_idx >= 0 && + candidate.mapping_quality() > surjected_group[best_idx].mapping_quality())) { + best_score = candidate.score(); + best_idx = (int64_t) i; + } + } + + if (best_idx < 0) { + // No mapped secondary could be promoted. Warn once so the user knows + // why an expected promotion didn't happen. + if (!warned_about_no_promotable_secondary.test_and_set()) { + #pragma omp critical (cerr) + { + cerr << "warning:[Surjector] --promote-secondary: cannot find a surjectable " + << "secondary alignment for read \"" + << surjected_group[primary_idx].name() + << "\" with an unsurjectable primary alignment."; + if (warn_about_input_collation) { + cerr << " Ensure the input is collated by read name so that secondary " + << "alignments are grouped with their primary alignment."; + } + cerr << " Suppressing further warnings." << endl; + } + } + return false; + } + + // Promote the chosen secondary in place of the unmapped primary. + Alignment promoted = std::move(surjected_group[best_idx]); + promoted.set_is_secondary(false); + // Record provenance so downstream consumers can tell this happened. + set_annotation(promoted, "promoted_from_secondary", true); + + // Remove the old unmapped primary and the now-redundant secondary copy. + // Erase the higher index first so the lower index stays valid. + size_t hi = (size_t) max(primary_idx, best_idx); + size_t lo = (size_t) min(primary_idx, best_idx); + surjected_group.erase(surjected_group.begin() + hi); + surjected_group.erase(surjected_group.begin() + lo); + + // Place the promoted alignment at the front as the new primary. + surjected_group.insert(surjected_group.begin(), std::move(promoted)); + + return true; + } + + bool Surjector::promote_secondary_pair_if_primary_unmapped(vector& surjected1, + vector& surjected2) const { + + if (!promote_secondary_on_failed_surjection) { + return false; + } + + // The two mate vectors must be index-aligned (mate 1 and mate 2 of the + // same alignment pair share an index). If they somehow aren't, do + // nothing rather than risk mispairing mates. + if (surjected1.size() != surjected2.size() || surjected1.empty()) { + return false; + } + + auto is_mapped = [](const Alignment& aln) { + return aln.path().mapping_size() != 0; + }; + + // Only act when BOTH mates of the primary pair (index 0) failed to + // surject. A half-mapped primary pair is a valid SAM state and is left + // alone. + if (is_mapped(surjected1[0]) || is_mapped(surjected2[0])) { + return false; + } + + // Find the best secondary pair to promote. A candidate is eligible if at + // least one of its mates surjected. Rank by summed mapped-mate score; + // break ties toward fully-mapped pairs, then toward higher summed + // mapping quality, for determinism. + int64_t best_idx = -1; + int32_t best_score = numeric_limits::min(); + bool best_fully_mapped = false; + int32_t best_mapq = numeric_limits::min(); + for (size_t k = 1; k < surjected1.size(); ++k) { + bool m1 = is_mapped(surjected1[k]); + bool m2 = is_mapped(surjected2[k]); + if (!m1 && !m2) { + // Neither mate surjected; not a usable candidate. + continue; + } + int32_t score = (m1 ? surjected1[k].score() : 0) + (m2 ? surjected2[k].score() : 0); + bool fully_mapped = m1 && m2; + int32_t mapq = (m1 ? surjected1[k].mapping_quality() : 0) + + (m2 ? surjected2[k].mapping_quality() : 0); + bool better = false; + if (score != best_score) { + better = score > best_score; + } + else if (fully_mapped != best_fully_mapped) { + // Same score: prefer the fully-mapped pair. + better = fully_mapped; + } + else { + better = mapq > best_mapq; + } + if (best_idx < 0 || better) { + best_idx = (int64_t) k; + best_score = score; + best_fully_mapped = fully_mapped; + best_mapq = mapq; + } + } + + if (best_idx < 0) { + // No secondary pair could be promoted. Warn once (shared with the + // single-end path's flag). + if (!warned_about_no_promotable_secondary.test_and_set()) { + #pragma omp critical (cerr) + { + cerr << "warning:[Surjector] --promote-secondary: cannot find a surjectable " + << "secondary alignment for read pair \"" + << surjected1[0].name() + << "\" with an unsurjectable primary alignment."; + if (warn_about_input_collation) { + cerr << " Ensure the input is collated by read name so that secondary " + << "alignments are grouped with their primary alignment."; + } + cerr << " Suppressing further warnings." << endl; + } + } + return false; + } + + // Promote: swap the chosen secondary pair into the primary slot for both + // mates, keeping the two vectors index-aligned. Update is_secondary + // flags: the promoted pair becomes primary, the demoted pair becomes + // secondary. Record provenance on the promoted mates. + surjected1[best_idx].set_is_secondary(false); + surjected2[best_idx].set_is_secondary(false); + surjected1[0].set_is_secondary(true); + surjected2[0].set_is_secondary(true); + set_annotation(surjected1[best_idx], "promoted_from_secondary", true); + set_annotation(surjected2[best_idx], "promoted_from_secondary", true); + std::swap(surjected1[0], surjected1[best_idx]); + std::swap(surjected2[0], surjected2[best_idx]); + + return true; + } + Alignment Surjector::make_null_alignment(const Alignment& source) { Alignment null; null.set_name(source.name()); diff --git a/src/surjector.hpp b/src/surjector.hpp index 70fe220089b..ed8768d0950 100644 --- a/src/surjector.hpp +++ b/src/surjector.hpp @@ -137,6 +137,11 @@ using namespace std; /// And have we complained about hitting it? mutable atomic_flag warned_about_subgraph_size = ATOMIC_FLAG_INIT; + + /// Have we already warned that a read's primary failed to surject and no + /// mapped secondary was available to promote? Used to emit that warning + /// only once across all threads. + mutable atomic_flag warned_about_no_promotable_secondary = ATOMIC_FLAG_INIT; bool prune_suspicious_anchors = false; int64_t max_tail_anchor_prune = 4; @@ -185,8 +190,36 @@ using namespace std; bool annotate_off_reference_pos = false; /// How far we will traverse the graph in search of a reference position? size_t off_reference_pos_search_limit = 20000; - - + + /// If a read's primary alignment fails to surject (becomes unmapped), + /// promote the best-scoring secondary alignment that does successfully + /// surject to be the new primary in its place. This avoids emitting + /// secondary alignments alongside an unmapped primary, which violates + /// the SAM/BAM spec. Only secondary alignments already present in the + /// input are considered; no realignment is ever performed. + bool promote_secondary_on_failed_surjection = false; + + /// If true, the one-time "no promotable secondary" warning includes a + /// note that the input may not be collated by read name. Set this only + /// for standalone vg surject; alignment subcommands (giraffe, map, + /// mpmap) always deliver a read's full multimapping together. + bool warn_about_input_collation = false; + + /// If the primary in surjected_group is unmapped, promote the highest-scoring + /// mapped, non-supplementary secondary to primary. No realignment is performed. + /// Requires promote_secondary_on_failed_surjection. Emits a one-time warning + /// if no promotable secondary exists. Returns true iff a promotion occurred. + bool promote_secondary_if_primary_unmapped(vector& surjected_group) const; + + /// Paired-end counterpart of promote_secondary_if_primary_unmapped. + /// surjected1[k] and surjected2[k] are the two mates of the same alignment + /// pair (index 0 = primary, k > 0 = secondaries). Acts only when both mates + /// of the primary pair are unmapped; promotes the secondary pair with the + /// highest summed score. The two vectors remain the same length and + /// index-aligned. Returns true iff a promotion occurred. + bool promote_secondary_pair_if_primary_unmapped(vector& surjected1, + vector& surjected2) const; + protected: /// Do the extra score setup for the DP-only Aligner. diff --git a/src/unittest/surject.cpp b/src/unittest/surject.cpp index 4de83ba1574..26497e61628 100644 --- a/src/unittest/surject.cpp +++ b/src/unittest/surject.cpp @@ -7,6 +7,7 @@ #include "catch.hpp" #include "surjector.hpp" #include "aligner.hpp" +#include "annotation.hpp" #include "bdsg/hash_graph.hpp" #include "bdsg/overlays/path_position_overlays.hpp" @@ -972,5 +973,282 @@ TEST_CASE("Supplementary alignments can be generated", "[surject]") { } } +TEST_CASE("Surjector promotes a mapped secondary when the primary fails to surject", "[surject][promote]") { + + // A minimal graph is enough; promote_secondary_if_primary_unmapped only + // inspects the already-surjected alignments, it does not touch the graph. + bdsg::HashGraph graph; + handle_t h1 = graph.create_handle("ACGT"); + path_handle_t p = graph.create_path_handle("ref"); + graph.append_step(p, h1); + bdsg::PositionOverlay pos_graph(&graph); + + Surjector surjector(&pos_graph); + + // Helper to build a surjected alignment. If mapped, give it one mapping so + // path().mapping_size() > 0; if unmapped, leave the path empty. + auto make_aln = [&](const string& name, bool secondary, bool mapped, int32_t score, + int32_t mapq, bool supplementary) { + Alignment aln; + aln.set_name(name); + aln.set_sequence("ACGT"); + aln.set_is_secondary(secondary); + aln.set_score(score); + aln.set_mapping_quality(mapq); + if (mapped) { + auto* m = aln.mutable_path()->add_mapping(); + m->mutable_position()->set_node_id(graph.get_id(h1)); + } + if (supplementary) { + set_annotation(aln, "supplementary", true); + } + return aln; + }; + + SECTION("Disabled by default: no promotion even if primary is unmapped") { + vector group; + group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary + group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary + REQUIRE(surjector.promote_secondary_on_failed_surjection == false); + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == false); + // Group is unchanged: primary still unmapped, secondary still present. + REQUIRE(group.size() == 2); + REQUIRE(group[0].path().mapping_size() == 0); + REQUIRE(group[0].is_secondary() == false); + } + + surjector.promote_secondary_on_failed_surjection = true; + + SECTION("Promotes the single mapped secondary into the primary slot") { + vector group; + group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary + group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == true); + // Exactly one alignment remains: the promoted secondary as new primary. + REQUIRE(group.size() == 1); + REQUIRE(group[0].is_secondary() == false); + REQUIRE(group[0].path().mapping_size() > 0); + REQUIRE(group[0].score() == 30); + REQUIRE(get_annotation(group[0], "promoted_from_secondary") == true); + } + + SECTION("Promotes the highest-scoring mapped secondary") { + vector group; + group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary + group.push_back(make_aln("r", true, true, 20, 30, false)); // lower-scoring secondary + group.push_back(make_aln("r", true, true, 45, 35, false)); // higher-scoring secondary + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == true); + REQUIRE(group.front().is_secondary() == false); + REQUIRE(group.front().score() == 45); + // The other mapped secondary is retained (still marked secondary). + REQUIRE(group.size() == 2); + REQUIRE(group[1].is_secondary() == true); + REQUIRE(group[1].score() == 20); + } + + SECTION("No-op when the primary surjected successfully") { + vector group; + group.push_back(make_aln("r", false, true, 50, 60, false)); // mapped primary + group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == false); + REQUIRE(group.size() == 2); + REQUIRE(group[0].is_secondary() == false); + REQUIRE(group[0].score() == 50); + } + + SECTION("No promotion when no secondary could be surjected") { + vector group; + group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary + group.push_back(make_aln("r", true, false, 0, 0, false)); // unmapped secondary + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == false); + REQUIRE(group.size() == 2); + REQUIRE(group[0].is_secondary() == false); + REQUIRE(group[0].path().mapping_size() == 0); + } + + SECTION("Supplementary alignments are never promoted") { + vector group; + group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary + group.push_back(make_aln("r", true, true, 90, 60, true)); // mapped but supplementary + bool promoted = surjector.promote_secondary_if_primary_unmapped(group); + REQUIRE(promoted == false); + REQUIRE(group.size() == 2); + REQUIRE(group[0].is_secondary() == false); + REQUIRE(group[0].path().mapping_size() == 0); + } +} + +TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to surject", + "[surject][promote]") { + + // As with the single-end promotion test, the graph is only needed to + // construct a Surjector; promote_secondary_pair_if_primary_unmapped inspects + // the already-surjected alignments and never touches the graph. + bdsg::HashGraph graph; + handle_t h1 = graph.create_handle("ACGT"); + path_handle_t p = graph.create_path_handle("ref"); + graph.append_step(p, h1); + bdsg::PositionOverlay pos_graph(&graph); + + Surjector surjector(&pos_graph); + + // Build one mate of a surjected pair. Mapped mates get a single mapping so + // path().mapping_size() > 0; unmapped mates leave the path empty. + auto make_aln = [&](const string& name, bool secondary, bool mapped, int32_t score, + int32_t mapq) { + Alignment aln; + aln.set_name(name); + aln.set_sequence("ACGT"); + aln.set_is_secondary(secondary); + aln.set_score(score); + aln.set_mapping_quality(mapq); + if (mapped) { + auto* m = aln.mutable_path()->add_mapping(); + m->mutable_position()->set_node_id(graph.get_id(h1)); + } + return aln; + }; + + SECTION("Disabled by default: no promotion even if the primary pair is unmapped") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary mate 1 + mate2.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary mate 2 + mate1.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary mate 1 + mate2.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary mate 2 + REQUIRE(surjector.promote_secondary_on_failed_surjection == false); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == false); + REQUIRE(mate1.size() == 2); + REQUIRE(mate2.size() == 2); + REQUIRE(mate1[0].path().mapping_size() == 0); + REQUIRE(mate2[0].path().mapping_size() == 0); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + } + + surjector.promote_secondary_on_failed_surjection = true; + + SECTION("Promotes a fully-mapped secondary pair into the primary slot") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + mate1.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary pair + mate2.push_back(make_aln("r", true, true, 25, 35)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == true); + // Vectors stay the same length and index-aligned. + REQUIRE(mate1.size() == 2); + REQUIRE(mate2.size() == 2); + // The promoted pair is now the primary (index 0) for both mates. + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + REQUIRE(mate1[0].path().mapping_size() > 0); + REQUIRE(mate2[0].path().mapping_size() > 0); + REQUIRE(mate1[0].score() == 30); + REQUIRE(mate2[0].score() == 25); + REQUIRE(get_annotation(mate1[0], "promoted_from_secondary") == true); + REQUIRE(get_annotation(mate2[0], "promoted_from_secondary") == true); + // The demoted former-primary pair is now secondary. + REQUIRE(mate1[1].is_secondary() == true); + REQUIRE(mate2[1].is_secondary() == true); + } + + SECTION("Promotes the secondary pair with the highest summed mapped-mate score") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + mate1.push_back(make_aln("r", true, true, 20, 30)); // secondary pair, summed 40 + mate2.push_back(make_aln("r", true, true, 20, 30)); + mate1.push_back(make_aln("r", true, true, 45, 35)); // secondary pair, summed 80 + mate2.push_back(make_aln("r", true, true, 35, 25)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == true); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + // The higher summed-score pair (45 + 35 = 80) wins. + REQUIRE(mate1[0].score() == 45); + REQUIRE(mate2[0].score() == 35); + } + + SECTION("Prefers a fully-mapped pair over a half-mapped pair with equal summed score") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + // Half-mapped secondary: only mate 1 mapped, summed score 50. + mate1.push_back(make_aln("r", true, true, 50, 40)); + mate2.push_back(make_aln("r", true, false, 0, 0)); + // Fully-mapped secondary: both mapped, summed score 50 (30 + 20). + mate1.push_back(make_aln("r", true, true, 30, 20)); + mate2.push_back(make_aln("r", true, true, 20, 20)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == true); + // The fully-mapped pair is chosen despite the tie in summed score. + REQUIRE(mate1[0].path().mapping_size() > 0); + REQUIRE(mate2[0].path().mapping_size() > 0); + REQUIRE(mate1[0].score() == 30); + REQUIRE(mate2[0].score() == 20); + } + + SECTION("Promotes a half-mapped secondary pair when no fully-mapped pair exists") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + // Only mate 2 of the secondary surjected. + mate1.push_back(make_aln("r", true, false, 0, 0)); + mate2.push_back(make_aln("r", true, true, 40, 30)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == true); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + // The mapped mate carries the alignment; its partner stays unmapped. + REQUIRE(mate1[0].path().mapping_size() == 0); + REQUIRE(mate2[0].path().mapping_size() > 0); + REQUIRE(mate2[0].score() == 40); + } + + SECTION("No-op when either mate of the primary pair surjected") { + vector mate1, mate2; + // Primary pair is half-mapped: a valid SAM state, left untouched. + mate1.push_back(make_aln("r", false, true, 50, 60)); + mate2.push_back(make_aln("r", false, false, 0, 0)); + mate1.push_back(make_aln("r", true, true, 30, 40)); + mate2.push_back(make_aln("r", true, true, 30, 40)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == false); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate1[0].score() == 50); + } + + SECTION("No promotion when no secondary pair has a mapped mate") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + mate1.push_back(make_aln("r", true, false, 0, 0)); // unmapped secondary pair + mate2.push_back(make_aln("r", true, false, 0, 0)); + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == false); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + REQUIRE(mate1[0].path().mapping_size() == 0); + REQUIRE(mate2[0].path().mapping_size() == 0); + } + + SECTION("No-op when the mate vectors are not index-aligned") { + vector mate1, mate2; + mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair + mate2.push_back(make_aln("r", false, false, 0, 0)); + mate1.push_back(make_aln("r", true, true, 30, 40)); // secondary only present for mate 1 + bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(promoted == false); + REQUIRE(mate1.size() == 2); + REQUIRE(mate2.size() == 1); + } +} + } } diff --git a/test/t/07_vg_map.t b/test/t/07_vg_map.t index ddc4147284e..1a3db3f1934 100644 --- a/test/t/07_vg_map.t +++ b/test/t/07_vg_map.t @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap PATH=../bin:$PATH # for vg -plan tests 60 +plan tests 61 vg construct -m 1000 -r small/x.fa -v small/x.vcf.gz >x.vg vg index -x x.xg -g x.gcsa -k 11 x.vg @@ -215,3 +215,5 @@ is "$(samtools view t1.sam | grep T1 | grep T2 | grep T3 | wc -l | sed 's/^[[:sp rm tagged1.fq tagged2.fq t1.bam t2.bam t3.bam t1.sam rm -f x.vg x.xg x.gcsa x.gcsa.lcp + +is $(vg map --help 2>&1 | grep -c "promote-secondary") 1 "vg map advertises --promote-secondary" diff --git a/test/t/15_vg_surject.t b/test/t/15_vg_surject.t index 6d25d8cab39..b996b39a609 100644 --- a/test/t/15_vg_surject.t +++ b/test/t/15_vg_surject.t @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap PATH=../bin:$PATH # for vg -plan tests 78 +plan tests 85 vg construct -r small/x.fa >j.vg vg index -x j.xg j.vg @@ -278,3 +278,33 @@ vg map -d g -f reads/ts.fq | vg surject -x g.xg -b --off-ref-position - > g.bam is $(samtools view g.bam | grep "NR:Z:x:8+" | wc -l | sed 's/^[[:space:]]*//') "1" "off reference reads can be annotated with the nearest reference position" rm g.xg g.gcsa g.gcsa.lcp g.bam + +# --promote-secondary is advertised and accepted +is $(vg surject --help 2>&1 | grep -c "promote-secondary") 1 "vg surject advertises --promote-secondary" + +# --promote-secondary is a safe no-op on reads whose primaries surject fine: +# every input read still yields a primary SAM record, and no reads are dropped. +vg surject -p x -x x.xg -t 1 --promote-secondary -s j.gam > promote.sam +is $(grep -v "^@" promote.sam | wc -l | sed 's/^[[:space:]]*//') 100 "surject with --promote-secondary keeps all reads when primaries surject" +is $(grep -v "^@" promote.sam | awk '{if(and($2,256)==0) print}' | wc -l | sed 's/^[[:space:]]*//') 100 "surject with --promote-secondary emits one primary per read when primaries surject" + +rm promote.sam + +# --promote-secondary on interleaved paired input is a safe no-op when the +# primary pair surjects fine: both mates are still emitted as mapped primaries +# and the pairing relationship is preserved. +echo '{"name": "read/2", "sequence": "CAAATAA", "path": {"mapping": [{"position": {"node_id": 1}, "edit": [{"from_length": 7, "to_length": 7}]}]}, "fragment_prev": {"name": "read/1"}}{"name": "read/1", "sequence": "CTTATTT", "path": {"mapping": [{"position": {"node_id": 1, "is_reverse": true}, "edit": [{"from_length": 7, "to_length": 7}]}]}, "fragment_next": {"name": "read/2"}}' | vg view -JGa - > pairpromote.gam +vg surject -p x -x x.xg -i -t 1 --promote-secondary -s pairpromote.gam > pairpromote.sam +is $(grep -v "^@" pairpromote.sam | wc -l | sed 's/^[[:space:]]*//') 2 "surject -i with --promote-secondary keeps both mates when the primary pair surjects" +is $(grep -v "^@" pairpromote.sam | awk '{if(and($2,256)==0 && and($2,4)==0) print}' | wc -l | sed 's/^[[:space:]]*//') 2 "surject -i with --promote-secondary emits both mates as mapped primaries when the primary pair surjects" +is $(grep -v "^@" pairpromote.sam | cut -f1 | sort -u | wc -l | sed 's/^[[:space:]]*//') 1 "surject -i with --promote-secondary preserves the shared QNAME across the pair" + +rm pairpromote.gam pairpromote.sam + +# --promote-secondary with interleaved GAF input is rejected with a clear message +# (paired promotion needs GAM so pairs can be grouped by read name). +echo '{"name": "read/2", "sequence": "CAAATAA", "path": {"mapping": [{"position": {"node_id": 1}, "edit": [{"from_length": 7, "to_length": 7}]}]}, "fragment_prev": {"name": "read/1"}}{"name": "read/1", "sequence": "CTTATTT", "path": {"mapping": [{"position": {"node_id": 1, "is_reverse": true}, "edit": [{"from_length": 7, "to_length": 7}]}]}, "fragment_next": {"name": "read/2"}}' | vg view -JGa - > pairgaf.gam +vg convert x.xg -G pairgaf.gam -t 1 > pairgaf.gaf +is $(vg surject -p x -x x.xg -i -G --promote-secondary -s pairgaf.gaf 2>&1 >/dev/null | grep -c "interleaved GAF input is not supported") 1 "surject -i -G with --promote-secondary reports that interleaved GAF is unsupported" + +rm -f pairgaf.gam pairgaf.gaf diff --git a/test/t/33_vg_mpmap.t b/test/t/33_vg_mpmap.t index 2210251c130..edf4b08896c 100644 --- a/test/t/33_vg_mpmap.t +++ b/test/t/33_vg_mpmap.t @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap PATH=../bin:$PATH # for vg -plan tests 25 +plan tests 26 # Exercise the GBWT @@ -174,4 +174,6 @@ is "$(samtools view t3.bam | grep T4 | grep T5 | grep T6 | grep read2 | wc -l | rm tagged1.fq tagged2.fq t1.bam t2.bam t3.bam rm x.vg x.gam xy.vg xy.xg xy.gcsa xy.snarls xy.dist xy.sam +is $(vg mpmap --help 2>&1 | grep -c "promote-secondary") 1 "vg mpmap advertises --promote-secondary" + diff --git a/test/t/50_vg_giraffe.t b/test/t/50_vg_giraffe.t index 0035ec4052e..c4f0e1c1b3c 100644 --- a/test/t/50_vg_giraffe.t +++ b/test/t/50_vg_giraffe.t @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap PATH=../bin:$PATH # for vg -plan tests 89 +plan tests 90 vg construct -a -r small/x.fa -v small/x.vcf.gz >x.vg vg index -x x.xg x.vg @@ -335,3 +335,5 @@ is "$(vg view --extract-tag PARAMS_JSON longread.gam | jq '.["track-provenance"] rm -f longread.gam 1mb1kgp.vg 1mb1kgp.dist 1mb1kgp.giraffe.gbz 1mb1kgp.shortread.withzip.min 1mb1kgp.shortread.zipcodes log.txt +is $(vg giraffe --help 2>&1 | grep -c "promote-secondary") 1 "vg giraffe advertises --promote-secondary" + From fa6545a494bf514bc9aaf2a40a32dd80b1dc9a04 Mon Sep 17 00:00:00 2001 From: gaoj66-roche Date: Mon, 20 Jul 2026 17:15:16 +0000 Subject: [PATCH 2/5] Cleanup --- src/mapping_quality_calculator.cpp | 16 ---------------- src/minimizer_mapper.cpp | 10 +++++----- src/minimizer_mapper_from_chains.cpp | 6 +++--- src/subcommand/mpmap_main.cpp | 15 --------------- 4 files changed, 8 insertions(+), 39 deletions(-) diff --git a/src/mapping_quality_calculator.cpp b/src/mapping_quality_calculator.cpp index b21f058c918..29f12819d56 100644 --- a/src/mapping_quality_calculator.cpp +++ b/src/mapping_quality_calculator.cpp @@ -271,22 +271,6 @@ void MappingQualityCalculator::compute_mapping_quality(vector& alignm for (size_t i = 1; i < alignments.size(); ++i) { alignments[0].add_secondary_score(alignments[i].score()); } - - // Compute meaningful MAPQs for all non-primary alignments using the same - // score vector. The primary-specific adjustments (identity scaling, cluster - // blending, mq_estimate cap) are not applied to secondaries. - if (alignments.size() > 1) { - vector raw_scores(alignments.size()); - for (size_t i = 0; i < alignments.size(); ++i) { - raw_scores[i] = alignments[i].score(); - } - vector all_mapqs = compute_all_mapping_qualities(raw_scores); - for (size_t i = 0; i < alignments.size(); ++i) { - if (i == max_idx) continue; - int32_t mq = (i < all_mapqs.size()) ? all_mapqs[i] : 0; - alignments[i].set_mapping_quality(min(mq, max_mapping_quality)); - } - } } void MappingQualityCalculator::compute_paired_mapping_quality(pair, vector>& alignment_pairs, diff --git a/src/minimizer_mapper.cpp b/src/minimizer_mapper.cpp index 7024c43d07c..5f0248d348c 100644 --- a/src/minimizer_mapper.cpp +++ b/src/minimizer_mapper.cpp @@ -1136,9 +1136,9 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { crash_unless(!mappings.empty()); // Compute MAPQ if not unmapped. Otherwise use 0 instead of the 50% this would give us. - // Use exact mapping quality - double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : - get_regular_aligner()->mapq_calc->compute_max_mapping_quality(scores, false); + // Use exact mapping quality + double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : + get_regular_aligner()->mapq_calc->compute_max_mapping_quality(scores, false) ; #ifdef print_minimizer_table double uncapped_mapq = mapq; @@ -1180,7 +1180,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { // Make sure to clamp 0-60. mappings.front().set_mapping_quality(max(min(mapq, 60.0), 0.0)); - + if (!supplementaries.empty()) { // Estimate a mapping quality for the supplementaries // TODO: only count the score of the overlapping portion of other alignments @@ -1204,7 +1204,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { for (size_t i = 0; i < mappings.size(); i++) { // For each output alignment in score order auto& out = mappings[i]; - + // Assign primary and secondary status out.set_is_secondary(i > 0); } diff --git a/src/minimizer_mapper_from_chains.cpp b/src/minimizer_mapper_from_chains.cpp index a4190713604..27325ab379b 100644 --- a/src/minimizer_mapper_from_chains.cpp +++ b/src/minimizer_mapper_from_chains.cpp @@ -987,8 +987,8 @@ vector MinimizerMapper::map_from_chains(Alignment& aln) { // Because the winning alignment won't necessarily *always* have the // maximum score, we need to use compute_first_mapping_quality and not // compute_max_mapping_quality. - double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : - get_regular_aligner()->mapq_calc->compute_first_mapping_quality(scaled_scores, false, &multiplicity_by_alignment); + double mapq = (mappings.front().path().mapping_size() == 0) ? 0 : + get_regular_aligner()->mapq_calc->compute_first_mapping_quality(scaled_scores, false, &multiplicity_by_alignment) ; #ifdef debug_write_minimizers #pragma omp critical @@ -1106,7 +1106,7 @@ vector MinimizerMapper::map_from_chains(Alignment& aln) { for (size_t i = 0; i < mappings.size(); i++) { // For each output alignment in score order auto& out = mappings[i]; - + // Assign primary and secondary status out.set_is_secondary(i > 0); } diff --git a/src/subcommand/mpmap_main.cpp b/src/subcommand/mpmap_main.cpp index 41c4da1f7ef..b83f9e8c704 100644 --- a/src/subcommand/mpmap_main.cpp +++ b/src/subcommand/mpmap_main.cpp @@ -98,8 +98,6 @@ static void error_if_negative(const Logger& logger, double value, const string& } } -/// Returns true if the given surjected multipath alignment is mapped (has at -/// least one aligned base). static bool mp_aln_is_mapped(const multipath_alignment_t& mp_aln) { for (size_t i = 0; i < mp_aln.subpath_size(); ++i) { if (mp_aln.subpath(i).path().mapping_size() > 0) { @@ -109,19 +107,6 @@ static bool mp_aln_is_mapped(const multipath_alignment_t& mp_aln) { return false; } -/// Apply "promote secondary on failed surjection" to a single read's surjected -/// multipath alignments, keeping the parallel path_positions array in sync. -/// -/// mp_alns holds the read's surjected primary followed by its secondaries (and -/// possibly supplementaries appended at the end); path_positions holds the -/// matching (path name, is_reverse, offset) for each. If the primary (the first -/// non-supplementary entry) is unmapped, the best-scoring mapped, -/// non-supplementary secondary is promoted into its place (both the alignment -/// and its position). Only alignments already present are considered; no -/// realignment is performed. Supplementaries are never promoted. -/// -/// warned is a shared flag used to emit the "no promotable secondary" warning -/// only once. Returns true iff a promotion occurred. static bool promote_secondary_mp(vector& mp_alns, vector>& path_positions, const string& read_name, From 96c9e2cc2d3146da58e9e033981c3662b9af6630 Mon Sep 17 00:00:00 2001 From: gaoj66-roche Date: Tue, 21 Jul 2026 17:09:46 +0000 Subject: [PATCH 3/5] Simplify surject_main code --- src/subcommand/surject_main.cpp | 420 +++++++++++---------------- src/surjecting_alignment_emitter.cpp | 2 +- 2 files changed, 163 insertions(+), 259 deletions(-) diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 4b7ddc1ac95..2bffd5d74d8 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -22,7 +22,6 @@ #include #include #include -#include #include "../utility.hpp" #include "../surjector.hpp" #include "../hts_alignment_emitter.hpp" @@ -533,7 +532,6 @@ int main_surject(int argc, char** argv) { surjector.promote_secondary_on_failed_surjection = false; } - // Count our threads int thread_count = vg::get_thread_count(); @@ -570,77 +568,17 @@ int main_surject(int argc, char** argv) { unique_ptr alignment_emitter = get_alignment_emitter("-", output_format, sequence_dictionary, thread_count, xgidx, ALIGNMENT_EMITTER_FLAG_HTS_RAW | (spliced * ALIGNMENT_EMITTER_FLAG_HTS_SPLICED)); - - // Emit one already-surjected read pair's multimappings, pairing mates by - // strand and handling supplementaries/unpaired mates. surjected1[k] and - // surjected2[k] need NOT be index-aligned here: mates are matched by - // reference strand, exactly as in the default interleaved path. - auto emit_surjected_pair = [&](vector& surjected1, vector& surjected2) { - // pair up non-supplementary alignments - unordered_map, size_t> strand_idx1, strand_idx2; - for (size_t i = 0; i < surjected1.size(); ++i) { - if (!is_supplementary(surjected1[i])) { - const auto& pos = surjected1[i].refpos(0); - strand_idx1[make_pair(pos.name(), pos.is_reverse())] = i; - } - } - for (size_t i = 0; i < surjected2.size(); ++i) { - if (!is_supplementary(surjected2[i])) { - const auto& pos = surjected2[i].refpos(0); - strand_idx2[make_pair(pos.name(), pos.is_reverse())] = i; - } - } - for (size_t i = 0; i < surjected1.size(); ++i) { - const auto& pos = surjected1[i].refpos(0); - auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse())); - if (!is_supplementary(surjected1[i]) && it != strand_idx2.end()) { - alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len); - } - else { - if (is_supplementary(surjected1[i]) && !has_annotation(surjected1[i], "mate_info")) { - string annotation; - if (!strand_idx2.empty()) { - const auto& mate = it != strand_idx2.end() ? surjected2[it->second] : surjected2[strand_idx2.begin()->second]; - annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), false)); - } - else { - annotation = std::move(mate_info("", -1, false, false)); - } - set_annotation(surjected1[i], "mate_info", annotation); - } - alignment_emitter->emit_single(std::move(surjected1[i])); - } - } - for (size_t i = 0; i < surjected2.size(); ++i) { - const auto& pos = surjected2[i].refpos(0); - auto it = strand_idx1.find(make_pair(pos.name(), !pos.is_reverse())); - if (is_supplementary(surjected2[i]) || it == strand_idx1.end()) { - if (is_supplementary(surjected2[i]) && !has_annotation(surjected2[i], "mate_info")) { - string annotation; - if (!strand_idx1.empty()) { - const auto& mate = it != strand_idx1.end() ? surjected1[it->second] : surjected1[strand_idx1.begin()->second]; - annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), true)); - } - else { - annotation = std::move(mate_info("", -1, false, true)); - } - set_annotation(surjected2[i], "mate_info", annotation); - } - alignment_emitter->emit_single(std::move(surjected2[i])); - } - } + + // Callback to determine if alignments should be grouped: if promote_secondary is enabled, we want to group by read name + // so we can access both primary and secondary alignments for a read. If promote_secondary is off, then we can just + // process each alignment independently and each 'group' will just be a single alignment. + // NOTE: In order for this to work, the input file must be collated by read name. + function same_group = [promote_secondary](const Alignment& a, const Alignment& b) { + return promote_secondary && a.name() == b.name(); }; - - if (interleaved && promote_secondary) { - // Paired secondary promotion needs a read pair's whole multimapping - // (its primary pair followed by its secondary pairs) together. We - // read grouped pairs (consecutive pair-records sharing a read name) - // so a primary pair and its secondary pairs are never split across - // worker threads. This requires the input to be collated by read - // name. If both mates of the primary pair fail to surject, the best - // surjectable secondary pair is promoted into its place. + if (interleaved) { using AlnPair = pair; - function&)> process_pair_group = [&](vector& group) { + function&)> process_group = [&](vector& group) { if (group.empty()) { return; } @@ -650,12 +588,39 @@ int main_surject(int argc, char** argv) { if (watchdog) { watchdog->check_in(thread_num, group.front().first.name()); } - // Surject every pair-record, building index-aligned mate - // vectors: promoted1[k] and promoted2[k] are the two mates of - // the same alignment pair, index 0 being the primary pair. - // Supplementaries are collected separately and emitted after. - vector promoted1, promoted2; - vector extra1, extra2; + // Make sure that the alignments are actually paired with each other + // (proper fragment_prev/fragment_next). We want to catch people giving us + // un-interleaved GAMs as interleaved. + // TODO: Integrate into for_each_interleaved_pair_parallel when running on Alignments. + const Alignment& src1 = group.front().first; + const Alignment& src2 = group.front().second; + if (src1.has_fragment_next()) { + // Alignment 1 comes first in fragment + if (src1.fragment_next().name() != src2.name() || + !src2.has_fragment_prev() || + src2.fragment_prev().name() != src1.name()) { +#pragma omp critical (cerr) + adjacent_but_not_paired_error(logger, src1.name(), src2.name()); + } + } else if (src2.has_fragment_next()) { + // Alignment 2 comes first in fragment + if (src2.fragment_next().name() != src1.name() || + !src1.has_fragment_prev() || + src1.fragment_prev().name() != src2.name()) { +#pragma omp critical (cerr) + adjacent_but_not_paired_error(logger, src1.name(), src2.name()); + } + } else { + // Alignments aren't paired up at all +#pragma omp critical (cerr) + adjacent_but_not_paired_error(logger, src1.name(), src2.name()); + } + + vector alns1, alns2; + vector extras1, extras2; + // Surject each pair in the group and run secondary promotion (if enabled) + // on the non-supplementary results. Append any supplementaries to the end + // after promotion is done for (auto& rec : group) { Alignment& src1 = rec.first; Alignment& src2 = rec.second; @@ -663,143 +628,119 @@ int main_surject(int argc, char** argv) { ensure_alignment_is_for_graph(logger, src1, *xgidx); ensure_alignment_is_for_graph(logger, src2, *xgidx); } + + // Preprocess read to set metadata before surjection set_metadata(src1); set_metadata(src2); - auto s1 = surjector.surject(src1, paths, subpath_global, spliced); - auto s2 = surjector.surject(src2, paths, subpath_global, spliced); - // Keep the non-supplementary (primary of this record) mate - // for pair-level promotion; route supplementaries to extra. - int64_t keep1 = -1, keep2 = -1; - for (size_t j = 0; j < s1.size(); ++j) { - if (keep1 < 0 && !is_supplementary(s1[j])) { - keep1 = (int64_t) j; + + // Surject + auto surjected1 = surjector.surject(src1, paths, subpath_global, spliced); + auto surjected2 = surjector.surject(src2, paths, subpath_global, spliced); + + // Split into non-supplementary and supplementary + int64_t k1 = -1, k2 = -1; + for (size_t j = 0; j < surjected1.size(); ++j) { + if (k1 < 0 && !is_supplementary(surjected1[j])) { + k1 = (int64_t) j; } else { - extra1.emplace_back(std::move(s1[j])); + extras1.emplace_back(std::move(surjected1[j])); } } - for (size_t j = 0; j < s2.size(); ++j) { - if (keep2 < 0 && !is_supplementary(s2[j])) { - keep2 = (int64_t) j; + for (size_t j = 0; j < surjected2.size(); ++j) { + if (k2 < 0 && !is_supplementary(surjected2[j])) { + k2 = (int64_t) j; } else { - extra2.emplace_back(std::move(s2[j])); + extras2.emplace_back(std::move(surjected2[j])); } } - // Fall back to a null (unmapped) alignment if somehow - // all results were supplementary, so downstream refpos - // access stays valid. - auto null_with_refpos = [&](const Alignment& src) { + auto make_placeholder = [](const Alignment& src) { Alignment a; a.set_name(src.name()); a.set_sequence(src.sequence()); a.set_quality(src.quality()); - a.add_refpos(); + a.add_refpos(); // empty refpos return a; }; - promoted1.emplace_back(keep1 >= 0 ? std::move(s1[keep1]) : null_with_refpos(src1)); - promoted2.emplace_back(keep2 >= 0 ? std::move(s2[keep2]) : null_with_refpos(src2)); + // Add placeholders if all surjected alignments were supplementary + // so the two vectors are still index-aligned + // Keep exactly one non-supplementary alignment per group so `alns1` and `alns2` + // are index-aligned + alns1.emplace_back(k1 >= 0 ? std::move(surjected1[k1]) : make_placeholder(src1)); + alns2.emplace_back(k2 >= 0 ? std::move(surjected2[k2]) : make_placeholder(src2)); } - // Promote the best surjectable secondary pair if both mates of - // the primary pair failed to surject. Warns once if none. - surjector.promote_secondary_pair_if_primary_unmapped(promoted1, promoted2); - // Emit the (possibly reordered) pairs plus any supplementaries. - for (auto& a : extra1) { - promoted1.emplace_back(std::move(a)); - } - for (auto& a : extra2) { - promoted2.emplace_back(std::move(a)); + + // No-op when promote_secondary_on_failed_surjection is false. + surjector.promote_secondary_pair_if_primary_unmapped(alns1, alns2); + + // Append back the supplementary alignments to the end of each vector + // now that secondary promotion is done. + for (auto& a : extras1) { + alns1.emplace_back(std::move(a)); } - emit_surjected_pair(promoted1, promoted2); - total_reads_surjected += 2 * group.size(); - if (watchdog) { - watchdog->check_out(thread_num); + for (auto& a : extras2) { + alns2.emplace_back(std::move(a)); } - clear_crash_context(); - } catch (const std::exception& ex) { - report_exception(ex); - } - }; - // Two pair-records belong to the same group if they share a read name. - function pairs_in_same_group = - [](const AlnPair& a, const AlnPair& b) { - return a.first.name() == b.first.name(); - }; - if (input_format == "GAM") { - get_input_file(file_name, [&](istream& in) { - vg::io::ProtobufIterator cursor(in); - function get_pair = [&](AlnPair& dest) { - if (!cursor.has_current()) { - return false; - } - dest.first = std::move(cursor.take()); - if (!cursor.has_current()) { - // Odd number of records in an interleaved GAM. - adjacent_but_not_paired_error(logger, dest.first.name(), ""); - return false; + + // pair up non-supplementary alignments + unordered_map, size_t> strand_idx1, strand_idx2; + for (size_t i = 0; i < alns1.size(); ++i) { + if (!is_supplementary(alns1[i])) { + const auto& pos = alns1[i].refpos(0); + strand_idx1[make_pair(pos.name(), pos.is_reverse())] = i; } - dest.second = std::move(cursor.take()); - return true; - }; - vg::io::grouped_unpaired_for_each_parallel(get_pair, process_pair_group, pairs_in_same_group); - }); - } - } else if (interleaved) { - // GAM input is paired, and for HTS output reads need to know their pair partners' mapping locations. - // TODO: We don't preserve order relationships (like primary/secondary) beyond the interleaving. - function lambda = [&](Alignment& src1, Alignment& src2) { - try { - set_crash_context(src1.name() + ", " + src2.name()); - size_t thread_num = omp_get_thread_num(); - if (watchdog) { - watchdog->check_in(thread_num, src1.name() + ", " + src2.name()); } - // Make sure that the alignments are actually paired with each other - // (proper fragment_prev/fragment_next). We want to catch people giving us - // un-interleaved GAMs as interleaved. - // TODO: Integrate into for_each_interleaved_pair_parallel when running on Alignments. - if (src1.has_fragment_next()) { - // Alignment 1 comes first in fragment - if (src1.fragment_next().name() != src2.name() || - !src2.has_fragment_prev() || - src2.fragment_prev().name() != src1.name()) { - -#pragma omp critical (cerr) - adjacent_but_not_paired_error(logger, src1.name(), src2.name()); - + for (size_t i = 0; i < alns2.size(); ++i) { + if (!is_supplementary(alns2[i])) { + const auto& pos = alns2[i].refpos(0); + strand_idx2[make_pair(pos.name(), pos.is_reverse())] = i; } - } else if (src2.has_fragment_next()) { - // Alignment 2 comes first in fragment - if (src2.fragment_next().name() != src1.name() || - !src1.has_fragment_prev() || - src1.fragment_prev().name() != src2.name()) { - -#pragma omp critical (cerr) - adjacent_but_not_paired_error(logger, src1.name(), src2.name()); - + } + for (size_t i = 0; i < alns1.size(); ++i) { + const auto& pos = alns1[i].refpos(0); + auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse())); + if (!is_supplementary(alns1[i]) && it != strand_idx2.end()) { + alignment_emitter->emit_pair(std::move(alns1[i]), std::move(alns2[it->second]), max_frag_len); + } else { + // supplementary or unpaired + if (is_supplementary(alns1[i]) && !has_annotation(alns1[i], "mate_info")) { + // we need to annotate this supplementary with mate info for SAM/BAM conversion + string annotation; + if (!strand_idx2.empty()) { + // there is a non-supplementary alignment available (prefer the one consistent with this path strand) + const auto& mate = it != strand_idx2.end() ? alns2[it->second] : alns2[strand_idx2.begin()->second]; + annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), false)); + } else { + // we don't have access to the primary, but we can still record the read 1/2 identity + annotation = std::move(mate_info("", -1, false, false)); + } + set_annotation(alns1[i], "mate_info", annotation); + } + alignment_emitter->emit_single(std::move(alns1[i])); } - } else { - // Alignments aren't paired up at all -#pragma omp critical (cerr) - adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } - - if (validate) { - ensure_alignment_is_for_graph(logger, src1, *xgidx); - ensure_alignment_is_for_graph(logger, src2, *xgidx); + for (size_t i = 0; i < alns2.size(); ++i) { + const auto& pos = alns2[i].refpos(0); + auto it = strand_idx1.find(make_pair(pos.name(), !pos.is_reverse())); + if (is_supplementary(alns2[i]) || it == strand_idx1.end()) { + // this strand's surjection is unpaired or supplementary + if (is_supplementary(alns2[i]) && !has_annotation(alns2[i], "mate_info")) { + // we need to annotate this supplementary with mate info for SAM/BAM conversion + string annotation; + if (!strand_idx1.empty()) { + // there is a non-supplementary alignment available (prefer the one consistent with this path strand) + const auto& mate = it != strand_idx1.end() ? alns1[it->second] : alns1[strand_idx1.begin()->second]; + annotation = std::move(mate_info(mate.refpos(0).name(), mate.refpos(0).offset(), mate.refpos(0).is_reverse(), true)); + } else { + // we don't have access to the primary, but we can still record the read 1/2 identity + annotation = std::move(mate_info("", -1, false, true)); + } + set_annotation(alns2[i], "mate_info", annotation); + } + alignment_emitter->emit_single(std::move(alns2[i])); + } } - - // Preprocess read to set metadata before surjection - set_metadata(src1); - set_metadata(src2); - - // Surject - auto surjected1 = surjector.surject(src1, paths, subpath_global, spliced); - auto surjected2 = surjector.surject(src2, paths, subpath_global, spliced); - // Pair up mates by strand and emit (shared with the - // promotion path). - emit_surjected_pair(surjected1, surjected2); - - total_reads_surjected += 2; + total_reads_surjected += 2 * group.size(); if (watchdog) { watchdog->check_out(thread_num); } @@ -810,24 +751,32 @@ int main_surject(int argc, char** argv) { }; if (input_format == "GAM") { get_input_file(file_name, [&](istream& in) { - vg::io::for_each_interleaved_pair_parallel(in, lambda); + function&)> process_interleaved = [&](vector& alns) { + vector pairs; + for (size_t i = 0; i + 1 < alns.size(); i += 2) { + pairs.emplace_back(std::move(alns[i]), std::move(alns[i + 1])); + } + if (alns.size() % 2 != 0) { +#pragma omp critical (cerr) + adjacent_but_not_paired_error(logger, alns.back().name(), ""); + } + process_group(pairs); + }; + vg::io::grouped_unpaired_for_each_parallel(in, process_interleaved, same_group); }); } else { + // promote_secondary is disabled for GAF input so this effectively just surjects each pair by itself. auto gaf_checking_lambda = [&](Alignment& src1, Alignment& src2) { check_gaf_aln(src1); check_gaf_aln(src2); - return lambda(src1, src2); + vector group; + group.emplace_back(std::move(src1), std::move(src2)); + process_group(group); }; vg::io::gaf_paired_interleaved_for_each_parallel(*xgidx, file_name, gaf_checking_lambda); } - } else if (promote_secondary) { - // Secondary promotion needs a read's whole multimapping (primary - // plus its secondaries) together so it can promote a mapped - // secondary if the primary fails to surject. We use a grouped - // parallel reader that keeps consecutive same-named reads in one - // group and never splits a group across worker threads. This - // requires the input to be collated by read name (as produced by - // vg giraffe / map / mpmap with --max-multimaps > 1). + } else { + // TODO: We don't preserve order relationships (like primary/secondary). function&)> process_group = [&](vector& group) { if (group.empty()) { return; @@ -838,7 +787,6 @@ int main_surject(int argc, char** argv) { if (watchdog) { watchdog->check_in(thread_num, group.front().name()); } - // Surject every member of the group, collecting all results. vector surjected_group; for (auto& src : group) { if (validate) { @@ -850,8 +798,7 @@ int main_surject(int argc, char** argv) { surjected_group.emplace_back(std::move(s)); } } - // Promote the best mapped secondary if the primary failed to - // surject. Emits a one-time warning if none is available. + // No-op when promote_secondary_on_failed_surjection is false. surjector.promote_secondary_if_primary_unmapped(surjected_group); alignment_emitter->emit_singles(std::move(surjected_group)); total_reads_surjected += group.size(); @@ -863,58 +810,15 @@ int main_surject(int argc, char** argv) { report_exception(ex); } }; - get_input_file(file_name, [&](istream& in) { - // Single-threaded reader; group boundaries are on read name. - vg::io::ProtobufIterator cursor(in); - function get_read = [&](Alignment& dest) { - if (!cursor.has_current()) { - return false; - } - dest = std::move(cursor.take()); - return true; - }; - function in_same_group = - [](const Alignment& a, const Alignment& b) { - return a.name() == b.name(); - }; - vg::io::grouped_unpaired_for_each_parallel(get_read, process_group, in_same_group); - }); - } else { - // We can just surject each Alignment by itself. - // TODO: We don't preserve order relationships (like primary/secondary). - function lambda = [&](Alignment& src) { - try { - set_crash_context(src.name()); - size_t thread_num = omp_get_thread_num(); - if (watchdog) { - watchdog->check_in(thread_num, src.name()); - } - if (validate) { - ensure_alignment_is_for_graph(logger, src, *xgidx); - } - - // Preprocess read to set metadata before surjection - set_metadata(src); - - // Surject and emit the single read. - alignment_emitter->emit_singles(surjector.surject(src, paths, subpath_global, spliced)); - total_reads_surjected++; - if (watchdog) { - watchdog->check_out(thread_num); - } - clear_crash_context(); - } catch (const std::exception& ex) { - report_exception(ex); - } - }; if (input_format == "GAM") { get_input_file(file_name, [&](istream& in) { - vg::io::for_each_parallel(in,lambda); + vg::io::grouped_unpaired_for_each_parallel(in, process_group, same_group); }); } else { auto gaf_checking_lambda = [&](Alignment& src) { check_gaf_aln(src); - return lambda(src); + vector group{std::move(src)}; + process_group(group); }; vg::io::gaf_unpaired_for_each_parallel(*xgidx, file_name, gaf_checking_lambda); } @@ -978,13 +882,13 @@ int main_surject(int argc, char** argv) { // TODO: highly repetitive with the version above for Alignments // surject and record path positions vector> positions1, positions2; - auto surjected1 = surjector.surject(mp_src1, paths, positions1, subpath_global, spliced); + auto alns1 = surjector.surject(mp_src1, paths, positions1, subpath_global, spliced); auto surjected2 = surjector.surject(mp_src2, paths, positions2, subpath_global, spliced); // pair up non-supplementary alignments unordered_map, size_t> strand_idx1, strand_idx2; - for (size_t i = 0; i < surjected1.size(); ++i) { - if (!is_supplementary(surjected1[i])) { + for (size_t i = 0; i < alns1.size(); ++i) { + if (!is_supplementary(alns1[i])) { strand_idx1[make_pair(get<0>(positions1[i]), get<2>(positions1[i]))] = i; } } @@ -994,11 +898,11 @@ int main_surject(int argc, char** argv) { } } - for (size_t i = 0; i < surjected1.size(); ++i) { + for (size_t i = 0; i < alns1.size(); ++i) { auto it = strand_idx2.find(make_pair(get<0>(positions1[i]), !get<2>(positions1[i]))); - if (!is_supplementary(surjected1[i]) && it != strand_idx2.end() && get<1>(positions1[i]) >= 0) { + if (!is_supplementary(alns1[i]) && it != strand_idx2.end() && get<1>(positions1[i]) >= 0) { // the alignments are paired on this strand - surjected.emplace_back(std::move(surjected1[i]), std::move(surjected2[it->second])); + surjected.emplace_back(std::move(alns1[i]), std::move(surjected2[it->second])); // reorder the positions to deal with the mismatch in the interfaces // note: we don't move() path names so we can check against them later @@ -1012,7 +916,7 @@ int main_surject(int argc, char** argv) { } else { // supplementary or unpaired - if (is_supplementary(surjected1[i]) && !surjected1[i].has_annotation("mate_info")) { + if (is_supplementary(alns1[i]) && !alns1[i].has_annotation("mate_info")) { string annotation; if (!strand_idx2.empty()) { size_t idx = it != strand_idx2.end() ? it->second : strand_idx2.begin()->second; @@ -1021,9 +925,9 @@ int main_surject(int argc, char** argv) { else { annotation = std::move(mate_info("", -1, false, false)); } - surjected1[i].set_annotation("mate_info", annotation); + alns1[i].set_annotation("mate_info", annotation); } - surjected_unpaired1.emplace_back(std::move(surjected1[i])); + surjected_unpaired1.emplace_back(std::move(alns1[i])); // reorder the position to deal with the mismatch in the interfaces positions_unpaired1.emplace_back(); diff --git a/src/surjecting_alignment_emitter.cpp b/src/surjecting_alignment_emitter.cpp index 5f564d118a9..c2f2bc794af 100644 --- a/src/surjecting_alignment_emitter.cpp +++ b/src/surjecting_alignment_emitter.cpp @@ -72,7 +72,7 @@ void SurjectingAlignmentEmitter::surject_paired_alignments_in_place(vector Date: Wed, 22 Jul 2026 03:38:09 +0000 Subject: [PATCH 4/5] Add interleaved grouped parallel iterator --- src/subcommand/surject_main.cpp | 43 ++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 2bffd5d74d8..9494a04ad9f 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -569,13 +569,6 @@ int main_surject(int argc, char** argv) { output_format, sequence_dictionary, thread_count, xgidx, ALIGNMENT_EMITTER_FLAG_HTS_RAW | (spliced * ALIGNMENT_EMITTER_FLAG_HTS_SPLICED)); - // Callback to determine if alignments should be grouped: if promote_secondary is enabled, we want to group by read name - // so we can access both primary and secondary alignments for a read. If promote_secondary is off, then we can just - // process each alignment independently and each 'group' will just be a single alignment. - // NOTE: In order for this to work, the input file must be collated by read name. - function same_group = [promote_secondary](const Alignment& a, const Alignment& b) { - return promote_secondary && a.name() == b.name(); - }; if (interleaved) { using AlnPair = pair; function&)> process_group = [&](vector& group) { @@ -751,18 +744,23 @@ int main_surject(int argc, char** argv) { }; if (input_format == "GAM") { get_input_file(file_name, [&](istream& in) { - function&)> process_interleaved = [&](vector& alns) { - vector pairs; - for (size_t i = 0; i + 1 < alns.size(); i += 2) { - pairs.emplace_back(std::move(alns[i]), std::move(alns[i + 1])); - } - if (alns.size() % 2 != 0) { -#pragma omp critical (cerr) - adjacent_but_not_paired_error(logger, alns.back().name(), ""); - } - process_group(pairs); - }; - vg::io::grouped_unpaired_for_each_parallel(in, process_interleaved, same_group); + if (promote_secondary) { + // Group consecutive pairs by first-mate name so promotion has access to + // every multimap for both mates at once. Primary and secondary alignments + // of the same mate carry the same name, so exact equality suffices. + // Input must be collated by read name. + auto same_group_pair = [](const AlnPair& a, const AlnPair& b) { + return a.first.name() == b.first.name(); + }; + vg::io::grouped_interleaved_for_each_parallel(in, process_group, same_group_pair); + } else { + // No promotion: use the standard consecutive-pair reader. + vg::io::for_each_interleaved_pair_parallel(in, [&](Alignment& src1, Alignment& src2) { + vector group; + group.emplace_back(std::move(src1), std::move(src2)); + process_group(group); + }); + } }); } else { // promote_secondary is disabled for GAF input so this effectively just surjects each pair by itself. @@ -811,6 +809,13 @@ int main_surject(int argc, char** argv) { } }; if (input_format == "GAM") { + // Callback to determine if alignments should be grouped: if promote_secondary is enabled, we want to group by read name + // so we can access both primary and secondary alignments for a read. If promote_secondary is off, then we can just + // process each alignment independently and each 'group' will just be a single alignment. + // NOTE: In order for this to work, the input file must be collated by read name. + function same_group = [promote_secondary](const Alignment& a, const Alignment& b) { + return promote_secondary && a.name() == b.name(); + }; get_input_file(file_name, [&](istream& in) { vg::io::grouped_unpaired_for_each_parallel(in, process_group, same_group); }); From 4b30af1ca2a4feaca72a0114ffd5b0e4ac41a90b Mon Sep 17 00:00:00 2001 From: gaoj66-roche Date: Thu, 30 Jul 2026 21:45:06 +0000 Subject: [PATCH 5/5] Replace secondary promotion with rescue tagging via YF:i:1 Instead of actually promoting a secondary alignment to primary when the primary fails to surject, tag the best surjectable secondary with YF:i:1 so downstream tools may treat it as primary if they choose. The secondary flag and alignment ordering are left unchanged, avoiding the introduction of low-quality primary alignments that could skew variant calling. Rename --promote-secondary to --rescue-secondary throughout (surject, map, mpmap, giraffe), rename all related identifiers and the internal annotation key (promoted_from_secondary -> rescued_secondary), and change the BAM/SAM tag from ps:i:1 to YF:i:1. Update unit test assertions to match tag-only semantics. Co-Authored-By: Claude Sonnet 4.6 (1M context) --- src/alignment.cpp | 14 +-- src/hts_alignment_emitter.cpp | 2 +- src/hts_alignment_emitter.hpp | 8 +- src/subcommand/giraffe_main.cpp | 38 +++--- src/subcommand/map_main.cpp | 42 ++++--- src/subcommand/mpmap_main.cpp | 109 ++++++++--------- src/subcommand/surject_main.cpp | 56 +++++---- src/surjecting_alignment_emitter.cpp | 14 +-- src/surjecting_alignment_emitter.hpp | 2 +- src/surjector.cpp | 67 ++++------ src/surjector.hpp | 38 +++--- src/unittest/surject.cpp | 176 +++++++++++++++------------ 12 files changed, 277 insertions(+), 289 deletions(-) diff --git a/src/alignment.cpp b/src/alignment.cpp index f5e0e6e16da..8932343378e 100644 --- a/src/alignment.cpp +++ b/src/alignment.cpp @@ -689,9 +689,9 @@ string alignment_to_sam_internal(const Alignment& alignment, if (has_annotation(alignment, "nearest_ref_pos")) { sam << "\tNR:Z:" << get_annotation(alignment, "nearest_ref_pos"); } - if (has_annotation(alignment, "promoted_from_secondary")) { - if (get_annotation(alignment, "promoted_from_secondary")) { - sam << "\tps:i:1"; + if (has_annotation(alignment, "rescued_secondary")) { + if (get_annotation(alignment, "rescued_secondary")) { + sam << "\tYF:i:1"; } } @@ -1046,10 +1046,10 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header, string pos = get_annotation(alignment, "nearest_ref_pos"); bam_aux_append(bam, "NR", 'Z', pos.size() + 1, (uint8_t*) pos.c_str()); } - if (has_annotation(alignment, "promoted_from_secondary")) { - if (get_annotation(alignment, "promoted_from_secondary")) { + if (has_annotation(alignment, "rescued_secondary")) { + if (get_annotation(alignment, "rescued_secondary")) { int32_t val = 1; - bam_aux_append(bam, "ps", 'i', sizeof(int32_t), (uint8_t*) &val); + bam_aux_append(bam, "YF", 'i', sizeof(int32_t), (uint8_t*) &val); } } @@ -1059,7 +1059,7 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header, auto parsed_tags = parse_sam_tags(get_annotation(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" || get<0>(tag) == "ps") { + 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; } diff --git a/src/hts_alignment_emitter.cpp b/src/hts_alignment_emitter.cpp index 233d86cf209..09fbaebcc32 100644 --- a/src/hts_alignment_emitter.cpp +++ b/src/hts_alignment_emitter.cpp @@ -65,7 +65,7 @@ unique_ptr get_alignment_emitter(const string& filename, const 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_PROMOTE_SECONDARY); + flags & ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY); } } else { diff --git a/src/hts_alignment_emitter.hpp b/src/hts_alignment_emitter.hpp index a8973cd634d..b27230ac3e5 100644 --- a/src/hts_alignment_emitter.hpp +++ b/src/hts_alignment_emitter.hpp @@ -54,10 +54,10 @@ enum alignment_emitter_flags_t { ALIGNMENT_EMITTER_FLAG_HTS_OFF_REF_POSITION = 64, /// When surjecting, attempt to left align ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN = 128, - /// When surjecting, if a read's primary alignment fails to surject, promote - /// the best-scoring secondary that does surject to be the new primary - /// instead of emitting an unmapped primary alongside mapped secondaries. - ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_SECONDARY = 256 + /// 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 diff --git a/src/subcommand/giraffe_main.cpp b/src/subcommand/giraffe_main.cpp index 7198f23758f..f36bf097cf1 100644 --- a/src/subcommand/giraffe_main.cpp +++ b/src/subcommand/giraffe_main.cpp @@ -732,9 +732,11 @@ void help_giraffe(char** argv, const BaseOptionGroup& parser, const std::map 1)" << 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 @@ -800,7 +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_PROMOTE_SECONDARY = 1017; + constexpr int OPT_RESCUE_SECONDARY = 1017; constexpr int OPT_HAPLOTYPE_NAME = 1100; constexpr int OPT_KFF_NAME = 1101; @@ -912,8 +914,8 @@ 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 promote its best surjectable secondary? - bool promote_secondary = 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; @@ -1173,7 +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}, - {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, + {"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'}, @@ -1366,8 +1368,8 @@ int main_giraffe(int argc, char** argv) { annotate_off_ref_position = true; break; - case OPT_PROMOTE_SECONDARY: - promote_secondary = true; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; break; case OPT_LEFT_ALIGN: @@ -2069,7 +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("promote-secondary", promote_secondary); + report_flag("rescue-secondary", rescue_secondary); report_flag("left-align", left_align); report_flag("set-refpos", set_refpos); minimizer_mapper.set_refpos = set_refpos; @@ -2215,22 +2217,22 @@ int main_giraffe(int argc, char** argv) { // When surjecting, attempt to left align flags |= ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN; } - if (promote_secondary && minimizer_mapper.max_multimaps < 2) { - logger.warn() << "--promote-secondary requires --max-multimaps > 1; " + 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 promote in case of an " + << "so there are no secondary alignments to rescue in case of an " << "unsurjectable primary alignment. Ignoring." << endl; - promote_secondary = false; + rescue_secondary = false; } - if (promote_secondary) { + if (rescue_secondary) { if (paired && !interleaved) { - logger.warn() << "--promote-secondary with two-file paired input uses paired-end " + 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, promote a mapped secondary if the primary fails to surject - flags |= ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_SECONDARY; + // 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. diff --git a/src/subcommand/map_main.cpp b/src/subcommand/map_main.cpp index c36b4123d86..0f2b8f4f746 100644 --- a/src/subcommand/map_main.cpp +++ b/src/subcommand/map_main.cpp @@ -132,9 +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 - << " --promote-secondary in HTSlib output, if a read's primary fails to" << endl - << " surject, promote its best surjectable secondary" << endl - << " to primary (needs --max-multimaps > 1)" << 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 @@ -168,7 +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_PROMOTE_SECONDARY = 1008; + constexpr int OPT_RESCUE_SECONDARY = 1008; string matrix_file_name; string seq; string qual; @@ -247,7 +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 promote_secondary = false; + bool rescue_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -327,7 +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}, - {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, + {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY}, {"help", no_argument, 0, 'h'}, {0, 0, 0, 0} }; @@ -656,8 +658,8 @@ int main_map(int argc, char** argv) { comments_as_tags = true; break; - case OPT_PROMOTE_SECONDARY: - promote_secondary = true; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; break; case 'h': @@ -810,22 +812,22 @@ int main_map(int argc, char** argv) { paths = get_sequence_dictionary(ref_paths_name, {}, reference_assembly_names, *xgidx); } - if (promote_secondary && !hts_output) { - logger.warn() << "--promote-secondary has no effect unless surjecting to SAM, BAM, or CRAM " + if (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless surjecting to SAM, BAM, or CRAM " << "(--surject-to); ignoring." << endl; - promote_secondary = false; + rescue_secondary = false; } - if (promote_secondary && max_multimaps < 2) { - logger.warn() << "--promote-secondary requires --max-multimaps > 1; " + 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 promote in case of an " + << "so there are no secondary alignments to rescue in case of an " << "unsurjectable primary alignment. Ignoring." << endl; - promote_secondary = false; + rescue_secondary = false; } - if (promote_secondary && !interleaved_input && !fastq2.empty()) { - logger.warn() << "--promote-secondary with two-file paired input uses paired-end " + 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; @@ -833,9 +835,9 @@ int main_map(int argc, char** argv) { // Set up output to an emitter that will handle serialization and surjection int emitter_flags = ALIGNMENT_EMITTER_FLAG_NONE; - if (promote_secondary) { - // When surjecting, promote a mapped secondary if the primary fails to surject. - emitter_flags |= ALIGNMENT_EMITTER_FLAG_HTS_PROMOTE_SECONDARY; + if (rescue_secondary) { + // When surjecting, rescue a secondary if the primary fails to surject. + emitter_flags |= ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY; } unique_ptr alignment_emitter = get_alignment_emitter("-", output_format, paths, thread_count, xgidx, emitter_flags); diff --git a/src/subcommand/mpmap_main.cpp b/src/subcommand/mpmap_main.cpp index b83f9e8c704..8a4e2d2753c 100644 --- a/src/subcommand/mpmap_main.cpp +++ b/src/subcommand/mpmap_main.cpp @@ -107,7 +107,7 @@ static bool mp_aln_is_mapped(const multipath_alignment_t& mp_aln) { return false; } -static bool promote_secondary_mp(vector& mp_alns, +static bool rescue_secondary_mp(vector& mp_alns, vector>& path_positions, const string& read_name, std::atomic_flag& warned) { @@ -128,7 +128,7 @@ static bool promote_secondary_mp(vector& mp_alns, return false; } - // Find the best-scoring mapped, non-supplementary secondary to promote. + // Find the best-scoring mapped, non-supplementary secondary to rescue. int64_t best_idx = -1; int32_t best_score = numeric_limits::min(); for (size_t i = 0; i < mp_alns.size(); ++i) { @@ -148,13 +148,13 @@ static bool promote_secondary_mp(vector& mp_alns, } if (best_idx < 0) { - // No mapped secondary to promote. Warn once. + // No mapped secondary to tag. Warn once. if (!warned.test_and_set()) { #pragma omp critical (cerr) { - cerr << "warning:[vg mpmap] --promote-secondary was requested, but a read whose " + cerr << "warning:[vg mpmap] --rescue-secondary was requested, but a read whose " << "primary alignment failed to surject had no surjectable secondary alignment " - << "to promote (first seen for read \"" << read_name << "\"). This can happen " + << "to tag (first seen for read \"" << read_name << "\"). This can happen " << "because the reads were not generated with secondaries (run mpmap with " << "--max-multimaps > 1). Affected reads are left with an unmapped primary. " << "This warning is shown only once." << endl; @@ -163,30 +163,25 @@ static bool promote_secondary_mp(vector& mp_alns, return false; } - // Promote: swap the chosen secondary into the primary slot, mark it primary, - // and mark the demoted (unmapped) alignment secondary. Keep positions in sync. - mp_alns[best_idx].set_annotation("secondary", false); - mp_alns[primary_idx].set_annotation("secondary", true); - mp_alns[best_idx].set_annotation("promoted_from_secondary", true); - std::swap(mp_alns[primary_idx], mp_alns[best_idx]); - std::swap(path_positions[primary_idx], path_positions[best_idx]); + // Tag the best secondary so downstream tools can treat it as primary if + // they choose; the secondary annotation and ordering are left unchanged. + mp_alns[best_idx].set_annotation("rescued_secondary", true); return true; } -/// Paired counterpart of promote_secondary_mp for mpmap's surjected output. +/// Paired counterpart of rescue_secondary_mp for mpmap's surjected output. /// /// output_mp_aln_pairs holds one surjected pair per multimapping (index 0 is the /// primary pair, later indices are secondary pairs); path_positions is the /// matching parallel array of (mate1 pos, mate2 pos). If BOTH mates of the /// primary pair are unmapped, the secondary pair with the highest summed -/// mapped-mate score is promoted into index 0 (ties broken toward fully-mapped -/// pairs), keeping path_positions in sync and updating the "secondary" -/// annotations. A candidate is eligible if at least one of its mates is mapped. +/// mapped-mate score is tagged with YF:i:1 (ties broken toward fully-mapped +/// pairs). A candidate is eligible if at least one of its mates is mapped. /// Only alignments already present are considered; no realignment is performed. /// -/// warned is shared so the "no promotable secondary" warning is emitted once. -/// Returns true iff a promotion occurred. -static bool promote_secondary_pair_mp(vector>& output_mp_aln_pairs, +/// warned is shared so the "no rescuable secondary" warning is emitted once. +/// Returns true iff a secondary pair was tagged. +static bool rescue_secondary_pair_mp(vector>& output_mp_aln_pairs, vector, tuple>>& path_positions, const string& read_name, std::atomic_flag& warned) { @@ -237,9 +232,9 @@ static bool promote_secondary_pair_mp(vector 1). Affected read pairs are left with an unmapped primary. " << "This warning is shown only once." << endl; @@ -248,14 +243,10 @@ static bool promote_secondary_pair_mp(vector 1)" << 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 primary" << endl + << " (needs --alt-paths > 1)" << endl << " -N, --sample NAME add this sample name to output" << endl << " -R, --read-group NAME add this read group to output" << endl << " -p, --suppress-progress do not report progress to stderr" << endl @@ -435,7 +428,7 @@ int main_mpmap(int argc, char** argv) { constexpr int OPT_REF_NAME = 1039; constexpr int OPT_LINEAR_PATH = 1040; constexpr int OPT_LINEAR_INDEX = 1041; - constexpr int OPT_PROMOTE_SECONDARY = 1042; + constexpr int OPT_RESCUE_SECONDARY = 1042; string matrix_file_name; string graph_name; string gcsa_name; @@ -522,7 +515,7 @@ int main_mpmap(int argc, char** argv) { int default_num_alt_alns = 16; int num_alt_alns = default_num_alt_alns; bool agglomerate_multipath_alns = false; - bool promote_secondary = false; + bool rescue_secondary = false; double suboptimal_path_exponent = 1.25; double likelihood_approx_exp = 10.0; double likelihood_approx_exp_arg = numeric_limits::lowest(); @@ -622,7 +615,7 @@ int main_mpmap(int argc, char** argv) { {"same-strand", no_argument, 0, 'T'}, {"ref-paths", required_argument, 0, 'S'}, {"ref-name", required_argument, 0, OPT_REF_NAME}, - {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, + {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY}, {"output-fmt", required_argument, 0, 'F'}, {"snarls", required_argument, 0, 's'}, {"synth-tail-anchors", no_argument, 0, OPT_SUPPRESS_TAIL_ANCHORS}, @@ -797,8 +790,8 @@ int main_mpmap(int argc, char** argv) { reference_assembly_names.insert(optarg); break; - case OPT_PROMOTE_SECONDARY: - promote_secondary = true; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; break; case 's': @@ -1362,18 +1355,18 @@ int main_mpmap(int argc, char** argv) { ref_paths_name = ""; } - if (promote_secondary && !hts_output) { - logger.warn() << "--promote-secondary has no effect unless output format (-F) is " + if (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless output format (-F) is " << "SAM, BAM, or CRAM; ignoring." << endl; - promote_secondary = false; + rescue_secondary = false; } - if (promote_secondary && num_alt_alns < 2) { - logger.warn() << "--promote-secondary requires --alt-paths > 1; " + if (rescue_secondary && num_alt_alns < 2) { + logger.warn() << "--rescue-secondary requires --alt-paths > 1; " << "with the current setting only one alignment is produced per read, " - << "so there are no secondary alignments to promote in case of an " + << "so there are no secondary alignments to rescue in case of an " << "unsurjectable primary alignment. Ignoring." << endl; - promote_secondary = false; + rescue_secondary = false; } if (!reference_assembly_names.empty() && !hts_output) { @@ -2142,8 +2135,8 @@ int main_mpmap(int argc, char** argv) { // during distribution estimation vector> ambiguous_pair_buffer; - // Shared flag so the "no promotable secondary" warning is emitted only once. - std::atomic_flag warned_no_promotable_secondary = ATOMIC_FLAG_INIT; + // Shared flag so the "no rescuable secondary" warning is emitted only once. + std::atomic_flag warned_no_rescuable_secondary = ATOMIC_FLAG_INIT; // do unpaired multipath alignment and write to buffer function do_unpaired_alignments = [&](Alignment& alignment) { @@ -2190,12 +2183,12 @@ int main_mpmap(int argc, char** argv) { } } - if (promote_secondary) { - // If the primary failed to surject, promote the best mapped - // secondary in its place (positions kept in sync). Warns once if - // no surjectable secondary is available. - promote_secondary_mp(mp_alns, path_positions, alignment.name(), - warned_no_promotable_secondary); + if (rescue_secondary) { + // If the primary failed to surject, rescue and tag the best mapped + // secondary with YF:i:1. Warns once if no surjectable secondary + // is available. + rescue_secondary_mp(mp_alns, path_positions, alignment.name(), + warned_no_rescuable_secondary); } } @@ -2344,12 +2337,12 @@ int main_mpmap(int argc, char** argv) { } } - if (promote_secondary) { - // If both mates of the primary pair failed to surject, promote - // the best surjectable secondary pair into its place (positions - // kept in sync). Warns once if no surjectable secondary exists. - promote_secondary_pair_mp(output_mp_aln_pairs, path_positions, - alignment_1.name(), warned_no_promotable_secondary); + if (rescue_secondary) { + // If both mates of the primary pair failed to surject, rescue and + // tag the best surjectable secondary pair with YF:i:1. Warns once + // if no surjectable secondary pair is available. + rescue_secondary_pair_mp(output_mp_aln_pairs, path_positions, + alignment_1.name(), warned_no_rescuable_secondary); } } else { diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 9494a04ad9f..0c5fb171eb2 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -87,9 +87,11 @@ void help_surject(char** argv) { << " of the pre-surjected graph alignment in GR tag" << endl << " --off-ref-position annotate SAM records that become unmapped during" << endl << " surject with the nearest ref. position in the NR tag" << endl - << " --promote-secondary if a read's primary fails to surject, promote its best" << endl - << " surjectable secondary to primary instead of emitting an" << endl - << " unmapped primary (input must be collated by read name)" << 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 primary" << endl + << " (input must be collated by read name)" << endl << " -C, --compression N level for compression [0-9]" << endl << " -V, --no-validate skip checking whether alignments plausibly are" << endl << " against the provided graph" << endl @@ -153,7 +155,7 @@ int main_surject(int argc, char** argv) { constexpr int OPT_NO_PRUNE_LOW_CPLX = 1000; constexpr int OPT_OFF_REF_POS = 1001; - constexpr int OPT_PROMOTE_SECONDARY = 1002; + constexpr int OPT_RESCUE_SECONDARY = 1002; if (argc == 2) { help_surject(argv); @@ -192,7 +194,7 @@ int main_surject(int argc, char** argv) { bool validate = true; bool show_progress = false; bool left_align = false; - bool promote_secondary = false; + bool rescue_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -219,7 +221,7 @@ int main_surject(int argc, char** argv) { {"sam-output", no_argument, 0, 's'}, {"supplementary", no_argument, 0, 'u'}, {"off-ref-position", no_argument, 0, OPT_OFF_REF_POS}, - {"promote-secondary", no_argument, 0, OPT_PROMOTE_SECONDARY}, + {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY}, {"left-align", no_argument, 0, 'B'}, {"read-length", required_argument, 0, 'D'}, {"spliced", no_argument, 0, 'S'}, @@ -393,8 +395,8 @@ int main_surject(int argc, char** argv) { annotate_off_reference_pos = true; break; - case OPT_PROMOTE_SECONDARY: - promote_secondary = true; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; break; case 'h': @@ -515,21 +517,21 @@ int main_surject(int argc, char** argv) { surjector.report_supplementary = report_supplementary; surjector.left_align = left_align; surjector.multimap_to_all_paths = multimap; - surjector.promote_secondary_on_failed_surjection = promote_secondary; + surjector.rescue_secondary_on_failed_surjection = rescue_secondary; surjector.warn_about_input_collation = true; bool hts_output = (output_format == "SAM" || output_format == "BAM" || output_format == "CRAM"); - if (promote_secondary && !hts_output) { - logger.warn() << "--promote-secondary has no effect unless output is SAM, BAM, or CRAM; " + if (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless output is SAM, BAM, or CRAM; " << "ignoring." << endl; - promote_secondary = false; - surjector.promote_secondary_on_failed_surjection = false; + rescue_secondary = false; + surjector.rescue_secondary_on_failed_surjection = false; } - if (promote_secondary && input_format == "GAF") { - logger.warn() << "--promote-secondary is not supported with GAF input because GAF does not " + if (rescue_secondary && input_format == "GAF") { + logger.warn() << "--rescue-secondary is not supported with GAF input because GAF does not " << "distinguish primary and secondary alignments; ignoring." << endl; - promote_secondary = false; - surjector.promote_secondary_on_failed_surjection = false; + rescue_secondary = false; + surjector.rescue_secondary_on_failed_surjection = false; } // Count our threads @@ -662,8 +664,8 @@ int main_surject(int argc, char** argv) { alns2.emplace_back(k2 >= 0 ? std::move(surjected2[k2]) : make_placeholder(src2)); } - // No-op when promote_secondary_on_failed_surjection is false. - surjector.promote_secondary_pair_if_primary_unmapped(alns1, alns2); + // No-op when rescue_secondary_on_failed_surjection is false. + surjector.rescue_secondary_pair_if_primary_unmapped(alns1, alns2); // Append back the supplementary alignments to the end of each vector // now that secondary promotion is done. @@ -744,7 +746,7 @@ int main_surject(int argc, char** argv) { }; if (input_format == "GAM") { get_input_file(file_name, [&](istream& in) { - if (promote_secondary) { + if (rescue_secondary) { // Group consecutive pairs by first-mate name so promotion has access to // every multimap for both mates at once. Primary and secondary alignments // of the same mate carry the same name, so exact equality suffices. @@ -763,7 +765,7 @@ int main_surject(int argc, char** argv) { } }); } else { - // promote_secondary is disabled for GAF input so this effectively just surjects each pair by itself. + // rescue_secondary is disabled for GAF input so this effectively just surjects each pair by itself. auto gaf_checking_lambda = [&](Alignment& src1, Alignment& src2) { check_gaf_aln(src1); check_gaf_aln(src2); @@ -796,8 +798,8 @@ int main_surject(int argc, char** argv) { surjected_group.emplace_back(std::move(s)); } } - // No-op when promote_secondary_on_failed_surjection is false. - surjector.promote_secondary_if_primary_unmapped(surjected_group); + // No-op when rescue_secondary_on_failed_surjection is false. + surjector.rescue_secondary_if_primary_unmapped(surjected_group); alignment_emitter->emit_singles(std::move(surjected_group)); total_reads_surjected += group.size(); if (watchdog) { @@ -809,12 +811,12 @@ int main_surject(int argc, char** argv) { } }; if (input_format == "GAM") { - // Callback to determine if alignments should be grouped: if promote_secondary is enabled, we want to group by read name - // so we can access both primary and secondary alignments for a read. If promote_secondary is off, then we can just + // Callback to determine if alignments should be grouped: if rescue_secondary is enabled, we want to group by read name + // so we can access both primary and secondary alignments for a read. If rescue_secondary is off, then we can just // process each alignment independently and each 'group' will just be a single alignment. // NOTE: In order for this to work, the input file must be collated by read name. - function same_group = [promote_secondary](const Alignment& a, const Alignment& b) { - return promote_secondary && a.name() == b.name(); + function same_group = [rescue_secondary](const Alignment& a, const Alignment& b) { + return rescue_secondary && a.name() == b.name(); }; get_input_file(file_name, [&](istream& in) { vg::io::grouped_unpaired_for_each_parallel(in, process_group, same_group); diff --git a/src/surjecting_alignment_emitter.cpp b/src/surjecting_alignment_emitter.cpp index c2f2bc794af..400800924de 100644 --- a/src/surjecting_alignment_emitter.cpp +++ b/src/surjecting_alignment_emitter.cpp @@ -15,7 +15,7 @@ using namespace std; SurjectingAlignmentEmitter::SurjectingAlignmentEmitter(const PathPositionHandleGraph* graph, unordered_set paths, unique_ptr&& backing, bool prune_suspicious_anchors, bool add_graph_alignment_tag, bool report_supplementary, - bool add_off_ref_position_tag, bool left_align, bool promote_secondary) : surjector(graph), paths(paths), backing(std::move(backing)) { + bool add_off_ref_position_tag, bool left_align, bool rescue_secondary) : surjector(graph), paths(paths), backing(std::move(backing)) { // Configure the surjector surjector.prune_suspicious_anchors = prune_suspicious_anchors; @@ -23,7 +23,7 @@ SurjectingAlignmentEmitter::SurjectingAlignmentEmitter(const PathPositionHandleG surjector.report_supplementary = report_supplementary; surjector.annotate_off_reference_pos = add_off_ref_position_tag; surjector.left_align = left_align; - surjector.promote_secondary_on_failed_surjection = promote_secondary; + surjector.rescue_secondary_on_failed_surjection = rescue_secondary; } void SurjectingAlignmentEmitter::surject_alignments_in_place(vector& alns) const { @@ -107,8 +107,8 @@ void SurjectingAlignmentEmitter::emit_mapped_singles(vector>&& surject_alignments_in_place(mappings); // Each inner vector holds one read's whole multimapping (its primary // followed by its secondaries), so this is the right granularity to - // promote a mapped secondary if the primary failed to surject. - surjector.promote_secondary_if_primary_unmapped(mappings); + // rescue a secondary if the primary failed to surject. + surjector.rescue_secondary_if_primary_unmapped(mappings); } // Forward it along backing->emit_mapped_singles(std::move(alns_batch_caught)); @@ -136,9 +136,9 @@ void SurjectingAlignmentEmitter::emit_mapped_pairs(vector>&& a surject_paired_alignments_in_place(alns1_batch_caught[i], alns2_batch_caught[i], supplementary_batch.back()); // After surjection, alns1_batch_caught[i][k] and alns2_batch_caught[i][k] // are the two mates of the same alignment pair (index 0 is the primary - // pair). If both mates of the primary pair failed to surject, promote - // the best surjectable secondary pair into its place. - surjector.promote_secondary_pair_if_primary_unmapped(alns1_batch_caught[i], alns2_batch_caught[i]); + // pair). If both mates of the primary pair failed to surject, rescue + // the best surjectable secondary pair by tagging it with YF:i:1. + surjector.rescue_secondary_pair_if_primary_unmapped(alns1_batch_caught[i], alns2_batch_caught[i]); } // Forward it along backing->emit_mapped_pairs(std::move(alns1_batch_caught), std::move(alns2_batch_caught), std::move(tlen_limit_batch)); diff --git a/src/surjecting_alignment_emitter.hpp b/src/surjecting_alignment_emitter.hpp index c65287ea2fe..8a6392cfac7 100644 --- a/src/surjecting_alignment_emitter.hpp +++ b/src/surjecting_alignment_emitter.hpp @@ -37,7 +37,7 @@ class SurjectingAlignmentEmitter : public vg::io::AlignmentEmitter { unordered_set paths, unique_ptr&& backing, bool prune_suspicious_anchors = false, bool add_graph_alignment_tag = false, bool report_supplementary = false, bool add_off_ref_position_tag = false, - bool left_align = false, bool promote_secondary = false); + bool left_align = false, bool rescue_secondary = false); /// Force full length alignment in surjection resolution bool surject_subpath_global = true; diff --git a/src/surjector.cpp b/src/surjector.cpp index 0380276915b..91b4ab36195 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -588,9 +588,7 @@ using namespace std; if (source_aln->is_secondary() || (i != 0 && !is_supplementary(alns_out->back()))) { alns_out->back().set_is_secondary(true); - // Non-promoted secondaries get mapq 0 in surjected output; - // promote_secondary_if_primary_unmapped restores mapq on the - // winner if promotion occurs. + // Secondaries get mapq 0 in surjected output. alns_out->back().set_mapping_quality(0); } @@ -5315,9 +5313,9 @@ using namespace std; } } - bool Surjector::promote_secondary_if_primary_unmapped(vector& surjected_group) const { + bool Surjector::rescue_secondary_if_primary_unmapped(vector& surjected_group) const { - if (!promote_secondary_on_failed_surjection) { + if (!rescue_secondary_on_failed_surjection) { return false; } @@ -5345,7 +5343,7 @@ using namespace std; } // The primary is unmapped: look for the best-scoring mapped secondary to - // promote. Supplementary alignments are never eligible to become primary. + // rescue. Supplementary alignments are never eligible. int64_t best_idx = -1; int32_t best_score = numeric_limits::min(); for (size_t i = 0; i < surjected_group.size(); ++i) { @@ -5371,13 +5369,12 @@ using namespace std; } if (best_idx < 0) { - // No mapped secondary could be promoted. Warn once so the user knows - // why an expected promotion didn't happen. - if (!warned_about_no_promotable_secondary.test_and_set()) { + // No mapped secondary could be tagged. Warn once. + if (!warned_about_no_rescuable_secondary.test_and_set()) { #pragma omp critical (cerr) { - cerr << "warning:[Surjector] --promote-secondary: cannot find a surjectable " - << "secondary alignment for read \"" + cerr << "warning:[Surjector] --rescue-secondary: cannot find a surjectable " + << "secondary alignment to tag for read \"" << surjected_group[primary_idx].name() << "\" with an unsurjectable primary alignment."; if (warn_about_input_collation) { @@ -5390,29 +5387,17 @@ using namespace std; return false; } - // Promote the chosen secondary in place of the unmapped primary. - Alignment promoted = std::move(surjected_group[best_idx]); - promoted.set_is_secondary(false); - // Record provenance so downstream consumers can tell this happened. - set_annotation(promoted, "promoted_from_secondary", true); - - // Remove the old unmapped primary and the now-redundant secondary copy. - // Erase the higher index first so the lower index stays valid. - size_t hi = (size_t) max(primary_idx, best_idx); - size_t lo = (size_t) min(primary_idx, best_idx); - surjected_group.erase(surjected_group.begin() + hi); - surjected_group.erase(surjected_group.begin() + lo); - - // Place the promoted alignment at the front as the new primary. - surjected_group.insert(surjected_group.begin(), std::move(promoted)); + // Tag the best secondary so downstream tools can treat it as primary + // if they choose; the secondary flag itself is left unchanged. + set_annotation(surjected_group[best_idx], "rescued_secondary", true); return true; } - bool Surjector::promote_secondary_pair_if_primary_unmapped(vector& surjected1, + bool Surjector::rescue_secondary_pair_if_primary_unmapped(vector& surjected1, vector& surjected2) const { - if (!promote_secondary_on_failed_surjection) { + if (!rescue_secondary_on_failed_surjection) { return false; } @@ -5434,7 +5419,7 @@ using namespace std; return false; } - // Find the best secondary pair to promote. A candidate is eligible if at + // Find the best secondary pair to rescue. A candidate is eligible if at // least one of its mates surjected. Rank by summed mapped-mate score; // break ties toward fully-mapped pairs, then toward higher summed // mapping quality, for determinism. @@ -5473,13 +5458,13 @@ using namespace std; } if (best_idx < 0) { - // No secondary pair could be promoted. Warn once (shared with the + // No secondary pair could be tagged. Warn once (shared with the // single-end path's flag). - if (!warned_about_no_promotable_secondary.test_and_set()) { + if (!warned_about_no_rescuable_secondary.test_and_set()) { #pragma omp critical (cerr) { - cerr << "warning:[Surjector] --promote-secondary: cannot find a surjectable " - << "secondary alignment for read pair \"" + cerr << "warning:[Surjector] --rescue-secondary: cannot find a surjectable " + << "secondary alignment to tag for read pair \"" << surjected1[0].name() << "\" with an unsurjectable primary alignment."; if (warn_about_input_collation) { @@ -5492,18 +5477,10 @@ using namespace std; return false; } - // Promote: swap the chosen secondary pair into the primary slot for both - // mates, keeping the two vectors index-aligned. Update is_secondary - // flags: the promoted pair becomes primary, the demoted pair becomes - // secondary. Record provenance on the promoted mates. - surjected1[best_idx].set_is_secondary(false); - surjected2[best_idx].set_is_secondary(false); - surjected1[0].set_is_secondary(true); - surjected2[0].set_is_secondary(true); - set_annotation(surjected1[best_idx], "promoted_from_secondary", true); - set_annotation(surjected2[best_idx], "promoted_from_secondary", true); - std::swap(surjected1[0], surjected1[best_idx]); - std::swap(surjected2[0], surjected2[best_idx]); + // Tag the best secondary pair so downstream tools can treat them as + // primary if they choose; secondary flags and ordering are unchanged. + set_annotation(surjected1[best_idx], "rescued_secondary", true); + set_annotation(surjected2[best_idx], "rescued_secondary", true); return true; } diff --git a/src/surjector.hpp b/src/surjector.hpp index ed8768d0950..5b8a75dc631 100644 --- a/src/surjector.hpp +++ b/src/surjector.hpp @@ -139,9 +139,9 @@ using namespace std; mutable atomic_flag warned_about_subgraph_size = ATOMIC_FLAG_INIT; /// Have we already warned that a read's primary failed to surject and no - /// mapped secondary was available to promote? Used to emit that warning + /// mapped secondary was available to rescue? Used to emit that warning /// only once across all threads. - mutable atomic_flag warned_about_no_promotable_secondary = ATOMIC_FLAG_INIT; + mutable atomic_flag warned_about_no_rescuable_secondary = ATOMIC_FLAG_INIT; bool prune_suspicious_anchors = false; int64_t max_tail_anchor_prune = 4; @@ -192,32 +192,30 @@ using namespace std; size_t off_reference_pos_search_limit = 20000; /// If a read's primary alignment fails to surject (becomes unmapped), - /// promote the best-scoring secondary alignment that does successfully - /// surject to be the new primary in its place. This avoids emitting - /// secondary alignments alongside an unmapped primary, which violates - /// the SAM/BAM spec. Only secondary alignments already present in the - /// input are considered; no realignment is ever performed. - bool promote_secondary_on_failed_surjection = false; - - /// If true, the one-time "no promotable secondary" warning includes a + /// tag the best-scoring secondary alignment that does successfully + /// surject with YF:i:1, so downstream tools may treat it as primary. + /// The secondary flag is not changed. Only secondary alignments already + /// present in the input are considered; no realignment is ever performed. + bool rescue_secondary_on_failed_surjection = false; + + /// If true, the one-time "no rescuable secondary" warning includes a /// note that the input may not be collated by read name. Set this only /// for standalone vg surject; alignment subcommands (giraffe, map, /// mpmap) always deliver a read's full multimapping together. bool warn_about_input_collation = false; - /// If the primary in surjected_group is unmapped, promote the highest-scoring - /// mapped, non-supplementary secondary to primary. No realignment is performed. - /// Requires promote_secondary_on_failed_surjection. Emits a one-time warning - /// if no promotable secondary exists. Returns true iff a promotion occurred. - bool promote_secondary_if_primary_unmapped(vector& surjected_group) const; + /// If the primary in surjected_group is unmapped, tag the highest-scoring + /// mapped, non-supplementary secondary with YF:i:1. No realignment is performed. + /// Requires rescue_secondary_on_failed_surjection. Emits a one-time warning + /// if no rescuable secondary exists. Returns true iff a secondary was tagged. + bool rescue_secondary_if_primary_unmapped(vector& surjected_group) const; - /// Paired-end counterpart of promote_secondary_if_primary_unmapped. + /// Paired-end counterpart of rescue_secondary_if_primary_unmapped. /// surjected1[k] and surjected2[k] are the two mates of the same alignment /// pair (index 0 = primary, k > 0 = secondaries). Acts only when both mates - /// of the primary pair are unmapped; promotes the secondary pair with the - /// highest summed score. The two vectors remain the same length and - /// index-aligned. Returns true iff a promotion occurred. - bool promote_secondary_pair_if_primary_unmapped(vector& surjected1, + /// of the primary pair are unmapped; tags the secondary pair with the highest + /// summed score with YF:i:1. Returns true iff a secondary pair was tagged. + bool rescue_secondary_pair_if_primary_unmapped(vector& surjected1, vector& surjected2) const; protected: diff --git a/src/unittest/surject.cpp b/src/unittest/surject.cpp index 26497e61628..b139a2a054a 100644 --- a/src/unittest/surject.cpp +++ b/src/unittest/surject.cpp @@ -973,9 +973,9 @@ TEST_CASE("Supplementary alignments can be generated", "[surject]") { } } -TEST_CASE("Surjector promotes a mapped secondary when the primary fails to surject", "[surject][promote]") { +TEST_CASE("Surjector rescues a mapped secondary when the primary fails to surject", "[surject][rescue]") { - // A minimal graph is enough; promote_secondary_if_primary_unmapped only + // A minimal graph is enough; rescue_secondary_if_primary_unmapped only // inspects the already-surjected alignments, it does not touch the graph. bdsg::HashGraph graph; handle_t h1 = graph.create_handle("ACGT"); @@ -1009,42 +1009,48 @@ TEST_CASE("Surjector promotes a mapped secondary when the primary fails to surje vector group; group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary - REQUIRE(surjector.promote_secondary_on_failed_surjection == false); - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == false); + REQUIRE(surjector.rescue_secondary_on_failed_surjection == false); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == false); // Group is unchanged: primary still unmapped, secondary still present. REQUIRE(group.size() == 2); REQUIRE(group[0].path().mapping_size() == 0); REQUIRE(group[0].is_secondary() == false); } - surjector.promote_secondary_on_failed_surjection = true; + surjector.rescue_secondary_on_failed_surjection = true; - SECTION("Promotes the single mapped secondary into the primary slot") { + SECTION("Rescues the single mapped secondary by tagging it") { vector group; group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == true); - // Exactly one alignment remains: the promoted secondary as new primary. - REQUIRE(group.size() == 1); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == true); + // Group is unchanged in size; primary stays unmapped, secondary stays secondary. + REQUIRE(group.size() == 2); REQUIRE(group[0].is_secondary() == false); - REQUIRE(group[0].path().mapping_size() > 0); - REQUIRE(group[0].score() == 30); - REQUIRE(get_annotation(group[0], "promoted_from_secondary") == true); + REQUIRE(group[0].path().mapping_size() == 0); + REQUIRE(group[1].is_secondary() == true); + REQUIRE(group[1].path().mapping_size() > 0); + REQUIRE(group[1].score() == 30); + REQUIRE(get_annotation(group[1], "rescued_secondary") == true); } - SECTION("Promotes the highest-scoring mapped secondary") { + SECTION("Rescues the highest-scoring mapped secondary") { vector group; group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary group.push_back(make_aln("r", true, true, 20, 30, false)); // lower-scoring secondary group.push_back(make_aln("r", true, true, 45, 35, false)); // higher-scoring secondary - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == true); - REQUIRE(group.front().is_secondary() == false); - REQUIRE(group.front().score() == 45); - // The other mapped secondary is retained (still marked secondary). - REQUIRE(group.size() == 2); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == true); + // Group is unchanged in size and ordering. + REQUIRE(group.size() == 3); + REQUIRE(group[0].is_secondary() == false); + REQUIRE(group[0].path().mapping_size() == 0); + // The higher-scoring secondary (index 2) is tagged; the other is not. + REQUIRE(group[2].is_secondary() == true); + REQUIRE(group[2].score() == 45); + REQUIRE(get_annotation(group[2], "rescued_secondary") == true); REQUIRE(group[1].is_secondary() == true); REQUIRE(group[1].score() == 20); } @@ -1053,41 +1059,41 @@ TEST_CASE("Surjector promotes a mapped secondary when the primary fails to surje vector group; group.push_back(make_aln("r", false, true, 50, 60, false)); // mapped primary group.push_back(make_aln("r", true, true, 30, 40, false)); // mapped secondary - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == false); REQUIRE(group.size() == 2); REQUIRE(group[0].is_secondary() == false); REQUIRE(group[0].score() == 50); } - SECTION("No promotion when no secondary could be surjected") { + SECTION("No rescue when no secondary could be surjected") { vector group; group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary group.push_back(make_aln("r", true, false, 0, 0, false)); // unmapped secondary - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == false); REQUIRE(group.size() == 2); REQUIRE(group[0].is_secondary() == false); REQUIRE(group[0].path().mapping_size() == 0); } - SECTION("Supplementary alignments are never promoted") { + SECTION("Supplementary alignments are never rescued") { vector group; group.push_back(make_aln("r", false, false, 0, 0, false)); // unmapped primary group.push_back(make_aln("r", true, true, 90, 60, true)); // mapped but supplementary - bool promoted = surjector.promote_secondary_if_primary_unmapped(group); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_if_primary_unmapped(group); + REQUIRE(rescued == false); REQUIRE(group.size() == 2); REQUIRE(group[0].is_secondary() == false); REQUIRE(group[0].path().mapping_size() == 0); } } -TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to surject", - "[surject][promote]") { +TEST_CASE("Surjector rescues a secondary pair when the primary pair fails to surject", + "[surject][rescue]") { - // As with the single-end promotion test, the graph is only needed to - // construct a Surjector; promote_secondary_pair_if_primary_unmapped inspects + // As with the single-end rescue test, the graph is only needed to + // construct a Surjector; rescue_secondary_pair_if_primary_unmapped inspects // the already-surjected alignments and never touches the graph. bdsg::HashGraph graph; handle_t h1 = graph.create_handle("ACGT"); @@ -1114,15 +1120,15 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su return aln; }; - SECTION("Disabled by default: no promotion even if the primary pair is unmapped") { + SECTION("Disabled by default: no rescue even if the primary pair is unmapped") { vector mate1, mate2; mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary mate 1 mate2.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary mate 2 mate1.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary mate 1 mate2.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary mate 2 - REQUIRE(surjector.promote_secondary_on_failed_surjection == false); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == false); + REQUIRE(surjector.rescue_secondary_on_failed_surjection == false); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == false); REQUIRE(mate1.size() == 2); REQUIRE(mate2.size() == 2); REQUIRE(mate1[0].path().mapping_size() == 0); @@ -1131,34 +1137,36 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su REQUIRE(mate2[0].is_secondary() == false); } - surjector.promote_secondary_on_failed_surjection = true; + surjector.rescue_secondary_on_failed_surjection = true; - SECTION("Promotes a fully-mapped secondary pair into the primary slot") { + SECTION("Rescues a fully-mapped secondary pair by tagging it") { vector mate1, mate2; mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair mate2.push_back(make_aln("r", false, false, 0, 0)); mate1.push_back(make_aln("r", true, true, 30, 40)); // mapped secondary pair mate2.push_back(make_aln("r", true, true, 25, 35)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == true); - // Vectors stay the same length and index-aligned. + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == true); + // Vectors stay the same length and index-aligned; no reordering. REQUIRE(mate1.size() == 2); REQUIRE(mate2.size() == 2); - // The promoted pair is now the primary (index 0) for both mates. + // Primary pair (index 0) is unchanged: still unmapped, still primary. REQUIRE(mate1[0].is_secondary() == false); REQUIRE(mate2[0].is_secondary() == false); - REQUIRE(mate1[0].path().mapping_size() > 0); - REQUIRE(mate2[0].path().mapping_size() > 0); - REQUIRE(mate1[0].score() == 30); - REQUIRE(mate2[0].score() == 25); - REQUIRE(get_annotation(mate1[0], "promoted_from_secondary") == true); - REQUIRE(get_annotation(mate2[0], "promoted_from_secondary") == true); - // The demoted former-primary pair is now secondary. + REQUIRE(mate1[0].path().mapping_size() == 0); + REQUIRE(mate2[0].path().mapping_size() == 0); + // Secondary pair (index 1) is tagged but stays secondary. REQUIRE(mate1[1].is_secondary() == true); REQUIRE(mate2[1].is_secondary() == true); + REQUIRE(mate1[1].path().mapping_size() > 0); + REQUIRE(mate2[1].path().mapping_size() > 0); + REQUIRE(mate1[1].score() == 30); + REQUIRE(mate2[1].score() == 25); + REQUIRE(get_annotation(mate1[1], "rescued_secondary") == true); + REQUIRE(get_annotation(mate2[1], "rescued_secondary") == true); } - SECTION("Promotes the secondary pair with the highest summed mapped-mate score") { + SECTION("Rescues the secondary pair with the highest summed mapped-mate score") { vector mate1, mate2; mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair mate2.push_back(make_aln("r", false, false, 0, 0)); @@ -1166,13 +1174,15 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su mate2.push_back(make_aln("r", true, true, 20, 30)); mate1.push_back(make_aln("r", true, true, 45, 35)); // secondary pair, summed 80 mate2.push_back(make_aln("r", true, true, 35, 25)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == true); - REQUIRE(mate1[0].is_secondary() == false); - REQUIRE(mate2[0].is_secondary() == false); - // The higher summed-score pair (45 + 35 = 80) wins. - REQUIRE(mate1[0].score() == 45); - REQUIRE(mate2[0].score() == 35); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == true); + // The higher summed-score pair (45 + 35 = 80, index 2) is tagged; ordering unchanged. + REQUIRE(mate1[2].is_secondary() == true); + REQUIRE(mate2[2].is_secondary() == true); + REQUIRE(mate1[2].score() == 45); + REQUIRE(mate2[2].score() == 35); + REQUIRE(get_annotation(mate1[2], "rescued_secondary") == true); + REQUIRE(get_annotation(mate2[2], "rescued_secondary") == true); } SECTION("Prefers a fully-mapped pair over a half-mapped pair with equal summed score") { @@ -1185,30 +1195,34 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su // Fully-mapped secondary: both mapped, summed score 50 (30 + 20). mate1.push_back(make_aln("r", true, true, 30, 20)); mate2.push_back(make_aln("r", true, true, 20, 20)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == true); - // The fully-mapped pair is chosen despite the tie in summed score. - REQUIRE(mate1[0].path().mapping_size() > 0); - REQUIRE(mate2[0].path().mapping_size() > 0); - REQUIRE(mate1[0].score() == 30); - REQUIRE(mate2[0].score() == 20); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == true); + // The fully-mapped pair (index 2) is tagged; the half-mapped pair (index 1) is not. + REQUIRE(mate1[2].path().mapping_size() > 0); + REQUIRE(mate2[2].path().mapping_size() > 0); + REQUIRE(mate1[2].score() == 30); + REQUIRE(mate2[2].score() == 20); + REQUIRE(get_annotation(mate1[2], "rescued_secondary") == true); + REQUIRE(get_annotation(mate2[2], "rescued_secondary") == true); } - SECTION("Promotes a half-mapped secondary pair when no fully-mapped pair exists") { + SECTION("Rescues a half-mapped secondary pair when no fully-mapped pair exists") { vector mate1, mate2; mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair mate2.push_back(make_aln("r", false, false, 0, 0)); // Only mate 2 of the secondary surjected. mate1.push_back(make_aln("r", true, false, 0, 0)); mate2.push_back(make_aln("r", true, true, 40, 30)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == true); - REQUIRE(mate1[0].is_secondary() == false); - REQUIRE(mate2[0].is_secondary() == false); - // The mapped mate carries the alignment; its partner stays unmapped. - REQUIRE(mate1[0].path().mapping_size() == 0); - REQUIRE(mate2[0].path().mapping_size() > 0); - REQUIRE(mate2[0].score() == 40); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == true); + // Secondary pair (index 1) is tagged; both mates stay secondary and in place. + REQUIRE(mate1[1].is_secondary() == true); + REQUIRE(mate2[1].is_secondary() == true); + REQUIRE(mate1[1].path().mapping_size() == 0); + REQUIRE(mate2[1].path().mapping_size() > 0); + REQUIRE(mate2[1].score() == 40); + REQUIRE(get_annotation(mate1[1], "rescued_secondary") == true); + REQUIRE(get_annotation(mate2[1], "rescued_secondary") == true); } SECTION("No-op when either mate of the primary pair surjected") { @@ -1218,20 +1232,20 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su mate2.push_back(make_aln("r", false, false, 0, 0)); mate1.push_back(make_aln("r", true, true, 30, 40)); mate2.push_back(make_aln("r", true, true, 30, 40)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == false); REQUIRE(mate1[0].is_secondary() == false); REQUIRE(mate1[0].score() == 50); } - SECTION("No promotion when no secondary pair has a mapped mate") { + SECTION("No rescue when no secondary pair has a mapped mate") { vector mate1, mate2; mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair mate2.push_back(make_aln("r", false, false, 0, 0)); mate1.push_back(make_aln("r", true, false, 0, 0)); // unmapped secondary pair mate2.push_back(make_aln("r", true, false, 0, 0)); - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == false); REQUIRE(mate1[0].is_secondary() == false); REQUIRE(mate2[0].is_secondary() == false); REQUIRE(mate1[0].path().mapping_size() == 0); @@ -1243,8 +1257,8 @@ TEST_CASE("Surjector promotes a secondary pair when the primary pair fails to su mate1.push_back(make_aln("r", false, false, 0, 0)); // unmapped primary pair mate2.push_back(make_aln("r", false, false, 0, 0)); mate1.push_back(make_aln("r", true, true, 30, 40)); // secondary only present for mate 1 - bool promoted = surjector.promote_secondary_pair_if_primary_unmapped(mate1, mate2); - REQUIRE(promoted == false); + bool rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == false); REQUIRE(mate1.size() == 2); REQUIRE(mate2.size() == 1); }