Skip to content

Repository navigation

Support symmetric per-input positional weights #10

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

Closed
mikessh wants to merge 1 commit into from
Closed

Support symmetric per-input positional weights #10

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
14 changes: 14 additions & 0 deletions docs/gapblock.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -400,3 +400,17 @@ or offset, including when a matching junction has a positive group penalty.
Omitting all three options preserves the existing scorer. Group offsets are not
exposed by top-k or paired reducers; the paired maximum-score contract is
unchanged. These additions supply geometry, not statistical calibration.

Per-input positional geometry
-----------------------------

``score_matrix``, ``ball_batch`` and ``count_batch`` accept
``query_position_weights`` and ``reference_position_weights``: nonnegative integer
rows matching each full input sequence. Each aligned substitution uses the larger
of its query and reference weights. This rule is symmetric when the substitution
matrix is symmetric; gap costs and exact identity exclusion remain unchanged.
Supply both axes and do not combine them with ``position_weights_by_length``.

``paired_sum_count_batch`` accepts separate ``_alpha`` and ``_beta`` versions of
these arguments. Each lane retains its own coordinates; distances sum against the
same linked reference row. These are geometry options, not confidence estimates.
25 changes: 20 additions & 5 deletions include/seqtree/seqtree.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,9 @@ std::vector<int32_t> align_dist_matrix(const std::vector<std::string>& queries,
// the diagonal zero. `matrix` may be null for unit cost. threads <= 0 => hardware_concurrency.
// Optional integer weights are indexed [longer length][longer position], symmetric in q/r.
// All input frames must be covered; empty tables retain the existing unweighted scorer.
// Alternatively, per-input rows index each sequence's full coordinates. Aligned
// substitution costs use max(query weight, reference weight); gaps are unchanged.
// Both axes are required; these rows cannot combine with length profiles.
// Optional group IDs index a nonnegative rectangular penalty table added per pair.
// Exact exclusion remains encoded junction identity, independent of these groups.
std::vector<int32_t> gapblock_matrix(const std::vector<std::string>& queries,
Expand All @@ -233,7 +236,9 @@ std::vector<int32_t> gapblock_matrix(const std::vector<std::string>& queries,
const std::vector<std::vector<int32_t>>& weights = {},
const std::vector<int32_t>& query_groups = {},
const std::vector<int32_t>& reference_groups = {},
const std::vector<std::vector<int32_t>>& group_distances = {});
const std::vector<std::vector<int32_t>>& group_distances = {},
const std::vector<std::vector<int32_t>>& query_weights = {},
const std::vector<std::vector<int32_t>>& reference_weights = {});

// Exhaustive ball rows in reference order; hit edit fields are unset (zero).
std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>& queries,
Expand All @@ -243,7 +248,9 @@ std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>
const std::vector<std::vector<int32_t>>& weights = {},
const std::vector<int32_t>& query_groups = {},
const std::vector<int32_t>& reference_groups = {},
const std::vector<std::vector<int32_t>>& group_distances = {});
const std::vector<std::vector<int32_t>>& group_distances = {},
const std::vector<std::vector<int32_t>>& query_weights = {},
const std::vector<std::vector<int32_t>>& reference_weights = {});

// Exhaustive bounded reductions using the same gapblock cell. Top-k ties use
// (penalty, ref_id); paired penalty is max(lane scores). Exact exclusion uses
Expand All @@ -264,15 +271,19 @@ std::vector<std::vector<uint64_t>> gapblock_count_batch(const std::vector<std::s
const std::vector<std::vector<int32_t>>& weights = {},
const std::vector<int32_t>& query_groups = {},
const std::vector<int32_t>& reference_groups = {},
const std::vector<std::vector<int32_t>>& group_distances = {});
const std::vector<std::vector<int32_t>>& group_distances = {},
const std::vector<std::vector<int32_t>>& query_weights = {},
const std::vector<std::vector<int32_t>>& reference_weights = {});
std::pair<std::vector<std::vector<uint64_t>>, std::vector<std::vector<uint64_t>>> gapblock_count_mass_batch(const std::vector<std::string>& queries,
const std::vector<std::string>& refs, const std::vector<std::vector<int32_t>>& thresholds,
Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge,
const std::vector<int32_t>& prior, uint32_t width, int threads, bool exclude_exact,
const std::vector<std::vector<int32_t>>& weights = {},
const std::vector<int32_t>& query_groups = {},
const std::vector<int32_t>& reference_groups = {},
const std::vector<std::vector<int32_t>>& group_distances = {});
const std::vector<std::vector<int32_t>>& group_distances = {},
const std::vector<std::vector<int32_t>>& query_weights = {},
const std::vector<std::vector<int32_t>>& reference_weights = {});
std::vector<std::vector<uint64_t>> gapblock_paired_count_batch(
const std::vector<std::string>& qa, const std::vector<std::string>& qb,
const std::vector<std::string>& ra, const std::vector<std::string>& rb,
Expand All @@ -290,7 +301,11 @@ std::pair<std::vector<std::vector<uint64_t>>, std::vector<std::vector<uint64_t>>
const std::vector<int32_t>& query_groups_alpha, const std::vector<int32_t>& query_groups_beta,
const std::vector<int32_t>& reference_groups_alpha, const std::vector<int32_t>& reference_groups_beta,
const std::vector<std::vector<int32_t>>& group_distances_alpha,
const std::vector<std::vector<int32_t>>& group_distances_beta);
const std::vector<std::vector<int32_t>>& group_distances_beta,
const std::vector<std::vector<int32_t>>& query_weights_alpha = {},
const std::vector<std::vector<int32_t>>& query_weights_beta = {},
const std::vector<std::vector<int32_t>>& reference_weights_alpha = {},
const std::vector<std::vector<int32_t>>& reference_weights_beta = {});

// ---------------------------------------------------------------------------------------
// Plain string edit distances on raw characters, unit costs. Deliberately alphabet-agnostic:
Expand Down
46 changes: 39 additions & 7 deletions python/seqtree/gapblock.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -527,6 +527,7 @@ def score_matrix(
threads: int = 0,
position_weights_by_length=None,
query_group_ids=None, reference_group_ids=None, group_distances=None,
query_position_weights=None, reference_position_weights=None,
) -> ScoreMatrix:
"""Gap-block penalty of every query against every reference, in C++ with the GIL released.

Expand All @@ -541,6 +542,11 @@ def score_matrix(
position_weights_by_length: Optional mapping of full sequence lengths to nonnegative
integer weights. Weights index the longer sequence; both axes must be covered.
None or an empty mapping preserves unweighted defaults.
query_position_weights: Optional nonnegative integer rows, one per query and residue.
Supply with reference_position_weights; incompatible with length profiles.
Substitutions use max(query_weight, reference_weight) at aligned coordinates.
Gap costs remain unchanged.
reference_position_weights: Per-reference rows with the same contract.
query_group_ids: Optional integer IDs, one per query; supply all group options together.
reference_group_ids: Optional integer IDs, one per reference.
group_distances: Nonempty rectangular nonnegative integer penalty table.
Expand Down Expand Up @@ -579,7 +585,8 @@ def score_matrix(
cube = _prior_cube(gap_prior, width) if gap_prior is not None else []
return _gapblock_matrix(q, r, alphabet, matrix, gap_open, gap_extend, cube, width, threads,
_position_weights(position_weights_by_length, q, r),
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r))
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r),
*_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r))


def _batch_integer(value, name, lower=0, upper=(1 << 31) - 1):
Expand Down Expand Up @@ -644,6 +651,22 @@ def _position_weights(weights, *axes):
return [tables.get(length, []) for length in range(max(tables) + 1)]


def _sequence_weights(query_weights, reference_weights, by_length, queries, refs):
"""Integer full-coordinate weights; aligned substitutions use the larger weight."""
if query_weights is None and reference_weights is None:
return [], []
if query_weights is None or reference_weights is None:
raise ValueError("both query and reference positional weights are required")
if by_length:
raise ValueError("length and per-input weights cannot be combined")
def validate(rows, axis):
rows = [list(row) for row in rows]
if len(rows) != len(axis) or any(len(row) != len(seq) for row, seq in zip(rows, axis)):
raise ValueError("per-input weights must match their sequence axis and lengths")
return [[_batch_integer(value, "position weight") for value in row] for row in rows]
return validate(query_weights, queries), validate(reference_weights, refs)


def _threshold_rows(thresholds, n_queries):
rows = [list(row) for row in thresholds]
if len(rows) != n_queries:
Expand Down Expand Up @@ -677,7 +700,8 @@ def topk_batch(queries, refs, k=10, matrix=None, gap_open=None, gap_extend=1,

def count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1,
gap_prior=None, alphabet="aa", threads=0, exclude_exact=False, position_weights_by_length=None,
query_group_ids=None, reference_group_ids=None, group_distances=None, linear_mass=False):
query_group_ids=None, reference_group_ids=None, group_distances=None, linear_mass=False,
query_position_weights=None, reference_position_weights=None):
"""Count all reference scores <= each query's requested integer thresholds.

Optional positional weights and additive group distances match :func:`score_matrix`.
Expand All @@ -702,13 +726,15 @@ def count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_exten
reduction = _gapblock_count_mass_batch if linear_mass else _gapblock_count_batch
return reduction(q, r, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact,
_position_weights(position_weights_by_length, q, r),
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r))
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r),
*_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r))


def ball_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1,
gap_prior=None, alphabet="aa", threads=0, exclude_exact=False,
position_weights_by_length=None,
query_group_ids=None, reference_group_ids=None, group_distances=None):
query_group_ids=None, reference_group_ids=None, group_distances=None,
query_position_weights=None, reference_position_weights=None):
"""Exhaustive bounded ball rows, reference order, unset edit fields; no Q*N matrix.

One integer cutoff per query. Negative cutoffs yield empty rows. Optional integer
Expand All @@ -724,7 +750,8 @@ def ball_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend
go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads)
return _gapblock_ball_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width,
threads, exclude_exact, _position_weights(position_weights_by_length, q, r),
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r))
*_group_options(query_group_ids, reference_group_ids, group_distances, q, r),
*_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r))


def paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10, matrix=None,
Expand Down Expand Up @@ -769,7 +796,9 @@ def paired_sum_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, mat
position_weights_by_length_alpha=None, position_weights_by_length_beta=None,
query_group_ids_alpha=None, query_group_ids_beta=None,
reference_group_ids_alpha=None, reference_group_ids_beta=None,
group_distances_alpha=None, group_distances_beta=None):
group_distances_alpha=None, group_distances_beta=None,
query_position_weights_alpha=None, query_position_weights_beta=None,
reference_position_weights_alpha=None, reference_position_weights_beta=None):
"""Return (counts, linear masses) for linked rows under summed lane distances.

Shared matrix/gap costs/prior and separate full-coordinate lane weights and
Expand All @@ -788,11 +817,14 @@ def paired_sum_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, mat
go, ge, cube, width, threads = _batch_options((qa, qb, ra, rb), matrix, gap_open, gap_extend, gap_prior, threads)
ag = _group_options(query_group_ids_alpha, reference_group_ids_alpha, group_distances_alpha, qa, ra)
bg = _group_options(query_group_ids_beta, reference_group_ids_beta, group_distances_beta, qb, rb)
aw = _sequence_weights(query_position_weights_alpha, reference_position_weights_alpha, position_weights_by_length_alpha, qa, ra)
bw = _sequence_weights(query_position_weights_beta, reference_position_weights_beta, position_weights_by_length_beta, qb, rb)
return _gapblock_paired_sum_count_batch(qa, qb, ra, rb, thresholds, alphabet, matrix, go, ge,
cube, width, threads, exclude_exact,
_position_weights(position_weights_by_length_alpha, qa, ra),
_position_weights(position_weights_by_length_beta, qb, rb),
ag[0], bg[0], ag[1], bg[1], ag[2], bg[2])
ag[0], bg[0], ag[1], bg[1], ag[2], bg[2],
aw[0], bw[0], aw[1], bw[1])


class GapBlockIndex:
Expand Down
Loading
Loading