ddmc.clustering

The core DDMC model and its supporting functions.

ddmc.clustering

Dual data and motif clustering (DDMC).

Contains the DDMC model itself — a sklearn.mixture.GaussianMixture subclass that jointly clusters peptides on their phosphorylation signal and their sequence motif — and get_pspl_pssm_distances, the helper it uses to compare cluster motifs against kinase specificity profiles.

DDMC

DDMC(
    n_components: int,
    seq_weight: float,
    distance_method: Literal[
        "PAM250", "Binomial"
    ] = "Binomial",
    random_state: int | RandomState | None = None,
    max_iter: int = 200,
    tol: float = 0.0001,
)

Bases: GaussianMixture

Cluster peptides by both sequence similarity and condition-wise phosphorylation following an expectation-maximization algorithm.

DDMC subclasses sklearn.mixture.GaussianMixture and reuses its EM loop, but scores each peptide against each cluster using both the usual Gaussian mixture log-probability over its phosphorylation signal and a sequence-motif term (weighted by seq_weight), and refits both the Gaussian mixture parameters and the per-cluster sequence motif at every M step. See ddmc.binomial.Binomial and ddmc.pam250.PAM250 for the two available motif models.

Attributes set by fit: p_signal: The p_signal DataFrame passed to fit. sequences: p_signal.index, as an upper-cased numpy array. seq_dist: The fitted Binomial or PAM250 sequence-distance model. scores_: Per-peptide, per-cluster responsibilities (soft cluster assignments) of shape (n_peptides, n_components). seq_scores_: Per-peptide, per-cluster weighted sequence log-probabilities (seq_weight * seq_dist.logWeights) from the last E step, of shape (n_peptides, n_components).

Parameters:

Name	Type	Description	Default
n_components	int	The number of clusters to fit.	required
seq_weight	float	Weight applied to the sequence-motif log-probability relative to the Gaussian mixture log-probability when scoring each peptide against each cluster. 0 reduces DDMC to an ordinary Gaussian mixture model.	required
distance_method	Literal['PAM250', 'Binomial']	Which sequence-distance model to use for the motif term: "Binomial" (ddmc.binomial.Binomial) or "PAM250" (ddmc.pam250.PAM250).	'Binomial'
random_state	int | RandomState | None	Seed or numpy.random.RandomState controlling the random initialization of the underlying Gaussian mixture, for reproducibility.	None
max_iter	int	Maximum number of EM iterations to run.	200
tol	float	Convergence threshold on the change in per-sample average log-likelihood between EM iterations.	0.0001

Source code in ddmc/clustering.py

def __init__(
    self,
    n_components: int,
    seq_weight: float,
    distance_method: Literal["PAM250", "Binomial"] = "Binomial",
    random_state: int | np.random.RandomState | None = None,
    max_iter: int = 200,
    tol: float = 1e-4,
):
    """
    Args:
        n_components: The number of clusters to fit.
        seq_weight: Weight applied to the sequence-motif log-probability
            relative to the Gaussian mixture log-probability when
            scoring each peptide against each cluster. `0` reduces
            `DDMC` to an ordinary Gaussian mixture model.
        distance_method: Which sequence-distance model to use for the
            motif term: `"Binomial"` (`ddmc.binomial.Binomial`) or
            `"PAM250"` (`ddmc.pam250.PAM250`).
        random_state: Seed or `numpy.random.RandomState` controlling the
            random initialization of the underlying Gaussian mixture,
            for reproducibility.
        max_iter: Maximum number of EM iterations to run.
        tol: Convergence threshold on the change in per-sample average
            log-likelihood between EM iterations.
    """
    super().__init__(
        n_components=n_components,
        covariance_type="diag",
        n_init=2,
        max_iter=max_iter,
        tol=tol,
        random_state=random_state,
    )
    self.distance_method = distance_method
    self.seq_weight = seq_weight

fit

fit(p_signal: DataFrame) -> DDMC

Compute EM clustering.

Parameters:

Name	Type	Description	Default
p_signal	DataFrame	Dataframe of shape (number of peptides, number of samples) containing the phosphorylation signal. p_signal.index contains the length-11 AA sequence of each peptide, containing the phosphoacceptor in the middle and five AAs flanking it.	required

Returns:

Type	Description
DDMC	self, fit to p_signal.

Source code in ddmc/clustering.py

def fit(self, p_signal: pd.DataFrame) -> "DDMC":  # ty: ignore[invalid-method-override]
    """
    Compute EM clustering.

    Args:
        p_signal: Dataframe of shape (number of peptides, number of samples)
            containing the phosphorylation signal. `p_signal.index` contains
            the length-11 AA sequence of each peptide, containing the
            phosphoacceptor in the middle and five AAs flanking it.

    Returns:
        self, fit to `p_signal`.
    """
    assert isinstance(p_signal, pd.DataFrame), (
        "`p_signal` must be a pandas dataframe."
    )
    sequences = p_signal.index.values

    for i, seq in enumerate(sequences):
        assert isinstance(seq, str), (
            f"Sequence {seq} at index {i} is not a string. All sequences must be strings."
        )
        assert len(seq) == 11, (
            f"Sequence {seq} at index {i} is of length {len(seq)}. All sequences must be of length 11."
        )
        assert all([token.upper() in AAlist for token in seq]), (
            f"Sequence {seq} at index {i} contains invalid characters."
        )

    assert (
        p_signal.select_dtypes(include=[np.number]).shape[1] == p_signal.shape[1]
    ), "All values in `p_signal` should be numerical"

    self.p_signal = p_signal
    self._gen_peptide_distances(sequences, self.distance_method)

    if np.any(np.isnan(p_signal)):
        self._missing = True
        self.missing_d = np.isnan(p_signal)

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            p_signal = SoftImpute(verbose=False).fit_transform(p_signal)
    else:
        self._missing = False

    super().fit(p_signal)
    self.scores_ = self.predict_proba(p_signal)

    assert np.all(np.isfinite(self.scores_))
    assert np.all(np.isfinite(self.seq_scores_))
    return self

get_nonempty_clusters

get_nonempty_clusters() -> np.ndarray

List the clusters that at least one peptide is assigned to.

Returns:

Type	Description
ndarray	Sorted array of the distinct cluster indices present in
ndarray	self.labels(); shorter than n_components if any clusters
ndarray	are empty.

Source code in ddmc/clustering.py

def get_nonempty_clusters(self) -> np.ndarray:
    """List the clusters that at least one peptide is assigned to.

    Returns:
        Sorted array of the distinct cluster indices present in
        `self.labels()`; shorter than `n_components` if any clusters
        are empty.
    """
    return np.unique(self.labels())

get_pssms

get_pssms(
    PsP_background: bool = False, clusters: None = None
) -> tuple[np.ndarray, np.ndarray]

get_pssms(
    PsP_background: bool = False, *, clusters: list[int]
) -> np.ndarray

get_pssms(
    PsP_background: bool = False,
    clusters: list[int] | None = None,
) -> tuple[np.ndarray, np.ndarray] | np.ndarray

Compute position-specific scoring matrix of each cluster. Note, to normalize by amino acid frequency this uses either all the sequences in the data set or a collection of random MS phosphosites in PhosphoSitePlus.

Parameters:

Name	Type	Description	Default
PsP_background	bool	Whether or not PhosphoSitePlus should be used for background frequency.	False
clusters	list[int] | None	cluster indices to get pssms for	None

Returns:

Type	Description
tuple[ndarray, ndarray] | ndarray	If the clusters argument is used, an array of shape (len(clusters), 20, 11),
tuple[ndarray, ndarray] | ndarray	else two arrays, where the first (of shape (n_pssms,))
tuple[ndarray, ndarray] | ndarray	contains the clusters of the pssms in the second
tuple[ndarray, ndarray] | ndarray	(of shape (n_pssms, 20, 11)).

Source code in ddmc/clustering.py

def get_pssms(
    self, PsP_background: bool = False, clusters: list[int] | None = None
) -> tuple[np.ndarray, np.ndarray] | np.ndarray:
    """
    Compute position-specific scoring matrix of each cluster.
    Note, to normalize by amino acid frequency this uses either
    all the sequences in the data set or a collection of random MS phosphosites in PhosphoSitePlus.

    Args:
        PsP_background: Whether or not PhosphoSitePlus should be used for background frequency.
        clusters: cluster indices to get pssms for

    Returns:
        If the clusters argument is used, an array of shape (len(clusters), 20, 11),
        else two arrays, where the first (of shape (n_pssms,))
        contains the clusters of the pssms in the second
        (of shape (n_pssms, 20, 11)).
    """
    pssm_names, pssms = [], []
    if PsP_background:
        bg_seqs = BackgroundSeqs(self.sequences)
        back_pssm = compute_control_pssm(bg_seqs)
    else:
        back_pssm = np.zeros((len(AAlist), 11), dtype=float)

    l1 = list(np.arange(self.n_components))
    l2 = list(set(self.labels()))
    ec = [i for i in l1 + l2 if i not in l1 or i not in l2]
    for ii in range(self.n_components):
        # Check for empty clusters and ignore them, if there are
        if ii in ec:
            continue

        # Compute PSSM
        pssm = np.zeros((len(AAlist), 11), dtype=float)
        for jj, seq in enumerate(self.sequences):
            seq = seq.upper()
            for kk, aa in enumerate(seq):
                pssm[AAlist.index(aa), kk] += self.scores_[jj, ii - 1]
                if ii == 1 and not PsP_background:
                    back_pssm[AAlist.index(aa), kk] += 1.0

        # Normalize by position across residues
        for pos in range(pssm.shape[1]):
            if pos == 5:
                continue
            pssm[:, pos] /= np.mean(pssm[:, pos])
            if ii == 1 and not PsP_background:
                back_pssm[:, pos] /= np.mean(back_pssm[:, pos])

        # Normalize to background PSSM to account for AA frequencies per position
        old_settings = np.seterr(divide="ignore", invalid="ignore")
        pssm /= back_pssm.copy()
        np.seterr(**old_settings)

        # Log2 transform
        pssm = np.ma.log2(pssm)
        pssm = pssm.filled(0)
        pssm = np.nan_to_num(pssm)
        pssm = pd.DataFrame(pssm)
        pssm.index = AAlist

        # Normalize phosphoacceptor position to frequency
        df = pd.DataFrame({"Sequence": self.sequences})
        df["Cluster"] = self.labels()
        clSeq = df[df["Cluster"] == ii]["Sequence"]
        clSeq = pd.DataFrame(frequencies(clSeq)).T
        tm = np.mean([clSeq.loc["S", 5], clSeq.loc["T", 5], clSeq.loc["Y", 5]])
        for p_site in ["S", "T", "Y"]:
            pssm.loc[p_site, 5] = np.log2(clSeq.loc[p_site, 5] / tm)

        pssms.append(np.clip(pssm, a_min=0, a_max=3))
        pssm_names.append(ii)

    pssm_names, pssms = np.array(pssm_names), np.array(pssms)

    if clusters is not None:
        return pssms[
            [np.where(pssm_names == cluster)[0][0] for cluster in clusters]
        ]

    return pssm_names, pssms

has_empty_clusters

has_empty_clusters() -> bool

Checks whether the most recent call to fit() resulted in empty clusters.

Returns:

Type	Description
bool	True if any of the n_components clusters has no peptides
bool	assigned to it.

Source code in ddmc/clustering.py

def has_empty_clusters(self) -> bool:
    """
    Checks whether the most recent call to fit() resulted in empty clusters.

    Returns:
        True if any of the `n_components` clusters has no peptides
        assigned to it.
    """
    check_is_fitted(self, ["scores_"])
    return self.get_nonempty_clusters().size != self.n_components

impute

impute() -> pd.DataFrame

Imputes missing values in the dataset passed in fit() and returns the imputed dataset.

Returns:

Type	Description
DataFrame	A copy of the p_signal passed to fit, with each peptide's
DataFrame	missing samples filled in from its assigned cluster's center.

Source code in ddmc/clustering.py

def impute(self) -> pd.DataFrame:
    """
    Imputes missing values in the dataset passed in fit() and returns the
    imputed dataset.

    Returns:
        A copy of the `p_signal` passed to `fit`, with each peptide's
        missing samples filled in from its assigned cluster's center.
    """
    p_signal = self.p_signal.copy()
    labels = self.labels()  # cluster assignments
    centers = self.transform()  # samples x clusters
    for ii in range(p_signal.shape[0]):
        p_signal.iloc[ii, np.isnan(p_signal.iloc[ii, :])] = centers[
            np.isnan(p_signal.iloc[ii, :]), labels[ii] - 1
        ]
    assert np.all(np.isfinite(p_signal))
    return p_signal

labels

labels() -> np.ndarray

Find cluster assignment with highest likelihood for each peptide.

Returns:

Type	Description
ndarray	Array of shape (n_peptides,) giving each peptide's cluster
ndarray	index. Equivalent to predict().

Source code in ddmc/clustering.py

def labels(self) -> np.ndarray:
    """Find cluster assignment with highest likelihood for each peptide.

    Returns:
        Array of shape (n_peptides,) giving each peptide's cluster
        index. Equivalent to `predict()`.
    """
    return self.predict()

predict

predict() -> np.ndarray

Provided the current model parameters, predict the cluster each peptide belongs to.

Returns:

Type	Description
ndarray	Array of shape (n_peptides,) giving the index of the
ndarray	highest-likelihood cluster for each peptide in self.p_signal.

Source code in ddmc/clustering.py

def predict(self) -> np.ndarray:  # ty: ignore[invalid-method-override]
    """Provided the current model parameters, predict the cluster each peptide belongs to.

    Returns:
        Array of shape (n_peptides,) giving the index of the
        highest-likelihood cluster for each peptide in `self.p_signal`.
    """
    check_is_fitted(self, ["scores_"])
    return np.argmax(self.scores_, axis=1)

predict_upstream_kinases

predict_upstream_kinases(
    PsP_background: bool = True,
) -> pd.DataFrame

Compute matrix-matrix similarity between kinase specificity profiles and cluster PSSMs to identify upstream kinases regulating clusters.

Parameters:

Name	Type	Description	Default
PsP_background	bool	Whether or not PhosphoSitePlus should be used for the background amino acid frequency when building each cluster's PSSM (see get_pssms).	True

Returns:

Type	Description
DataFrame	DataFrame of shape (n_kinases, n_nonempty_clusters) with a
DataFrame	Frobenius distance between each kinase's specificity profile and
DataFrame	each cluster's PSSM; smaller values indicate a better match.

Source code in ddmc/clustering.py

def predict_upstream_kinases(
    self,
    PsP_background: bool = True,
) -> pd.DataFrame:
    """Compute matrix-matrix similarity between kinase specificity profiles
    and cluster PSSMs to identify upstream kinases regulating clusters.

    Args:
        PsP_background: Whether or not PhosphoSitePlus should be used
            for the background amino acid frequency when building each
            cluster's PSSM (see `get_pssms`).

    Returns:
        DataFrame of shape (n_kinases, n_nonempty_clusters) with a
        Frobenius distance between each kinase's specificity profile and
        each cluster's PSSM; smaller values indicate a better match.
    """
    kinases, pspls = get_pspls()
    clusters, pssms = self.get_pssms(PsP_background=PsP_background)
    distances = get_pspl_pssm_distances(
        pspls,
        pssms,
        as_df=True,
        pssm_names=clusters,
        kinases=kinases,
    )
    return distances

score

score() -> float

Generate score of the fitting.

Returns:

Type	Description
float	The lower bound on the log-likelihood of the fitted model
float	(self.lower_bound_, set by GaussianMixture.fit).

Source code in ddmc/clustering.py

def score(self) -> float:  # ty: ignore[invalid-method-override]
    """Generate score of the fitting.

    Returns:
        The lower bound on the log-likelihood of the fitted model
        (`self.lower_bound_`, set by `GaussianMixture.fit`).
    """
    check_is_fitted(self, ["lower_bound_"])
    return self.lower_bound_

transform

transform(as_df: Literal[False] = False) -> np.ndarray

transform(as_df: Literal[True]) -> pd.DataFrame

transform(as_df: bool = False) -> np.ndarray | pd.DataFrame

Return cluster centers.

Parameters:

Name	Type	Description	Default
as_df	bool	Whether or not the result should be wrapped in a dataframe with labeled axes.	False

Returns:

Type	Description
ndarray | DataFrame	The cluster centers, either a np array or pd df of shape (n_samples, n_components).

Source code in ddmc/clustering.py

def transform(self, as_df: bool = False) -> np.ndarray | pd.DataFrame:
    """
    Return cluster centers.

    Args:
        as_df: Whether or not the result should be wrapped in a dataframe with labeled axes.

    Returns:
        The cluster centers, either a np array or pd df of shape (n_samples, n_components).
    """
    check_is_fitted(self, ["means_"])
    assert self.means_ is not None
    centers = self.means_.T
    if as_df:
        centers = pd.DataFrame(
            centers,
            index=self.p_signal.columns,
            columns=np.arange(self.n_components),
        )
    return centers

get_pspl_pssm_distances

get_pspl_pssm_distances(
    pspls: ndarray,
    pssms: ndarray,
    as_df: Literal[False] = False,
    pssm_names: Sequence | ndarray | None = None,
    kinases: Sequence | ndarray | None = None,
) -> np.ndarray

get_pspl_pssm_distances(
    pspls: ndarray,
    pssms: ndarray,
    as_df: Literal[True],
    pssm_names: Sequence | ndarray | None = None,
    kinases: Sequence | ndarray | None = None,
) -> pd.DataFrame

get_pspl_pssm_distances(
    pspls: ndarray,
    pssms: ndarray,
    as_df: bool = False,
    pssm_names: Sequence | ndarray | None = None,
    kinases: Sequence | ndarray | None = None,
) -> np.ndarray | pd.DataFrame

Computes a distance matrix between PSPLs and PSSMs.

Parameters:

Name	Type	Description	Default
pspls	ndarray	kinase specificity profiles of shape (n_kinase, 20, 9)	required
pssms	ndarray	position-specific scoring matrices of shape (n_pssms, 20, 11)	required
as_df	bool	Whether or not the returned matrix should be returned as a dataframe. Requires pssm_names and kinases.	False
pssm_names	Sequence | ndarray | None	list of names for the pssms of shape (n_pssms,)	None
kinases	Sequence | ndarray | None	list of names for the pspls of shape (n_kinase,)	None

Returns:

Type	Description
ndarray | DataFrame	Distance matrix of shape (n_kinase, n_pssms).

Source code in ddmc/clustering.py

def get_pspl_pssm_distances(
    pspls: np.ndarray,
    pssms: np.ndarray,
    as_df: bool = False,
    pssm_names: Sequence | np.ndarray | None = None,
    kinases: Sequence | np.ndarray | None = None,
) -> np.ndarray | pd.DataFrame:
    """
    Computes a distance matrix between PSPLs and PSSMs.

    Args:
        pspls: kinase specificity profiles of shape (n_kinase, 20, 9)
        pssms: position-specific scoring matrices of shape (n_pssms, 20, 11)
        as_df: Whether or not the returned matrix should be returned as a
            dataframe. Requires pssm_names and kinases.
        pssm_names: list of names for the pssms of shape (n_pssms,)
        kinases: list of names for the pspls of shape (n_kinase,)

    Returns:
        Distance matrix of shape (n_kinase, n_pssms).
    """
    assert pssms.shape[1:3] == (20, 11)
    assert pspls.shape[1:3] == (20, 9)
    pssms = np.delete(pssms, [5, 10], axis=2)
    dists = np.linalg.norm(pspls[:, None, :, :] - pssms[None, :, :, :], axis=(2, 3))
    if as_df:
        kinases = list(kinases) if kinases is not None else None
        pssm_names = list(pssm_names) if pssm_names is not None else None
        dists = pd.DataFrame(dists, index=kinases, columns=pssm_names)
    return dists