mkt.databases.pssm

Position-specific residue statistics: background-relative information content with a substitution-aware (Henikoff) pseudocount prior.

General-purpose MSA-column conservation scoring, independent of any particular alignment. A column’s information content is the KL divergence (in bits) of its smoothed residue distribution against an amino-acid background, where the smoothing imputes unseen residues toward the substitution neighbors of the residues actually observed (Henikoff & Henikoff 1996, as in PSI-BLAST PSSM construction). This credits conservative substitution (e.g. K/R) and behaves sensibly at low effective-count columns, unlike a bare consensus-frequency threshold or plain Shannon entropy.

SubstitutionPseudocounts wraps a Biopython substitution matrix (default BLOSUM62) and provides the prior only – the conditional target probabilities, the per-column pseudocount / smoothed frequencies, and the information content against a supplied background. Estimating a dataset-specific background and analyzing particular columns is the caller’s responsibility.

Module Attributes

`STR_AA20`	The 20 standard amino acids, in the order used for all count/frequency vectors.
`DICT_AA20_INDEX`	Amino acid -> index into a length-20 vector (`STR_AA20` order).
`STR_GAP_CHARS`	Characters treated as gap/unknown and excluded from column conservation.
`ARR_SWISSPROT_BG`	Swiss-Prot background frequencies q_a; the default KL background.

Functions

`column_conservation`(sequences[, gap_chars, ...])	Per-column consensus residue and its frequency for a set of aligned sequences.
`column_counts`(sequences[, gap_chars, weights])	Per-column (weighted) residue counts for a set of aligned sequences.
`column_information_content`(sequences[, ...])	Per-column background-relative information content (bits) for an MSA.
`consurf_grade_boundaries`(scores[, n_bins])	Interior quantile boundaries splitting `scores` into `n_bins` grades.
`consurf_grades`(scores[, n_bins, bool_ascending])	Assign ConSurf-style discrete conservation grades to continuous scores.
`henikoff_weights`(sequences[, gap_chars])	Henikoff & Henikoff (1994) position-based sequence weights.

Classes

SubstitutionPseudocounts([matrix, beta, ...])

Substitution-aware pseudocount / information-content scorer for MSA columns.

mkt.databases.pssm.ARR_SWISSPROT_BG = array([0.08256605, 0.05534428, 0.04063251, 0.05464371, 0.01371097, 0.03933147, 0.06745396, 0.07085669, 0.02271817, 0.05964772, 0.09667734, 0.05854684, 0.02421938, 0.0386309 , 0.04703763, 0.0657526 , 0.05344275, 0.01080865, 0.02922338, 0.068755 ])

Swiss-Prot background frequencies q_a; the default KL background.

Type:: np.ndarray

mkt.databases.pssm.DICT_AA20_INDEX = {'A': 0, 'C': 4, 'D': 3, 'E': 6, 'F': 13, 'G': 7, 'H': 8, 'I': 9, 'K': 11, 'L': 10, 'M': 12, 'N': 2, 'P': 14, 'Q': 5, 'R': 1, 'S': 15, 'T': 16, 'V': 19, 'W': 17, 'Y': 18}

Amino acid -> index into a length-20 vector (STR_AA20 order).

Type:: dict[str, int]

mkt.databases.pssm.STR_AA20 = 'ARNDCQEGHILKMFPSTWYV'

The 20 standard amino acids, in the order used for all count/frequency vectors.

Type:: str

mkt.databases.pssm.STR_GAP_CHARS = '-X'

Characters treated as gap/unknown and excluded from column conservation.

Type:: str

class mkt.databases.pssm.SubstitutionPseudocounts(matrix: str = 'BLOSUM62', beta: float = 10.0, background: ndarray | None = None)[source]

Bases: object

Substitution-aware pseudocount / information-content scorer for MSA columns.

Parameters:

matrix (str) – Biopython substitution matrix name (default "BLOSUM62"); its implied background and conditional target probabilities P(a | b) are recovered from the matrix itself (no hardcoded frequency tables).
beta (float) – Henikoff/PSI-BLAST pseudocount strength; the total pseudocount mass is beta * (distinct residue types observed) per column (default 10).
background (np.ndarray | None) – KL background q_a (length-20, STR_AA20 order) for column_information(); defaults to ARR_SWISSPROT_BG. Supply a dataset-specific composition for “surprising relative to this dataset”.

__init__(matrix: str = 'BLOSUM62', beta: float = 10.0, background: ndarray | None = None) → None

_build_target_probs() → ndarray[source]

Conditional target probabilities P(a | b) recovered from the matrix.

From the log-odds scores S_ab (half-bits) form the odds Q_ab = 2^(S_ab/2). The implied background f is the distribution whose joint f_a f_b Q_ab has marginal f – i.e. Q f = 1 – and P(a | b) = f_a Q_ab, column-normalized to absorb the log-odds rounding.

background: ndarray | None = None

beta: float = 10.0

column_information(counts: ndarray) → float | None[source]

Background-relative information content of a column, in bits.

IC = sum_a p_a log2(p_a / q_a) (KL of the smoothed column against background); >= 0, higher = more conserved-and-surprising. An empty column returns None rather than scoring the bare pseudocount against the background.

Parameters:: counts (np.ndarray) – Length-20 (weighted) residue counts n_a.
Returns:: Information content in bits, or None for an empty column.
Return type:: float | None

matrix: str = 'BLOSUM62'

pseudocounts(counts: ndarray) → ndarray[source]

Substitution-aware pseudocount vector B * g_a for a column.

g_a = sum_b (n_b / N) P(a | b) spreads each observed residue’s mass onto its substitution neighbors; B = beta * (# distinct residue types observed).

Parameters:: counts (np.ndarray) – Length-20 (weighted) residue counts n_a in STR_AA20 order.
Returns:: Length-20 pseudocount vector B * g_a (zeros for an empty column).
Return type:: np.ndarray

smoothed_freqs(counts: ndarray) → ndarray | None[source]

Posterior-mean residue frequencies p_a = (n_a + B g_a) / (N + B).

Returns:: Length-20 smoothed distribution, or None for an empty column (N = 0).
Return type:: np.ndarray | None

property target_probs: ndarray: 20 x 20 conditional target probabilities P[a, b] = P(a | b) (columns sum to 1).

mkt.databases.pssm.column_conservation(sequences: list[str], gap_chars: str = '-X', weights: list[float] | None = None) → list[tuple[str | None, float]][source]

Per-column consensus residue and its frequency for a set of aligned sequences.

The column-level conservation primitive: any subset of aligned sequences (e.g. a whole MSA, or the members of a clade in a hierarchical tree) can be scored without instantiating a stateful analyzer. Gap/unknown characters in gap_chars are excluded; a column that is all gaps yields (None, 0.0).

Optional per-sequence weights generalize the unweighted residue counts to a weighted consensus (e.g. Henikoff sequence weights to down-weight redundant subfamilies); None reproduces the unweighted counts exactly. Because the reported fraction is a ratio of summed weights it is invariant to the overall weight scale.

Parameters:

sequences (list[str]) – Aligned sequence strings of equal length.
gap_chars (str) – Characters treated as gap/unknown and excluded from the consensus (default: STR_GAP_CHARS).
weights (list[float] | None) – Optional per-sequence weights (same length as sequences). None weights every sequence equally (default: None).

Returns:

One (consensus_aa, fraction) tuple per column, where fraction is the consensus residue’s (weighted) frequency among the observed (non-gap) residues.

Return type:

list[tuple[str | None, float]]

mkt.databases.pssm.column_counts(sequences: list[str], gap_chars: str = '-X', weights: list[float] | None = None) → ndarray[source]

Per-column (weighted) residue counts for a set of aligned sequences.

The counts primitive underlying background-relative scoring: builds one length-20 count vector per alignment column (in STR_AA20 order) from any subset of aligned sequences, without instantiating a stateful analyzer. Gap/unknown characters in gap_chars (and any residue outside STR_AA20) are excluded. Optional per-sequence weights (e.g. henikoff_weights()) generalize the raw counts to weighted counts; None reproduces the unweighted counts.

Parameters:

sequences (list[str]) – Aligned sequence strings of equal length.
gap_chars (str) – Characters treated as gap/unknown and excluded (default: STR_GAP_CHARS).
weights (list[float] | None) – Optional per-sequence weights (same length as sequences); None weights every sequence equally (default: None).

Returns:

(n_columns, 20) array of (weighted) residue counts in STR_AA20 order.

Return type:

np.ndarray

mkt.databases.pssm.column_information_content(sequences: list[str], scorer: SubstitutionPseudocounts | None = None, gap_chars: str = '-X', weights: list[float] | None = None) → list[float | None][source]

Per-column background-relative information content (bits) for an MSA.

Convenience wrapper mirroring column_conservation(): builds the weighted per-column counts (column_counts()) and scores each column with SubstitutionPseudocounts.column_information(), so any subset of aligned sequences (a whole MSA, a clade, a bootstrap resample) is scored through the single shared implementation. Pass henikoff_weights() as weights to down-weight redundant subfamilies, and a scorer with a dataset-specific background for “surprising relative to this dataset”.

Parameters:

sequences (list[str]) – Aligned sequence strings of equal length.
scorer (SubstitutionPseudocounts | None) – Pseudocount/information scorer; a default SubstitutionPseudocounts (BLOSUM62, Swiss-Prot background) is built when None.
gap_chars (str) – Characters treated as gap/unknown and excluded (default: STR_GAP_CHARS).
weights (list[float] | None) – Optional per-sequence weights (default: None -> unweighted).

Returns:

One information-content value (bits) per column; None for an all-gap column.

Return type:

list[float | None]

mkt.databases.pssm.consurf_grade_boundaries(scores: list[float | None], n_bins: int = 9) → ndarray[source]

Interior quantile boundaries splitting scores into n_bins grades.

The ConSurf grading is equal-frequency (quantile) binning of a continuous conservation score. This returns the n_bins - 1 interior score cut points (ascending), so a caller can map a specific grade cutoff back to a score threshold (e.g. the score at the grade 6|7 boundary). None/NaN scores are ignored when estimating the quantiles.

Parameters:

scores (list[float | None]) – Per-position conservation scores (any continuous scale).
n_bins (int) – Number of grades (default 9, the ConSurf convention).

Returns:

Ascending array of the n_bins - 1 interior quantile boundaries; empty if no finite scores are supplied.

Return type:

numpy.ndarray

mkt.databases.pssm.consurf_grades(scores: list[float | None], n_bins: int = 9, bool_ascending: bool = True) → list[int | None][source]

Assign ConSurf-style discrete conservation grades to continuous scores.

Bins the finite scores into n_bins equal-frequency (quantile) grades – the scheme ConSurf uses to turn a continuous conservation score into a small ordinal scale. With bool_ascending (higher score = more conserved, e.g. information content or percent identity) the most-conserved bin is grade n_bins and the most-variable is grade 1; pass bool_ascending=False for rate-type scores where lower = more conserved (e.g. Rate4Site). Grade n_bins is always the most conserved end regardless.

Parameters:

scores (list[float | None]) – Per-position conservation scores; None/NaN entries get grade None.
n_bins (int) – Number of grades (default 9, the ConSurf convention).
bool_ascending (bool) – Whether a higher score means more conserved (default True).

Returns:

One grade in 1..n_bins per input score (None where the score was missing).

Return type:

list[int | None]

mkt.databases.pssm.henikoff_weights(sequences: list[str], gap_chars: str = '-X') → ndarray[source]

Henikoff & Henikoff (1994) position-based sequence weights.

Down-weights redundant (near-duplicate) sequences so a few over-represented subfamilies do not dominate a column’s composition – the standard correction for phylogenetic non-independence in a large MSA. Each sequence’s weight is sum_columns 1 / (k_col * n_{col, residue}) over its non-gap positions, where k_col is the number of distinct residue types in the column and n_{col, residue} the count of this sequence’s residue there. Returned unnormalized (the scale cancels in any weighted frequency).

Parameters:

sequences (list[str]) – Aligned sequence strings of equal length.
gap_chars (str) – Characters treated as gap/unknown and skipped (default: STR_GAP_CHARS).

Returns:

Length-len(sequences) array of unnormalized sequence weights.

Return type:

np.ndarray