Lab 2 — Hidden populations: mixtures on real single-cell data

Author

Hicham Janati — FSML

This lab works on real single cell gene expression data (Muraro et al., 2016). Each row is a cell from a human pancreas, each column a gene: the number measures how strongly that gene is expressed in that cell (already normalized and log-transformed). Different pancreatic cell populations express different genes, but for now the cells’ biological identities are hidden. Your job: find the populations with a Gaussian mixture, then work out what they are.

Boxes as always: 🔮 predict before running (no AI) · 🧠 explain in your own words · ✍️ short task, boilerplate given · 🤖 a prompt that sends you to find something the lab did not hand you.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

BLUE, RED, GREEN, ORANGE, GRAY = "#1f6fb2", "#c92a1a", "#1c7a43", "#e8871a", "#8a8a8a"
plt.rcParams["figure.figsize"] = (6.5, 4.2)
plt.rcParams["axes.spines.top"] = plt.rcParams["axes.spines.right"] = False

def scatter_groups(xy, groups, title="", ax=None):
    """Scatter of 2D points colored by a group label (any dtype)."""
    if ax is None:
        _, ax = plt.subplots()
    for g in np.unique(groups):
        m = groups == g
        ax.scatter(xy[m, 0], xy[m, 1], s=6, label=str(g), alpha=0.7)
    ax.set_title(title); ax.legend(fontsize=8, markerscale=2)
    return ax

Part 1 — Intro to pandas

# the data live on the course website: no file to download
DATA = "https://hichamjanati.github.io/media/teaching/SMV/"
df = pd.read_csv(DATA + "muraro_endocrine_student.csv")
print(df.shape)
df.head()

One row per cell, identified by cell_id; every other column is a gene, named by its symbol. So df.shape reads: cells × (1 + genes).

genes = df.columns[1:]            # everything except cell_id
print(len(genes), "genes; first five:", list(genes[:5]))
print("INS in the file?", "INS" in genes)

Selecting one column gives a Series:

df["INS"].head()

Selecting several columns takes a Python list inside the brackets and gives a DataFrame:

df[["INS", "GCG"]].head()

✍️ Try it. Select the four hormone genes INS, GCG, SST, PPY at once and call .describe() on the result.

hormones = ["INS", "GCG", "SST", "PPY"]
# TODO: describe the four hormone columns at once

Filtering rows is a boolean condition inside the brackets, here the cells whose INS expression exceeds 2:

df[df["INS"] > 2].shape

✍️ Try it. How many cells have GCG above 6?

# TODO: count the cells with GCG above 6

The opposite of a condition is ~condition:

df[~(df["INS"] > 2)].shape

✍️ Try it. Check that the two counts (INS above 2, and its opposite) add up to the number of cells.

# TODO: the two counts, added up

Rows and columns together use .loc[condition, columns]:

df.loc[df["INS"] > 2, ["INS", "GCG"]].head()

🔮 Predict. INS is the insulin gene. If every cell in the file were the same kind of cell, what would a histogram of INS look like, and what would change if there were several kinds?

✍️ Your answer:

…

plt.hist(df["INS"], bins=30, color=GRAY)
plt.xlabel("INS expression (log scale)"); plt.ylabel("number of cells")
plt.show()

Two bumps, with a valley near 2. Three short tasks, one per cell:

✍️ Your turn.
  • the mean expression of INS across all cells;
  • the number of cells with INS above 2, the right-hand bump;
  • among those cells, the mean of GCG; then the same mean among the other cells. 
# TODO: mean INS
# TODO: how many cells have INS above 2
# TODO: mean GCG among INS-high cells, then among the others

🧠 Explain. Cells rich in INS have clearly lower GCG than the rest. Before any model: what does that already suggest about the population of cells?

✍️ Your answer:

…

One last idiom: models want a plain numpy array, df[genes].values (shape cells × genes), with no cell_id and no column names.

X = df[genes].values
print(X.shape, X.dtype)

Part 2 — Look at the geometry before modeling

500 genes = 500 dimensions. Before asking a mixture model to find populations in there, ask a cheaper question: does the cloud even look structured?

2.1 The cloud in two dimensions

✍️ Your turn. Fit PCA(n_components=2, random_state=0) on X, store the projected data in Z_2d, and print the explained-variance ratio of the two components.

from sklearn.decomposition import PCA

# TODO: pca_2d = ...
# TODO: Z_2d = ...
# TODO: print the explained-variance ratio

🔮 Predict. If distinct cell populations exist, must they be visible in a PC1-vs-PC2 scatter? Which populations could hide, and where?

✍️ Your answer:

…

plt.scatter(Z_2d[:, 0], Z_2d[:, 1], s=6, color=GRAY, alpha=0.6)
plt.xlabel("PC1"); plt.ylabel("PC2"); plt.title("cells in the first two PCs")
plt.show()

🧠 Explain. What structure do you see in the plot, and what do you expect it to mean for the modeling decisions ahead?

✍️ Your answer:

…

2.2 How many dimensions to model?

🧠 Explain. Why reduce the dimension at all before fitting a Gaussian mixture? Count: how many parameters does a single Gaussian have on all 500 genes, and how many cells do we have to estimate them?

✍️ Your answer:

…

So we reduce first: we compute 100 components and then decide how many to keep, using the cumulative explained variance.

✍️ Your turn. Fit PCA(n_components=100, random_state=0) on X, store the projection in Z_all, and plot the cumulative explained-variance ratio against the number of components.

# TODO: pca = ...
# TODO: Z_all = ...
# TODO: plot the cumulative explained variance against the number of components
plt.xlabel("components"); plt.ylabel("cumulative explained variance"); plt.show()

✍️ Your turn. Choose your pca_dim. There is no single right number; pick one, and write one sentence justifying it from the curve.

# TODO: pca_dim to be defined as your choice, an integer

✍️ Your answer:

…

# TODO: pca = ...
# TODO: Z_all = ...

Z = Z_all[:, :pca_dim]                          # the representation we will model
print(f"modeling on {pca_dim} PCs, keeping {pca.explained_variance_ratio_[:pca_dim].sum():.0%} of the variance")

Part 3 — Choosing the mixture

We model Z with a Gaussian mixture whose components have full covariances: any ellipsoid, a different one per component. One parameter is left to choose, the number of components K.

🔮 Predict. How many parameters does a K-component mixture with full covariances have in dimension d = pca_dim? Write the formula.

✍️ Your answer:

…

✍️ Your turn. Evaluate it for K = 2 … 6 with your pca_dim, and print the number of cells per parameter.

d, n_cells = pca_dim, len(df)
for K in [2, 3, 4, 5, 6]:
    # TODO: npar = ...
    # TODO: print K, npar and the cells per parameter
    pass

🧠 Explain. Does the cells-per-parameter figure make you nervous about K = 6? About the pca_dim you chose?

✍️ Your answer:

…

✍️ Your turn. For each K in Ks, fit GaussianMixture(K, covariance_type=“full”, n_init=5, random_state=0) on Z and store its bic(Z) in the dictionary.

from sklearn.mixture import GaussianMixture

Ks = [2, 3, 4, 5, 6, 7, 8]
bic = {}
for K in Ks:
    # TODO: bic[K] = the BIC on Z of a K-component full-covariance mixture
    pass

for K in Ks:
    print(f"K = {K}: BIC = {bic[K]:.0f}")
plt.plot(Ks, [bic[K] for K in Ks], "o-", color=BLUE)
plt.xlabel("K"); plt.ylabel("BIC (lower = better)"); plt.show()

🧠 Explain. What is the role of n_init=5 in the code above ?

✍️ Your answer:

…

🧠 Explain. What K does the BIC prefer? Compare with one of your friends who chose a different pca_dim: do you agree on the best K?

✍️ Your answer:

…

BIC weighs fit against parameter count. Look at the models it likes, not just the number. First on the PCs the models were fitted on:

✍️ Your turn. List the three K values you want to look at.

# TODO: selected_Ks to be defined as a list of three integers
def fit_gmm(K, n_init=5, random_state=0):
    return GaussianMixture(K, covariance_type="full", n_init=n_init, random_state=random_state).fit(Z)

labels_of = {K: fit_gmm(K).predict(Z) for K in selected_Ks}

fig, axes = plt.subplots(1, 3, figsize=(13, 3.8))
for ax, K in zip(axes, selected_Ks):
    scatter_groups(Z[:, :2], labels_of[K], f"K = {K}: sizes {np.bincount(labels_of[K]).tolist()}", ax)
    ax.set_xlabel("PC1"); ax.set_ylabel("PC2")
plt.tight_layout(); plt.show()

Then on a t-SNE map:

from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, perplexity=30, init="pca", random_state=0).fit_transform(Z)

fig, axes = plt.subplots(1, 3, figsize=(13, 3.8))
for ax, K in zip(axes, selected_Ks):
    scatter_groups(tsne, labels_of[K], f"K = {K} on the t-SNE map", ax)
    ax.set_xticks([]); ax.set_yticks([])
plt.tight_layout(); plt.show()

🧠 Explain. What do these visualizations tell you about the data and about the three models?

✍️ Your answer:

…

🧠 Explain. Choose one model to interpret in Part 4 and justify your choice.

✍️ Your answer:

…

K_chosen = selected_Ks[0]      # change if you argued for another one above
gm = fit_gmm(K_chosen)
clusters = gm.predict(Z)
print("K =", K_chosen, "| cluster sizes:", np.bincount(clusters).tolist())
print("weights:", gm.weights_.round(3), "| converged:", gm.converged_)

Part 4 — What are these clusters? Ask the genes

Back to the original gene columns. A gene characterizes a cluster if it is expressed much more inside than outside:

\[ \Delta_{kj} = \text{mean of gene } j \text{ in cluster } k \;-\; \text{mean of gene } j \text{ outside cluster } k . \]

✍️ Your turn. For each cluster, compute Δ for every gene: the mean over the rows inside the cluster minus the mean over the rows outside it (see the code in Part 1). Then store the three genes with the largest Δ in top[k].

top = {}
for k in range(K_chosen):
    inside = clusters == k        # boolean mask over cells: True where the cell was assigned to cluster k
    # TODO: the Series deltas to be defined as "mean inside minus mean outside", one value per gene
    top[k] = deltas.nlargest(3).to_dict()

for k, s in top.items():
    print(f"cluster {k} (n = {(clusters == k).sum()}):", ", ".join(f"{g} (+{v:.1f})" for g, v in s.items()))

🧠 Explain. These are the genes most characteristic of each inferred population. Most clusters have a top gene with a large Δ (above 3); if one cluster has nothing above 1, what might that tell you before you know any biology?

✍️ Your answer:

…

🤖 With AI. You now have gene names and no biology. Use an assistant to learn some, one cluster at a time: “In human pancreatic single-cell data, what cell type or biological role is associated with the genes X, Y, Z? Explain briefly; do not assume a cluster label.” Fill the table below: top-3 genes · what you learn the genes do · your proposed cell type · a confidence (high / medium / low) and why. (Naming: “PP cells” and “gamma cells” are the same population.)

cluster top-3 genes what the genes do your cell type confidence
0
1
2
3
…

✍️ Your answer:

…

Part 5 — The reveal

The cells were annotated by the original authors. Load the annotations now, and only now, and cross them with your clusters:

labels = pd.read_csv(DATA + "muraro_endocrine_labels.csv")
merged = df[["cell_id"]].assign(cluster=clusters).merge(labels, on="cell_id")
print(merged["muraro_celltype"].value_counts().to_dict())
pd.crosstab(merged["cluster"], merged["muraro_celltype"])

🧠 Explain, label switching. Muraro calls one population “beta”; the mixture calls it “cluster 1”. Did the GMM make a mistake, and what happens to the mixture model if you permute the component numbers?

✍️ Your answer:

…

✍️ Your turn. Read the correspondence off the table, one line per cluster (cluster → cell type), and mark any cluster that has no clean match.

✍️ Your answer:

…

Comparing with the experts, despite label switching. The adjusted Rand index (ARI) is a score that can be used to compare two different clusterings (of two models for e.g A and B) even if label names do not match. It looks at every pair of cells and asks whether the labelings of model A agree on “same group / not same group”, so it never needs cluster numbers to match. ARI = 1 for identical partitions and ARI approximately 0 for unrelated ones.

✍️ Your turn. Compute the ARI between your clusters and the Muraro annotations with adjusted_rand_score (from sklearn.metrics); then the same for the other two models in labels_of, a small scores table.

# TODO: import the score

for K in selected_Ks:
    # TODO: print the ARI of labels_of[K] against merged["muraro_celltype"]
    pass

🧠 Explain. The scores are high but not 1. Looking at the table, which cells are costing the points ?

✍️ Your answer:

…