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..8932343378e 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, "rescued_secondary")) { + if (get_annotation(alignment, "rescued_secondary")) { + sam << "\tYF:i:1"; + } + } sam << "\n"; return sam.str(); @@ -1041,14 +1046,20 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header, string pos = get_annotation(alignment, "nearest_ref_pos"); bam_aux_append(bam, "NR", 'Z', pos.size() + 1, (uint8_t*) pos.c_str()); } - + if (has_annotation(alignment, "rescued_secondary")) { + if (get_annotation(alignment, "rescued_secondary")) { + int32_t val = 1; + bam_aux_append(bam, "YF", 'i', sizeof(int32_t), (uint8_t*) &val); + } + } + // TODO: it would be nice wrap htslib and set the other tags this way as well if (has_annotation(alignment, "tags")) { // encode the alignments SAM tags auto parsed_tags = parse_sam_tags(get_annotation(alignment, "tags")); for (const auto& tag : parsed_tags) { - if (get<0>(tag) == "AS" || get<0>(tag) == "RG" || get<0>(tag) == "SS" || get<0>(tag) == "GR" || get<0>(tag) == "NR") { + if (get<0>(tag) == "AS" || get<0>(tag) == "RG" || get<0>(tag) == "SS" || get<0>(tag) == "GR" || get<0>(tag) == "NR" || get<0>(tag) == "YF") { // we handle these tags separately continue; } diff --git a/src/hts_alignment_emitter.cpp b/src/hts_alignment_emitter.cpp index c0b8f576ca8..09fbaebcc32 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_RESCUE_SECONDARY); } } else { diff --git a/src/hts_alignment_emitter.hpp b/src/hts_alignment_emitter.hpp index d047dc262fc..b27230ac3e5 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, 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/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..f36bf097cf1 100644 --- a/src/subcommand/giraffe_main.cpp +++ b/src/subcommand/giraffe_main.cpp @@ -732,6 +732,11 @@ 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 +802,7 @@ int main_giraffe(int argc, char** argv) { constexpr int OPT_OFF_REF_POSITION = 1014; constexpr int OPT_LEFT_ALIGN = 1015; constexpr int OPT_NO_REC_MODE = 1016; + constexpr int OPT_RESCUE_SECONDARY = 1017; constexpr int OPT_HAPLOTYPE_NAME = 1100; constexpr int OPT_KFF_NAME = 1101; @@ -908,6 +914,9 @@ int main_giraffe(int argc, char** argv) { // When surjecting, should we annotate the off-reference reads with the nearest reference position? bool annotate_off_ref_position = false; + // When surjecting, if a read's primary fails to surject, should we rescue its best surjectable secondary? + bool rescue_secondary = false; + // For GAM format, should we report in named-segment space instead of node ID space? bool named_coordinates = false; @@ -1166,6 +1175,7 @@ int main_giraffe(int argc, char** argv) { {"ref-name", required_argument, 0, OPT_REF_NAME}, {"add-graph-aln", no_argument, 0, OPT_ADD_GRAPH_ALIGNMENT}, {"off-ref-position", no_argument, 0, OPT_OFF_REF_POSITION}, + {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY}, {"left-align", no_argument, 0, OPT_LEFT_ALIGN}, {"named-coordinates", no_argument, 0, OPT_NAMED_COORDINATES}, {"discard", no_argument, 0, 'n'}, @@ -1357,6 +1367,10 @@ int main_giraffe(int argc, char** argv) { case OPT_OFF_REF_POSITION: annotate_off_ref_position = true; break; + + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; + break; case OPT_LEFT_ALIGN: left_align = true; @@ -2057,6 +2071,7 @@ int main_giraffe(int argc, char** argv) { report_flag("interleaved", interleaved); report_flag("add-graph-aln", add_graph_alignment); report_flag("off-ref-position", annotate_off_ref_position); + report_flag("rescue-secondary", rescue_secondary); report_flag("left-align", left_align); report_flag("set-refpos", set_refpos); minimizer_mapper.set_refpos = set_refpos; @@ -2202,6 +2217,23 @@ int main_giraffe(int argc, char** argv) { // When surjecting, attempt to left align flags |= ALIGNMENT_EMITTER_FLAG_HTS_LEFT_ALIGN; } + if (rescue_secondary && minimizer_mapper.max_multimaps < 2) { + logger.warn() << "--rescue-secondary requires --max-multimaps > 1; " + << "with the current setting only one alignment is produced per read, " + << "so there are no secondary alignments to rescue in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + rescue_secondary = false; + } + if (rescue_secondary) { + if (paired && !interleaved) { + logger.warn() << "--rescue-secondary with two-file paired input uses paired-end " + << "promotion semantics: promotion fires only when both mates of the " + << "primary pair fail to surject, and the demoted pair is kept in the " + << "output as secondary records. Use -i if your input is interleaved." << endl; + } + // When surjecting, rescue a secondary if the primary fails to surject + flags |= ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY; + } // We send along the positional graph when we have it, and otherwise we send the GBWTGraph which is sufficient for GAF output. // TODO: What if we need both a positional graph and a NamedNodeBackTranslation??? diff --git a/src/subcommand/map_main.cpp b/src/subcommand/map_main.cpp index cc1114fc2fb..0f2b8f4f746 100644 --- a/src/subcommand/map_main.cpp +++ b/src/subcommand/map_main.cpp @@ -132,6 +132,11 @@ void help_map(char** argv) { << " --ref-paths FILE ordered list of paths in graph, one per line" << endl << " or HTSlib .dict, for HTSLib @SQ headers" << endl << " --ref-name NAME reference assembly in graph for HTSlib output" << endl + << " --rescue-secondary in HTSlib output, tag the best secondary alignment" << endl + << " (tag with YF:i:1) when primary alignment is not" << endl + << " surjectable; downstream tools may treat secondary" << endl + << " alignments tagged with YF:i:1 as if they were" << endl + << " primary (needs --max-multimaps > 1)" << endl << " -X, --compare realign -G GAM input, writing alignment with" << endl << " \"correct\" field set to overlap with input" << endl << " -v, --refpos-table for efficient testing output a table of" << endl @@ -165,6 +170,7 @@ int main_map(int argc, char** argv) { constexpr int OPT_COMMENTS_AS_TAGS = 1005; constexpr int OPT_MAX_GAP_LENGTH = 1006; constexpr int OPT_XDROP_ALIGNMENT = 1007; + constexpr int OPT_RESCUE_SECONDARY = 1008; string matrix_file_name; string seq; string qual; @@ -243,6 +249,7 @@ int main_map(int argc, char** argv) { uint32_t max_gap_length = 40; bool log_time = false; bool comments_as_tags = false; + bool rescue_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -322,6 +329,7 @@ int main_map(int argc, char** argv) { {"gaf", no_argument, 0, '%'}, {"log-time", no_argument, 0, '^'}, {"comments-as-tags", no_argument, 0, OPT_COMMENTS_AS_TAGS}, + {"rescue-secondary", no_argument, 0, OPT_RESCUE_SECONDARY}, {"help", no_argument, 0, 'h'}, {0, 0, 0, 0} }; @@ -650,6 +658,10 @@ int main_map(int argc, char** argv) { comments_as_tags = true; break; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; + break; + case 'h': case '?': /* getopt_long already printed an error message. */ @@ -800,9 +812,35 @@ int main_map(int argc, char** argv) { paths = get_sequence_dictionary(ref_paths_name, {}, reference_assembly_names, *xgidx); } + if (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless surjecting to SAM, BAM, or CRAM " + << "(--surject-to); ignoring." << endl; + rescue_secondary = false; + } + + if (rescue_secondary && max_multimaps < 2) { + logger.warn() << "--rescue-secondary requires --max-multimaps > 1; " + << "with the current setting only one alignment is produced per read, " + << "so there are no secondary alignments to rescue in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + rescue_secondary = false; + } + + if (rescue_secondary && !interleaved_input && !fastq2.empty()) { + logger.warn() << "--rescue-secondary with two-file paired input uses paired-end " + << "promotion semantics: promotion fires only when both mates of the " + << "primary pair fail to surject, and the demoted pair is kept in the " + << "output as secondary records. Use -i if your input is interleaved." << endl; + } + // Set up output to an emitter that will handle serialization and surjection + int emitter_flags = ALIGNMENT_EMITTER_FLAG_NONE; + if (rescue_secondary) { + // When surjecting, rescue a secondary if the primary fails to surject. + emitter_flags |= ALIGNMENT_EMITTER_FLAG_HTS_RESCUE_SECONDARY; + } unique_ptr 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..8a4e2d2753c 100644 --- a/src/subcommand/mpmap_main.cpp +++ b/src/subcommand/mpmap_main.cpp @@ -98,6 +98,158 @@ static void error_if_negative(const Logger& logger, double value, const string& } } +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; +} + +static bool rescue_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 rescue. + 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 tag. Warn once. + if (!warned.test_and_set()) { + #pragma omp critical (cerr) + { + cerr << "warning:[vg mpmap] --rescue-secondary was requested, but a read whose " + << "primary alignment failed to surject had no surjectable secondary alignment " + << "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; + } + } + return false; + } + + // 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 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 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 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) { + 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] --rescue-secondary was requested, but a read pair whose " + << "primary alignment failed to surject had no surjectable secondary alignment " + << "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 read pairs are left with an unmapped primary. " + << "This warning is shown only once." << endl; + } + } + return false; + } + + // Tag the best secondary pair so downstream tools can treat them as primary + // if they choose; secondary annotations and ordering are left unchanged. + output_mp_aln_pairs[best_idx].first.set_annotation("rescued_secondary", true); + output_mp_aln_pairs[best_idx].second.set_annotation("rescued_secondary", true); + 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 +294,11 @@ 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 + << " --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 @@ -271,6 +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_RESCUE_SECONDARY = 1042; string matrix_file_name; string graph_name; string gcsa_name; @@ -357,6 +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 rescue_secondary = false; double suboptimal_path_exponent = 1.25; double likelihood_approx_exp = 10.0; double likelihood_approx_exp_arg = numeric_limits::lowest(); @@ -456,6 +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}, + {"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}, @@ -629,6 +789,10 @@ int main_mpmap(int argc, char** argv) { case OPT_REF_NAME: reference_assembly_names.insert(optarg); break; + + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; + break; case 's': snarls_name = require_exists(logger, optarg); @@ -1190,6 +1354,20 @@ int main_mpmap(int argc, char** argv) { << "when output format (-F) is SAM, BAM, or CRAM." << endl; ref_paths_name = ""; } + + if (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless output format (-F) is " + << "SAM, BAM, or CRAM; ignoring." << endl; + rescue_secondary = false; + } + + 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 rescue in case of an " + << "unsurjectable primary alignment. Ignoring." << endl; + rescue_secondary = false; + } if (!reference_assembly_names.empty() && !hts_output) { logger.warn() << "Reference assembly names (--ref-name) are only used " @@ -1957,6 +2135,9 @@ int main_mpmap(int argc, char** argv) { // during distribution estimation vector> ambiguous_pair_buffer; + // 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) { #ifdef record_read_run_times @@ -2001,6 +2182,14 @@ int main_mpmap(int argc, char** argv) { get<2>(suppl_positions[j]), get<1>(suppl_positions[j])); } } + + 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); + } } if (is_rna) { @@ -2147,6 +2336,14 @@ int main_mpmap(int argc, char** argv) { } } } + + 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 { 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..0c5fb171eb2 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -21,6 +21,7 @@ #include "../xg.hpp" #include #include +#include #include "../utility.hpp" #include "../surjector.hpp" #include "../hts_alignment_emitter.hpp" @@ -86,6 +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 + << " --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 @@ -149,6 +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_RESCUE_SECONDARY = 1002; if (argc == 2) { help_surject(argv); @@ -187,6 +194,7 @@ int main_surject(int argc, char** argv) { bool validate = true; bool show_progress = false; bool left_align = false; + bool rescue_secondary = false; int c; optind = 2; // force optind past command positional argument @@ -213,6 +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}, + {"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'}, @@ -386,6 +395,10 @@ int main_surject(int argc, char** argv) { annotate_off_reference_pos = true; break; + case OPT_RESCUE_SECONDARY: + rescue_secondary = true; + break; + case 'h': case '?': help_surject(argv); @@ -504,6 +517,22 @@ int main_surject(int argc, char** argv) { surjector.report_supplementary = report_supplementary; surjector.left_align = left_align; surjector.multimap_to_all_paths = multimap; + 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 (rescue_secondary && !hts_output) { + logger.warn() << "--rescue-secondary has no effect unless output is SAM, BAM, or CRAM; " + << "ignoring." << endl; + rescue_secondary = false; + surjector.rescue_secondary_on_failed_surjection = false; + } + 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; + rescue_secondary = false; + surjector.rescue_secondary_on_failed_surjection = false; + } // Count our threads int thread_count = vg::get_thread_count(); @@ -541,126 +570,172 @@ 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)); - + 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) { + using AlnPair = pair; + function&)> process_group = [&](vector& group) { + if (group.empty()) { + return; + } try { - set_crash_context(src1.name() + ", " + src2.name()); + set_crash_context(group.front().first.name()); size_t thread_num = omp_get_thread_num(); if (watchdog) { - watchdog->check_in(thread_num, src1.name() + ", " + src2.name()); + watchdog->check_in(thread_num, group.front().first.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. + 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()); - + 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()); } - - if (validate) { - ensure_alignment_is_for_graph(logger, src1, *xgidx); - ensure_alignment_is_for_graph(logger, src2, *xgidx); + + 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; + if (validate) { + 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); + + // 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 { + extras1.emplace_back(std::move(surjected1[j])); + } + } + for (size_t j = 0; j < surjected2.size(); ++j) { + if (k2 < 0 && !is_supplementary(surjected2[j])) { + k2 = (int64_t) j; + } else { + extras2.emplace_back(std::move(surjected2[j])); + } + } + 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(); // empty refpos + return a; + }; + // 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)); } - - // 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); - + + // 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. + for (auto& a : extras1) { + alns1.emplace_back(std::move(a)); + } + for (auto& a : extras2) { + alns2.emplace_back(std::move(a)); + } + // 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); + 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; } } - for (size_t i = 0; i < surjected2.size(); ++i) { - if (!is_supplementary(surjected2[i])) { - const auto& pos = surjected2[i].refpos(0); + 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; } } - - - for (size_t i = 0; i < surjected1.size(); ++i) { - const auto& pos = surjected1[i].refpos(0); + 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(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 { + 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(surjected1[i]) && !has_annotation(surjected1[i], "mate_info")) { + 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() ? surjected2[it->second] : surjected2[strand_idx2.begin()->second]; + 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 { + } 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); + set_annotation(alns1[i], "mate_info", annotation); } - alignment_emitter->emit_single(std::move(surjected1[i])); + alignment_emitter->emit_single(std::move(alns1[i])); } } - for (size_t i = 0; i < surjected2.size(); ++i) { - const auto& pos = surjected2[i].refpos(0); + 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(surjected2[i]) || it == strand_idx1.end()) { + if (is_supplementary(alns2[i]) || it == strand_idx1.end()) { // this strand's surjection is unpaired or supplementary - if (is_supplementary(surjected2[i]) && !has_annotation(surjected2[i], "mate_info")) { + 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() ? surjected1[it->second] : surjected1[strand_idx1.begin()->second]; + 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 { + } 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); + set_annotation(alns2[i], "mate_info", annotation); } - alignment_emitter->emit_single(std::move(surjected2[i])); + alignment_emitter->emit_single(std::move(alns2[i])); } } - - total_reads_surjected += 2; + + total_reads_surjected += 2 * group.size(); if (watchdog) { watchdog->check_out(thread_num); } @@ -671,36 +746,62 @@ 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); + 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. + // 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 { + // 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); - 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 { - // We can just surject each Alignment by itself. // TODO: We don't preserve order relationships (like primary/secondary). - function lambda = [&](Alignment& src) { + function&)> process_group = [&](vector& group) { + if (group.empty()) { + return; + } try { - set_crash_context(src.name()); + set_crash_context(group.front().name()); size_t thread_num = omp_get_thread_num(); if (watchdog) { - watchdog->check_in(thread_num, src.name()); + watchdog->check_in(thread_num, group.front().name()); } - if (validate) { - ensure_alignment_is_for_graph(logger, src, *xgidx); + 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)); + } } - - // 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++; + // 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) { watchdog->check_out(thread_num); } @@ -710,13 +811,21 @@ int main_surject(int argc, char** argv) { } }; if (input_format == "GAM") { + // 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 = [rescue_secondary](const Alignment& a, const Alignment& b) { + return rescue_secondary && a.name() == b.name(); + }; 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); } @@ -780,13 +889,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; } } @@ -796,11 +905,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 @@ -814,7 +923,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; @@ -823,9 +932,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 af5a789e4a3..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) : 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,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.rescue_secondary_on_failed_surjection = rescue_secondary; } void SurjectingAlignmentEmitter::surject_alignments_in_place(vector& alns) const { @@ -71,7 +72,7 @@ void SurjectingAlignmentEmitter::surject_paired_alignments_in_place(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 + // 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)); @@ -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, 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 f876aaa62ed..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 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 dc19d9f76fc..91b4ab36195 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -588,6 +588,8 @@ using namespace std; if (source_aln->is_secondary() || (i != 0 && !is_supplementary(alns_out->back()))) { alns_out->back().set_is_secondary(true); + // Secondaries get mapq 0 in surjected output. + alns_out->back().set_mapping_quality(0); } if (annotate_with_all_path_scores) { @@ -5311,6 +5313,178 @@ using namespace std; } } + bool Surjector::rescue_secondary_if_primary_unmapped(vector& surjected_group) const { + + if (!rescue_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 + // 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) { + 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 tagged. Warn once. + if (!warned_about_no_rescuable_secondary.test_and_set()) { + #pragma omp critical (cerr) + { + 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) { + 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; + } + + // 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::rescue_secondary_pair_if_primary_unmapped(vector& surjected1, + vector& surjected2) const { + + if (!rescue_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 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. + 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 tagged. Warn once (shared with the + // single-end path's flag). + if (!warned_about_no_rescuable_secondary.test_and_set()) { + #pragma omp critical (cerr) + { + 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) { + 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; + } + + // 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; + } + 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..5b8a75dc631 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 rescue? Used to emit that warning + /// only once across all threads. + mutable atomic_flag warned_about_no_rescuable_secondary = ATOMIC_FLAG_INIT; bool prune_suspicious_anchors = false; int64_t max_tail_anchor_prune = 4; @@ -185,8 +190,34 @@ 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), + /// 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, 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 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; 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: /// Do the extra score setup for the DP-only Aligner. diff --git a/src/unittest/surject.cpp b/src/unittest/surject.cpp index 4de83ba1574..b139a2a054a 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,296 @@ TEST_CASE("Supplementary alignments can be generated", "[surject]") { } } +TEST_CASE("Surjector rescues a mapped secondary when the primary fails to surject", "[surject][rescue]") { + + // 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"); + 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.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.rescue_secondary_on_failed_surjection = true; + + 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 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[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("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 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); + } + + 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 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 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 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 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 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 rescues a secondary pair when the primary pair fails to surject", + "[surject][rescue]") { + + // 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"); + 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 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.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); + REQUIRE(mate2[0].path().mapping_size() == 0); + REQUIRE(mate1[0].is_secondary() == false); + REQUIRE(mate2[0].is_secondary() == false); + } + + surjector.rescue_secondary_on_failed_surjection = true; + + 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 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); + // 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); + // 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("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)); + 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 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") { + 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 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("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 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") { + 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 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 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 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); + 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 rescued = surjector.rescue_secondary_pair_if_primary_unmapped(mate1, mate2); + REQUIRE(rescued == 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" +