plotPCA: write and plot the per-sample loadings in the untransposed layout

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

galaxy/wrapper/plotPCA.xml

@@ -79,16 +79,17 @@
             <param name="outFileNameData" value="True" />
             <output name="outFileName" file="plotPCA_result2.png" ftype="png" compare="sim_size" delta="12000" />
             <output name="output_outFileNameData" ftype="tabular">
-                <!-- The duplicate-sample column is mathematically zero, so it comes
-                     out as tiny BLAS-dependent noise (~1e-18) that varies across
-                     platforms. Assert structure, the stable PC1 coordinates (sign
-                     may flip across BLAS, hence -?), and the leading eigenvalue;
-                     do not pin the degenerate zero column or eigenvalue. -->
+                <!-- The two samples are the same file, so PC1 loads both equally
+                     (+-1/sqrt(2)) with eigenvalue 3 and PC2 is degenerate: its
+                     loadings and eigenvalue (~1e-35) are BLAS-dependent noise.
+                     Assert structure, the stable PC1 loadings (sign may flip
+                     across BLAS, hence -?), and the leading eigenvalue; do not pin
+                     the degenerate second component. -->
                 <assert_contents>
                     <has_n_columns n="4" />
                     <has_line_matching expression="Component\tbowtie2-test1.bam\tbowtie2-test1.bam\tEigenvalue" />
-                    <has_line_matching expression="1\t-?1\.99949\d*\t.*\t2\.9999\d*" />
-                    <has_line_matching expression="2\t-?0\.96079\d*\t.*" />
+                    <has_line_matching expression="1\t-?0\.70710\d*\t-?0\.70710\d*\t2\.9999\d*" />
+                    <has_line_matching expression="2\t-?\d.*\t-?\d.*\t.*" />
                 </assert_contents>
             </output>
         </test>

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

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)

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."""

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

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