Skip to content
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

Jump to bottom

Rework single-end MAPQ computation #328

Merged
merged 4 commits into from
Nov 15, 2023
Merged

Rework single-end MAPQ computation #328

merged 4 commits into from
Nov 15, 2023

Commits on Nov 15, 2023

  1. Rework single-end MAPQ computation 

    First, there was a bug in MAPQ computation due to initializing
best_score to -1000. In the first iteration, min_mapq_diff is updated thus:

min_mapq_diff = std::max(0, alignment.score - best_score);

This is just alignment.score + 1000, so most of the time, min_mapq_diff
ends up being something like 1600 or so for typical read lengths. (This is
clamped to 60 in the output, but still.)

Here, we keep track of the highest and second-highest alignment score in
best_score and second_best_score. Then, mapq is computed such that it is
* 60 for second_best_score=0 and
* 0 for second_best_score=best_score,
where anything in between is interpolated.

See #25
    @marcelm
    marcelm committed Nov 15, 2023
    Configuration menu
    Copy the full SHA
    0376572 View commit details
    Browse the repository at this point in the history

  2. Remove Alignment::mapq field 

    Keep mapping quality separate from the alignment itself appears to make more
sense:
- Mapping quality for an alignment is not a property of an alignment per se,
  but is derived from its relationship to other alignments.
- get_alignment returns an Alignment, but leaves mapq uninitialized, that
  is, we must remember to fill it in afterwards, which can easily be
  forgotten.
- Sam::add_pair is passed two mapq values and also two Alignment instances,
  which also contain mapq values. It is unclear for the caller which of the
  values is chosen.
- Adding a mapq parameter to Sam::add makes it mirror Sam::add_pair.
    @marcelm
    marcelm committed Nov 15, 2023
    Configuration menu
    Copy the full SHA
    aeb7308 View commit details
    Browse the repository at this point in the history

  3. Changelog entry

    @marcelm
    marcelm committed Nov 15, 2023
    Configuration menu
    Copy the full SHA
    c9c685b View commit details
    Browse the repository at this point in the history

  4. Ensure we never report MAPQ=0 if best and 2nd best score differs

    @marcelm
    marcelm committed Nov 15, 2023
    Configuration menu
    Copy the full SHA
    9120362 View commit details
    Browse the repository at this point in the history