Subgroup detection in genotype data using invariant coordinate selection
 Daniel Fischer^{1}Email authorView ORCID ID profile,
 Mervi Honkatukia^{1},
 Maria TuiskulaHaavisto^{1},
 Klaus Nordhausen^{2, 4},
 David Cavero^{3},
 Rudolf Preisinger^{3} and
 Johanna Vilkki^{1}
DOI: 10.1186/s1285901715899
© The Author(s) 2017
Received: 9 August 2016
Accepted: 9 March 2017
Published: 16 March 2017
Abstract
Background
The current gold standard in dimension reduction methods for highthroughput genotype data is the Principle Component Analysis (PCA). The presence of PCA is so dominant, that other methods usually cannot be found in the analyst’s toolbox and hence are only rarely applied.
Results
We present a modern dimension reduction method called ’Invariant Coordinate Selection’ (ICS) and its application to highthroughput genotype data. The more commonly known Independent Component Analysis (ICA) is in this framework just a special case of ICS. We use ICS on both, a simulated and a real dataset to demonstrate first some deficiencies of PCA and how ICS is capable to recover the correct subgroups within the simulated data. Second, we apply the ICS method on a chicken dataset and also detect there two subgroups. These subgroups are then further investigated with respect to their genotype to provide further evidence of the biological relevance of the detected subgroup division. Further, we compare the performance of ICS also to five other popular dimension reduction methods.
Conclusion
The ICS method was able to detect subgroups in data where the PCA fails to detect anything. Hence, we promote the application of ICS to highthroughput genotype data in addition to the established PCA. Especially in statistical programming environments like e.g. R, its application does not add any computational burden to the analysis pipeline.
Keywords
ICS PCA Genotype data Classification Dimension reductionBackground
The fast progress in analyzing variations in the genome by deep sequencing has led to a plethora of high density genotyping arrays in many livestock species. Thereby also, the amount of single nucleotide polymorphism (SNP) data that is available for analyzing the genetic relationships between different samples is constantly growing. One common approach to handle this type of data and to identify e.g. subpopulations, is the application of dimension reduction methods such as Principle Component Analysis (PCA). Currently PCA is established to be the standard approach in clustering genotype data, see e.g. [1]. However, as we will demonstrate with an simulation example, there are drawbacks and pitfalls in the PCA approach. In a PCA, the principle components are ordered according to the variance they explain, but there is no theoretical justification that the component with the largest explained variance also contains the information required, e.g. to separate subgroups within the data. A vivid counterexample is a large hamburger. If it is large enough, the component with the largest variation goes through the diameter of the burger, but to separate the subgroups, one would need a direction from bottom to top. That means, even in this simple threedimensional case the interesting component would be only the second one. Hence, the interesting components might explain only a small fraction of the variance and consequently are easily missed by checking only the few first or last components. For an overview and a more theoretical background on the application of PCA in genotype data, see [2] or [3].
Other dimension reduction methods, such as Invariant Coordinate Selection (ICS), are not commonly applied to genomic data. ICS is a modern multivariate method originally introduced as generalized PCA in [4], but then established as ICS in the seminal paper [5] to avoid a name mismatch with a different generalized PCA approach, see e.g. [6].
The basic idea of ICS is to use two different scatter matrices and to compare how they differ. Different choices of scatter matrices lead then to different applications of ICS. Currently, ICS has been e.g. applied to nearreal time retrieval of low stratiform cloud coverage [7]. Further, ICS was used to enhance the discrimination between snow and ice clouds and detection of broken, thin clouds [8] and also for studies of developmental canalization and the identification of divergent and stabilizing selection [9].
To discuss possible problems with PCA, we will first present a basic counterexample to show that PCA does not necessarily identify cluster structures in a dataset and after that apply PCA and ICS on a real example genotype dataset. Further, for both datasets we will compare the two methods also with other methods used by the Bioinformatics community. For that we apply also tdistributed Stochastic Neighbor Embedding (tSNE) [10], Isomap [11], Locally Linear Embedding (LLE) [12], kernel PCA (kPCA) [13] and Diffusion Maps (DM) [14] to the simulated and the real data. For completeness, we also check the performance of a Linear Discriminant Analysis (LDA) for the simulation example.
Methods
Simulated data
First we simulated a dataset as an example that PCA is not always capable of detecting clusters in highdimensional data. Consider three 10variate normal populations with N _{10}(μ _{ i },Σ), where
μ _{1}=(−μ ^{∗},μ ^{∗},0,0,0,0,0,0,0,0)^{⊤},
μ _{2}=(μ ^{∗},0,0,0,0,0,0,0,0,0)^{⊤} and
μ _{3}=(0,−μ ^{∗},0,0,0,0,0,0,0,0)^{⊤} and
\(\boldsymbol {\Sigma }=\text {diag}\left (\frac {1}{5},\frac {1}{5},1,1,1,2,2,2,2,2\right)\) with μ ^{∗}=2. From each population we simulated then 100 samples. In order to hide the clearly visible subpopulations, we further rotated the simulated observations with a random orthogonal matrix. Note that the rotation has no impact on the performance of PCA as the method is rotation invariant.
Additional file 1: Figure S1 shows the simulated data before the rotation and Additional file 1: Figure S2 the data after it. In the latter one, the groups are clearly not visible anymore.
Chicken data
The highdensity genotype data consists of 749 chicken from 4 generations. The last generation is the largest group with 603 samples. The other generations contain 50, 46 and 50 samples. The data consists of sequence based variation data from 7 genomic regions, covering approx. 35% of the genome. The regions have been preselected based on previous studies as containing loci affecting eggquality traits, see [15] and [16]. As reference genome we used galgal4.
In total there were 157,528 genotypes measured in those regions. See Additional file 1: Figure S3 for the locations of the used regions on the chicken reference genome.
In addition to the genotype data, also a set of 15 different breeding values was available for all chicken. These were, besides others, egg production in period 3 to 7, egg production in period 9 to 12 and feed intake. We use this data as real data example and will follow up on the biological findings only for one detected subgroup in order to keep the focus on the method.
Invariant coordinate selection
An interpretation can then be given as follows. First S _{1} is used to whiten the data, i.e. uncorrelate the variables and standardize the scales. Then perform on the whitened data PCA using S _{2}. Therefore the idea is to see if S _{2} finds still some interesting structure after removing second order information as measured by S _{1}.
The univariate concept of kurtosis can be seen as the ratio of two (standardized) scale measures and similarly \(\mathbf {S}_{1}^{1}\mathbf {S}_{2}\) can hence be seen as a multivariate extension of this concept. Therefore the eigenvalues contained in D can be interpreted as generalized kurtosis values as measured by S _{1} and S _{2}. In the special case of S _{1}=COV and S _{2}=COV_{4} it can be shown that the diagonal elements in D are a linear function of the classical measures of kurtosis of the components in z [18].
And for example when searching clusters it is wellknown that large clusters can be found often in directions with small kurtosis and outliers and small clusters in directions with large kurtosis. This means that invariant coordinates are very suitable for searching for groups as the components are ordered according to their (generalized) kurtosis. As actually [5] show, in the context of mixtures of elliptical distributions with proportional scatter matrices, ICS finds Fisher’s linear discriminant subspace without knowing the group memberships. Hence, when using ICS for exploratory data analysis usually most attention is paid to the components with extreme generalized kurtosis values, like for example the first 3–5 and last 3–5 components. For more details about ICS see [4, 5, 18, 19].
As practical considerations we would however like to point out that there is no general best combination of scatter matrices and the performance might depend on the choice of S _{1} and S _{2}. The choice S _{1}=COV and S _{2}=COV_{4} is however wellestablished and for example also a solution to the independent component problem (ICA) if x follows it (see eg [20] for further details). ICA has been applied in the context of genetic data e.g. in [21].
Furthermore, ICS is however currently limited to the case when p<n−1 as otherwise scatter matrices are always proportional to each other, see [22] for details. Therefore if p≥n−1, then one can for example first perform dimension reduction using PCA, resulting in a n×n matrix where the nth eigenvalue is zero. Then ICS is only applied to a subspace which is known to have variation and is of smaller dimension than n−1. This is for example standard practice in many multivariate methods which are limited to the p≥n−1 case, like for example the highdimensional noisy ICA approaches [23].
Distance measure, distance groups and statistical testing
For the simulated data the classification decision based on the scatterplot matrices from PCA and ICS was done by applying a kmeans algorithm to the desired components. The classification results of the different dimension reduction methods were then evaluated using the adjusted Rand index [24]. In the real data example, the classification decision was done by visual inspection of the figures.
In order to calculate the genetic distance of two different groups in a region of interest, we followed a basic approach. Assuming two subpopulations A and B have been identified in the data, we determined first at each loci l=1,2,… the most common genotypes for both groups and denote these \(\tilde {G}_{A,l}\) respective \(\tilde {G}_{B,l}\). Then, we compared if these genotypes match between the two groups, by setting G _{ l }=1, if \(\tilde {G}_{A,l} = \tilde {G}_{B,l}\) and 0 else. Afterwards we calculated a moving average of length 1000 across the data and calculated in each window the average level of agreement. Let W=w _{1},w _{2},… be the set of all windows of length 1000 with w _{1}=l _{1},…,l _{1000},w _{2}=l _{2},…,l _{1001}, the average level of agreement in window i is then \(\bar {x}_{w_{i}}=\sum _{\forall l \in w_{i}} G_{l} / 1000\). For the sake of simplicity, we calculated the moving average also across chromosomal borders.
For all windows w _{ i } with level of agreement between two subpopulations \(\bar {x}_{w_{i}} \leq 0.4\), the individual distance of each individuum in the one group was calculated to the average of the other group. For that, we use again the most common genotype for each loci in the subpopulation coded as 0,1,2 and then we calculated the Manhattan distance of each individuum from the standard population to that.
Results
To evaluate the performance of the different dimension reduction methods to unravel the original cluster structure, we first clustered the plain simulated data using kmeans with the constrain of three classes (k=3). For the classification result, we then calculated the adjusted Rand Index for the 3×3 table between the original class labels and the result of the kmeans clustering. Next, we performed a PCA, followed again by a kmeans clustering using the first two components for classification. Also for this classification result table we calculated the adjusted Rand index. Then we applied ICS onto the dataset and calculated in the same way again the Rand index for the kmeans applied to the last two components. To compare the results to other popular dimension reduction methods, we applied also tSNE, Isomap, LLE, kPCA and DM to the simulation data and calculated the corresponding Rand indices. Further, we searched with each dimension reduction method the same dimensionality, that was d=2 for the simulated and d=7 for the real chicken data.
The results for the other dimension reduction methods were rather weak. Whereas the tSNE method was almost as good as the ICS (Randindex 0.93), the four others clearly were outperformed by these two mothods. Isomap had a Rand index of 0.71, LLE had a value of 0.48, DM had also only 0.50 and the kPCA method had with 0.42 even a value smaller than the PCA had. That means, none of these methods was able to fully recover the original data. The corresponding Figures S4–S11 can be found in the Additional file 1. To calculate the tSNE we used the RPackage tsne [27], Isomap is implemented in the RPackage RDRToolbox [28] and LLE in lle [29]. For kPCA we used the kernmap package [30] and for DM the destiny package [31].
The lda function applied to the simulation data resulted in an errorfree separation of the data and had consequently a Rand index of 1. However, the ICS method is with 0.94 not too far away from that optimum. In absolute numbers, 6 out of 300 observations were mislabeled using the ICS function. LDA cannot be applied to the real example data, as the identification of subgroups is done without any prior knowledge and as such supervised methods like LDA cannot be applied to the problem.
Before analyzing the phenotypical particularities of the identified subgroup, we also test the performance of the other dimension reduction methods on the real chicken data. Here, kPCA and LLE are able to identify the same clusters as ICS does, but tSNE and Isomap fail to identify any clear cluster structures. In case of tSNE we tried both k=7 and k=2, but in neither case any obvious subgroup could be identifed. Diffusion map, however, apparently identifies another subgroup. The corresponding scatterplot matrices can also be found in the Additional file 1. We used the default settings and protocols as provided by the different packages. That means, e.g. for LLE we calculated the optimal number of neighbors as 17.
Members of the red subgroup, identified by the ICS method are all offsprings from the same father and mainly from the same mother. From the 20 members of the subgroup only one individual (indicated by green) has a different mother. The subpopulation indicated in blue is also formed by a family. Seven chicken from this population have the same father and mother. Further, the father of those 7 chicken can also be found in this group.
A region of approximate length 4Mb (Chr2:70,348,41374,448,870), containing 1340 SNPs was identified by calculating the genetic similarity between the deviating (red) family and the remaining population. The genetic similarity was calculated with a moving average using windows of size 1 kb. Areas, where the average level of agreement drops below 0.4 are considered to be the major cause for the difference between the divergent red family and the main population. Additional file 1: Figure S12 shows the level of agreement across the considered chromosomal regions. Also for the blue subpopulation we could identify in a candidate region a similar way.
Next, we calculated for each chicken within the main population the Manhattan distance between the mode genotype values of the deviating red family in the region of interest and the individual genotypes. There we could clearly identify three subgroups within the main population, see Additional file 1: Figure S13. We denote those subgroups as close, intermediate and far.
When breeding values of 15 production values were compared between the red subpopulation and the main population, significant differences were seen in 10 traits (The twosided MannWhitney test was significant at level α=0.05). These were then tested further using a generalized MannWhitney test for directional alternatives. This means, we tested for a directional trend of the phenotypes with respect to the close, the intermediate and the far group.
Discussion
We applied the modern dimension reduction method ICS to a simulation example and compared it to the commonly used PCA method to visualize some deficiency of the PCA approach. Further, we applied the other, modern dimension reduction methods tSNE, Isomap, LLE, kPCA and DM to the simulation data. Here, in the controlled environment we could clearly see that the PCA method was not able to identify all three true groups in the simulated data, but the ICS method, however, was. From the other tested methods, only tSNE was able to recover all three subgroups, but all other four tested methods failed doing so. Some of them separated a single subgroup, but mixed the remaining two groups into a single large cluster. When the methods were then applied to a highdensity genotype chicken data, the PCA method could not identify any subgroups. The ICS method clearly identified two subgroups consisting of 20, respective 10 samples, that share the same family background. Two (kPCA and LLE) of the five other methods, however, also detected the same subgroups in the real chicken data. The other three methods failed to identify any clear cluster structures.
In the scatterplot matrices some outlying observations could be identified by tSNE (see Additional file 1), but not as evident as in the ICS case. A closer look at component 3 showed e.g., that some of the chicken with a value larger 25 are related but the most of them are unrelated. In terms of calculation times, the ICS needed around 0.2 s, whereas the tSNE run took around three minutes. The other used methods needed at most only a few seconds for the calculation.
We considered also the red subgroup identified by ICS closer. It was superior in more than half of the available breeding values compared to the standard chicken population. Within the standard chicken population we could identify three subgroups that were either genetically close, intermediate or far away from this subgroup based on the most deviating chromosomal region. In addition, these three groups of the main population showed a directional trend in many traits, especially in the important production values P2 and P3. Also the blue subgroup is deviating in five breeding values from the main population, including the production values P3.
There were no other combinations with those parents in the data available so that no further investigations could be conducted to identify the reason for the subgroups to behave in such a different way. The biological explanation for the difference is beyond the scope of this paper.
The identification of three groups (main group and two subgroups) within the data is remarkable. As all the chicken originate from the same line, one would not assume any subpopulation structures and by applying a PCA to it, we did not identify any. ICS identified two subpopulations that were thereafter also seen to differ from the main population for some of the phenotypes for production traits.
Further, we could identify strongly deviating genetic regions between the subpopulations and the main group and followed exemplary up on the one that corresponds to the red subgroup. Within that, we calculated then the genetic distance of the remaining chicken to the identified subpopulation and could see that chicken genetically more similar with regard to the deviating region to the subpopulation also have better production values. Moreover, we could identify a directional relationship between the genetic similarity in that region and certain production values.
Conclusion
We presented here an alternative dimension reduction method that is already used in other scientific fields, but that has not yet made its way to the genomic community. However, although ICS is superior over PCA in the current scenario, its purpose is not to replace PCA or any other dimension reduction method, but it is rather considered to be another tool in the dense genotype data analysis toolbox.
Its good results for both, the simulation and the real dataset encourage its use also for other genomic datasets to further evaluate its performance in a larger scale. Compared to other, modern dimension reduction methods, we saw that there is a large variation in the performance of each method, depending on the dataset. For our data, only ICS showed good results in the simulation as well as in the real data set, Isomap and Diffusion map had the weakest results for both setups. tSNE only performed well in the simulation setup and LLE only for the real data.
Abbreviations
 DM:

Diffusion maps
 ICA:

Independent component analysis
 ICS:

Invariant coordinate selection
 kPCA:

kernel PCA
 LDA:

Linear discriminant analysis
 LLE:

Locally linear embedding
 PCA:

Principle component analysis
 SNP:

Single nucleotide polymorphism
 tSNE:

tDistributed stochastic neighbor embedding
Declarations
Acknowledgements
The authors would like to thank everyone at Steffen Weigend’s group at the Friedrich Loeffler Institute, Greifswald, Germany, for preparing the raw genotype data. Further, we would like to thank the anonymous reviewers for their valuable comments.
Funding
This research was financially supported by Lohmann Tierzucht GmbH, the design analysis and conclusions of the study were done by the authors, and the material was received from the Lohmann Tierzucht breeding programme.
Availability of data and materials
Data is not publicly shared, as it is only used as an example and no direct results are derived from it. The data is, however, a subset of a larger study that will be published later.
Authors’ contributions
DF, KN and JV wrote the manuscript. DF and KN implemented the Rscripts. MH, MTH and JV provided the genetical background and interpretation. DC and RP provided and processed the used data. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Notapplicable.
Ethics approval and consent to participate
Chicken lines are owned bythe breeding company, Lohmann Tierzucht (LTZ). The progeny of these birds are used for intracommunity trade and exports. EU Directive 158/2009 requires the examination of blood samples for certain pathogens. Further testing is required for multiple health certificates for third countries. Blood samples obtained from this routine health monitoring, are used to extract DNA to genotype animals. The mentioned DNA were extracted from routine blood samples in the inhouse laboratory at LTZ.
All birds are kept and supervised in strict compliance of the animal welfare laws of Lower Saxony, Germany. All staff taking the blood samples regularly participate in training programs for adecuate sampling in strict observation of current animal welfare regulations.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Solovieff N, Hartley SW, Baldwin CT, Perls TT, Steinberg MH, Sebastiani P. Clustering by genetic ancestry using genomewide snp data. BMC Genet. 2010; 11. doi:10.1186/1471215611108.
 Patterson N, Price AL, Reich D. Population structure and eigenanalysis. PLoS Genet. 2006; 2:190. doi:10.1371/journal.pgen.0020190.View ArticleGoogle Scholar
 Ma S, Dai Y. Principal component analysis based methods in bioinformatics studies. Brief Bioinforma. 2010; 12:714–22. doi:10.1093/bib/bbq090.View ArticleGoogle Scholar
 Caussinus H, Ruiz A. Interesting Projections of Multidimensional Data by Means of Generalized Principal Component Analyses In: Momirović K, Mildner V, editors. Compstat: Proceedings in Computational Statistics, 9th Symposium held at Dubrovnik, Yugoslavia, 1990. Heidelberg: PhysicaVerlag HD: 1990. p. 121–6. doi:10.1007/9783642500961_19.
 Tyler DE, Critchley F, Dümbgen L, Oja H. Invariant coordinate selection. J R Stat Soc Series B. 2009; 71:549–92. doi:10.1111/j.14679868.2009.00706.x.View ArticleGoogle Scholar
 Vidal R, Ma Y, Sastry SS. Generalized Principal Component Analysis. New York: Springer; 2016.View ArticleGoogle Scholar
 Musial JP, Hüsler F, Sütterlin M, Neuhaus C, Wunderle S. Daytime low stratiform cloud detection on avhrr imagery. Remote Sensing. 2014; 6(6):5124. doi:10.3390/rs6065124.View ArticleGoogle Scholar
 Musial JP, Hüsler F, Sütterlin M, Neuhaus C, Wunderle S. Probabilistic approach to cloud and snow detection on advanced very high resolution radiometer (avhrr) imagery. Atmos Meas Tech. 2014; 7(3):799–822. doi:10.5194/amt77992014.View ArticleGoogle Scholar
 Bookstein FL, Mitteroecker P. Comparing covariance matrices by relative eigenanalysis, with applications to organismal biology. Evol Biol. 2013; 41(2):336–50. doi:10.1007/s1169201392605.View ArticleGoogle Scholar
 van der Maaten LJP, Hinton GE. Visualizing highdimensional data using tsne. J Mach Learn Res. 2008; 9:2579–605.Google Scholar
 Tenenbaum JB, de Silva V, Langford JC. A global geometric framework for nonlinear dimensionality reduction. Science. 2000; 290:2319–23.View ArticlePubMedGoogle Scholar
 Roweis S, Saul L. Nonlinear dimensionality reduction by locally linear embedding. Science. 2000; 290:2323–6.View ArticlePubMedGoogle Scholar
 Schölkopf B, Smola A, Müller KR. Nonlinear component analysis as a kernel eigenvalue problem. Neural Comput. 1998; 10:1299–319.View ArticleGoogle Scholar
 Coifman RR, Lafon S, Lee AB, Maggioni M, Nadler B, Warner F, Zucker SW. Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps. PNAS. 2005; 102:7426–31.View ArticlePubMedPubMed CentralGoogle Scholar
 TuiskulaHaavisto M, Honkatukia M, Preisinger R, Schmutz M, de Koning DJ, Wei WH, Vilkki J. Quantitative trait loci affecting eggshell traits in an f2 population. Animal Genet. 2011; 42:293–9.View ArticlePubMedGoogle Scholar
 Honkatukia M, TuiskulaHaavisto M, Arango J, Tabell J, Schmutz M, Preisinger R, Vilkki J. Qtl mapping of egg albumen quality in egg layers. Genet Sel Evol. 2013; 45:31.View ArticlePubMedPubMed CentralGoogle Scholar
 Nordhausen K, Tyler DE. A cautionary note on robust covariance plugin methods. Biometrika. 2015. doi:10.1093/biomet/asv022.
 Nordhausen K, Oja H, Ollila E. Multivariate Models and the First Four Moments. Singapore: World Scientific; 2011, pp. 267–87. doi:10.1142/9789814340564_0016.
 Nordhausen K, Oja H, Tyler DE. Tools for exploring multivariate data: The package ICS. J Stat Softw. 2008; 28(6):1–31. doi:10.18637/jss.v028.i06.View ArticleGoogle Scholar
 Miettinen J, Taskinen S, Nordhausen K, Oja H. Fourth moments and independent component analysis. Statist Sci. 2015; 30(3):372–90. doi:10.1214/15STS520.View ArticleGoogle Scholar
 Tapio M, Tapio I, Grislis Z, Holm LE, Jeppsson S, Kantanen J, Miceikiene I, Olsaker I, Viinalass H, Eythorsdottir E. Native breeds demonstrate high contributions to the molecular variation in northern european sheep. Mol Ecol. 2005; 14(13):3951–63. doi:10.1111/j.1365294X.2005.02727.x.View ArticlePubMedGoogle Scholar
 Tyler DE. A note on multivariate location and scatter statistics for sparse data sets. Stat Probab Lett. 2010; 80(17–18):1409–13. doi:10.1016/j.spl.2010.05.006.View ArticleGoogle Scholar
 Oja H, Nordhausen K. Independent Component Analysis In: ElShaarawi AH, Piegorsch W, editors. Encyclopedia of Environmetrics. New Jersey: John Wiley & Sons: 2012. p. 1352–1360.Google Scholar
 Rand WM. Objective criteria for the evaluation of clustering methods. J Am Stat Assoc. 1971; 66(336):846–50.View ArticleGoogle Scholar
 Fischer D, Oja H, Schleutker J, Sen PK, Wahlfors T. Generalized MannWhitney type tests for microarray experiments. Scand J Stat. 2014; 41:672–92. doi:10.1111/sjos.12055.View ArticleGoogle Scholar
 Fischer D, Oja H. MannWhitney type tests for microarray experiments: The R package gMWT. J Stat Softw. 2015; 65(1):1–19. doi:10.18637/jss.v065.i09.
 Donaldson J. Tsne: TDistributed Stochastic Neighbor Embedding for R (tSNE). 2016. R package version 0.13. http://CRAN.Rproject.org/package=tsne. Accessed 30 Nov 2016.
 Bartenhagen C. RDRToolbox: A Package for Nonlinear Dimension Reduction with Isomap and LLE. 2014. R package version 1.20.0. https://www.bioconductor.org/packages/release/bioc/html/RDRToolbox.html. Accessed 30 Nov 2016.
 Diedrich H, Abel M. Lle: Locally Linear Embedding. 2012. R package version 1.1. http://CRAN.Rproject.org/package=lle. Accessed 30 Nov 2016.
 Karatzoglou A, Smola A, Hornik K, Zeileis A. kernlab – an S4 package for kernel methods in R. J Stat Softw. 2004; 11(9):1–20.View ArticleGoogle Scholar
 Angerer P, Haghverdi L, Büttner M, Theis FJ, Marr C, Buettner F. destiny – diffusion maps for largescale singlecell data in R. Bioinformatics. 2015. doi:10.1093/bioinformatics/btv715. http://bioinformatics.oxfordjournals.org/content/early/2015/12/13/bioinformatics.btv715.full.pdf+html.
 Zheng X, Levine D, Shen J, Gogarten S, Laurie C, Weir B. A highperformance computing toolset for relatedness and principal component analysis of snp data. Bioinformatics. 2012; 28:3326–8. doi:10.1093/bioinformatics/bts606.View ArticlePubMedPubMed CentralGoogle Scholar