Sequence distance models

DDMC supports two interchangeable ways of scoring how well a peptide sequence matches a cluster's motif, selected via DDMC(..., distance_method=...).

ddmc.binomial

ddmc.binomial

Binomial sequence-distance model used by ddmc.clustering.DDMC.

Contains

• AAfreq / AAlist: reference amino acid frequencies and the fixed amino acid ordering used throughout the package.
• Position weight matrix helpers (position_weight_matrix, fast_position_weight_matrix, frequencies, GenerateBinarySeqID).
• Background phosphosite sequence sampling from PhosphoSitePlus (BackgroundSeqs, BackgProportions, CountPsiteTypes, and their cached loaders).
• The Binomial class: for each cluster, models how enriched each amino acid is at each position (relative to the background) using the binomial-probability approach of Schwartz & Gygi, Nat Biotechnol 2005 (doi:10.1038/nbt1146), and scores every peptide sequence against each cluster's model.

Binomial

Binomial(seqs: ndarray)

Binomial sequence-distance model, used by ddmc.clustering.DDMC when distance_method="Binomial".

For each cluster, scores how enriched each amino acid is at each position of a peptide's sequence relative to a background distribution of phosphosites, following Schwartz & Gygi, Nat Biotechnol 2005 (doi:10.1038/nbt1146).

Attributes:

Name	Type	Description
background		Background PWM (amino acid frequency per position) of shape (len(AAlist), n_pos), built from BackgroundSeqs(seqs).
n_aa		Number of amino acids (len(AAlist)).
n_pos		Number of sequence positions (11).
foreground_flat		Flattened one-hot encoding of seqs, of shape (n_seqs, n_aa * n_pos).
logWeights		Log-probability of each sequence under each cluster's current binomial model, of shape (n_seqs, n_clusters). Set to the scalar 0.0 until from_summaries is first called.

Parameters:

Name	Type	Description	Default
seqs	ndarray	The length-11 peptide sequences being clustered.	required

Source code in ddmc/binomial.py

def __init__(self, seqs: np.ndarray):
    """
    Args:
        seqs: The length-11 peptide sequences being clustered.
    """
    # Background sequences
    self.background = fast_position_weight_matrix(BackgroundSeqs(seqs))
    foreground: np.ndarray = GenerateBinarySeqID(seqs)
    self.n_aa, self.n_pos = foreground.shape[1], foreground.shape[2]
    # Flattened, float view of the one-hot foreground used for fast
    # matrix multiplication in from_summaries (replacing einsum, which
    # is much slower on the boolean input and is called every EM step).
    self.foreground_flat = foreground.reshape(foreground.shape[0], -1).astype(
        np.float32
    )

    self.logWeights = 0.0
    assert np.all(np.isfinite(self.background))
    assert np.all(np.isfinite(self.foreground_flat))

from_summaries

from_summaries(weightsIn: ndarray) -> None

Refit each cluster's binomial model from the current soft cluster assignments, and update self.logWeights with each sequence's log-probability under its (updated) cluster model.

Parameters:

Name	Type	Description	Default
weightsIn	ndarray	Soft cluster assignments (responsibilities) of shape (n_seqs, n_clusters), i.e. exp(log_resp) from the EM E step.	required

Source code in ddmc/binomial.py

def from_summaries(self, weightsIn: np.ndarray) -> None:
    """Refit each cluster's binomial model from the current soft cluster
    assignments, and update `self.logWeights` with each sequence's
    log-probability under its (updated) cluster model.

    Args:
        weightsIn: Soft cluster assignments (responsibilities) of shape
            (n_seqs, n_clusters), i.e. `exp(log_resp)` from the EM E step.
    """
    k_flat = weightsIn.T.astype(np.float32) @ self.foreground_flat
    k = k_flat.reshape(-1, self.n_aa, self.n_pos)
    betaA = np.sum(weightsIn, axis=0)[:, None, None] - k
    betaA = np.clip(betaA, 0.001, np.inf)
    probmat = sc.betainc(betaA, k + 1, 1 - self.background)
    probmat_flat = probmat.reshape(probmat.shape[0], -1).astype(np.float32)
    tempp = self.foreground_flat @ probmat_flat.T
    self.logWeights = np.log(tempp)

BackgroundSeqs

BackgroundSeqs(forseqs: ndarray) -> list[str]

Build a background data set of length-11 phosphosite motifs sampled from PhosphoSitePlus, matching the proportion of pY, pT, and pS sites found in the foreground set of sequences.

Note this PsP data set contains 51976 pY, 226131 pS, 81321 pT Source: https://www.phosphosite.org/staticDownloads.action - Phosphorylation_site_dataset.gz - Last mod: Wed Dec 04 14:56:35 EST 2019 Cite: Hornbeck PV, Zhang B, Murray B, Kornhauser JM, Latham V, Skrzypek E PhosphoSitePlus, 2014: mutations, PTMs and recalibrations. Nucleic Acids Res. 2015 43:D512-20. PMID: 25514926

Parameters:

Name	Type	Description	Default
forseqs	ndarray	The foreground peptide sequences whose pY/pS/pT proportions the background set should match.	required

Returns:

Type	Description
list[str]	Length-11 background peptide sequences sampled from PhosphoSitePlus,
list[str]	with the phosphoacceptor lowercased, in pY/pS/pT order.

Source code in ddmc/binomial.py

def BackgroundSeqs(forseqs: np.ndarray) -> list[str]:
    """Build a background data set of length-11 phosphosite motifs sampled from
    PhosphoSitePlus, matching the proportion of pY, pT, and pS sites found in
    the foreground set of sequences.

    Note this PsP data set contains 51976 pY, 226131 pS, 81321 pT
    Source: https://www.phosphosite.org/staticDownloads.action -
    Phosphorylation_site_dataset.gz - Last mod: Wed Dec 04 14:56:35 EST 2019
    Cite: Hornbeck PV, Zhang B, Murray B, Kornhauser JM, Latham V, Skrzypek E PhosphoSitePlus, 2014: mutations,
    PTMs and recalibrations. Nucleic Acids Res. 2015 43:D512-20. PMID: 25514926

    Args:
        forseqs: The foreground peptide sequences whose pY/pS/pT proportions
            the background set should match.

    Returns:
        Length-11 background peptide sequences sampled from PhosphoSitePlus,
        with the phosphoacceptor lowercased, in pY/pS/pT order.
    """
    # Get porportion of psite types in foreground set
    forw_pYn, forw_pSn, forw_pTn = CountPsiteTypes(forseqs)
    forw_tot = forw_pYn + forw_pSn + forw_pTn

    pYf = forw_pYn / forw_tot
    pSf = forw_pSn / forw_tot
    pTf = forw_pTn / forw_tot

    refseqs, backg_pYn = _load_reference_seqs()
    len_bg = len(refseqs)

    # Make sure there are enough pY peptides to meet proportions
    if backg_pYn >= len_bg * pYf:
        pYn = int(len_bg * pYf)
        pSn = int(len_bg * pSf)
        pTn = int(len_bg * pTf)

    # Not enough pYs, adjust number of peptides based on maximum number of pY peptides
    else:
        tot_p = int(backg_pYn / pYf)
        pYn = backg_pYn
        pSn = int(tot_p * pSf)
        pTn = int(tot_p * pTf)

    # Build background sequences (cached, since for fixed reference data
    # this only depends on the pY/pS/pT proportions of the foreground set)
    return list(_cached_background_proportions(pYn, pSn, pTn))

BackgProportions

BackgProportions(
    refseqs: list[str], pYn: int, pSn: int, pTn: int
) -> list[str]

Slice length-11 motifs out of the +/-7 AA reference sequences, keeping up to the requested number of pY, pS, and pT sites.

Parameters:

Name	Type	Description	Default
refseqs	list[str]	Raw +/-7 AA PhosphoSitePlus reference sequences.	required
pYn	int	Maximum number of pY motifs to keep.	required
pSn	int	Maximum number of pS motifs to keep.	required
pTn	int	Maximum number of pT motifs to keep.	required

Returns:

Type	Description
list[str]	The length-11 background motifs (phosphoacceptor lowercased),
list[str]	concatenated in pY, pS, pT order.

Source code in ddmc/binomial.py

def BackgProportions(refseqs: list[str], pYn: int, pSn: int, pTn: int) -> list[str]:
    """Slice length-11 motifs out of the +/-7 AA reference sequences, keeping
    up to the requested number of pY, pS, and pT sites.

    Args:
        refseqs: Raw +/-7 AA PhosphoSitePlus reference sequences.
        pYn: Maximum number of pY motifs to keep.
        pSn: Maximum number of pS motifs to keep.
        pTn: Maximum number of pT motifs to keep.

    Returns:
        The length-11 background motifs (phosphoacceptor lowercased),
        concatenated in pY, pS, pT order.
    """
    y_seqs: list[str] = []
    s_seqs: list[str] = []
    t_seqs: list[str] = []

    pR = ["y", "t", "s"]
    for seq in refseqs:
        if seq[7] not in pR:
            continue

        motif = str(seq)[7 - 5 : 7 + 6].upper()
        assert len(motif) == 11, f"Wrong sequence length. Sliced: {motif}, Full: {seq}"
        assert motif[5].lower() in pR, (
            f"Wrong central AA in background set. Sliced: {motif}, Full: {seq}"
        )

        if motif[5] == "Y" and len(y_seqs) < pYn:
            y_seqs.append(motif)

        if motif[5] == "S" and len(s_seqs) < pSn:
            s_seqs.append(motif)

        if motif[5] == "T" and len(t_seqs) < pTn:
            t_seqs.append(motif)

    return y_seqs + s_seqs + t_seqs

CountPsiteTypes

CountPsiteTypes(X) -> tuple[int, int, int]

Count the number of different phosphorylation types in an MS data set.

Parameters:

Name	Type	Description	Default
X	list[str]	The list of peptide sequences.	required

Returns:

Type	Description
tuple[int, int, int]	tuple[int, int, int]: The number of pY, pS, and pT sites.

Source code in ddmc/binomial.py

def CountPsiteTypes(X) -> tuple[int, int, int]:
    """Count the number of different phosphorylation types in an MS data set.

    Args:
        X (list[str]): The list of peptide sequences.

    Returns:
        tuple[int, int, int]: The number of pY, pS, and pT sites.
    """
    X = np.char.upper(X)

    # Find the center amino acid
    cA = int((len(X[0]) - 1) / 2)

    phospho_aminos = [seq[cA] for seq in X]
    pS = phospho_aminos.count("S")
    pT = phospho_aminos.count("T")
    pY = phospho_aminos.count("Y")
    return pY, pS, pT

position_weight_matrix

position_weight_matrix(
    seqs: list[str],
    pseudoC: OrderedDict[str, float] = AAfreq,
) -> Any

Build a position weight matrix (PWM) of a given set of same-length sequences.

Parameters:

Name	Type	Description	Default
seqs	list[str]	Sequences (all the same length) to build the PWM from.	required
pseudoC	OrderedDict[str, float]	Per-amino-acid pseudocounts to add before normalizing, keyed by one-letter amino acid code. Defaults to AAfreq.	AAfreq

Returns:

Type	Description
Any	A Biopython PositionWeightMatrix (amino acid frequency per
Any	position, normalized to sum to 1 down each column) of shape
Any	(len(AAlist), sequence length).

Source code in ddmc/binomial.py

def position_weight_matrix(
    seqs: list[str], pseudoC: OrderedDict[str, float] = AAfreq
) -> Any:
    """Build a position weight matrix (PWM) of a given set of same-length sequences.

    Args:
        seqs: Sequences (all the same length) to build the PWM from.
        pseudoC: Per-amino-acid pseudocounts to add before normalizing,
            keyed by one-letter amino acid code. Defaults to `AAfreq`.

    Returns:
        A Biopython `PositionWeightMatrix` (amino acid frequency per
        position, normalized to sum to 1 down each column) of shape
        (len(AAlist), sequence length).
    """
    return frequencies(seqs).normalize(pseudocounts=pseudoC)

fast_position_weight_matrix

fast_position_weight_matrix(seqs: list[str]) -> np.ndarray

Build a (len(AAlist), seq_length) PWM of a given set of same-length sequences, equivalent to position_weight_matrix but without the overhead of Biopython's general-purpose alignment machinery.

Parameters:

Name	Type	Description	Default
seqs	list[str]	Sequences, all of the same length, to build the PWM from.	required

Returns:

Type	Description
ndarray	Array of shape (len(AAlist), sequence length) giving the
ndarray	pseudocount-smoothed frequency of each amino acid at each position.

Source code in ddmc/binomial.py

def fast_position_weight_matrix(seqs: list[str]) -> np.ndarray:
    """Build a (len(AAlist), seq_length) PWM of a given set of same-length
    sequences, equivalent to `position_weight_matrix` but without the
    overhead of Biopython's general-purpose alignment machinery.

    Args:
        seqs: Sequences, all of the same length, to build the PWM from.

    Returns:
        Array of shape (len(AAlist), sequence length) giving the
        pseudocount-smoothed frequency of each amino acid at each position.
    """
    seq_len = len(seqs[0])
    # Convert to fixed-width bytes and view as a 2D uint8 array so the
    # char->index lookup is a single vectorized gather instead of a nested
    # Python loop over every character of every sequence.
    seqs_bytes = np.asarray(seqs, dtype=f"S{seq_len}")
    byte_view = seqs_bytes.view(np.uint8).reshape(len(seqs), seq_len)
    idx = _AAbyteLookup[byte_view]

    counts = np.zeros((len(AAlist), seq_len))
    for pos in range(seq_len):
        counts[:, pos] = np.bincount(idx[:, pos], minlength=len(AAlist))

    return (counts + _pseudoCounts[:, None]) / (
        counts.sum(axis=0, keepdims=True) + _pseudoCounts.sum()
    )

frequencies

frequencies(seqs: list[str]) -> Any

Build a per-position amino acid counts matrix of a given set of same-length sequences.

Parameters:

Name	Type	Description	Default
seqs	list[str]	Sequences, all of the same length, to count.	required

Returns:

Type	Description
Any	A Biopython FrequencyPositionMatrix giving the raw count of each
Any	amino acid at each position across seqs.

Source code in ddmc/binomial.py

def frequencies(seqs: list[str]) -> Any:
    """Build a per-position amino acid counts matrix of a given set of same-length sequences.

    Args:
        seqs: Sequences, all of the same length, to count.

    Returns:
        A Biopython `FrequencyPositionMatrix` giving the raw count of each
        amino acid at each position across `seqs`.
    """
    return motifs.create(seqs, alphabet="".join(AAlist)).counts

GenerateBinarySeqID

GenerateBinarySeqID(
    seqs: list[str] | ndarray,
) -> np.ndarray

Build a one-hot encoding of amino acid identity at each position, for every sequence.

Parameters:

Name	Type	Description	Default
seqs	list[str] | ndarray	Length-11 peptide sequences to encode.	required

Returns:

Type	Description
ndarray	Boolean array of shape (len(seqs), len(AAlist), 11), where
ndarray	result[i, j, k] is True if sequence i has amino acid AAlist[j]
ndarray	at position k.

Source code in ddmc/binomial.py

def GenerateBinarySeqID(seqs: list[str] | np.ndarray) -> np.ndarray:
    """Build a one-hot encoding of amino acid identity at each position, for every sequence.

    Args:
        seqs: Length-11 peptide sequences to encode.

    Returns:
        Boolean array of shape (len(seqs), len(AAlist), 11), where
        `result[i, j, k]` is True if sequence `i` has amino acid `AAlist[j]`
        at position `k`.
    """
    res = np.zeros((len(seqs), len(AAlist), 11), dtype=bool)
    for ii, seq in enumerate(seqs):
        for pos, aa in enumerate(seq):
            res[ii, AAlist.index(aa.upper()), pos] = 1
    return res

ddmc.pam250

ddmc.pam250

PAM250 sequence-distance model used by ddmc.clustering.DDMC.

Contains the PAM250 class, which scores peptide sequences against each cluster by their average PAM250 substitution-matrix similarity to the other sequences currently assigned to that cluster, and get_pam250_scores, which precomputes the full pairwise PAM250 similarity matrix used to do so.

PAM250

PAM250(seqs: list[str])

PAM250 sequence-distance model, used by ddmc.clustering.DDMC when distance_method="PAM250".

Scores each peptide sequence against a cluster by its (responsibility weighted) average pairwise PAM250 substitution score against every other sequence, using the fixed set of pairwise scores computed once at construction time.

Attributes:

Name	Type	Description
background		Pairwise PAM250 similarity matrix between all input sequences, of shape (n_seqs, n_seqs).
logWeights		Log-probability (average PAM250 score) of each sequence under each cluster's current model, of shape (n_seqs, n_clusters). Set to the scalar 0.0 until from_summaries is first called.

Parameters:

Name	Type	Description	Default
seqs	list[str]	The length-11 peptide sequences being clustered.	required

Source code in ddmc/pam250.py

def __init__(self, seqs: list[str]):
    """
    Args:
        seqs: The length-11 peptide sequences being clustered.
    """
    # Compute all pairwise distances. Cast to float32 once here rather
    # than in from_summaries, which runs every EM iteration and would
    # otherwise re-convert this (potentially large) int8 matrix each time.
    self.background = get_pam250_scores(seqs).astype(np.float32)
    self.logWeights = 0.0

from_summaries

from_summaries(weightsIn: ndarray) -> None

Update self.logWeights with each sequence's responsibility weighted average PAM250 similarity to all sequences, per cluster.

Parameters:

Name	Type	Description	Default
weightsIn	ndarray	Soft cluster assignments (responsibilities) of shape (n_seqs, n_clusters), i.e. exp(log_resp) from the EM E step.	required

Source code in ddmc/pam250.py

def from_summaries(self, weightsIn: np.ndarray) -> None:
    """Update `self.logWeights` with each sequence's responsibility
    weighted average PAM250 similarity to all sequences, per cluster.

    Args:
        weightsIn: Soft cluster assignments (responsibilities) of shape
            (n_seqs, n_clusters), i.e. `exp(log_resp)` from the EM E step.
    """
    sums = np.sum(weightsIn, axis=0)
    sums = np.clip(sums, 0.00001, np.inf)  # Avoid empty cluster divide by 0
    self.logWeights = (self.background @ weightsIn) / sums

get_pam250_scores

get_pam250_scores(seqs: list[str]) -> np.ndarray

Compute the full pairwise PAM250 similarity matrix between sequences.

Parameters:

Name	Type	Description	Default
seqs	list[str]	Sequences (all the same length) to score pairwise.	required

Returns:

Type	Description
ndarray	Symmetric array of shape (len(seqs), len(seqs)), where entry
ndarray	[i, j] is the summed PAM250 substitution score between
ndarray	seqs[i] and seqs[j] (aligned position-by-position, no gaps).

Source code in ddmc/pam250.py

def get_pam250_scores(seqs: list[str]) -> np.ndarray:
    """Compute the full pairwise PAM250 similarity matrix between sequences.

    Args:
        seqs: Sequences (all the same length) to score pairwise.

    Returns:
        Symmetric array of shape (len(seqs), len(seqs)), where entry
        `[i, j]` is the summed PAM250 substitution score between
        `seqs[i]` and `seqs[j]` (aligned position-by-position, no gaps).
    """
    pam250 = substitution_matrices.load("PAM250")
    seq_idx = np.array(
        [[pam250.alphabet.find(aa) for aa in seq] for seq in seqs],
        dtype=np.int8,
    )

    # convert to np array
    pam250m = np.array(pam250.values(), dtype=np.int8).reshape(pam250.shape)

    out = np.zeros((seq_idx.shape[0], seq_idx.shape[0]), dtype=np.int8)
    i_idx, j_idx = np.tril_indices(seq_idx.shape[0])
    out[i_idx, j_idx] = np.sum(pam250m[seq_idx[i_idx], seq_idx[j_idx]], axis=1)

    i_upper = np.triu_indices_from(out, k=1)
    out[i_upper] = out.T[i_upper]  # pylint: disable=unsubscriptable-object
    return out