7. OTU Methods in Numerical Taxonomy

7.3.1.1 Aristotle’s Metaphysics of Essentialism

Essentialists believe true essences exist. Aristotle defined essence as “the what it was to be” for a thing (Aristotle 2016) or sometimes used the shorter term “the what it is” for approximating the same idea. Aristotle defined essence as analyzed entities. For example, in his logical works (Aristotle 2012), Aristotle links essence to definition—“a definition is an account (logos) that signifies an essence” or “there will be an essence only of those things whose logos is a definition,” “the essence of each thing is what it is said to be intrinsically.” Aristotle’s essence is related to a thing. In his classic definition of essence, a thing’s essence is that which it is said to be per se, or is the property without which the thing would not exist as itself. Therefore, essence is constitutive of a thing, which is eternal, universal, most irreducible, unchanging. The classical example of essence definition is that Aristotle defined the human essence: to be human is to be a “rational animal.” The “rational” is the species or specific difference, which separates humans from everything else in the class of animals.

Another important essentialist concept is forms. For Plato, forms are abstract objects, existing completely outside space and time. Aristotle criticized and rejected Plato’s theory of forms but not the notion of form itself. For Aristotle, forms inform the thing’s particular matter. Thus, forms do not exist independently of things—every form is the form of some things. He called “the form” of a particular thing (Albritton 1957).

The naturalists in seventeenth and eighteenth centuries who attempted to classify living organisms were actually inspired by the terms and definitions from the Aristotelian system of logic. Darwin traced evolutionary ideas (the principle of natural selection) back to Aristotle in his later editions of the book of The Origin of Species (Darwin 1872), and commented on Aristotle on forms.

Aristotelian system attempted to discover and define the essence of a taxonomic group (“real nature” or “what makes the thing what it is”). Aristotelian logic is a system of analyzed entities rather than a system of unanalyzed entities.

7.3.1.2 Natural Classification: Numerical-Based Empiricism

Numerical taxonomy is a largely empirical science (Sokal 1985; Sneath 1995).

Natural classification was developed based on the philosophy of empiricism and was affected by the philosophy of the inductive sciences. Some key concepts of information content and predictivity provided by John Gilmour (Sneath 1995) who was one of the early contributors to the founding of numerical taxonomy were adapted from the Victorian philosophers of science John Stuart Mill (1886) and William Whewell (1840).

Sokal and Sneath criticized that Aristotelian logic cannot lead to biological taxonomy because the properties of taxonomy cannot be inferred from the definitions [see (Sokal and Sneath 1963, p. 12; Sokal 1962; Sneath and Sokal 1973, p. 20)]. Thus, they emphasized that natural taxa do not necessarily possess any single specified feature, which is the Aristotelian concept of an essence of a taxon, because natural or “general” classification or natural groups are in logic unanalyzed entities (Sokal and Sneath 1963, pp. 19–20).

Sneath and Sneath (1962) differentiated the terms “polythetic” and “monothetic” in taxonomic groups. The monothetic groups define a unique set of features that is both sufficient and necessary for membership in the group thus defined. The monothetic groups have the risk to misclassify natural phenetic groups and hence do not yield “natural” taxa. In contrast, in a polythetic group, organisms that have the greatest number of shared character states are placed together, and hence no single state is either essential to group membership or is sufficient to make an organism a member of the group.

For Sneath and Sokal, natural taxa are usually not fully polythetic, and it is possible that they are never fully polythetic; however, in practice, they suggest considering the possibility of a taxon being fully polythetic due to limited knowledge of whether any characters are common to all members (Sneath and Sokal 1973, pp. 20–21).

Numerical taxonomists emphasize that biological taxonomy is a system of unanalyzed entities. Their properties cannot be inferred from the definitions, instead through classification. In his “Taxonomy and Philosophy” (Gilmour 1940), Gilmour stated that classification enables us to make generalizations and predictions. The more information obtained by classification enables us to make more predictions. The more predictions the classification is more “natural.”

For Gilmour, the most “natural” system/classification contains the most information, which is highly predictive and serves the most general purpose (Vernon 1988), while less “natural” classification is very useful for specific purposes. Generally Sokal and Sneath held with Gilmour’s review of “natural” classification. However, classification serves as special purpose. Such classification is a special classification and is often “arbitrary.” Thus Gilmour’s classification conveys less information than a general or “natural” one. For Sokal and Sneath, a natural taxonomic group or natural taxa of any rank can be classified. For example, the genetic natural classification can be supported (Vernon 1988). They defined a natural classification in phenetics as one in which the members of each taxon at each level are on the average more similar to each other than they are to members of other taxa at corresponding levels (Sokal 1986).

The philosophical attitude of numerical taxonomy in systematics is empiricism and consequently is not committed to the existence of biological species (Sokal and Crovello 1970). Put together, numerical taxonomy is both operational and empirical. Numerical-based natural classification is an empiricism. Empiricism in taxonomy is the only reasonable approach for arranging organized nature (Sokal 1963, p. 55). Sneath and Sokal (1973, p. 362) emphasized that numerical taxonomy is an empiricism and operationism and thought that the biological species definition was difficult to apply and in practice was nonoperational.

7.4.2.1 Data Matrix/OTU Table

In numerical taxonomy data matrix (or OTU table) (Table 7.1), rows represent characters (attributes/descriptors/variables), while columns represent OTUs (taxa/analysis units/observations).

The setting of ecological data matrix is different from that in a numerical taxonomic matrix, in which rows represent objects (subjects/stands/quadrats), while columns represent descriptors (species). Although data matrices in ecology and numerical taxonomy have different structures (Legendre and Legendre 2012, p. 60), however, knowing the distinctions between attributes and individuals is very important because it is often confused between attributes and individuals in ecology. In ecology, attributes are species and individuals are objects (Table 7.2).

The structures of microbiome abundance data matrix (or OUT/ASV table) are relatively simple. It is a hybrid of numerical data matrix and numerical taxonomic matrix, in which rows represent OTUs (ASVs/sub-OTUs/variables), while columns represent samples (subjects) or reverse (Table 7.3). Here samples are analysis units; OTUs or species are variables. We do not need to distinct attributes and individuals. However, when we use microbiome software, we need to specify which data matrix structure is used. For example, the phyloseq package has been developed based on ecological software vegan package. When creating data object, it asks the user to specify whether taxa are_rows = TRUE or FALSE, which is confirming the structure of data matrix: whether rows represent taxa/OTUs/ASVs. In its initial version, the phyloseq package even used species to stand for taxa.

7.4.2.2 R-Technique and Q-Technique

Data matrices in Tables 7.1, 7.2 and 7.3 can be examined using two kinds of techniques (Cattell 1952): the R-technique and the Q-technique (originated by Stephenson under the name of inverted factor technique (Stephenson 1936, 1952, 1953)). In general, in an association matrix, a so-called Q matrix is formed to compare all pairs of individuals, while a so-called R matrix (dispersion or correlation matrix) is formed for comparisons between the variates. Numerical taxonomy starts logically with the choice of the t entities to be classified (OTUs) and n attributes or variables (generally characters in terms of biological work), which gives the n x t matrix. The R-technique examines the association of pairs of characters (rows) over all OTUs (columns), while the Q-technique, in contrast, examines the association of pairs of OTUs (columns) over all characters (rows). In numerical taxonomy, the R-technique refers to correlations among characters based on OTUs, considering the correlation of pairs of characters present in a population, while the Q-technique refers to the qualifications of relations between pairs of taxa, frequently species, based on a preferably large number of characters (Sokal 1963, pp. 125–126), considering the correlations of pairs of individuals present in a population in terms of their characters. In other words, Q-technique refers to the study of similarity between pairs of OTUs, and R-technique refers to the study of similarity between pairs of characters (Sneath and Sokal 1973, p. 256).

The adoption of using the R-technique and the Q-technique in numerical taxonomy has been criticized at the beginning after the publication of Sokal and Sneath’s 1963 book. The duality of the R or Q problem is one of some fundamental (Williams and Dale 1966) and statistical (Sneath 1967) problems in numerical taxonomy. For example, factor analysis compares characters but finally results in a grouping of OTUs. To overcome the confusion over the use of Q and R technique, Williams and Dale (1966) suggested that (1) Q and R should refer to the purpose of the analysis, i.e., refer to the matrix, and (2) A-space and I-space rather than R-space or Q-space should be used that is operated upon when the relationships are represented in a hyperspace.

The organisms can be represented as points in a multivariate character space (A-space, or attribute space), and conversely the characters can be represented in an I-space (individual space) (Sneath 1967). A-space and I-space refer to the model: the former is used for a model in attribute-space and the latter is used for a model in the individual-space. A matrix of inter-individual distances implies relationships between points in an A-space, while an inter-individual correlation matrix implies angles in an I-space. They both are Q-technique. An attribute-correlation matrix is R-technique which concerns with angles, while an individual-distance matrix is Q-technique which concerns with point.

Sneath (1967; Sneath and Sokal 1973) accepted the two suggestions of Williams and Dale (1966). In numerical taxonomy, R-technique leads to a classification of characters, Q-technique to one of OTUs. But they emphasized formally both R-technique and Q-technique have the same main mathematical steps because by transposing the data matrix in an R study, the characters (rows) become the individuals comparable to the former OTUs, and the actual OTUs or taxa (columns) become the characteristics (attributes) over which association is computed (Sneath and Sokal 1973, p. 116). In either space the character values are the coordinates of the points in the hyperspace. Based on the distinction between A-space and I-space in Williams and Dale (1966), in numerical taxonomy, A-space (attribute space) has formally n dimensions (one for each attribute or character) with t points that represent the OTUs. I-space (individual space) has formally t dimensions (one for each OTU) with n points representing the attributes or characters (Sneath and Sokal 1973) although less than n (in A-space) or t (in I-space) dimensions respectively could be after attempting the dimensionality reduction (p. 116).

Cattell (1966) restricted Q- and R-techniques to factor analysis, and applied Q-technique as a customary practice to numerical taxonomy in terms of cluster analysis of OTUs. Both Q- and R-techniques are employed in cluster as well as factor analysis although numerical taxonomy mainly emphasized the Q-technique. By using the notations of A-space and I-space, most numerical taxonomic methods are operated by Q-technique on an A-space, but some cluster and factor analyses are Q-techniques operating on an I-space.

Overall, Sneath and Sokal emphasized that Q- and R-techniques have been adapted in the general sense in numerical taxonomy. Conveniently a Q-study is used for the quantification of relations between organisms according to their characters aiming to produce classifications of organisms; while an R-study refers to classifying characters into groups according to the organisms that possess them, in which the characters become in effect the OTUs, and the procedures are similar to Q-study, but operate on the transpose of the n x t matrix [See (Sneath 1967; Sneath and Sokal 1973, pp. 115–116)].

Basically, the uses of Q- and R-techniques (or analyses) in numerical taxonomy have been adopted in numerical ecology (Legendre and Legendre 2012, pp. 266–267), whereas Q measures (Q mode study) are used for the association coefficients (the relationships) among objects given the set of descriptors and R measures (R mode analysis) are used for the association coefficients (the relationships) among descriptors (variables) for the set of objects in the study.

In summary, Q analyses refer to the studies that start with computing an association matrix among objects, whereas R analyses refer to the studies that start with computing an association matrix among descriptors (variables). Numerical taxonomy defines a character as a variable (feature) and OTU as an analysis unit, which leads to the different understanding of R-technique and Q-technique in some other taxonomists and ecology. Later microbiome study adopted the methods of numerical taxonomy, but not the definition of character, where samples are observations and OTUs are variables (features/characteristics). The way data matrix is used in microbiome study is consistent with ecology. In microbiome study, Q- and R-techniques (or analyses) are not explicitly differentiated but in some software that are adopted from ecology the terms and according data matrices are often used. In terms of R-technique for relationship among variables and Q-technique for relationship among individuals, the microbiome study mainly uses R-technique to find the relationship among variables (OTUs/ASVs/taxa) over the samples.

7.4.2.3 Similarity Coefficients

In numerical taxonomy, the term “similarity coefficient” is most commonly used, although in the strict sense the term “coefficient of resemblance” would be more appropriate. The term “coefficients of similarity” or “coefficients of dissimilarity” is also used. Numerical taxonomy roughly groups four types of coefficients—distance coefficients, association coefficients, correlation coefficients, and probabilistic similarity coefficients. They collectively refer to as coefficients of resemblance or similarity and are used to compute resemblances between taxa.

The space that measures similarity does not need to be Euclidean in the strict sense, but it should be metric, in other words, determined by a metric function (Williams and Dale 1966; Johnson 1970).

Based on the requirements for similarity coefficients, Sneath and Sokal (1973) described the four axioms over the entire set of OTUs that the measures of dissimilarity or functions (converting similarity coefficients into measures of dissimilarity) should be satisfied (p. 120). This is very important when we apply these similarity coefficients or develop the new measures of similarity. However, in microbiome literature, this topic has been rarely mentioned. Thus, here we briefly describe these four axioms before we formally introduce similarity coefficients.

All these four axioms are true for the Euclidean distances (which additionally obeys Pythagoras’ theorem: a² + b² = c², i.e., the lengths of the square whose side is the hypotenuse (the side opposite the right angle) is equal to the sum of the lengths of the squares on the triangle’s other two sides), and thus they are metrics. However, not all measures of dissimilarity are metrics in numerical taxonomy. For example, there exist the cases that the difference in character states between two OTUs a and b is known; however, the function between a and b may still be zero. Thus, the axiom 4 is not fulfilled for all OTUs in the taxon (which are called as pseudometric or semimetric systems) (Sneath and Sokal 1973, p. 121). For ease of understanding of the relationships, numerical taxonomy mainly focuses on the coefficients whose properties lead to metric spaces.

Here the pair-function is called an ultrametric pair-function for any pair of OTUs. A phenogram is one example of this function (Sneath and Sokal 1973, p. 121).

Distance Coefficients

They are used to measure the distance between OTUs in a space defined in various ways, including Euclidean distance in a character space of one or more dimensions and non-Euclidean metrices. In the numerical taxonomy, the two commonly used distances are Manhattan distance (mean character difference) (Cain and Harrison 1958) and Euclidean distance (taxonomic distance) (Sokal 1961). Distance coefficients are the converse (or complement) of similarity coefficients. The distance is a measure that is complement of similarity, which are, in fact, measures of dissimilarity. Actually, we can consider complementary functions of association or correlation coefficients as measures of dissimilarity in an analogous way. It was shown that the maximum distance possible between OTUs would be in two-dimensional model, would be in three-dimensional model, and in an n-space, a hyperspace of n dimensions (Sokal 1963, pp. 143–144). Table 7.4 summarizes some measures of taxonomic distance compared by Sokal and Sneath.

Table 7.4

Selected measures of taxonomic distance in numerical taxonomy

Definition	Reference and comments
The mean character difference of (Cain and Harrison 1958):	It is defined as the absolute (positive) values of the differences between the OTUs for each character It is a special case of Manhattan or city-block metric It has the advantages of its simplicity and it is a metric It also has some disadvantages, such as always underestimating the Euclidean distance between the taxa in space, lacking of the desirable attributes of the taxonomic distance or its square (Sneath and Sokal 1973)
The taxonomic distance of (Sokal 1961): And the average distance is defined as:	The taxonomic distance between two OTUs j and k is defined in (Sokal 1961; Rohlf 1965; Sokal 1963) This is a generalized version of the Euclidean distance between two points in an n-dimensional space Like to M. C. D., it is a special case of Manhattan or city-block metric It was first used by Heincke as early as 1898 (Sokal and Sneath 1963) was used in (Sokal 1961) and suggested to standardize character state codes for each character The expected value of d very closely approaches and the expected variance of d ≈ 1/n which approaches zero as n tends to infinity (Sneath and Sokal 1973)
The coefficient of diverge of (Clark 1952):	The ratio varies between zero and unity Originally each OTU represents a number of specimens, while single values are used in numerical taxonomy (Sokal 1963)
The coefficient of racial likeness (Sokal and Sneath 1963):	This formula was developed based on Karl Pearson (1926)’s coefficient of racial likeness (Pearson 1926), in which the rows of the Q-type matrix are standardized The OTUs are represented by only single values (with a variance of one) rather than means of the ith character for entity j (Sokal 1963) The strengths of this coefficient (Sneath and Sokal 1973, p. 185) are (1) taking the variance of the estimates into consideration and (2) permitting estimates of the degree of misclassification of individuals in two populations
The generalized distance developed by Mahalanobis (1936) and Rao (1948): , Where W−1 is the inverse of the pooled variance-covariance (dispersion) matrix within samples with dimension n × n δJK is a vector of differences between means of samples J and K for all characters	The generalized distance is equivalent to distances between mean discriminant values in a generalized discriminant function, and therefore measuring the distance as a function of the overlap between pairs of populations and transforming the original distances so as to maximize the power of discrimination between individual specimens in the constructed groups (i.e., the OTUs) (Sneath and Sokal 1973, p. 128) It is computed by first obtaining the maximal variance between pairs of groups relative to the pooled variance within groups and then maximizing the difference between pairs of means of characters (Sneath and Sokal 1973, p. 127) The generalized distance has the same strengths as the coefficient of racial likeness (Sneath and Sokal 1973, p. 185): (1) accounting for the variance of the estimates into consideration and (2) permitting estimates of the degree of misclassification of individuals in two populations It is particularly useful in numerical taxonomy (Sneath and Sokal 1973, p. 127) However, it also has some weaknesses (Sneath and Sokal 1973, p. 128): (1) it may distort the original distances considerably; (2) small eigenvalues may make a large haphazard contribution to the distance when a large number of characters is employed; (3) its value is heavily dependent on assumptions of multivariate normal distributions within the OTUs, which will usually not hold above the population or species level

In summary, generally the association coefficients include both positive and negative matches in their calculation of coefficients. Both similarity and dissimilarity coefficients are used for measuring resemblance. For binary variables, Jaccard (similarity) coefficient was used in Sneath (1957b) and single-matching coefficient was used in Sokal and Michener (1958). Sneath and Sokal (1973) advocated negative matches and hence favored the modified Gower’s coefficient which allows negative matches in two-state characters. They thought that the modified Gower’s coefficient appeared to be a very attractive index for expressing phenetic similarity between two OTUs based on mixed types of characters (p. 136).

For continuous variables, the Pearson product-moment correlation coefficient was most often used in both psychology and ecology Q-type studies previously (Stephenson 1936). At the end of 1950s, it was first proposed to use in numerical taxonomy by (Michener and Sokal 1957; Sokal and Michener 1958), and has been used since then and become the most frequently employed similarity coefficients in numerical taxonomy for calculating between pairs of OTUs (Sneath and Sokal 1973, p. 137). The coefficient, computed between taxa j and k, is given as:

(7.1)

where X_ij denotes the character state value of character i in OTU j, is the mean of all state values for OTU j, and n is the number of character sampled. The values of correlation coefficients range from −1 to +1.

However, overall Sneath and Sokal (1973, p. 140) believe that correlation coefficients are usually the most suitable measure in numerical taxonomy to evaluate OTUs. They do not think above undesirable properties would not manifest themselves in large data matrices involving many OTUs with the large number of characters. Additionally, when the interpretation of taxonomic structure is based on phenograms, it has been observed that the correlation coefficients are the most appropriate measure.

In summary, due to heterogeneity of OTUs, the basic assumptions of the bivariate normal frequency distribution cannot be met, so in psychology and numerical taxonomy literature, the character scores are often standardized or a percentage scale or a ratio of the variable are used before applying this coefficient (Sokal and Sneath 1963, p. 142).

One important kind of similarity coefficients is the information statistics. It has the advantage of additive property. Thus information statistics permits the partitioning of indices for larger groups into information statistic for subgroups (Sneath and Sokal 1973).

In ecology, information statistics is used to measure the frequencies or probabilities of the occurrence of the various species over the collection of stands. Sneath and Sokal introduced one probabilistic similarity coefficient based on information-theoretical concepts into numerical taxonomy to measure the homogeneity of the system by partitioning or subpartitioning sets of OTUs (Sneath and Sokal 1973, p. 120).

(7.2)

where m_i is the number of different states in character i, and p_ig is the observed proportion of the t OTUs exhibiting state g for character i. . When the n characters are not correlated, the total information of the group can be obtained by summing the separate values of .

The constructed information or entropy can also be applied to ecological data. In an ecological setting characters (rows) are species and OTUs (columns) are stands or quadrats, and hence for data matrices in ecology, this formula is to measure the frequencies of any one species that are recorded over a series of stands (Sneath and Sokal 1973, p. 142). To conduct hypothesis testing of “states” m_i, these frequencies must be grouped into classes to provide m_i “states” (e.g., “absent,” “rare,” “common,” and “abundant”) of a species (or character) i.

An approximately x² distribution with 2n degrees of freedom can be described by this formula in terms of frequencies for taxon J containing t_J OTUs. Where the quantity I measures the total information (entropy, surprisal) in the entire sample with the lower value of I indicating the more homogeneous a taxon will be. Thus, one goal of clustering OTUs using information statistics is to obtain taxa with low values of I (Sneath and Sokal 1973, p. 142).

Compared to most resemblance coefficients that only apply to pairs of single OTUs or pairs of taxa, another advantage of the information statistics is it’s applicable to sets of OTUs. Sneath and Sokal (1973, pp. 143–145) introduced an index based on two concepts: joint information and mutual information discussed by (Estabrook 1967; Orloci 1969). Given the correlation of characters, joint information I(h, i) is the union of the information content of two characters over a set of OTUs; mutual information I(h; i) is their intersection. The relation between these quantities is given below:

(7.3)

The mutual information I(h; i) is the ratio of the sums of the information content exclusive to characters h and i divided by the total information possessed jointly by h and i.

Overall, information statistics have important properties (Sneath and Sokal 1973, p. 120), such as are usually additive, are distributed as chi-square, and are probabilistic, and therefore a statistical hypothesis testing can be conducted through a standard statistical means. Additionally, the special cases of these coefficients can be represented as distance coefficients.

The denoising methods that are used to classify ASVs or sub-OTUs in current microbiome research actually take the same approaches as in probabilistic similarity coefficients (see Chap. 8 for details).

The choice of which resemblance coefficient not only is determined by the nature of the original data matrix, such as association coefficients for binary and categorical data matrices and distance and correlation for continuous variables, but also is based on the type of resemblance. Sokal (1985) reviewed that classifications have been based on correlation coefficients than based on similarity using distance-based measures. For binary data, correlations are not especially helpful and standardization does not help either.

In summary, regarding these three coefficients of association, distance measures, and correlation, Sokal and Sneath (1963, p. 167) thought that association and distance are easier to interpret conceptually and relatively simpler than correlation. Among these three coefficients, they recommended using the simple matching coefficient and Jaccard coefficient in case of excluding negative matches, the coefficient of taxonomic distance d_jk based on standardized characters, and Pearson’s product-moment r_jk, respectively. In their 1973 book Sneath and Sokal favor more strongly than in their 1963 book on binary similarity coefficient because of the following: (1) its simplicity and (2) possible relationship to information theory; (3) the similarity between fundamental units of variation can be estimated; and (4) related to natural measures of similarity or distance between fundamental genetic units (amino acid or nucleotide sequences) (Sneath and Sokal 1973, p. 147).

7.4.3.1 Similarity/Resemblance Matrix

The main techniques used for constructing taxonomic structure are clustering techniques, which are applicable to all three types of similarity coefficients (distance, association, and correlation coefficients). The more complicated method factor analysis can be used to construct taxonomic structure.

Four sets of methods have been developed in numerical taxonomy: (1) clustering, (2) factor analysis, (3) ordination, and (4) discriminant analysis. We can collectively call them as “numerical taxonomic methods” or “OTU methods” since they are used to classify or identify the OTUs. The first two sets of methods were mainly proposed or developed in Sokal and Sneath’s early work, such as in their 1963 book,whereas the last two sets of methods were mainly adopted into numerical taxonomy in their later works (e.g., 1973 book). The first three sets of methods are used for classification, and the discriminant analysis is used for both classification and identification. We describe them separately below.

One goal of clustering in numerical taxonomy is to recognize the patterns of distribution of OTUs and groups of OTUs (taxa) in a space (commonly a hyperdimensional, phenetic A-space) (Sneath and Sokal 1973, p. 192).

7.4.3.2 Definitions of Cluster in Numerical Taxonomy

When the patterns are descripted, they are typically conducted at least within a three-space (mathematically a higher-dimensional space): multidimensional Euclidean, non-Euclidean, and nonmetric spaces. Generally, regular and random distributions are the two observed dispersion patterns that can be compared. The clusters in numerical taxonomy are defined as the sets of OTUs in phenetic hyperspace that exhibit neither random nor regular distribution patterns and that meet one or more of various criteria.

7.4.3.3 Commonly Clustering Methods in Defining Groups of OTUs

Numerical classification can be performed by arranging of organisms in one of three ways (Sneath and Sokal 1973, p. 200): (1) hierarchical clustering systems of partitions (nested, mutually exclusive systems), (2) nonhierarchic clustering systems of partitions, and (3) ordination-the positioning of OTUs in relation to multidimensional character axes in an A space. In this section, we introduce clustering including both hierarchical clustering and non-hierarchic clustering. In the next section, we will briefly review ordination techniques in numerical taxonomy, and the more detailed introduction will be presented in Chap. 10 when microbiome study is covered.

In numerical taxonomy, cluster analysis is referred to as a large class of numerical techniques for defining groups of related OTUs based on high similarity coefficients. Clustering methods typically represent the similarity matrix among taxonomic units by dendrograms (a treelike arrangement). The grouped OTUs by clustering may form a coarser partition from a finer one, or may separate an entire set of OTUs into increasingly finer partitions.

Clustering can be described from eight aspects: (1) agglomerative versus divisive; (2) hierarchic versus non-hierarchic; (3) non-overlapping versus overlapping; (4) sequential versus simultaneous; (5) local versus global criteria; (6) direct versus iterative solutions; (7) weighted versus unweighted; and (8) non-adaptive versus adaptive. Among these eight methods, the sequential, agglomerative, hierarchic, non-overlapping (the acronym SAHN) clustering methods are the most frequently employed strategies for finding clusters (Sneath and Sokal 1973, p. 214).

The grouping of OTUs in numerical taxonomy has been most often employed by sequential clustering and by agglomerative techniques and rarely by divisive techniques. Agglomerative techniques start with a set of t separate entities, group them into successively fewer than t sets, arriving eventually at a single set containing all t OTUs. In contrast, divisive techniques start with all t OTUs in one set, subdividing it into one or more subsets, until each of the subsets is subdivided. Numerical taxonomy defines hierarchies in terms of nestedness as well as mutual exclusiveness. Thus, hierarchic clustering is to classify organisms into nested, mutually exclusive taxa. The hierarchic classifications are referred to as stratified. In a non-overlapping method, taxa at any one rank are mutually exclusive; that is, if the OTUs have been contained within one taxon, they cannot also be members of a second taxon of the same rank. For the details, the interested readers are referred to Sneath and Sokal’s 1973 book (Sneath and Sokal 1973, pp. 202–214).

Here, we summarize the commonly SAHN clustering methods used for defining groups of OTUs. Among them, four mainly described cluster methods are (1) single linkage clustering, (2) complete linkage clustering, (3) average linkage clustering, and (4) central or nodal clustering.

7.4.3.4 Factor Analysis

Taxonomic structure can also be constructed through factor analysis. Based on (Sokal and Sneath 1963), it is Sokal (1958b) who first used factor analysis to analyze taxonomic relationships from a similarity matrix of correlation coefficients, in which factor analysis is used for describing the complex interrelationships among taxa in terms of the smaller number of factors. Here, OTU is used as an observed variable, while taxon is used as a latent factor. The factor loading represents the degree of similarity between an OTU and the average aspect of the taxon. The higher the factor loading, the more typical the OTU belongs to the taxon. Later multiple factor analysis (Harman 1976) was initially used to phenetic classification by (Rohlf and Sokal 1962), in which centroid factor extraction with subsequent rotation is performed to obtain “simple structure” for revealing the interrelationships among OTUs.

Employing factor analysis to classification actually is to define taxa (e.g., species) in terms of similarity of OTUs by using factor analysis. The later bioinformatic work using factor analysis to define species originated in this work of numerical taxonomy. However, the challenge is how many factors should be extracted to define the taxa. In factor analysis the number of factors is usually arbitrarily chosen. Additionally, factor analysis cannot draw straight line across a dendrogram to yield the groups, suggesting that the factor analysis only can approximately, but not exactly form the groups at the same hierarchical level (Rohlf and Sokal 1962; Sokal 1963).

However, factor analysis has both practical and theoretical problems. Practically, before a factor analysis is conducted, factor scores and correlations of the characters are needed; thus factor analysis is not really efficient. Additionally, similar as it is a challenge to find clusters of highly correlated characters in cluster analysis, it is even risk to find factors when the characters are highly correlated. Thus, Sneath and Sokal suggested that we might still use the whole set of characters to estimate phenetic similarity from factor scores (Sneath and Sokal 1973, p. 99). Theoretically, it is difficult to know how to weight the extracted factors in measuring the similarity between OTUs. Thus, we need more research to develop an independent criterion for weighting factors.

7.4.3.5 Ordination Methods

Ordinations also present classifications, in which the data are first required to be transformed before being turned into ordinations. The ordination methods in phenetic taxonomy are referred to the methods that describe the occupation of phenetic space by the OTUs being ordinated (Sokal 1986). The ordination methods in application to phenetic taxonomy differentiate from clustering methods, in which OTUs are arranged in a continuous hyperspace with its dimensions being defined by the characters of the OTUs, whereas clustering methods assign OTUs to groups (Sokal and Sneath 1963; Sneath and Sokal 1973; Clifford and Stephenson 1975; Dunn and Everitt 1982; de Queiroz and Good 1997). Both clustering and ordination methods are based on similarities. However, clustering methods aim to base nested hierarchical taxonomies on similarity, whereas ordination methods do not attempt to give hierarchical results, which was reviewed as more accurately representing similarities than do by clustering methods (de Queiroz and Good 1997).

Narrowly OTU-method in numerical taxonomy only includes commonly used clustering methods and factor analysis. Based on (de Queiroz and Good 1997), Sokal pioneered the biological application of ordination techniques, e.g., in (Rohlf and Sokal 1962). However, in their 1963 book, Sokal and Sneath did not include ordination techniques as an approach to phenetic taxonomy. Ordination techniques were described as a potential approach to phenetic taxonomy in their 1973 book (Sneath and Sokal 1973), including principal component analysis (PCA), principal coordinate analysis (PCoA), nonmetric multidimensional scaling (NMDS), and seriation. PCA, PCoA, and NMDS have been readily adopted in microbiome study. We will introduce and illustrate PCA, PCoA, and NMDS as well as other ordination methods with real microbiome data in Chap. 10 (Sect. 10.​3). Here, we briefly describe seriation method. For details, the interested readers are referred to Sneath and Sokal (1973) and Legendre and Legendre (2012) and the references therein.

Seriation is to permute the rows and columns of a similarity matrix in such a way as to arrange the highest similarities (or the lowest distances) near the principal diagonal of the similarity (resemblance) matrix and with the similarity values decreasing orderly away from the diagonal. Thus, seriation method essentially is how to operationally define order among the coefficients in a similarity matrix (Sneath and Sokal 1973, pp. 250–251).

7.4.3.6 Discriminant Analysis in Identifying OTUs

Numerical taxonomy makes distinction between classification and identification. As we described in Sect. 7.2.3, at the beginning, unequal weighting was not promoted in numerical taxonomy. However, Sokal and Sneath recognized that strictly equal weighting is impossible in their later publications, such as in (Sneath and Sokal 1973). They distinguished identification from classification, in which identification or assignment is used to assign a new entity (Dagnelie 1966).

Classification is used in the sense of making classes, clusters, or taxa. By contrast, identification is used for “identifying” new individuals (OTUs/taxa). However, when we identify an individual as new one, we actually know that an existing class is unacceptable. The cluster analyses can effectively work on the classification. But to identify new individuals, it requires that some strategies can combine the two procedures of classification and identification. The discriminant functions (analyses) are the main methods used in identification.

Along with the distinction between classification and identification, Sneath and Sokal believed that weighting is needed. The goal of weighting characters in numerical taxonomy is for identification. Sneath and Sokal (1973, p. 400) thought that weighting can maximize the probability of correctly identifying unknown specimens from a few close or overlapping taxa through a set of multivariate statistical methods.

Sneath and Sokal have specifically introduced three discriminant methods in describing or discriminate OTUs (or taxa): discriminant function analysis (DFA), and multiple discriminant analysis in Mahalanobis’ D-space and canonical variates analysis. We can view these multivariate statistical methods as extensions of taxon distance models in Sect. 7.4.2.3. The purpose of taxon distance models in Sect. 7.4.2.3 is to classify the OTUs, while the purpose of discriminant analyses is to perform both classification and identification, i.e., to discriminate the OTUs (or taxa).

7.4.3.6.1 Discriminant Function Analysis (DFA)

DFA or linear discriminant analysis (LDA) or linear discriminant function analysis (LDFA) is a generalization of Fisher’s linear discriminant method (Fisher 1936). LDFA is a member of the family of classification methods known as supervised learning, in which group membership of individuals is used to maximize group differences during the process. In this method, a linear combination of features that characterizes or separates two or more classes of objects is to be identified and to be used as a linear classifier or for dimensionality reduction before later classification. LDA requires that groups (i.e., the class label) are known a priori.

In numerical taxonomy, the linear discriminant function is a linear function z of characters describing OTUs that weight and separate the characters in such a way to identify two taxa: maximizing number of the OTUs in one taxon has high values for z and maximizing number of the OTUs in another taxon has low values, so that z serves as a much better discriminant of the two taxa than does any one character taken singly (Sneath and Sokal 1973, p. 400). By performing the linear discriminant function, (1) the n characters are almost always reduced to a smaller set, m; and (2) the function also has maximal variance between groups relative to the pooled variance within groups.

The dichotomous discriminant function that was originally developed by (Fisher 1936) was applied to two taxa, and was later generalized to multiple taxa. In LDA the pooled variances and covariances between the m characters within each taxon are calculated. They are calculated by a weighted average of the variances and covariances of the characters in taxa J and K in its original form, while in its generalized form, the variances and covariances are calculated from pooling all the taxa being considered. Additionally, the means of each character for each taxon are calculated to represent the centroids of the taxa. Suppose J, K, and L represent three taxa for two character axes, X₁ and X₂. LDA can be briefly summarized as the following eleven steps. For the details, the readers are referred to Sneath and Sokal (1973, pp. 401–405).

• Step 1: The linear discriminant function generates a set of taxa (i.e., three taxa J, K, and L) for all the character axes (i.e., X₁ and X₂). The generated taxa are represented by the clusters of individuals (OTUs).

• Step 2: After clustering individuals (OTUs), we now only concern the descriptive parameters of the clusters (i.e., three taxa J, K, and L), i.e., the centroids and the dispersions of the clusters and do not concern with the values for the individuals.

• Step 3: First calculate the variances and covariances between the characters, yielding a set of (here three) m × m matrices (here m = 2); then average these variances and covariances to give a pooled within-groups variance-covariance matrix W.

• Step 4: Calculate the discriminant function between J and K. Multiply the inverted W matrix by the vector to obtain the discriminant function as a vector z: z_JK = W⁻¹δ_JK. This vector consists of a series of weights, w_i, for characters 1, 2, …, m, it is symbolized as z₁, z₂, …, z_m.

• Step 5: Calculate a discriminant score DS (here the discriminant score scale for taxa J and K). Multiply these weights by the observed character values of the unknown individual, and sum them to give a discriminant score DS: DS_u = z₁X_1u + z₂X_2u + … + z_mX_mu.

• Step 6: Use this discriminant score to calculate three reference scores for the centroid of J, the centroid of K, and the point midway between them: DS_J, DS_K, and DS_0.5 respectively for discriminating between members of J and K. They are given by:

• 

• Step 7: Calculate the length of the line between the centroids of J and K and allocate the unknown to taxon. First, obtain the square root of Mahalanobis’ D², which is the length of the line between the centroids of J and K measured in discriminant function units (the absolute difference between the scores DS_J and DS_K is also equal to D²); then allocate the unknown to taxon J if the observed score for an unknown,DS_u, lies on the DS_J side, otherwise to taxon K if on the D DS_K side. A different discriminant function, and discriminant scores, is calculated for each pair of taxa. The square root of Mahalanobis’ D² (Mahalanobis 1936) is used to measure the length of the line between the centroids of two taxa J and K in discriminant function units, which equals to the absolute difference between the discriminant scores of these two taxa DS_J and DS_K.

• Step 8: Interpret discriminant function vectors. The vector angle between two discriminant function vectors can be interpreted as measuring contrasts of form of the taxa, where the small angle measures similar contrasts of form and large angles represent distinct contrasts.

• Step 9: Perform statistical hypothesis testing of the centroids using the distance scores. The distance between the scores DS_J and DS_K (Sneath and Sokal 1973, pp. 403–404) has two important uses. One is that we can perform a null statistical hypothesis through a F-test to test whether the centroids are significantly different with a testing ratio:

• , with m and (t_J + t_K − m − 1) degrees of freedom.

The test ratio is related to Hotellings’ T² (Rao 1952):

(7.4)

• Step 10: Perform variable (character) selection. We also can use the distance scores to determine each character contributing to D². If the characters with little discriminatory power can be identified, then drop them from the analysis. Alternatively, we can choose the best few characters from the set by the F test. If we do not consider correlations between characters, the percent contribution of character i is 100 × (z_iδ_i/D²), where z_i and δ_i are the ith elements of vectors z_JK and δ_JK.

• Step 11: Calculate the probability of misclassification. A primary purpose of a discriminant function is to minimize the probability of wrong assignment of unknown individuals. We can calculate the probability of misclassification of assignment of unknown individuals by assuming a multivariate normal distribution and also that the unknown does belong to J or K (and not to some distant cluster).

7.4.3.6.2 Multiple Discriminant Analysis (MDA)

The original dichotomous discriminant function has two limitations (Sneath and Sokal 1973, p. 404): (1) The algorithm used in the original dichotomous discriminant function implies that the unknown either belongs to one or the other of the two taxa being considered. Thus, when an unknown taxon is located far off in the space, belongs to a quite different cluster, and hence it really belongs to a third cluster. However, using the dichotomous discriminant function, this unknown taxon may have almost any discriminant score, and thus is assigned to one or the other of the two taxa under consideration. (2) It cannot test against a large number of discriminant functions when there are many taxa.

MDA in Mahalanobis’ D-space (Mahalanobis 1936) was developed to largely overcome the above two problems of the original dichotomous DFA. MDA is a multivariate dimensionality reduction technique. It is not directly used to perform classification, instead merely supports classification by compressing the signal down to a lower-dimensional space amenable to classification. MDA typically projects the multivariate signal down to a dimensional space with less one dimension of the number of categories. Mahalanobis (1936) generalized calculation of D² between the centroids from the discriminant scores to between any pair of points so that the square root of D² become simply the Euclidean distance in the D-space. The Mahalanobis’ D² between any pair of points f and g is given:

(7.5)

MDA method transforms the original space into a new space, where the original axes are stretched and also skewed, and results in not any more right angles between axes.

The length of a unit in dimension equals to the original units times by the square root of the element of the matrix w⁻¹ (Sneath and Sokal 1973, p. 405). The direction cosine between the dimensions has an interpretation of correlation coefficients (Gower 1967b).

7.5.1.1 Heterogeneity of the Column Vectors of the Data Matrix

In the OTU Data Matrix, the similarity between OTUs j and k are measured over a set of n characters. The problem is that these OTUs are not taken as random samples from a common population, but instead are really taken from a heterogeneous sample wherein each variate estimates a different character. In other words, the character states for any one OTU are not samples that estimate a single mean and variance. Thus, because the column vectors of the data matrix are heterogeneous, it is difficult to assume the sampling parametric distributions for justifying statistical test; and hence it is inappropriate to make ordinary tests of significance (Sokal and Sneath 1963, p. 313). Since a parametric model cannot be easily specified, in ecology and microbiome study, a permutation test is often used to perform a hypothesis testing.

7.5.1.2 Hypothesis Testing of Individual Similarity Coefficients

An appropriate parametric statistical test in multivariate statistics requires the variables (in this case, characters) are independent so that there are p variables (characters) in the data matrix, and then the data matrix has p dimension with each of the p variables representing one dimension of the data. The problem is that in OTU data matrix, the characters (rows of the data matrix) are correlated with each other, and thus a similarity coefficient based on p characters would not reflect p independent dimensions of variation (Sneath and Sokal 1973, p. 163). Thus, standardization of characters is often recommended in numerical taxonomy to alleviate the heterogeneity of column vectors. However, standardization of every character state to have a population mean of zero and variance of unity only mitigate some degree of heterogeneity of column vectors. Actually, put another way, statistical hypothesis tests in numerical taxonomy only have partially succeeded for entire classifications (similarity matrices), but not for individual coefficients (Sneath and Sokal 1973, pp. 163–164). In numerical taxonomy, scaling continuous characters (variables)(lengths, percentages, ratios) is recommended by either standardization or ranging scaling methods. Because of scaling characters, correlation or taxonomic distance matrices between OTUs hence are altered. Two states of characters should be coded as binary, while ordered discrete multistate characters (qualitative traits and counts) can be analyzed as continuous characters (Sokal 1986).

These two problems actually reveal the challenges of directly conducting the hypothesis testing of individual similarity coefficients in numerical taxonomy, and ecology as well as in later microbiome study. Performing hypothesis testing pairwise correlation or partial correlation also poses a challenge. In a similar reason, beta diversity measures such as Bray-Curtis, Jaccard, and Sørensen dissimilarities are usually used as input data for ordination and association analysis rather than directly being tested. This links up the concepts and topics in multivariate statistics, with which we will discuss in Chapters 10 and 11.

7.6.2.1 Summary of Characteristics of Numerical Taxonomy

In numerical taxonomy classifications are called nature or “general” classifications, and the constructed taxa have “natural” properties or are “natural” taxa/“natural taxa.” For Sneath and Sokal (1962), “natural taxonomic groups” are “polytypic” (Beckner 1959) or better called as “polythetic” (meaning “many-arrangement”), which is against “monothetic” (meaning “one-arrangement”).

Sokal and Sneath tried to avoid confusing phylogeny with classfication (Huxley 1869) and emphasized that taxonomic classifications are distinctive from phylogenetic classifications. Basically numerical taxonomy does not attempt to force phylogenetic criteria on to taxonomic groupings. Sokal and Sneath separated the taxonomic process (which is to be based on affinity or resemblance) from phylogenetic speculation (Sneath and Sokal 1962). They emphasized that taxonomic relationships are “static” (Michener 1957) or preferably are called “phenetic” (Sneath and Sokal 1962; Cain and Harrison 1960), which are to be evaluated purely based on the resemblance among the materials. In numerical taxonomy it is logically fallacious to make phylogenetics deductions to arriving at taxonomies. Sokal and Sneath descibed a fallacy of circular reasoning (Sokal and Sneath 1963, p. 21):

“Since we have only an infinitesimal portion of phylogenetic theory in the fossil record, it is almost impossible to establish nature taxa on a phylogenetic basis. Conversely, it is unsound to derive a definitive phylogeny from a tentative natural classification.”

In numerical taxonomy, taxa are based on correlations between features (Sneath and Sokal 1962), and taxonomic rank is defined by similarity (Sneath 1957a). The phenetic resemblance between the taxonomic entities is estimated using the similarity and dissimilarity coefficients (Sneath and Sokal 1962).

Three convenient types of coefficients (association, correlation, and distance) were suggested in numerical taxonomy as suitable measures of similarity between pairs of taxa (Sokal 1961). Specifically association coefficient is utilized to measure “present-absent” (binary) characters, which was first described by Sneath (Sneath 1957a). Binary data have attractive properties of simplicity and the potential to decompose variation among OTUs into fundamental units of information. The simple matching coefficient (including its one-complement coefficient, the average Manhattan distance coefficient) and the Jaccard coefficient were recommended for use in numerical taxonomy among the many binary association coefficients that most have been used in ecology (Sokal and Sneath 1963; Sokal 1986).

The ordinary Pearson product-moment correlation coefficient is still preferred to use as measures of similarity between OTUs because it is less sensitive to size differences, although it is non-metrics and non-invariant to the direction (Sokal 1986). In contrast, cosine coefficient is metrics and invariant to the direction of correlation coefficients. For continuous or ordered multistate characters (variables), both taxonomic distances and Manhattan distances can be used to calculate similarity and dissimilarity. Distance coefficients determine the similarities between two taxa as a function of their distance in a multiple dimensional space with the characters as coordinates. They are employed to measure continuous or more than two states (ordered multistate) of the characters (variables) (Sokal 1958a; Morishima and Oka 1960). Between taxonomic distances and Manhattan distances that can calculate similarity and dissimilarity, Sokal (1986) preferred taxonomic distances because of their geometric and heuristic properties. Based on Sokal (1961), the use of distance to measure taxonomic similarity can be located earliest by Heincke (1898) and others (Anderson and Abbe 1934; Pearson 1926; Rao 1952; Mahalanobis 1936). Both correlation and cosine coefficients are sensitive to data points close to the centroid. However, classifications in numerical taxonomy are preferred to be based on similarity using correlation (or cosine) coefficients than by using distance-based measures because the former usually provide better effects of classifications to exhibit size and shape differences in the data (Sokal and Rohlf 1980; Sokal 1986).

Taxonomic groups are numerically established by cluster analysis of values of similarity. To exploit the taxonomic structure, numerical taxonomy uses different clustering methods to group OTUs into clusters for achieving the taxonomic structure. Three commonly used clustering methods in numerical taxonomy are single linkage clustering, complete linkage clustering, and average linkage clustering methods. Among these three clustering methods, the average-linkage clustering method was promoted; and particularly the unweighted pair-group method using arithmetic averages (UPGMA) was also preferred for use because the average clustering method can present a compromise solution to form a cluster (Sokal 1986).

It is also because single-linkage and complete-linkage clustering methods represent the two extreme criteria and both could not give taxonomically satisfactory results (Hartigan 1985), although single linkage clustering has several elegant mathematical properties (Sokal 1986). Thus they both are used in special cases only (Sokal 1963; Sneath and Sokal 1973).

OTUs were introduced as a pragmatic definition to classify groups of organisms, allowing for taxonomy being classified quantitatively or numerically. Thus, the concept of OTUs and the field of numerical taxonomy were established against the subjective nature of classification in classical Linnaean taxonomy. Through creating the concept of OTUs, the term “taxon” can still be retained for its proper function role played in taxonomic groupings. Thus, this indicates that OTUs are not true taxa instead are operated as taxa. OTUs are used in constructing classifications without considering any phylogenetic information. Actually by using the OTUs, the classification was considered purely as a study of similarity, of “static” relationships in which phylogeny played no part (Vernon 1988). Thus, when constructing phylogenetic classifications, the static relationships and phylogeny should be separately assessed.

7.6.2.2 Remarks on Characteristics of Numerical Taxonomy

The approach of phonetic numerical taxonomy (or phenetics) and commonly clustering-based OTU methods have several challenges.

First, phonetic numerical taxonomy (Gilmour naturalness) assumes that high similarity (in fact largely clustered) will also lead to high predictivity. This has been identified as an unfinished problem (Farris 1979; Sokal 1985) because the linkage between similarity and predictivity has not been defined clearly. Actually, on one hand, it is challenging to show how to maximize similarity; on the other hand, the relationship between classification and predictivity has not yet been established although a close relation has obviously been observed between classifications in which objects within any given taxon resemble each other closely and predictivity of characters within that taxon (Sokal 1985). Predictivity is the most important criterion of a classification in terms of Gilmour naturalness. Thus, although it is difficult to achieve, generally high predictivity is the goal of a natural classification (Sokal 1985).

For example, phenetic clustering is commonly used in specifically and biologically estimating phylogeny or phylogeny reconstruction, and is also used to analyze relationships among potentially conspecific populations, as well as geographic variation and genetic continuity, and consequently, of species limits (Highton et al. 1989). Both above applications of phenetic clustering have been criticized having serious limitations because of their assumptions are questionable (de Queiroz and Good 1997). When phenetic clustering is used for estimating phylogeny, it assumes that similarity is correspondent to recency of common ancestry, while when it is used for analyzing patterns of geographic variation and genetic continuity among populations, it assumes that similarity is correspondent to degree of genetic continuity.

The fundamental weakness of numerical taxonomy and hence its clustering-based OTU methods is its questionable assumptions. Because similarity relationships do not strictly conform to the pattern of a nested hierarchy, and thus attempts to base nested hierarchical taxonomies on similarity will never be fully satisfactory (de Queiroz and Good 1997).

Second, the classifications in biological taxonomy have the property that the resulting classification should be hierarchical and nonoverlapping; and hence requires the resulting classification represents faithfully as possible the actual similarities among all pairs of OTUs resemblance matrix, and the measures (e.g., the cophenetic correlation coefficient) can express the goodness of a classification (Sokal and Rohlf 1962; Sokal 1985). However, it was shown that different clustering methods have different goodness of fits by cophenetic correlation coefficients to the same resemblance matrix (Sokal 1985). Both artificial data (Sneath 1966) and real data (Sokal 1985) have shown average-linkage clustering outperformed complete and single linkage clustering methods to measure the resemblances.

Overall, in one hand, importantly phenetic clustering OTUs expects to result in nested hierarchical structure. However, similarity that does not exhibit a strictly nested hierarchical structure has important consequences for its analysis using phenetic clustering (de Queiroz and Good 1997). On the other hand, even random distributions of OTUs will result in the formation of clusters [see(Sneath and Sokal 1973, p. 252; Dunn and Everitt 1982, p. 94; de Queiroz and Good 1997)].

Third, whether clustering methods or ordination analysis methods better represent taxonomic structure is still arguable. The phenetic criteria of goodness of classification are restricted to hierarchic Linnaean classifications. However, non-Linnaean modes of data analysis would provide better fits to the original resemblance matrix than dendrograms generated by hierarchic cluster analysis methods; thus, ordinations could summarize the taxonomic structure, and hence various ordination analysis methods have been employed in phenetic taxonomy; on the other hand, the hierarchic clustering analysis models have been used as standard methods for classification in the narrow sense (Sokal 1985). Given robust methods for similarity coefficients, Sokal (1985) considered that cluster analysis may ultimately result in obtaining stable and repeatable phenetic classification although it is not possible to achieve complete stability.

7.6.2.3 Remarks on Impact of Numerical Taxonomy on Microbiome Research

The philosophy, OTU concept, framework of statistical analysis of phonetic numerical taxonomy, and especially commonly clustering-based OTU methods have been adopted into microbiome research. All these have provided a solid foundation in the development of microbiome research, especially in the theoretical and conceptual frameworks. The philosophy and methodology of numerical taxonomic procedure have posed large impacts on early microbiome study. However, the limitations and challenges of phonetic numerical taxonomy both as a method for clustering OTUs, and as a statistical analysis method for assessing taxonomies and microbiomes still remain, and even more severe when using OTU methods to define species.

However, both philosophy and methodology of numerical taxonomy have their own fundamental weaknesses. When the philosophy and methodology of numerical taxonomy, and especially the clustering-based OTU methods, are adopted into microbiome study to define species, these weaknesses have not disappeared but somehow have been highlighted.

The methodology and especially the clustering based-OTU methods in numerical taxonomy have played an important role in defining OTUs in bioinformatic analysis of microbiome. In the meantime, we want to emphasize that both the philosophy and methodology of numerical taxonomic procedure have their limitations. When the clustering based-OTU methods are used in microbiome study to define species, tasks have been actually assigned: not only for construction of taxonomy but also implicitly requiring for a biological or functional role that the defined species should have. This is beyond the original means described in numerical taxonomy.

Numerical taxonomy along with its clustering based-OTU methods was a reform movement that arose in response to a dissatisfaction with the subjective and nonquantitative nature of taxonomy as it existed in the late 1950s (Vernon 1988). However, the applications of phenetic clustering methods within the context of taxonomy and in the biological sciences have never been entirely satisfactory (de Queiroz and Good 1997). The application of clustering based-OTU methods to bioinformatic classification of OTUs in the field of microbiome still neither is entirely satisfactory.

Actually, a reform movement in response to a dissatisfaction using clustering based-OTU methods to define taxonomy that has already started in the bioinformatic analysis of microbiome sequences. We will comprehensively review the bioinformatic analysis methods and the movement beyond OTUs in Chap. 8.

Numerical taxonomy in its original setting mainly is a purely descriptive science; its potential to be used for formulating hypotheses has not been fully developed in numerical taxonomy. The potential for formulating hypotheses is developed by microbiome researchers in their differential abundance analysis and association analysis. We will cover statistical analysis methods and models in Chaps. 9, 10, 11, 12, 13, 14, 15, 16, 17 and 18.