Update documentation

Author: ncaptier

docs/source/Examples.rst

@@ -8,4 +8,5 @@ Omics data analysis
     :name: thumbnail_gallery
 
     _examples/transcriptomic_ICA
-    _examples/stability_study
+    _examples/stability_study
+    _examples/cell_cycle_ICA.ipynb

docs/source/ICA.rst

@@ -0,0 +1,57 @@
+Icasso algorithm
+================
+
+The main algorithm of **stabilized-ica**, `sica.base.StabilizedICA <modules/generated/sica.base.StabilizedICA.html#sica.base.StabilizedICA>`_, solves the ICA problem 
+and stabilizes the resulting components through multiple runs of the ICA solver and a final aggregation step. Our implementation is mainly based on the Icasso method 
+which was developped by J. Himberg and A. Hyvarinen in 2003 [1]_. There are three main steps : 
+
+   #. An ICA solver is chosen and its hyperparameters are fixed. The user also sets the number of components to extract (we will call it `n_components`). Then, the solver is run `n_runs` times with different initial conditions and, optionally, with different bootstrap samples of the original data set. At the end, `n_components*n_runs` are extracted.   
+
+   #. The `n_components*n_runs` components are clustered into `n_components` clusters using hierarchical agglomerative clustering with average linkage criterion. The similarity metric is the absolute value of the pearson correlation coefficient between estimates of the ICA components.   
+
+   #. Finally, the centrotype of each cluster is computed and returned as a stabilized independent component. The centrotype is the point in the cluster which has the maximum sum of similarities to other points in the cluster. Besides, for each cluster, an index which quantifies its compactness and isolation is computed in order to assess the quality of the associated stabilized components (i.e its centrotype). It is computed as the difference between the average intra-cluster similarities and average extra-cluster similarities.
+
+.. figure:: images/icasso.png
+   :align: center
+
+   Illustration of the Icasso method with `n_runs = 4` and `n_components = 5`.
+
+Stability of the final components
+---------------------------------
+This algorithm gives us a way to measure the reliability of the ICA components [1] , [2]_. Intuitively, if an ICA component is reliable it should be estimated at every run of the Icasso method and should be associated with a tight cluster at the end. Thus we can assess the stability of a component generated by the Icasso algorithm by measuring the "tightness" of its associated cluster of `n_runs` components. For a given final ICA source :math:`m` and its associated cluster :math:`C_m` we compute its stability index :math:`I_q(C_m)` :
+
+.. math:: I_q (C_m) = \underbrace{\dfrac{1}{|C_m|^2} \sum_{i , j \in C_m} \sigma_{ij}}_{\text{mean intra-cluster similarity}} - \underbrace{\dfrac{1}{|C_m||C_{-m}|} \sum_{i \in C_m} \sum_{j \in C_{-m}} \sigma_{ij}}_{\text{mean extra-cluster similarity}} \quad \in \quad [0 ,1]
+
+where :math:`|C_m|` corresponds to the number of points in the clusters, :math:`|C_{-m}|` corresponds to the number of points outside the cluster and :math:`\sigma_{ij}` corresponds to the similarity beween components i and j (i.e absolute value of Pearson correlation).
+
+This stability index can measure the robustness of the ICA components with respect to two main uncertainty sources :
+
+   * The **algorithmic uncertainty** is measured through the multiple runs of the chosen ICA solver with different random initializations.     
+   * The **statistical uncertainty** is measured through the multiple runs of the chosen ICA solver with different bootstrap samples as input.
+
+This measure of stability allows us to rank the final ICA components extracted with Icasso. Besides, it brings the light on the ICA estimates that most likely relate to meaningful and interpretable latent sources. On the contrary, the most unstable estimates are likely to correspond to noisy directions that could originate from diverse sources such as numerical instabilities, the ICA model inadequacy for our real data, or statistical limitations due to a small sample size.
+
+Choice of the number of components
+----------------------------------
+The choice of the number of ICA components to extract is a crucial step but yet very difficult. Indeed, contrary to PCA, the choice of the order of decomposition may have an impact on all the extracted ICA components. When going from k components to k+1, there is no guarantee that the k first components will remain the same. Actually, it has been noticed that a choice of too few ICA components could lead to "fused components" that would not capture the heterogeneity of the data (under-decomposition) whereas a choice of too many components could lead to noisy "splitted components" (over-decomposition).
+
+The choice of the order of decomposition relates to the determination of the intrinsic data dimension. Many methods address this question (ex: using PCA and the Kaiser rule) but a few are specific to ICA. In [3]_ U. Kairov et al. empirically developped a method for choosing an optimal number of components. In particular they showed that, for gene expression data, this method led to more interpretable and more reproducible components. Their protocol relies on the natural ranking that stabilized ICA offers through the stability index associated to each component. It could be roughly summarized with two steps:
+
+   * For a number of extracted components varying between :math:`n_{min}` and :math:`n_{max}` (ex: :math:`n_{min} = 2 ,3 , ... , n_{max} = 100`), they run the stabilized ICA algorithm and plot the associated stability profile (i.e the stability indexes of the components ordered from the most stable to the least one). At the end, within a single plot they obtain several stability profiles, associated to different order of decomposition (cf figure below).
+   * Looking at these stability profiles, they try to select a sufficiently large number of components to capture the information contained in the original data set, but small enough to avoid producing an excessive amount of highly unstable components. With their experiments with gene expression data sets, they noticed that an inflection point in the distribution of the stability profiles often offers such a satisfying trade-off.
+
+.. figure:: images/MSTD.png
+   :align: center
+   :scale: 50 %
+
+The module `sica.base.MSTD <modules/generated/sica.base.MSTD.html#sica.base.MSTD>`_ allows the user to plot the stability profiles for several order of decomposition. It offers an insightful visualization tool to select a relevant number of ICA components to extract, using the reasoning described above. 
+
+**Note 1:** There is no guarantee to obtain a distribution of stability profiles similar to the figure above or to those plotted in [3]. This very much depends on the data and the choice of hyperparamters for the ICA solver. However the philosophy of the method developped in [3] remains the same : the plot of this distribution should help you finding a satisfying compromise between many components and a high overall stability.
+
+**Note 2:** Many other methods could be used to infer the intrinsic dimension of a data set and selecting a relevant number of components before using **stabilized-ica**. We advise the user to explore the `scikit-dimension package <https://scikit-dimension.readthedocs.io/en/latest/>`_ . 
+
+.. topic:: References:
+
+    .. [1] : J. Himberg and A. Hyvarinen, "Icasso: software for investigating the reliability of ICA estimates by clustering and visualization," 2003 IEEE XIII Workshop on Neural Networks for Signal Processing (IEEE Cat. No.03TH8718), 2003, pp. 259-268, doi: 10.1109/NNSP.2003.1318025.
+    .. [2] : F. Meinecke, A. Ziehe, M. Kawanabe and K.-R. Müller, “Estimating the Reliability of ICA Projections,” in Advances in Neural Information Processing Systems 14, MIT Press, 2002.
+    .. [3] : Kairov U, Cantini L, Greco A, Molkenov A, Czerwinska U, Barillot E, Zinovyev A. Determining the optimal number of independent components for reproducible transcriptomic data analysis. BMC Genomics. 2017 Sep 11;18(1):712. doi: 10.1186/s12864-017-4112-9. PMID: 28893186; PMCID: PMC5594474.

docs/source/Icasso.rst

@@ -0,0 +1,41 @@
+Mutual Nearest Neighbors method
+===============================
+In many applications, one interesting question is: Are the stabilized ICA components extracted from one data set reproducible over several data sets ? For instance, working on a gene expression data set with samples associated to a specific disease, it may be worth wondering whether the extracted ICA components are reproducible for heatlhy samples or not.
+
+A meta-analysis with several data sets is required to address such a question. The idea would be to extract the stabilized ICA components for each data set, and then identify those that are reproducible over several data sets. 
+
+The Reciprocal Best Hits method
+-------------------------------
+One of the most natural method to do that in practice is called Reciprocal Best Hits (RBH). It was first developped to find orthologous genes among different species (i.e genes in different species that originated by vertical descent from a single gene of the last common ancestor). 
+
+The RBH method relies on a simple rule. For two components :math:`x_i^{1}` and :math:`x_j^{2}` extracted from data sets :math:`\mathcal{D}_1` and :math:`\mathcal{D}_2`, we consider that they are linked if :math:`x_j^{2}` is the nearest neighbor of :math:`x_i^{1}` among the points of :math:`\mathcal{D}_2` and :math:`x_i^{1}` is the nearest neighbor of :math:`x_j^{2}` among the points of :math:`\mathcal{D}_1`. Applying this rule for every component of each data set, we obtain a graph whose nodes represent the different components and edges highlight their reproducibility accross data sets.
+
+**Note :** The RBH method can deal with any metric. However for stabilized ICA components it seems quite natural to use a correlation metric (i.e :math:`1 - |\rho_{ij}|`) since they can be interpreted as finite samples from latent distributions. By default, `sica.mutualknn.MNNgraph <modules/generated/sica.mutualknn.MNNgraph.html#sica.mutualknn.MNNgraph>`_ uses the Pearson correlation coefficients.
+
+An illustration of the RBH method with gene expression data
+-----------------------------------------------------------
+The following figure shows a `cytoscape <https://cytoscape.org/>`_ visualization of a RBH graph obtained for seven different lung cancer gene expression data sets. First, we applied `sica.base.StabilizedICA <modules/generated/sica.base.StabilizedICA.html#sica.base.StabilizedICA>`_ on each data set. Then the RBH graph was computed with `sica.mutualknn.MNNgraph <modules/generated/sica.mutualknn.MNNgraph.html#sica.mutualknn.MNNgraph>`_ . Each node corresponds to an extracted stabilized ICA component and each color is associated to a specific data set. 
+
+.. figure:: images/cytoscape.png
+   :align: center
+   :scale: 70 %
+
+We observe, on the left of the graph, several tight, isolated and fully connected communities with one component per data set. Such a topology may indicate reproducible biological signals, present in all the data sets.
+
+A natural extension of RBH: the mutual nearest neighbors method
+---------------------------------------------------------------
+The mutual k-nearest neighbors approach offers a straightforward extension to the RBH method described above. Indeed, it simply consists in considering the reciprocal members of the k-nearest neighborhood rather than only the reciprocal nearest neighbors. The module **sica.mutualknn** allows the user to set the metric as well as the parameter k. When k=1, this is equivalent to the RBH approach.
+
+The following figure illustrates the building of a mutual nearest neighbors graph with k=2 and two different data sets A and B. An edge will be drawn between two components :math:`x_A` and :math:`x_B` from A and B if and only if :math:`x_A` belongs to the 2-nearest neighborhood of :math:`x_B` (only considering points of A) and vice versa.
+
+.. figure:: images/MNN.png
+   :align: center 
+   :scale: 70 %
+
+   Comparison of two situations for the building of the mutual k-nearest neighbors graph (k=2)
+
+**Note :** The mutual k-nearest neighbors approach is and extension of the classical k-nearest neighbors method. It seems to date back to the early eighties in the context of clustering [1]_. Nowadays, it is growing more and more popular and emerges as a performing alternative to k-nearest neighbors in many machine learning algorithms. In particular, it is claimed to reduce the effect of hubness in high dimension. However, to our knowledge, the theoretical properties of such a method are not yet fully understood. Thus, you shoud choose a metric and a value for k wisely and interpret the resulting graph with care.
+
+.. topic:: References:
+
+    .. [1] : K. Chidananda Gowda and G. Krishna. Agglomerative clustering using the concept of mutual nearest neighbourhood. Pattern Recognition, 10(2):105–112, 1978.