Clear code smell flagged by vulture, jscpd, and radon/xenon

analysis/figures/figure2f_h.py

@@ -35,7 +35,7 @@ def makeFigure():
     return f
 
 
-def groupDrugs(X: anndata, label_name: str):
+def groupDrugs(X: anndata.AnnData, label_name: str):
     """Groups drugs of similar category"""
     labels = X.obs[label_name]
     names = np.unique(labels)

analysis/figures/figure3.py

@@ -104,20 +104,6 @@ def makeFigure():
     return f
 
 
-def plot_cell_perc_corr(cellDF: pd.DataFrame, pop1: str, pop2: str, ax: Axes):
-    """Plots correlation of cell percentages against each other"""
-    newDF = pd.DataFrame()
-    newDF2 = pd.DataFrame()
-    newDF[[pop1, "Condition"]] = cellDF.loc[cellDF["Cell Type"] == pop1][
-        ["Cell Type Percentage", "Condition"]
-    ]
-    newDF2[[pop2, "Condition"]] = cellDF.loc[cellDF["Cell Type"] == pop2][
-        ["Cell Type Percentage", "Condition"]
-    ]
-    newDF = newDF.merge(newDF2, on="Condition")
-    sns.scatterplot(newDF, x=pop1, y=pop2, hue="Condition", ax=ax)
-
-
 def plot_cell_perc_comp_corr(
     X: anndata.AnnData, cellDF: pd.DataFrame, pop: str, comp: int, ax: Axes, unique=None
 ):

analysis/figures/figure4e_k.py

@@ -53,7 +53,7 @@ def makeFigure():
 
 
 def plot_correlation_cmp_cell_count_perc(
-    X: anndata, cmp: int, cellcountDF: pd.DataFrame, ax: Axes, cellPerc=True
+    X: anndata.AnnData, cmp: int, cellcountDF: pd.DataFrame, ax: Axes, cellPerc=True
 ):
     """Plot component weights by cell type count or percentage for a cell type"""
     yt = np.unique(X.obs["Condition"])

analysis/figures/figureS9e_f.py

@@ -48,8 +48,8 @@ def makeFigure():
 
 def logreg_weights_scores(
     X: anndata.AnnData, y: pd.Series, prediction: str
-) -> pd.DataFrame:
-    """Fit logistic regression model, return coefficients of that model"""
+) -> tuple[pd.DataFrame, float]:
+    """Fit logistic regression, returning its coefficients and accuracy score."""
     status = y[prediction]
     cond_factors = np.array(X.uns["Pf2_A"])
     lr = logistic_regression("accuracy").fit(cond_factors, status)

analysis/gating.py

@@ -5,7 +5,11 @@
 
 
 def gateThomsonCellsLeiden(X) -> npt.ArrayLike:
-    """Manually gates cell types for Thomson PaCMAP"""
+    """Manually gates cell types for Thomson PaCMAP.
+
+    Run by hand to (re)build ThomsonCellTypes.parquet, which
+    :func:`gateThomsonCells` then reads; deliberately has no in-repo caller.
+    """
     sc.pp.neighbors(X, n_neighbors=15, use_rep="projections", random_state=0)
     sc.tl.leiden(X, resolution=3, random_state=0)
     X.obs["Cell Type"] = X.obs.leiden.replace(thomson_layer1).astype(str)
@@ -26,7 +30,11 @@ def gateThomsonCells(X) -> npt.ArrayLike:
 
 
 def Thomson_Doublet():
-    """Detects doublets in scRNA-seq"""
+    """Detects doublets in scRNA-seq.
+
+    Run by hand to (re)build ThomsonDoublets.csv; deliberately has no in-repo
+    caller.
+    """
     import hdf5plugin  # noqa: F401
     import vsparse
 
@@ -49,38 +57,6 @@ def Thomson_Doublet():
     X.obs["doublet"].to_csv("analysis/data/Thomson/ThomsonDoublets.csv")
 
 
-def getHiResOldLupus(X) -> npt.ArrayLike:
-    """Manually gates cell types for SLE PaCMAP"""
-    X.obs["Cell Type Old2"] = X.obs["Cell Type Old"].astype(str)
-    X.obs.cell_type_lympho = X.obs.cell_type_lympho.astype(str)
-    X.obs.loc[X.obs["cell_type_lympho"] != "nan", "Cell Type Old2"] = X.obs.loc[
-        X.obs["cell_type_lympho"] != "nan"
-    ].cell_type_lympho.values
-    X.obs = X.obs.replace({"Cell Type Old2": cell_type_conv})
-
-    return X
-
-
-cell_type_conv = {
-    "T4_naive": "T4 Naive",
-    "B_naive": "B Naive",
-    "CytoT_GZMH+": "T8 GZMH",
-    "T4_em": "T4 EM",
-    "NK_dim": "NK Dim",
-    "CytoT_GZMK+": "T8 GZMK",
-    "T8_naive": "T8 Naive",
-    "ncM": "nCM",
-    "B_mem": "B mem",
-    "T4_reg": "T4 Reg",
-    "cDC": "cDC2",
-    "T_mait": "T mait",
-    "B_plasma": "B plasma",
-    "NK_bright": "NK Bright",
-    "B_atypical": "B Atypical",
-    "cM": "CM",
-}
-
-
 thomson_layer1 = {
     "0": "NK Cells",
     "1": "T Cells",
@@ -179,7 +155,8 @@ def getHiResOldLupus(X) -> npt.ArrayLike:
 }
 
 
-# Taken directly from PopAlign and used to annotate level 1
+# Taken directly from PopAlign and used to annotate level 1. Reference data
+# consulted when curating the layer maps above by hand -- no code reads it.
 marker_genes_1 = {
     "Monocytes": [
         "CD14",
@@ -211,7 +188,8 @@ def getHiResOldLupus(X) -> npt.ArrayLike:
 }
 
 
-# Take from various resources to annotate level 2
+# Taken from various resources to annotate level 2. Reference data consulted
+# when curating the layer maps above by hand -- no code reads it.
 marker_genes_2 = {
     "B cells": ["PXK", "MS4A1", "CD19", "CD74", "CD79A", "BANK1", "PTPRC", "CR2"],
     "B Memory": ["NPIB15", "BACH2", "IL7", "NMBR", "MS4A1", "MBL2", "LY86", "CD27"],

scrise/annotation_alignment.py

@@ -294,6 +294,40 @@ def cell_type_alignment(
     )
 
 
+def _component_metrics(
+    loadings_matrix: np.ndarray,
+    codes: np.ndarray,
+    n_types: int,
+    signed: bool,
+    n_permutations: int,
+    rng: np.random.Generator,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    """Per-component AUROCs, permutation p-values, tau, eta^2 and epsilon^2."""
+    n_cells, n_comps = loadings_matrix.shape
+    aurocs_mat = np.zeros((n_comps, n_types), dtype=float)
+    p_vals_mat = np.ones((n_comps, n_types), dtype=float)
+    tau_vec = np.zeros(n_comps, dtype=float)
+    eta2_vec = np.zeros(n_comps, dtype=float)
+    eps2_vec = np.zeros(n_comps, dtype=float)
+
+    for comp_idx in range(n_comps):
+        y = loadings_matrix[:, comp_idx].astype(float)
+        if signed:
+            y = np.abs(y)
+
+        auc = compute_auroc_per_cell_type(y, codes, n_types)
+        aurocs_mat[comp_idx] = auc
+        tau_vec[comp_idx] = compute_tau(auc, baseline=0.0)
+        eta2_vec[comp_idx] = compute_eta_squared(y, codes, n_types)
+        eps2_vec[comp_idx] = compute_kruskal_epsilon_squared(y, codes, n_types)
+
+        p_vals_mat[comp_idx] = _component_p_values(
+            y, codes, n_types, n_cells, auc, n_permutations, rng
+        )
+
+    return aurocs_mat, p_vals_mat, tau_vec, eta2_vec, eps2_vec
+
+
 def score_cell_type_alignment(
     data: anndata.AnnData | np.ndarray | pd.DataFrame,
     cell_types: pd.Series | np.ndarray | str | None = None,
@@ -353,26 +387,9 @@ def score_cell_type_alignment(
     rng = _as_generator(random_state)
 
     component_labels = [i + 1 for i in range(n_comps)]
-    aurocs_mat = np.zeros((n_comps, n_types), dtype=float)
-    p_vals_mat = np.ones((n_comps, n_types), dtype=float)
-    tau_vec = np.zeros(n_comps, dtype=float)
-    eta2_vec = np.zeros(n_comps, dtype=float)
-    eps2_vec = np.zeros(n_comps, dtype=float)
-
-    for comp_idx in range(n_comps):
-        y = loadings_matrix[:, comp_idx].astype(float)
-        if signed:
-            y = np.abs(y)
-
-        auc = compute_auroc_per_cell_type(y, codes, n_types)
-        aurocs_mat[comp_idx] = auc
-        tau_vec[comp_idx] = compute_tau(auc, baseline=0.0)
-        eta2_vec[comp_idx] = compute_eta_squared(y, codes, n_types)
-        eps2_vec[comp_idx] = compute_kruskal_epsilon_squared(y, codes, n_types)
-
-        p_vals_mat[comp_idx] = _component_p_values(
-            y, codes, n_types, n_cells, auc, n_permutations, rng
-        )
+    aurocs_mat, p_vals_mat, tau_vec, eta2_vec, eps2_vec = _component_metrics(
+        loadings_matrix, codes, n_types, signed, n_permutations, rng
+    )
 
     # Joint BH FDR correction across all (component x cell_type) tests
     if aurocs_mat.size > 1:

scrise/factor_io.py

@@ -111,39 +111,45 @@ def _read_raw_dataset(raw_path: str) -> anndata.AnnData:
         return anndata.read_h5ad(raw_path)
 
 
-def _attach_raw_data(adata: anndata.AnnData, raw_path: str) -> None:
-    """Match the factors' cells and genes against raw data and attach ``X``."""
-    raw = _read_raw_dataset(raw_path)
-
+def _subset_raw_cells(adata: anndata.AnnData, raw: anndata.AnnData) -> anndata.AnnData:
+    """Align the raw cells to the factors' barcodes, retrying via ``cell_barcode``."""
     if not isinstance(raw.obs, pd.DataFrame):
         raise TypeError("raw.obs must be an in-memory pandas DataFrame.")
 
-    # Match cells by index or cell_barcode column
     if not np.all(adata.obs_names.isin(raw.obs_names)) and "cell_barcode" in raw.obs:
         raw.obs.index = pd.Index(raw.obs["cell_barcode"].astype(str))
 
-    common_cells = adata.obs_names[adata.obs_names.isin(raw.obs_names)]
-    if len(common_cells) == 0:
+    if not adata.obs_names.isin(raw.obs_names).any():
         raise ValueError(
             "No matching cell barcodes found between factors and raw data."
         )
-    raw_sub = raw[adata.obs_names, :].copy()
+    return raw[adata.obs_names, :].copy()
 
+
+def _subset_raw_genes(
+    adata: anndata.AnnData, raw_sub: anndata.AnnData
+) -> anndata.AnnData:
+    """Align the raw genes to the factors' names, retrying via ``gene_ids``."""
     if not isinstance(raw_sub.var, pd.DataFrame):
         raise TypeError("raw_sub.var must be an in-memory pandas DataFrame.")
 
-    # Match genes
     if (
         not np.all(adata.var_names.isin(raw_sub.var_names))
         and "gene_ids" in raw_sub.var
         and "gene_ids" in adata.var
     ):
         raw_sub.var.index = pd.Index(raw_sub.var["gene_ids"].astype(str))
 
-    common_genes = adata.var_names[adata.var_names.isin(raw_sub.var_names)]
-    if len(common_genes) == 0:
+    if not adata.var_names.isin(raw_sub.var_names).any():
         raise ValueError("No matching gene names found between factors and raw data.")
-    raw_sub = raw_sub[:, adata.var_names].copy()
+    return raw_sub[:, adata.var_names].copy()
+
+
+def _attach_raw_data(adata: anndata.AnnData, raw_path: str) -> None:
+    """Match the factors' cells and genes against raw data and attach ``X``."""
+    raw_sub = _subset_raw_genes(
+        adata, _subset_raw_cells(adata, _read_raw_dataset(raw_path))
+    )
 
     from parafac2.normalize import prepare_dataset
 

scrise/factorization.py

@@ -269,6 +269,12 @@ def pf2(
         Additional keyword arguments forwarded to ``parafac2_nd`` (e.g.
         ``n_inner``, ``callback``), for options not otherwise exposed here.
 
+    condition_key : str, optional (default: None)
+        Column in ``X.obs`` holding the condition labels, used to derive
+        ``condition_unique_idxs`` when that column is not already present.
+    adata : anndata.AnnData, optional (default: None)
+        Alias for ``X``; supply either one, not both.
+
     Returns
     -------
     anndata.AnnData

scrise/plotting/factors.py

@@ -111,6 +111,12 @@ def plot_condition_factors(
         Custom colors for condition group labels. If None, uses default palette.
     group_cond : bool, optional (default: False)
         If True and cond_group_labels provided, sorts conditions by group.
+    control_pattern : str, optional (default: None)
+        Substring identifying the control conditions whose median and spread
+        set the normalization. Overrides ``ThomsonNorm``'s implied 'CTRL'.
+    control_conditions : Sequence[str], optional (default: None)
+        Explicit control condition names, taking precedence over
+        ``control_pattern``. With neither given, every condition contributes.
     """
     yt = pd.Series(np.unique(data.obs[cond]))
     X = np.array(data.uns["Pf2_A"])
@@ -152,6 +158,22 @@ def plot_condition_factors(
     ax.set(xlabel="Component")
 
 
+def _draw_unit_scaled_heatmap(X: np.ndarray, yt, ax: Axes) -> None:
+    """Draw a factor heatmap rescaled so the largest magnitude sits at +/-1."""
+    X = X / np.max(np.abs(np.asarray(X)))
+    sns.heatmap(
+        data=X,
+        xticklabels=np.arange(1, X.shape[1] + 1),
+        yticklabels=yt,
+        ax=ax,
+        center=0,
+        cmap=cmap,
+        vmin=-1,
+        vmax=1,
+    )
+    ax.set(xlabel="Component")
+
+
 def plot_eigenstate_factors(data: anndata.AnnData, ax: Axes):
     """Plot eigen-state factors as a heatmap showing cell state patterns.
 
@@ -168,23 +190,8 @@ def plot_eigenstate_factors(data: anndata.AnnData, ax: Axes):
     ax : matplotlib.axes.Axes
         Matplotlib axes object to plot on.
     """
-    rank = data.uns["Pf2_B"].shape[1]
-    xticks = np.arange(1, rank + 1)
-    X = data.uns["Pf2_B"]
-    X = X / np.max(np.abs(np.array(X)))
-    yt = np.arange(1, rank + 1)
-
-    sns.heatmap(
-        data=X,
-        xticklabels=xticks,
-        yticklabels=yt,
-        ax=ax,
-        center=0,
-        cmap=cmap,
-        vmin=-1,
-        vmax=1,
-    )
-    ax.set(xlabel="Component")
+    X = np.asarray(data.uns["Pf2_B"])
+    _draw_unit_scaled_heatmap(X, np.arange(1, X.shape[1] + 1), ax)
 
 
 def plot_gene_factors(data: anndata.AnnData, ax: Axes, weight=0.08, trim=True):
@@ -209,7 +216,6 @@ def plot_gene_factors(data: anndata.AnnData, ax: Axes, weight=0.08, trim=True):
     trim : bool, optional (default: True)
         If True, filters genes based on the weight parameter. If False, shows all genes.
     """
-    rank = data.varm["Pf2_C"].shape[1]
     X = np.array(data.varm["Pf2_C"])
     yt = data.var.index.values
 
@@ -224,22 +230,7 @@ def plot_gene_factors(data: anndata.AnnData, ax: Axes, weight=0.08, trim=True):
         yt = yt[kept_idxs]
 
     ind = reorder_table(X)
-    X = X[ind]
-    X = X / np.max(np.abs(X))
-    yt = yt[ind]
-    xticks = np.arange(1, rank + 1)
-
-    sns.heatmap(
-        data=X,
-        xticklabels=xticks,
-        yticklabels=yt,
-        ax=ax,
-        center=0,
-        cmap=cmap,
-        vmin=-1,
-        vmax=1,
-    )
-    ax.set(xlabel="Component")
+    _draw_unit_scaled_heatmap(X[ind], yt[ind], ax)
 
 
 def reorder_table(projs: np.ndarray) -> np.ndarray: