diff --git a/CHANGES.txt b/CHANGES.txt index 56dfe7914d..552cf3488a 100755 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -16,6 +16,7 @@ * alignmentSieve output order matches input order exactly * --missingDataAsZero no longer takes bases exceeding chromosome bounds as 0 values but rather purges the bins * large scale values precision slightly altered with new backend (f32 vs f64) +* plotPCA (untransposed, the default) again writes and plots the loadings of each sample on the principal components (rows of V^T); the re-implementation wrote the PC scores of the first rows of the matrix under the sample names. --transpose was correct 3.5.6 * minimal supported python version raised to 3.9 (numpy >= 2 support); NaN handling switched to np.nan diff --git a/galaxy/wrapper/plotPCA.xml b/galaxy/wrapper/plotPCA.xml index f5fe600e88..3f62bf9da6 100644 --- a/galaxy/wrapper/plotPCA.xml +++ b/galaxy/wrapper/plotPCA.xml @@ -79,16 +79,17 @@ - + - - + + diff --git a/galaxy/wrapper/test-data/plotPCA_result1.png b/galaxy/wrapper/test-data/plotPCA_result1.png index b043d20f58..9909578eeb 100644 Binary files a/galaxy/wrapper/test-data/plotPCA_result1.png and b/galaxy/wrapper/test-data/plotPCA_result1.png differ diff --git a/galaxy/wrapper/test-data/plotPCA_result2.png b/galaxy/wrapper/test-data/plotPCA_result2.png index 1f4338467c..9909578eeb 100644 Binary files a/galaxy/wrapper/test-data/plotPCA_result2.png and b/galaxy/wrapper/test-data/plotPCA_result2.png differ diff --git a/galaxy/wrapper/test-data/plotPCA_result2.tabular b/galaxy/wrapper/test-data/plotPCA_result2.tabular index 3ae66a2321..c239deebdb 100644 --- a/galaxy/wrapper/test-data/plotPCA_result2.tabular +++ b/galaxy/wrapper/test-data/plotPCA_result2.tabular @@ -1,3 +1,3 @@ Component bowtie2-test1.bam bowtie2-test1.bam Eigenvalue -1 1.9994942044176225 0.0 2.9999999999999996 -2 -0.9607959164084681 0.0 0.0 +1 0.7071067811865476 0.7071067811865475 2.9999999999999996 +2 0.7071067811865475 -0.7071067811865476 1.558261167288188e-35 diff --git a/pydeeptools/deeptools/correlation.py b/pydeeptools/deeptools/correlation.py index cd7a35acb6..cb225eb2bd 100644 --- a/pydeeptools/deeptools/correlation.py +++ b/pydeeptools/deeptools/correlation.py @@ -497,19 +497,22 @@ def plot_pca(self, plot_filename=None, PCs=[1, 2], plot_title='', image_format=N 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() if self.transpose: - # With samples as observations, U * S already gives each sample's + # Samples are the observations: U * S == X @ V 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 + Wt = (U * S).T + else: + # Samples are the variables: what is written and plotted per sample + # is its loading on each component, i.e. the rows of V^T + # (components, samples). U * S would be the scores of the rows + # (bins) of the matrix, one per selected bin, not per sample. + Wt = Vt if plot_filename is not None: n = n_bars = len(self.labels) @@ -519,13 +522,6 @@ def plot_pca(self, plot_filename=None, PCs=[1, 2], plot_title='', image_format=N 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) diff --git a/pydeeptools/deeptools/test/test_plotPCA.py b/pydeeptools/deeptools/test/test_plotPCA.py index a284677764..dac41372ae 100644 --- a/pydeeptools/deeptools/test/test_plotPCA.py +++ b/pydeeptools/deeptools/test/test_plotPCA.py @@ -158,22 +158,21 @@ def test_plotPCA_ntop_zero_uses_all_rows(): 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) + number of retained components (rows); every row still holds one loading + per sample, so the plot of PC1 vs PC2 is well-defined.""" + data = _run_pca(["--ntop", "2"]) + assert data.shape == (2, 8) 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).""" +def test_plotPCA_requested_PC_beyond_ntop_errors_cleanly(): + """Asking for a component that the retained rows cannot provide must exit + with a clear message rather than crash with an IndexError.""" 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() + args = "-in {0}test_samples.npz -o {1} --ntop 2 --PCs 1 3".format(TEST_DATA, plotfile.name).split() try: with pytest.raises(SystemExit) as exc: deeptools.plotPCA.main(args) @@ -183,6 +182,25 @@ def test_plotPCA_ntop_below_samples_plot_errors_cleanly(): os.remove(plotfile.name) +def test_plotPCA_default_table_holds_sample_loadings(): + """In the default (untransposed) layout the table row for component i is + the loading of every sample on PC i, i.e. row i of V^T from the SVD of the + standardised (ntop rows x samples) matrix: the rows are orthonormal and + equal, up to the per-component sign, to a numpy PCA on the same rows. + (The bin scores U*S, one row per selected bin, are neither.)""" + data = _run_pca() + W = data[:, 1:7] + np.testing.assert_allclose(W @ W.T, np.eye(6), atol=1e-6) + + mat = np.load(TEST_DATA + "test_samples.npz")["matrix"].astype(float) + m = mat[np.argpartition(mat.var(axis=1), -500)[-500:]] + m = (m - m.mean(axis=0)) / m.std(axis=0) + _, _, Vt = np.linalg.svd(m - m.mean(axis=0), full_matrices=False) + # PC1 dominates and is well separated, so its loadings are stable + # across BLAS backends; compare it sign-invariantly. + np.testing.assert_allclose(_sign_fix(W[:1].T)[:, 0], _sign_fix(Vt[:1].T)[:, 0], rtol=1e-5, atol=1e-8) + + 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.""" diff --git a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.png b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.png index e27190e1c4..deb6ececb1 100644 Binary files a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.png and b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.png differ diff --git a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.tsv b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.tsv index 2825477a63..b568eac027 100644 --- a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.tsv +++ b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_default.tsv @@ -1,7 +1,7 @@ Component wt1 wt2 wt3 kd1 kd2 kd3 Eigenvalue -1 -1.9593953629718786 -0.0599263205226625 0.015514438980507314 0.3608726269539977 -0.06417725736658203 0.012713770441420292 5.807692278755936 -2 -3.359410695125286 0.1453930968029311 -0.13474864281257853 -0.012411165664075682 -0.12363598259527717 -0.060663076611993313 0.07423028883557825 -3 -2.335822964750958 -0.0998772259871388 -0.03413541974064263 -0.07165342296216205 0.05170385813600817 0.13732570368291072 0.04897177773493676 -4 -0.3078978392240127 -0.004185214180292667 -0.013424388619809201 0.051462455776524134 -0.0794287737098722 0.07458644503209262 0.03680941552538939 -5 -1.280108686154136 0.0779530163424985 -0.04768720563614561 -0.023214274278217935 -0.0324056770005561 0.0742066588868442 0.026706723301448194 -6 -3.1274153845376063 -0.17202014078313563 0.08493798982333454 0.014293535309782737 0.01850856292486998 0.1575339674343244 0.017613563942900697 +1 0.41015948328129415 0.40714083168340265 0.4051015001575071 0.4088383560914448 0.40895832745985133 0.4092708162212304 5.807692278755935 +2 -0.10564992452808587 -0.13903863717051682 -0.7036293245034178 -0.031273785139079135 0.5057293580802189 0.46655320303121317 0.07423028883557847 +3 -0.3072826224468924 0.8073633171644164 -0.366048315157121 0.2231968254778915 -0.12059197662562816 -0.23535345181275683 0.048971777734937216 +4 -0.16638072203492837 -0.3999296950298342 -0.0911203187723264 0.8821546210006602 -0.10280300047534902 -0.12371559972677873 0.03680941552539014 +5 0.8255039368543798 0.008574695725392222 -0.4455301735837803 0.057473412505828825 -0.3107944624876728 -0.1416905014440879 0.026706723301448895 +6 -0.13055798419418163 0.05524299298321737 0.0026704255973279928 0.023587232780211485 -0.6747494521159683 0.7239147136476167 0.017613563942900895 diff --git a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.png b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.png index 283f5e1845..4d9506455f 100644 Binary files a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.png and b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.png differ diff --git a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.tsv b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.tsv index 9822010ab4..b568eac027 100644 --- a/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.tsv +++ b/pydeeptools/deeptools/test/test_plotPCA/test_plotPCA_ggplot.tsv @@ -1,7 +1,7 @@ Component wt1 wt2 wt3 kd1 kd2 kd3 Eigenvalue -1 -1.9593953629718677 0.05992632052266172 0.015514438980501043 0.360872626954001 -0.06417725736657884 0.012713770441428606 5.807692278755935 -2 -3.3594106951252853 -0.1453930968029286 -0.134748642812576 -0.012411165664076072 -0.12363598259528208 -0.060663076611996596 0.07423028883557839 -3 -2.335822964750958 0.09987722598713754 -0.03413541974064124 -0.07165342296216681 0.05170385813600347 0.1373257036829112 0.0489717777349372 -4 -0.30789783922401226 0.004185214180292798 -0.013424388619808516 0.051462455776524196 -0.07942877370987256 0.0745864450320923 0.036809415525390125 -5 -1.2801086861541353 -0.07795301634249786 -0.04768720563614555 -0.02321427427821972 -0.032405677000559435 0.07420665888684398 0.026706723301448902 -6 -3.127415384537607 0.17202014078313335 0.08493798982333509 0.01429353530978131 0.018508562924871238 0.15753396743432724 0.017613563942900885 +1 0.41015948328129415 0.40714083168340265 0.4051015001575071 0.4088383560914448 0.40895832745985133 0.4092708162212304 5.807692278755935 +2 -0.10564992452808587 -0.13903863717051682 -0.7036293245034178 -0.031273785139079135 0.5057293580802189 0.46655320303121317 0.07423028883557847 +3 -0.3072826224468924 0.8073633171644164 -0.366048315157121 0.2231968254778915 -0.12059197662562816 -0.23535345181275683 0.048971777734937216 +4 -0.16638072203492837 -0.3999296950298342 -0.0911203187723264 0.8821546210006602 -0.10280300047534902 -0.12371559972677873 0.03680941552539014 +5 0.8255039368543798 0.008574695725392222 -0.4455301735837803 0.057473412505828825 -0.3107944624876728 -0.1416905014440879 0.026706723301448895 +6 -0.13055798419418163 0.05524299298321737 0.0026704255973279928 0.023587232780211485 -0.6747494521159683 0.7239147136476167 0.017613563942900895