from typing import Optional
import numpy as np
import pandas as pd
import scipy.stats
import umap
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
[docs]
def get_histogram_series(s: pd.Series, bins: np.ndarray) -> pd.Series:
hist_data, bin_edges = np.histogram(s, bins=bins)
ret = pd.Series(
hist_data,
index=pd.MultiIndex.from_arrays(
[bin_edges[:-1], bin_edges[1:], range(len(hist_data))],
names=["bin_number", "bin_start", "bin_end"],
),
)
return ret
[docs]
def calculate_coefficient_variation(df: pd.DataFrame) -> pd.Series:
"""
Compute the coefficient of variation (CV) for each column in a DataFrame.
The coefficient of variation is defined as the ratio of the standard deviation
to the mean, expressed as a percentage. This function uses the biased standard
deviation (normalization by N) as implemented in `scipy.stats.variation`.
Parameters
----------
df : pandas.DataFrame
Input DataFrame containing numeric values (e.g., log2-transformed data).
Each column will be processed independently.
Returns
-------
pandas.Series
Series containing the coefficient of variation (in percent) for each column.
The index corresponds to the columns of the input DataFrame.
See Also
--------
scipy.stats.variation : Function used to compute the coefficient of variation.
Examples
--------
>>> import pandas as pd
>>> import numpy as np
>>> df = pd.DataFrame({'A': [1, 2, 3], 'B': [4, 5, 6]})
>>> calculate_coefficient_variation(df)
A 40.825
B 16.330
Name: coef_of_var, dtype: float64
"""
cv = scipy.stats.variation(df, axis=0) * 100
cv = pd.Series(cv, index=df.columns).rename("coef_of_var")
return cv
[docs]
def calculate_coef_of_var_and_mean(
df: pd.DataFrame,
) -> pd.DataFrame:
"""Calculate coefficient of variation and mean for each column in the dataframe,
the mean calculated on both log2 and linear scale.
Parameters
----------
df : pd.DataFrame
The input dataframe containing the linear values (non-log transformed).
Returns
-------
pd.DataFrame
A dataframe with columns 'mean_log2', 'mean' and 'coef_of_var' for each column
in the input dataframe.
"""
cv = calculate_coefficient_variation(df)
means = df.mean().rename("mean")
means_logs = np.log2(df).mean().rename("mean_log2")
return pd.concat(
[
means_logs,
means,
cv,
],
axis=1,
).rename_axis("name")
[docs]
def get_coefficient_variation(
data: pd.DataFrame,
drop_columns: Optional[list[str]] = None,
group: str = "group",
):
"""
Extracts the coefficients of variation in each group.
:param data: pandas.DataFrame with samples as rows and protein identifiers as columns
(with additional columns 'group', 'sample' and 'subject'). The values
should be the original intensities for massspectrometry-based
measurements.
:param list drop_columns: column labels to be dropped from the dataframe
:param str group: column label containing group identifiers.
:return: Pandas dataframe with columns 'name' (protein identifier),
'x' (coefficient of variation), 'y' (mean) and 'group'.
Example::
result = get_coefficient_variation(data, drop_columns=['sample', 'subject'], group='group')
"""
formated_df = data
if drop_columns is not None:
formated_df = data.drop(drop_columns, axis=1)
cvs = formated_df.groupby(group).apply(func=calculate_coef_of_var_and_mean)
cvs = cvs.reset_index()[["name", "mean_log2", "mean", "coef_of_var", group]]
cvs_df = cvs
return cvs_df
[docs]
def run_pca(
data,
drop_cols=["sample", "subject"],
group="group",
annotation_cols=["sample"],
components=2,
dropna=True,
):
"""
Performs principal component analysis and returns the values of each component for each sample
and each protein, and the loadings for each protein.
For information visit
https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.
:param data: pandas dataframe with samples as rows and protein identifiers as columns
(with additional columns 'group', 'sample' and 'subject').
:param list drop_cols: column labels to be dropped from the dataframe.
:param str group: column label containing group identifiers.
:param list annotation_cols: list of columns to be added in the scatter plot annotation
:param int components: number of components to keep.
:param bool dropna: if True removes all columns with any missing values.
:return: tuple: 1) three pandas dataframes: components, loadings and variance; 2)
xaxis and yaxis titles with components loadings for plotly.
Example::
result = run_pca(data, drop_cols=['sample', 'subject'], group='group',
components=2, dropna=True)
"""
np.random.seed(112736)
var_exp = []
args = {}
if data.empty:
raise ValueError("Dataframe is empty.")
df = data.copy()
annotations = pd.DataFrame()
if annotation_cols is not None:
if len(list(set(annotation_cols).intersection(data.columns))) > 0:
annotations = data.set_index(group)[annotation_cols]
drop_cols_int = list(set(drop_cols).intersection(df.columns))
if len(drop_cols_int) > 0:
df = df.drop(drop_cols_int, axis=1)
y = df[group].tolist()
df = df.set_index(group)
df = df.select_dtypes(["number"])
if dropna:
df = df.dropna(axis=1)
X = df.values
if X.size > 0 and X.shape[1] > components:
pca = PCA(n_components=components)
X = pca.fit_transform(X)
var_exp = pca.explained_variance_ratio_
loadings = pd.DataFrame(pca.components_.transpose())
loadings.index = df.columns
# ? apply calculation before transposition on columns (for performance)
values = {
index: np.sqrt(np.power(row, 2).sum()) for index, row in loadings.iterrows()
}
loadings["value"] = loadings.index.map(values.get)
loadings = loadings.sort_values(by="value", ascending=False)
args = {
"x_title": "PC1" + " ({0:.2f})".format(var_exp[0]),
"y_title": "PC2" + " ({0:.2f})".format(var_exp[1]),
"group": "group",
}
if components == 2:
resultDf = pd.DataFrame(X, index=y, columns=["x", "y"])
resultDf = resultDf.assign(**annotations)
resultDf = resultDf.reset_index()
resultDf.columns = ["group", "x", "y"] + annotation_cols
loadings.columns = ["x", "y", "value"]
if components > 2:
args.update({"z_title": "PC3" + " ({0:.2f})".format(var_exp[2])})
resultDf = pd.DataFrame(X, index=y)
resultDf = resultDf.assign(**annotations)
resultDf = resultDf.reset_index()
pca_cols = []
loading_cols = []
if components > 3:
pca_cols = [str(i) for i in resultDf.columns[4:]]
loading_cols = [str(i) for i in loadings.columns[3:]]
resultDf.columns = ["group", "x", "y", "z"] + pca_cols
loadings.columns = ["x", "y", "z"] + loading_cols
return (resultDf, loadings, var_exp), args
[docs]
def run_tsne(
data,
drop_cols=["sample", "subject"],
group="group",
annotation_cols=["sample"],
components=2,
perplexity=40,
max_iter=1000,
init="pca",
dropna=True,
):
"""
Performs t-distributed Stochastic Neighbor Embedding analysis.
For more information visit
https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.html.
:param data: pandas dataframe with samples as rows and protein identifiers as columns
(with additional columns 'group', 'sample' and 'subject').
:param list drop_cols: column labels to be dropped from the dataframe.
:param str group: column label containing group identifiers.
:param int components: dimension of the embedded space.
:param list annotation_cols: list of columns to be added in the scatter plot annotation
:param int perplexity: related to the number of nearest neighbors that is used
in other manifold learning algorithms.
Consider selecting a value between 5 and 50.
:param int max_iter: maximum number of iterations for the optimization (at least 250).
:param str init: initialization of embedding ('random', 'pca' or
numpy array of shape n_samples x n_components).
:param bool dropna: if True removes all columns with any missing values.
:return: Two dictionaries:
1) pandas dataframe with embedding vectors,
2) xaxis and yaxis titles for plotly.
Example::
result = run_tsne(data,
drop_cols=['sample', 'subject'],
group='group',
components=2,
perplexity=40,
max_iter=1000,
init='pca',
dropna=True
)
"""
result = {}
args = {}
df = data.copy()
if len(set(drop_cols).intersection(df.columns)) == len(drop_cols):
df = df.drop(drop_cols, axis=1)
df = df.set_index(group)
if dropna:
df = df.dropna(axis=1)
df = df.select_dtypes(["number"])
X = df.values
y = df.index
annotations = pd.DataFrame()
if annotation_cols is not None:
if len(list(set(annotation_cols).intersection(data.columns))) > 0:
annotations = data[annotation_cols]
if X.size > 0:
tsne = TSNE(
n_components=components,
verbose=0,
perplexity=perplexity,
max_iter=max_iter,
init=init,
)
X = tsne.fit_transform(X)
args = {"x_title": "C1", "y_title": "C2"}
if components == 2:
resultDf = pd.DataFrame(X, index=y, columns=["x", "y"])
resultDf = resultDf.reset_index()
resultDf.columns = ["group", "x", "y"]
if components > 2:
args.update({"z_title": "C3"})
resultDf = pd.DataFrame(X, index=y)
resultDf = resultDf.reset_index()
cols = []
if len(components) > 4:
cols = resultDf.columns[4:]
resultDf.columns = ["group", "x", "y", "z"] + cols
resultDf = resultDf.join(annotations)
result["tsne"] = resultDf
return result, args
[docs]
def run_umap(
data,
drop_cols=["sample", "subject"],
group="group",
annotation_cols=["sample"],
n_neighbors=10,
min_dist=0.3,
metric="cosine",
dropna=True,
):
"""
Performs Uniform Manifold Approximation and Projection.
For more information vist https://umap-learn.readthedocs.io.
:param data: pandas dataframe with samples as rows and protein identifiers as columns
(with additional columns 'group', 'sample' and 'subject').
:param list drop_cols: column labels to be dropped from the dataframe.
:param str group: column label containing group identifiers.
:param list annotation_cols: list of columns to be added in the scatter plot annotation
:param int n_neighbors: number of neighboring points used
in local approximations of manifold structure.
:param float min_dist: controls how tightly the embedding is allowed compress points together.
:param str metric: metric used to measure distance in the input space.
:param bool dropna: if True removes all columns with any missing values.
:return: Two dictionaries:
1) pandas dataframe with embedding of the training data in low-dimensional space,
2) xaxis and yaxis titles for plotly.
Example::
result = run_umap(data,
drop_cols=['sample', 'subject'],
group='group',
n_neighbors=10,
min_dist=0.3,
metric='cosine',
dropna=True
)
"""
np.random.seed(1145536)
result = {}
args = {}
df = data.copy()
if len(set(drop_cols).intersection(df.columns)) == len(drop_cols):
df = df.drop(drop_cols, axis=1)
df = df.set_index(group)
if dropna:
df = df.dropna(axis=1)
df = df.select_dtypes(["number"])
X = df.values
y = df.index
annotations = pd.DataFrame()
if annotation_cols is not None:
if len(list(set(annotation_cols).intersection(data.columns))) > 0:
annotations = data[annotation_cols]
if not X.size:
return result, args
X = umap.UMAP(
n_neighbors=n_neighbors, min_dist=min_dist, metric=metric
).fit_transform(X)
args = {"x_title": "C1", "y_title": "C2"}
resultDf = pd.DataFrame(X, index=y)
resultDf = resultDf.reset_index()
cols = []
if len(resultDf.columns) > 3:
cols = resultDf.columns[3:]
resultDf.columns = ["group", "x", "y"] + cols
resultDf = resultDf.join(annotations)
result["umap"] = resultDf
return result, args
