feat: lightweight PCA backend (scipy/numpy SVD)

pydeeptools/deeptools/correlation.py

@@ -15,8 +15,7 @@
 import matplotlib.markers
 import matplotlib.colors as pltcolors
 from deeptools.utilities import toString, convertCmap
-from sklearn.decomposition import PCA
-from sklearn.preprocessing import StandardScaler
+from scipy.linalg import svd
 
 class Correlation:
     """
@@ -454,48 +453,80 @@ def plot_pca(self, plot_filename=None, PCs=[1, 2], plot_title='', image_format=N
 
         fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(plotWidth, plotHeight), layout="constrained")
 
-        # Filter
-        m = self.matrix
-        rvs = m.var(axis=1)
-        if self.transpose:
-            m = m[np.nonzero(rvs)[0], :]
-            rvs = rvs[np.nonzero(rvs)[0]]
-        if self.ntop > 0 and m.shape[0] > self.ntop:
-            m = m[np.argpartition(rvs, -self.ntop)[-self.ntop:], :]
-            rvs = rvs[np.argpartition(rvs, -self.ntop)[-self.ntop:]]
+        # Work on a float copy so the transforms below take effect and
+        # self.matrix (which callers may reuse) is left untouched.
+        m = self.matrix.astype(float, copy=True)
 
-        # log2 (if requested)
+        # log2 (if requested). Applied before variance filtering so the
+        # transform actually influences which rows are kept.
         if self.log2:
-            self.matrix = np.log2(self.matrix + 0.01)
+            m = np.log2(m + 0.01)
 
-        # Row center / transpose
+        # Row center (incompatible with --transpose, enforced by the CLI).
         if self.rowCenter and not self.transpose:
-            _ = self.matrix.mean(axis=1)
-            self.matrix -= _[:, None]
-        if self.transpose:
-            m = m.T
+            m -= m.mean(axis=1)[:, None]
 
-        # Center and scale
-        scaler = StandardScaler()
-        m2 = scaler.fit_transform(m)
+        # Filter to the most variable rows (variance computed post-transform).
+        rvs = m.var(axis=1)
+        if self.transpose:
+            nz = np.nonzero(rvs)[0]
+            m = m[nz, :]
+            rvs = rvs[nz]
+        if self.ntop > 0 and m.shape[0] > self.ntop:
+            idx = np.argpartition(rvs, -self.ntop)[-self.ntop:]
+            m = m[idx, :]
+            rvs = rvs[idx]
 
-        # PCA
-        pca = PCA()
-        Wt = pca.fit_transform(m2)
+        if self.transpose:
+            m = m.T
 
-        # % variance, eigenvalues
-        variance = pca.explained_variance_ratio_
+        # Center and scale each column to zero mean and unit variance
+        # (equivalent to sklearn's StandardScaler, using population std/ddof=0).
+        col_mean = m.mean(axis=0)
+        col_std = m.std(axis=0)
+        col_std[col_std == 0] = 1.0
+        m2 = (m - col_mean) / col_std
+
+        # PCA via SVD of the (re-)centered matrix, mirroring sklearn's PCA().
+        n_samples = m2.shape[0]
+        X = m2 - m2.mean(axis=0)
+        U, S, Vt = svd(X, full_matrices=False)
+
+        # Deterministic sign convention (sklearn's svd_flip, u_based_decision).
+        max_abs_cols = np.argmax(np.abs(U), axis=0)
+        signs = np.sign(U[max_abs_cols, range(U.shape[1])])
+        U *= signs
+        Vt *= signs[:, None]
+
+        # Projected coordinates: U * S == X @ V
+        Wt = U * S
+
+        # Eigenvalues and % variance explained.
+        eigenvalues = (S ** 2) / (n_samples - 1)
+        variance = eigenvalues / eigenvalues.sum()
         pvar = variance / variance.sum()
-        eigenvalues = pca.explained_variance_
 
         if self.transpose:
-            # Use the projected coordinates for the transposed matrix
-            Wt = np.dot(m, Wt.T).T
+            # With samples as observations, U * S already gives each sample's
+            # projection onto the PCs (rows=samples, cols=components). Orient as
+            # (components, samples) to match the indexing used below.
+            Wt = Wt.T
 
         if plot_filename is not None:
             n = n_bars = len(self.labels)
             if eigenvalues.size < n:
                 n_bars = eigenvalues.size
+            # The requested principal components must exist.
+            if max(PCs) > eigenvalues.size:
+                sys.exit("Cannot plot PC{}: only {} principal component(s) are "
+                         "available. Reduce --PCs or increase --ntop.\n".format(max(PCs), eigenvalues.size))
+            # In the untransposed layout each point is a sample indexed along
+            # the component axis, so there must be at least as many components
+            # as samples (i.e. enough usable rows / a large enough --ntop).
+            if not self.transpose and Wt.shape[1] < n:
+                sys.exit("Not enough principal components ({}) to plot {} "
+                         "samples; increase --ntop to at least the sample "
+                         "count.\n".format(Wt.shape[1], n))
             markers = itertools.cycle(matplotlib.markers.MarkerStyle.filled_markers)
             if cols is not None:
                 colors = itertools.cycle(cols)

pydeeptools/deeptools/test/test_plotPCA.py

@@ -1,6 +1,7 @@
 import os
 import filecmp
 import numpy as np
+import pytest
 from matplotlib.testing.compare import compare_images
 from tempfile import NamedTemporaryFile
 import deeptools.plotPCA
@@ -12,6 +13,63 @@
 print(ROOT)
 tolerance = 50
 
+
+def _run_pca(extra=None, plot=True):
+    """Run plotPCA over the shared test matrix and return the parsed
+    --outFileNameData table (header stripped). ``extra`` is a list of extra
+    CLI tokens. When ``plot`` is True a plot file is also requested so the
+    full plotting path runs; set it False to exercise only the numeric output
+    (e.g. small --ntop values whose plotting path is separately broken)."""
+    tsvfile = NamedTemporaryFile(suffix='.tsv', prefix='deeptools_testfile_', delete=False)
+    args = "-in {0}test_samples.npz --outFileNameData {1}".format(
+        TEST_DATA, tsvfile.name).split()
+    plotfile = None
+    if plot:
+        plotfile = NamedTemporaryFile(suffix='.png', prefix='deeptools_testfile_', delete=False)
+        args += ["-o", plotfile.name]
+    if extra:
+        args += extra
+    deeptools.plotPCA.main(args)
+    data = np.loadtxt(tsvfile.name, skiprows=1)
+    os.remove(tsvfile.name)
+    if plotfile is not None:
+        os.remove(plotfile.name)
+    return data
+
+
+def _sign_fix(coords, component_axis=1):
+    """PCA eigenvector signs are arbitrary: they flip across BLAS/platforms and
+    between implementations (e.g. sklearn vs a scipy/SVD rewrite). The sign
+    freedom is per principal component, so normalize each component's vector to
+    have a positive largest-magnitude entry.
+
+    ``component_axis`` says which axis indexes the principal components:
+    in the untransposed --outFileNameData table components are the columns
+    (axis=1); in the transposed table they are the rows (axis=0)."""
+    coords = np.array(coords, dtype=float)
+    if component_axis == 1:
+        for j in range(coords.shape[1]):
+            i = np.argmax(np.abs(coords[:, j]))
+            if coords[i, j] < 0:
+                coords[:, j] = -coords[:, j]
+    else:
+        for i in range(coords.shape[0]):
+            j = np.argmax(np.abs(coords[i, :]))
+            if coords[i, j] < 0:
+                coords[i, :] = -coords[i, :]
+    return coords
+
+
+# Golden eigenvalues captured from the sklearn-backed implementation over
+# test_samples.npz with the default --ntop 500. Eigenvalues are the portable
+# invariant (stable across BLAS backends and across the scipy/SVD rewrite);
+# untransposed per-feature coordinates are not (see test_plotPCA_default_eigenvalues).
+_GOLDEN_DEFAULT_EIGENVALUES = np.array([
+    5.807692278756, 0.074230288836, 0.048971777735,
+    0.036809415525, 0.026706723301, 0.017613563943,
+])
+
+
 def test_plotPCA_default():
     plotfile = NamedTemporaryFile(suffix='.png', prefix='deeptools_testfile_', delete=False)
     tsvfile  = NamedTemporaryFile(suffix='.tsv', prefix='deeptools_testfile_', delete=False)
@@ -51,4 +109,150 @@ def test_plotPCA_outFileNameData():
     np.testing.assert_allclose(eigenvalues, expected_eigenvalues, rtol=1e-5)
 
     os.remove(plotfile.name)
-    os.remove(tsvfile.name)
+    os.remove(tsvfile.name)
+
+
+def test_plotPCA_default_eigenvalues():
+    """Regression on the untransposed eigenvalues, the portable numeric
+    invariant of this path.
+
+    We deliberately do NOT assert the projected coordinates here. After PC1
+    the eigenvalues are tiny and near-degenerate (~0.07, 0.05, 0.04, ...), so
+    the corresponding eigenvectors are free to rotate within that subspace,
+    and the top-``ntop`` row selection (np.argpartition) breaks variance ties
+    differently across BLAS backends (Linux OpenBLAS vs macOS Accelerate).
+    The resulting per-feature coordinates are therefore not reproducible
+    across platforms/implementations. Coordinate-level regression is covered
+    by test_plotPCA_transpose, whose components are well separated and stable.
+    The default plot itself is still pinned by test_plotPCA_default (image
+    comparison)."""
+    data = _run_pca()
+    np.testing.assert_array_equal(data[:, 0], np.arange(1, 7))
+    np.testing.assert_allclose(data[:, -1], _GOLDEN_DEFAULT_EIGENVALUES, rtol=1e-5)
+
+
+def test_plotPCA_variance_matches_eigenvalues():
+    """The per-PC variance fraction shown on the axis labels / scree plot is the
+    eigenvalue proportion. Pin that relationship so the rewrite keeps the two in
+    sync (eigenvalues are monotonically non-increasing and normalize to 1)."""
+    eig = _run_pca()[:, -1]
+    assert np.all(np.diff(eig) <= 1e-9), "eigenvalues must be non-increasing"
+    pvar = eig / eig.sum()
+    np.testing.assert_allclose(pvar.sum(), 1.0, rtol=1e-9)
+    # PC1 dominates on this synthetic wt/kd matrix.
+    assert pvar[0] > 0.9
+
+
+def test_plotPCA_ntop_zero_uses_all_rows():
+    """--ntop 0 disables the top-variable-rows filter and therefore changes the
+    result relative to the default --ntop 500 (the test matrix has >500 rows)."""
+    default = _run_pca()
+    allrows = _run_pca(["--ntop", "0"])
+    assert allrows.shape == (6, 8)
+    # Different row selection -> different eigenvalues.
+    assert not np.allclose(default[:, -1], allrows[:, -1])
+    # Eigenvalues still normalize and stay ordered.
+    eig = allrows[:, -1]
+    assert np.all(np.diff(eig) <= 1e-9)
+
+
+def test_plotPCA_ntop_smaller_than_samples():
+    """When --ntop is below the sample count the table is truncated to the
+    number of retained components (rows)."""
+    # plot=False: the numeric table is well-defined even with 2 features.
+    data = _run_pca(["--ntop", "2"], plot=False)
+    assert data.shape == (2, 4)
+    np.testing.assert_array_equal(data[:, 0], np.arange(1, 3))
+    # First component carries all the variance for the 2-feature case.
+    np.testing.assert_allclose(data[0, -1], 12.0, rtol=1e-6)
+    assert abs(data[1, -1]) < 1e-6
+
+
+def test_plotPCA_ntop_below_samples_plot_errors_cleanly():
+    """Plotting with fewer retained components than samples cannot lay out the
+    scatter; the tool must exit with a clear message rather than crash with an
+    IndexError (previously a bug at correlation.py's scatter loop)."""
+    plotfile = NamedTemporaryFile(suffix='.png', prefix='deeptools_testfile_', delete=False)
+    args = "-in {0}test_samples.npz -o {1} --ntop 2".format(TEST_DATA, plotfile.name).split()
+    try:
+        with pytest.raises(SystemExit) as exc:
+            deeptools.plotPCA.main(args)
+        assert "principal component" in str(exc.value)
+    finally:
+        if os.path.exists(plotfile.name):
+            os.remove(plotfile.name)
+
+
+def test_plotPCA_PCs_selection_does_not_change_table():
+    """--PCs only selects which components are drawn; the numeric table always
+    contains every component, so it is independent of --PCs."""
+    default = _run_pca()
+    pcs13 = _run_pca(["--PCs", "1", "3"])
+    np.testing.assert_allclose(default, pcs13, rtol=1e-9, atol=1e-12)
+
+
+@pytest.mark.parametrize("extra, msg", [
+    (["--PCs", "2", "2"], "different principal components"),
+    (["--PCs", "0", "1"], "at least 1"),
+    (["--ntop", "-1"], "must be >= 0"),
+])
+def test_plotPCA_invalid_arguments_exit(extra, msg):
+    plotfile = NamedTemporaryFile(suffix='.png', prefix='deeptools_testfile_', delete=False)
+    args = "-in {0}test_samples.npz -o {1}".format(TEST_DATA, plotfile.name).split() + extra
+    try:
+        with pytest.raises(SystemExit) as exc:
+            deeptools.plotPCA.main(args)
+        assert msg in str(exc.value)
+    finally:
+        if os.path.exists(plotfile.name):
+            os.remove(plotfile.name)
+
+
+def test_plotPCA_requires_an_output():
+    with pytest.raises(SystemExit) as exc:
+        deeptools.plotPCA.main("-in {0}test_samples.npz".format(TEST_DATA).split())
+    assert "must be specified" in str(exc.value)
+
+
+# Golden values for the transposed PCA (samples as observations, so each row
+# of the table is a component's projection across the six samples). Captured
+# after fixing the projection bug; stored raw (components are rows -> axis=0).
+_GOLDEN_TRANSPOSE_COORDS = np.array([
+    [8.096369192617, 27.65422672360, -1.598082844166, -15.48892072797, -18.49767188707, -0.1659204570140],
+    [3.552671394141, -4.837722476763, 20.02992087542, 0.4882670876827, -7.713434681130, -11.51970219935],
+    [10.09925342625, -9.161879672887, 1.351881375625, -11.06368229223, -0.2099739792611, 8.984401142503],
+    [-11.75933735644, 2.772538322120, 7.637287613484, -8.588903203498, 5.490108421837, 4.448306202499],
+    [4.468893041249, 2.035996403779, -1.674415783401, -5.444924755301, 9.786060422166, -9.171609328491],
+    [3.400996688227e-15, 3.400996688227e-15, 3.400996688227e-15, 3.400996688227e-15, 3.400996688227e-15, 3.400996688227e-15],
+])
+_GOLDEN_TRANSPOSE_EIGENVALUES = np.array([
+    282.9918757435, 125.9323562327, 78.18623219731,
+    65.59902453902, 47.29051128753, 1.388013416800e-29,
+])
+
+
+def test_plotPCA_transpose():
+    """--transpose runs (previously crashed) and projects each sample onto the
+    PCs. Coordinates are compared sign-invariantly; the last component is a
+    numerical-zero residual so we skip its unstable sign."""
+    data = _run_pca(["--transpose"])
+    assert data.shape == (6, 8)
+    np.testing.assert_array_equal(data[:, 0], np.arange(1, 7))
+    # Transposed table: components are rows -> sign-fix per row (axis=0).
+    coords = _sign_fix(data[:, 1:7], component_axis=0)
+    golden = _sign_fix(_GOLDEN_TRANSPOSE_COORDS, component_axis=0)
+    # Compare the informative components; the final ~1e-15 residual row is noise.
+    np.testing.assert_allclose(coords[:-1], golden[:-1], rtol=1e-4, atol=1e-6)
+    np.testing.assert_allclose(data[:, -1], _GOLDEN_TRANSPOSE_EIGENVALUES, rtol=1e-4, atol=1e-6)
+    # Transposed eigenvalues differ from the untransposed layout.
+    assert not np.allclose(data[:, -1], _GOLDEN_DEFAULT_EIGENVALUES)
+
+
+def test_plotPCA_log2_and_rowCenter_affect_output():
+    """--log2 and --rowCenter now actually transform the data before the PCA,
+    so each changes the result relative to the default."""
+    default = _run_pca()
+    log2 = _run_pca(["--log2"])
+    rowcenter = _run_pca(["--rowCenter"])
+    assert not np.allclose(default[:, -1], log2[:, -1]), "--log2 was a no-op"
+    assert not np.allclose(default[:, -1], rowcenter[:, -1]), "--rowCenter was a no-op"

pyproject.toml

@@ -34,8 +34,6 @@ dependencies = [
     "pysam >= 0.23",
     "pyBigWig >= 0.3",
     "py2bit >= 0.3",
-    "pandas >= 2.2",
-    "scikit-learn >= 1.6",
     "deeptoolsintervals >= 0.1",
     "maturin"
 ]