Replace MDS solver; match variable naming to paper (#4)

Pandora now uses scikit-allel MDS instead of scikit-learn MDS and the confidence_level variable/CLI flag was renamed to convergence tolerance to match the terminology of the paper

Author: tschuelia

docs/cli_config.rst

@@ -12,7 +12,7 @@ Configuration options:
 - ``file_format``, default = ``EIGENSTRAT``, Name of the file format your dataset is in. Supported formats are ``ANCESTRYMAP``, ``EIGENSTRAT``, ``PED``, ``PACKEDPED``, ``PACKEDANCESTRYMAP``. For more information see Section `Input data`_ below.
 - ``convertf``, default = ``convertf``, File path pointing to an executable of Eigensoft's ``convertf`` tool. ``convertf`` is used if the provided dataset is not in ``EIGENSTRAT`` format. Default is ``convertf``. This will only work if ``convertf`` is installed systemwide.
 - ``bootstrap_convergence_check``, default = ``True``, If true, instead of computing ``n_replicates`` bootstraps and embeddings, Pandora will check for convergence once every ``max(10, threads)`` bootstrap embeddings are computed. If according to our heuristic (see TODO for more details) the bootstrap procedure converged, all remaining tasks are cancelled and the stability is determined uisng only the number of replicates computed when convergence is determined. Due to the runtime overhead of the convergence check compared to the runtime of MDS computations, we only advice using this convergence check for PCA analyses. Note that this parameter is only relevant if ``analysis_mode`` is ``AnalysisMode.BOOTSTRAP``.
-- ``bootstrap_convergence_confidence_level``, default=0.05, Determines the level of confidence when checking for bootstrap convergence. A value of :math:`X` means that we allow deviations of up to :math:`X * 100\%` between pairwise bootstrap comparisons and still assume convergence.
+- ``bootstrap_convergence_tolerance``, default=0.05, Determines the level of deviation tolerance when checking for bootstrap convergence. A value of :math:`X` means that we allow deviations of up to :math:`X * 100\%` between pairwise bootstrap comparisons and still assume convergence.
 - ``n_replicates``, default = 100, Number of bootstrap replicates or sliding windows to compute
 - ``keep_replicates``, default = ``false``, Whether to store all intermediate datasets files (``.geno``, ``.snp``, ``.ind``). Note that this will result in a substantial storage consumption. Note that in case of bootstrapping, the bootstrapped indices are stored as checkpoints for full reproducibility in any case.
 - ``n_components``, default = 10, Number of components to compute and compare for PCA or MDS analyses. We recommend 10 for PCA analyses and 2 for MDS analyses. The default is 10 since the default for ``embedding_algorithm`` is ``PCA``.

docs/getting_started.rst

@@ -46,12 +46,12 @@ have to add the path to ``smartpca`` to the ``config-example.yaml``.
 
 You should then see an output similar to this:::
 
-    Pandora version 1.0.3 released by The Exelixis Lab
+    Pandora version 2.0.0 released by The Exelixis Lab
     Developed by: Julia Haag
     Latest version: https://github.com/tschuelia/Pandora
     Questions/problems/suggestions? Please open an issue on GitHub.
 
-    Pandora was called at 06-Nov-2023 16:26:50 as follows:
+    Pandora was called at 18-Apr-2024 16:26:50 as follows:
 
     /Users/julia/micromamba/envs/pandora/bin/pandora -c config_example.yaml
 
@@ -63,7 +63,7 @@ You should then see an output similar to this:::
     n_replicates: 10
     keep_replicates: False
     bootstrap_convergence_check: True
-    bootstrap_convergence_confidence_level: 0.05
+    bootstrap_convergence_tolerance: 0.05
     n_components: 10
     embedding_algorithm: PCA
     smartpca: smartpca
@@ -85,7 +85,7 @@ You should then see an output similar to this:::
     [00:00:02] Running SmartPCA on the input dataset.
     [00:00:02] Plotting embedding results for the input dataset.
     [00:00:18] Drawing 10 bootstrapped datasets and running PCA.
-    [00:00:18] NOTE: Bootstrap convergence check is enabled. Will terminate bootstrap computation once convergence is determined. Convergence confidence level: 0.05
+    [00:00:18] NOTE: Bootstrap convergence check is enabled. Will terminate bootstrap computation once convergence is determined. Convergence tolerance: 0.05
     [00:00:27] Bootstrapping done. Number of replicates computed: 10
     [00:00:27] Comparing bootstrapping embedding results.
     [00:00:34] Plotting bootstrapping embedding results.

docs/usage.rst

@@ -40,11 +40,11 @@ The Command line interface, as well as when using the Eigen-based Pandora Python
 Smartpca is a powerful PCA tool that implements a lot of genotype-data specific routines and optimizations and provides a lot of useful options for meaningful PCA analyses such as outlier detection.
 Pandora supports all custom configuration settings of smartpca. See the Section :ref:`SmartPCA` for more information. For MDS analyses, Pandora will use
 smartpca to generate the Fst-distance matrix as input for MDS. Note that this distance matrix computes the distances between population and not between samples.
-The subsequent MDS analysis is performed using the scikit-learn MDS implementation.
+The subsequent MDS analysis is performed using the scikit-allel MDS implementation.
 If you have genotype data in Eigenfiles but want to be able to do a more flexible analysis, consider using the alternative NumPy interface. Pandora provides a method
 to load your genotype data in EIGENSTRAT format as numpy array.
 
-If you are using the NumPy-based Pandora interface, PCA and MDS is performed using the scikit-learn implementations. For both analyses, Pandora supports different types of data imputation, see the API documentation for more information.
+If you are using the NumPy-based Pandora interface, PCA and MDS is performed using the scikit-learn and scikit-allel implementations respectively. For both analyses, Pandora supports different types of data imputation, see the API documentation for more information.
 Per default, Pandora will apply SNP-wise mean imputation. The default distance metric for MDS analysis is the pairwise euclidean distance between all samples in your data. However, Pandora provides alternative distance metrics
 and allows you to define your own distance metric as well. Again, see the API documentation for further information.
 

pandora/bootstrap.py

@@ -29,7 +29,7 @@
 def _bootstrap_convergence_check(
     bootstraps: List[Union[NumpyDataset, EigenDataset]],
     embedding: EmbeddingAlgorithm,
-    bootstrap_convergence_confidence_level: float,
+    bootstrap_convergence_tolerance: float,
     threads: int,
     logger: Optional[loguru.Logger] = None,
 ):
@@ -44,14 +44,12 @@ def _bootstrap_convergence_check(
             f"Unrecognized embedding option {embedding}. Supported are 'pca' and 'mds'."
         )
 
-    return _bootstrap_converged(
-        embeddings, bootstrap_convergence_confidence_level, threads
-    )
+    return _bootstrap_converged(embeddings, bootstrap_convergence_tolerance, threads)
 
 
 def _bootstrap_converged(
     bootstraps: List[Embedding],
-    bootstrap_convergence_confidence_level: float,
+    bootstrap_convergence_tolerance: float,
     threads: int,
 ):
     """Checks for convergence by comparing the Pandora Stabilities for 10 subsets of the given list of bootstraps."""
@@ -76,7 +74,7 @@ def _bootstrap_converged(
             stabilities[j] = stability_s2
 
         relative_difference = abs(stability_s2 - stability_s1) / (stability_s1 + 1e-6)
-        if round(relative_difference, 3) > bootstrap_convergence_confidence_level:
+        if round(relative_difference, 3) > bootstrap_convergence_tolerance:
             return False
     return True
 
@@ -252,7 +250,7 @@ def run(
         self,
         threads: int,
         bootstrap_convergence_check: bool,
-        bootstrap_convergence_confidence_level: float,
+        bootstrap_convergence_tolerance: float,
         embedding: EmbeddingAlgorithm,
         logger: Optional[loguru.Logger] = None,
     ):
@@ -296,7 +294,7 @@ def run(
                         converged = _bootstrap_convergence_check(
                             bootstraps,
                             embedding,
-                            bootstrap_convergence_confidence_level,
+                            bootstrap_convergence_tolerance,
                             threads,
                             logger,
                         )
@@ -390,7 +388,7 @@ def bootstrap_and_embed_multiple(
     redo: bool = False,
     keep_bootstraps: bool = False,
     bootstrap_convergence_check: bool = True,
-    bootstrap_convergence_confidence_level: float = 0.05,
+    bootstrap_convergence_tolerance: float = 0.05,
     smartpca_optional_settings: Optional[Dict] = None,
     logger: Optional[loguru.Logger] = None,
 ) -> List[EigenDataset]:
@@ -434,8 +432,8 @@ def bootstrap_and_embed_multiple(
     bootstrap_convergence_check : bool, default=True
         Whether to automatically determine bootstrap convergence. If ``True``, will only compute as many replicates as
         required for convergence according to our heuristic (see Notes below).
-    bootstrap_convergence_confidence_level : float, default=0.05
-        Determines the level of confidence when checking for bootstrap convergence. A value of X means that we allow
+    bootstrap_convergence_tolerance : float, default=0.05
+        Determines the deviation tolerance when checking for bootstrap convergence. A value of X means that we allow
         deviations of up to :math:`X * 100\\%` between pairwise bootstrap comparisons and still assume convergence.
     smartpca_optional_settings : Dict, default=None
         Additional smartpca settings.
@@ -462,7 +460,7 @@ def bootstrap_and_embed_multiple(
     We first create 10 random subsets of size :math:`int(N/2)` by sampling from all :math:`N` replicates.
     We then compute the Pandora Stability (PS) for each of the 10 subsets and compute the relative difference of PS
     values between all possible pairs of subsets :math:`(PS_1, PS_2)` by computing :math:`\\frac{\\left|PS_1 - PS_2\\right|}{PS_2}`.
-    We assume convergence if all pairwise relative differences are below X * 100% were X is the set confidence level.
+    We assume convergence if all pairwise relative differences are below X * 100% were X is the set tolerance.
     If we determine that the bootstrap has converged, all remaining bootstrap computations are cancelled.
 
     (*) The reasoning for checking every ``max(10, threads)`` is the following: if Pandora runs on a machine with e.g. 48
@@ -517,7 +515,7 @@ def bootstrap_and_embed_multiple(
     bootstraps, finished_indices = parallel_bootstrap_process_manager.run(
         threads=threads,
         bootstrap_convergence_check=bootstrap_convergence_check,
-        bootstrap_convergence_confidence_level=bootstrap_convergence_confidence_level,
+        bootstrap_convergence_tolerance=bootstrap_convergence_tolerance,
         embedding=embedding,
         logger=logger,
     )
@@ -568,7 +566,7 @@ def bootstrap_and_embed_multiple_numpy(
     ] = euclidean_sample_distance,
     imputation: Optional[str] = "mean",
     bootstrap_convergence_check: bool = True,
-    bootstrap_convergence_confidence_level: float = 0.05,
+    bootstrap_convergence_tolerance: float = 0.05,
 ) -> List[NumpyDataset]:
     """Draws ``n_replicates`` bootstrap datasets of the provided NumpyDataset and performs PCA/MDS analysis (as
     specified by ``embedding``) for each bootstrap.
@@ -611,8 +609,8 @@ def bootstrap_and_embed_multiple_numpy(
     bootstrap_convergence_check : bool, default=True
         Whether to automatically determine bootstrap convergence. If ``True``, will only compute as many replicates as
         required for convergence according to our heuristic (see Notes below).
-    bootstrap_convergence_confidence_level : float, default=0.05
-        Determines the level of confidence when checking for bootstrap convergence. A value of X means that we allow
+    bootstrap_convergence_tolerance : float, default=0.05
+        Determines the level of deviation tolerance when checking for bootstrap convergence. A value of X means that we allow
         deviations of up to :math:`X * 100\\%` between pairwise bootstrap comparisons and still assume convergence.
 
     Returns
@@ -633,7 +631,7 @@ def bootstrap_and_embed_multiple_numpy(
     We first create 10 random subsets of size :math:`int(N/2)` by sampling from all :math:`N` replicates.
     We then compute the Pandora Stability (PS) for each of the 10 subsets and compute the relative difference of PS
     values between all possible pairs of subsets :math:`(PS_1, PS_2)` by computing :math:`\\frac{\\left|PS_1 - PS_2\\right|}{PS_2}`.
-    We assume convergence if all pairwise relative differences are below X * 100% were X is the set confidence level.
+    We assume convergence if all pairwise relative differences are below X * 100% were X is the set tolerance.
     If we determine that the bootstrap has converged, all remaining bootstrap computations are cancelled.
 
     (*) The reasoning for checking every ``max(10, threads)`` is the following: if Pandora runs on a machine with e.g. 48
@@ -666,7 +664,7 @@ def bootstrap_and_embed_multiple_numpy(
     bootstraps, _ = parallel_bootstrap_process_manager.run(
         threads=threads,
         bootstrap_convergence_check=bootstrap_convergence_check,
-        bootstrap_convergence_confidence_level=bootstrap_convergence_confidence_level,
+        bootstrap_convergence_tolerance=bootstrap_convergence_tolerance,
         embedding=embedding,
     )