Biscotti Substitution Matrix Calculation

This document describes the calculation of substitution scores from aligned sequences using Henikoff weighting, weighted amino acid pair counting, and log-odds scoring (BLOSUM-style).

1. Henikoff Sequence Weighting

To correct for unequal sequence representation while preserving overall diversity, Henikoff sequence weights are applied. For each sequence $i$, the weight is:

$$ w_i = \sum_p \frac{1}{r_p \cdot n_{a,p}} $$

Where:

$r_p$ = number of unique residues at alignment position $p$
$n_{a,p}$ = count of amino acid $a$ at position $p$

These weights ensure that sequences from overrepresented families do not dominate the calculation.

2. Weighted Amino Acid Pair Counting

For each cluster of sequences, weighted amino acid pair counts are collected across alignment positions:

$$ C(aa_1, aa_2)^+ = \begin{cases} 2 \cdot w_1 \cdot w_2 & \text{if } aa_1 \neq aa_2 \\ w_2 & \text{if } aa_1 = aa_2 \end{cases} $$

Counts from all clusters are then aggregated:

$$ C_\text{total}(aa_1, aa_2) = \sum_\text{clusters} C_\text{cluster}(aa_1, aa_2) $$

This produces a symmetric matrix of observed pair frequencies, weighted by sequence representation.

3. Observed and Expected Probabilities

The observed probability of aligning amino acids $\alpha_1$ and $\alpha_2$ is:

$$ P_\text{obs}(\alpha_1, \alpha_2) = C_\text{total}(\alpha_1, \alpha_2) + \alpha $$

where $\alpha = 0.1$ is a pseudocount to avoid zero probabilities.

The expected probability under the assumption of independent amino acid occurrence is:

$$ P_\text{exp}(\alpha_1, \alpha_2) = \begin{cases} P_\text{obs}(\alpha_1)^2 & \text{if } \alpha_1 = \alpha_2 \\ 2 \cdot P_\text{obs}(\alpha_1) \cdot P_\text{obs}(\alpha_2) & \text{if } \alpha_1 \neq \alpha_2 \end{cases} $$

4. Log-Odds Substitution Scores

Finally, substitution scores are computed as log-odds ratios:

$$ S(\alpha_1, \alpha_2) = \log_2 \frac{P_\text{obs}(\alpha_1, \alpha_2)}{P_\text{exp}(\alpha_1, \alpha_2)} $$

These scores are symmetric:

$$ S(\alpha_1, \alpha_2) = S(\alpha_2, \alpha_1) $$

The result is a $20 \times 20$ substitution matrix capturing amino acid preferences in the input sequences, weighted to account for redundancy.

Summary

Biscotti takes aligned sequences, applies Henikoff weighting to correct for overrepresentation, counts weighted amino acid pairs, and computes log-odds substitution scores with pseudocounts. The output is a robust substitution matrix suitable for downstream bioinformatics analyses.