diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 54adfb0c17..77207f5d93 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -67,6 +67,8 @@ void help_surject(char** argv) { << " -P, --prune-low-cplx prune short/low complexity anchors in realignment" << endl << " (on by default for long reads)" << endl << " --no-prune-low-cplx disable anchor pruning" << endl + << " -j, --prune-tail-region prune anchors lying fully inside mapper-declared" << endl + << " tail regions" << endl << " -I, --max-slide N look for offset duplicates of anchors up to N bp" << endl << " away when pruning " << "(default: " << Surjector::DEFAULT_MAX_SLIDE << ")" << endl @@ -187,6 +189,7 @@ int main_surject(int argc, char** argv) { bool validate = true; bool show_progress = false; bool left_align = false; + bool prune_tail_region = false; int c; optind = 2; // force optind past command positional argument @@ -218,6 +221,7 @@ int main_surject(int argc, char** argv) { {"spliced", no_argument, 0, 'S'}, {"prune-low-cplx", no_argument, 0, 'P'}, {"no-prune-low-cplx", no_argument, 0, OPT_NO_PRUNE_LOW_CPLX}, + {"prune-tail-region", no_argument, 0, 'j'}, {"max-slide", required_argument, 0, 'I'}, {"max-anchors", required_argument, 0, 'a'}, {"qual-adj", no_argument, 0, 'A'}, @@ -235,7 +239,7 @@ int main_surject(int argc, char** argv) { }; int option_index = 0; - c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:iGmcbsuBN:R:f:C:t:D:SPI:a:AE:LHMVw:r", + c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:iGmcbsuBN:R:f:C:t:D:SPjI:a:AE:LHMVw:r", long_options, &option_index); // Detect the end of the options. @@ -326,6 +330,10 @@ int main_surject(int argc, char** argv) { prune_anchors = false; break; + case 'j': + prune_tail_region = true; //remove surject anchors from tails + break; + case 'I': max_slide = parse(optarg); break; @@ -482,6 +490,7 @@ int main_surject(int argc, char** argv) { default_full_length_bonus); } surjector.prune_suspicious_anchors = *prune_anchors; + surjector.prune_tail_region_anchors = prune_tail_region; surjector.max_slide = max_slide; surjector.max_anchors = max_anchors; if (spliced) { diff --git a/src/surjector.cpp b/src/surjector.cpp index dc19d9f76f..8ef0053a81 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -336,11 +336,30 @@ using namespace std; } #endif - // we want to remove anchors that can be error-prone: short anchors in the tails and anchors in - // low complexity sequences + // Read mapper-declared tail lengths. These coordinates are + // relative to the stored read sequence and independent of path orientation. + // Missing annotations remain zero, making tail-region pruning a no-op. + size_t left_tail_length = 0, right_tail_length = 0; + if (source_aln) { + if (has_annotation(*source_aln, "left_tail_length")) { + left_tail_length = static_cast(get_annotation(*source_aln, "left_tail_length")); + } + if (has_annotation(*source_aln, "right_tail_length")) { + right_tail_length = static_cast(get_annotation(*source_aln, "right_tail_length")); + } + } + else { + if (source_mp_aln->has_annotation("left_tail_length")) { + left_tail_length = static_cast(*((const double*) source_mp_aln->get_annotation("left_tail_length").second)); + } + if (source_mp_aln->has_annotation("right_tail_length")) { + right_tail_length = static_cast(*((const double*) source_mp_aln->get_annotation("right_tail_length").second)); + } + } for (auto it = path_overlapping_anchors.begin(); it != path_overlapping_anchors.end(); ++it) { prune_and_trim_anchors(source_aln ? source_aln->sequence() : source_mp_aln->sequence(), - it->second.first, it->second.second); + it->second.first, it->second.second, + left_tail_length, right_tail_length); } // the surjected alignment for each path we overlapped @@ -4643,10 +4662,11 @@ using namespace std; } void Surjector::prune_and_trim_anchors(const string& sequence, vector& path_chunks, - vector>& step_ranges) const { - - if (!prune_suspicious_anchors && max_anchors > path_chunks.size()) { - // the setting don't require us to prune anything here + vector>& step_ranges, + size_t left_tail_length, size_t right_tail_length) const { + + if (!prune_suspicious_anchors && !prune_tail_region_anchors && max_anchors > path_chunks.size()) { + // the settings don't require us to prune anything here return; } @@ -4670,7 +4690,37 @@ using namespace std; } vector keep(path_chunks.size(), true); - + + // Clamp tail lengths to the read length + size_t read_length = sequence.size(); + left_tail_length = std::min(left_tail_length, read_length); + right_tail_length = std::min(right_tail_length, read_length); + size_t right_tail_begin = read_length - right_tail_length; + + if (prune_tail_region_anchors) { + // Drop anchors whose entire read interval lies inside a mapper-declared tail region. + for (int i = 0; i < path_chunks.size(); ++i) { + auto& chunk = path_chunks[i]; + size_t anchor_read_start = chunk.first.first - sequence.begin(); + size_t anchor_read_end = chunk.first.second - sequence.begin(); + + bool fully_in_left_tail = (left_tail_length > 0 && anchor_read_end <= left_tail_length); + bool fully_in_right_tail = (right_tail_length > 0 && anchor_read_start >= right_tail_begin); + + if (fully_in_left_tail || fully_in_right_tail) { +#ifdef debug_anchored_surject + cerr << "anchor " << i << " (read[" << anchor_read_start << ":" << anchor_read_end + << "]) pruned for lying fully inside " + << (fully_in_left_tail ? "left" : "right") + << " tail region" + << " (left_tail_length=" << left_tail_length + << ", right_tail_length=" << right_tail_length << ")" << endl; +#endif + keep[i] = false; + } + } + } + if (prune_suspicious_anchors) { #ifdef debug_anchored_surject cerr << "pruning suspicious anchors"; @@ -4682,7 +4732,12 @@ using namespace std; for (int i = 0; i < path_chunks.size(); ++i) { auto& chunk = path_chunks[i]; // Mark anchors that are themselves suspicious as not to be kept. - + + if (!keep[i]) { + // Already pruned by the tail-region check above; skip remaining checks. + continue; + } + // Short tails if ((chunk.first.first == path_chunks.front().first.first || chunk.first.second == path_chunks.back().first.second) // Is at either tail && (anchor_lengths[i] <= max_tail_anchor_prune || chunk.first.second - chunk.first.first <= max_tail_anchor_prune)) { // And is too short diff --git a/src/surjector.hpp b/src/surjector.hpp index 70fe220089..0ef4b8bb61 100644 --- a/src/surjector.hpp +++ b/src/surjector.hpp @@ -139,6 +139,9 @@ using namespace std; mutable atomic_flag warned_about_subgraph_size = ATOMIC_FLAG_INIT; bool prune_suspicious_anchors = false; + /// Remove anchors contained entirely within mapper-declared read tails. + /// Tail coordinates are relative to the stored read sequence. + bool prune_tail_region_anchors = false; int64_t max_tail_anchor_prune = 4; static constexpr int64_t DEFAULT_MAX_SLIDE = 6; /// Declare an anchor suspicious if it appears again at any offset up @@ -243,7 +246,8 @@ using namespace std; vector>& connections) const; void prune_and_trim_anchors(const string& sequence, vector& path_chunks, - vector>& step_ranges) const; + vector>& step_ranges, + size_t left_tail_length = 0, size_t right_tail_length = 0) const; /// Compute the widest end-inclusive interval of path positions that /// the realigned sequence could align to, or an interval where start >