Skip to content

vg surject: prune anchors in mapper-declared tail regions #4986

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Sagorikanag wants to merge 1 commit into
base: master
Choose a base branch
from
+79 −11
Open

vg surject: prune anchors in mapper-declared tail regions #4986

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion src/subcommand/surject_main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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'},
Expand All @@ -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.
Expand Down Expand Up @@ -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<int64_t>(optarg);
break;
Expand Down Expand Up @@ -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) {
Expand Down
73 changes: 64 additions & 9 deletions src/surjector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<size_t>(get_annotation<double>(*source_aln, "left_tail_length"));
}
if (has_annotation(*source_aln, "right_tail_length")) {
right_tail_length = static_cast<size_t>(get_annotation<double>(*source_aln, "right_tail_length"));
}
}
else {
if (source_mp_aln->has_annotation("left_tail_length")) {
left_tail_length = static_cast<size_t>(*((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<size_t>(*((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
Expand Down Expand Up @@ -4643,10 +4662,11 @@ using namespace std;
}

void Surjector::prune_and_trim_anchors(const string& sequence, vector<path_chunk_t>& path_chunks,
vector<pair<step_handle_t, step_handle_t>>& step_ranges) const {

if (!prune_suspicious_anchors && max_anchors > path_chunks.size()) {
// the setting don't require us to prune anything here
vector<pair<step_handle_t, step_handle_t>>& 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;
}

Expand All @@ -4670,7 +4690,37 @@ using namespace std;
}

vector<bool> 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";
Expand All @@ -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
Expand Down
6 changes: 5 additions & 1 deletion src/surjector.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -243,7 +246,8 @@ using namespace std;
vector<tuple<size_t, size_t, int32_t>>& connections) const;

void prune_and_trim_anchors(const string& sequence, vector<path_chunk_t>& path_chunks,
vector<pair<step_handle_t, step_handle_t>>& step_ranges) const;
vector<pair<step_handle_t, step_handle_t>>& 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 >
Expand Down
Loading