Skip to content

Commit b94030a

Browse files
adRn-sclaude
andcommitted
Compute PCA with scipy instead of scikit-learn
Replace the sklearn PCA and StandardScaler in Correlation.plot_pca with a scipy/numpy SVD implementation, reproducing sklearn's output (column standardization with population std, deterministic svd_flip sign convention, explained_variance_ from singular values). Drop the now-unused pandas and scikit-learn dependencies. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
1 parent 045a821 commit b94030a

2 files changed

Lines changed: 25 additions & 15 deletions

File tree

pydeeptools/deeptools/correlation.py

Lines changed: 25 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,7 @@
1515
import matplotlib.markers
1616
import matplotlib.colors as pltcolors
1717
from deeptools.utilities import toString, convertCmap
18-
from sklearn.decomposition import PCA
19-
from sklearn.preprocessing import StandardScaler
18+
from scipy.linalg import svd
2019

2120
class Correlation:
2221
"""
@@ -475,18 +474,31 @@ def plot_pca(self, plot_filename=None, PCs=[1, 2], plot_title='', image_format=N
475474
if self.transpose:
476475
m = m.T
477476

478-
# Center and scale
479-
scaler = StandardScaler()
480-
m2 = scaler.fit_transform(m)
481-
482-
# PCA
483-
pca = PCA()
484-
Wt = pca.fit_transform(m2)
485-
486-
# % variance, eigenvalues
487-
variance = pca.explained_variance_ratio_
477+
# Center and scale each column to zero mean and unit variance
478+
# (equivalent to sklearn's StandardScaler, using population std/ddof=0).
479+
col_mean = m.mean(axis=0)
480+
col_std = m.std(axis=0)
481+
col_std[col_std == 0] = 1.0
482+
m2 = (m - col_mean) / col_std
483+
484+
# PCA via SVD of the (re-)centered matrix, mirroring sklearn's PCA().
485+
n_samples = m2.shape[0]
486+
X = m2 - m2.mean(axis=0)
487+
U, S, Vt = svd(X, full_matrices=False)
488+
489+
# Deterministic sign convention (sklearn's svd_flip, u_based_decision).
490+
max_abs_cols = np.argmax(np.abs(U), axis=0)
491+
signs = np.sign(U[max_abs_cols, range(U.shape[1])])
492+
U *= signs
493+
Vt *= signs[:, None]
494+
495+
# Projected coordinates: U * S == X @ V
496+
Wt = U * S
497+
498+
# Eigenvalues and % variance explained.
499+
eigenvalues = (S ** 2) / (n_samples - 1)
500+
variance = eigenvalues / eigenvalues.sum()
488501
pvar = variance / variance.sum()
489-
eigenvalues = pca.explained_variance_
490502

491503
if self.transpose:
492504
# Use the projected coordinates for the transposed matrix

pyproject.toml

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,8 +34,6 @@ dependencies = [
3434
"pysam >= 0.23",
3535
"pyBigWig >= 0.3",
3636
"py2bit >= 0.3",
37-
"pandas >= 2.2",
38-
"scikit-learn >= 1.6",
3937
"deeptoolsintervals >= 0.1",
4038
"maturin"
4139
]

0 commit comments

Comments
 (0)