fasta_lake.blosum

fasta_lake.blosum

BLOSUM-guided amino acid substitution scoring.

BLOSUM95 is optimized for closely related sequences (>95% identity), appropriate for single amino acid variant (SAV) discovery within a species.

For de novo error correction, we also support: - BLOSUM62: general-purpose, good for cross-species comparisons - De novo error matrix: empirical substitution rates from de novo prediction errors (mass-equivalent substitutions, terminal errors, I/L confusion)

The de novo error matrix captures what de novo sequencing engines actually get wrong — which is different from evolutionary substitution patterns. Common de novo errors: - I/L confusion (isobaric, 100% confusion rate) - N/D confusion (deamidation, +0.98 Da) - Q/E confusion (deamidation, +0.98 Da) - K/Q confusion (near-isobaric, 0.036 Da difference) - Adjacent AA swaps (mass-preserving) - F/oxidized-M (near-isobaric)

Reference: - Henikoff & Henikoff (1992) BLOSUM matrices - Miura et al. (2025) ICS-guided refinement - Lebedev, Treit, Mann (2026) de novo error taxonomy

score

score(aa1: str, aa2: str, matrix: str = 'blosum95') -> int

Substitution score for aa1 -> aa2.

Parameters:

Name	Type	Description	Default
matrix	str	'blosum95' or 'denovo'.	'blosum95'

Source code in fasta_lake/blosum.py

def score(aa1: str, aa2: str, matrix: str = "blosum95") -> int:
    """
    Substitution score for aa1 -> aa2.

    Parameters
    ----------
    matrix : str
        ``'blosum95'`` or ``'denovo'``.
    """
    if matrix == "denovo":
        return _DENOVO_CONFUSIONS.get((aa1.upper(), aa2.upper()), -5)
    return BLOSUM95.get((aa1.upper(), aa2.upper()), -4)

substitutions

substitutions(aa: str, min_score: int = 0, matrix: str = 'blosum95') -> List[Tuple[str, int]]

Possible substitutions for an amino acid, sorted by score descending.

Parameters:

Name	Type	Description	Default
aa	str	Source amino acid.	required
min_score	int	Minimum score to include.	0
matrix	str	'blosum95' or 'denovo'.	'blosum95'

Returns:

Type	Description
list of (target_aa, score)

Source code in fasta_lake/blosum.py

def substitutions(
    aa: str,
    min_score: int = 0,
    matrix: str = "blosum95",
) -> List[Tuple[str, int]]:
    """
    Possible substitutions for an amino acid, sorted by score descending.

    Parameters
    ----------
    aa : str
        Source amino acid.
    min_score : int
        Minimum score to include.
    matrix : str
        ``'blosum95'`` or ``'denovo'``.

    Returns
    -------
    list of (target_aa, score)
    """
    aa = aa.upper()
    if aa not in AMINO_ACIDS:
        return []

    m = _DENOVO_CONFUSIONS if matrix == "denovo" else BLOSUM95
    subs = []
    for target in AMINO_ACIDS:
        if target != aa:
            s = m.get((aa, target), -5 if matrix == "denovo" else -4)
            if s >= min_score:
                subs.append((target, s))
    subs.sort(key=lambda x: -x[1])
    return subs

generate_variants

generate_variants(sequence: str, position: int, min_score: int = -1, matrix: str = 'blosum95', max_per_position: int = 3) -> List[Tuple[str, str, int]]

Generate single-AA variants at a position.

Returns:

Type	Description
list of (variant_sequence, mutation_notation, score)

Source code in fasta_lake/blosum.py

def generate_variants(
    sequence: str,
    position: int,
    min_score: int = -1,
    matrix: str = "blosum95",
    max_per_position: int = 3,
) -> List[Tuple[str, str, int]]:
    """
    Generate single-AA variants at a position.

    Returns
    -------
    list of (variant_sequence, mutation_notation, score)
    """
    if position < 0 or position >= len(sequence):
        return []
    aa = sequence[position].upper()
    subs = substitutions(aa, min_score=min_score, matrix=matrix)
    variants = []
    for target, s in subs[:max_per_position]:
        var_seq = sequence[:position] + target + sequence[position + 1 :]
        notation = f"{aa}{position + 1}{target}"
        variants.append((var_seq, notation, s))
    return variants

generate_all_variants

generate_all_variants(sequence: str, min_score: int = -1, matrix: str = 'blosum95', max_per_position: int = 3) -> List[Tuple[str, str, int]]

Generate all single-substitution variants for a sequence.

Parameters:

Name	Type	Description	Default
matrix	str	'blosum95' for evolutionary variants, 'denovo' for de novo error correction.	'blosum95'

Source code in fasta_lake/blosum.py

def generate_all_variants(
    sequence: str,
    min_score: int = -1,
    matrix: str = "blosum95",
    max_per_position: int = 3,
) -> List[Tuple[str, str, int]]:
    """
    Generate all single-substitution variants for a sequence.

    Parameters
    ----------
    matrix : str
        ``'blosum95'`` for evolutionary variants, ``'denovo'`` for
        de novo error correction.
    """
    all_vars = []
    for pos in range(len(sequence)):
        all_vars.extend(generate_variants(sequence, pos, min_score, matrix, max_per_position))
    return all_vars