# 2-Wayk-Means as a Model for Microbiome Samples

Authors: Weston J. Jackson; Ipsita Agarwal; Itsik Pe’er
Journal: Journal of Healthcare Engineering
 (2017)

---

## Abstract

Motivation. Microbiome sequencing allows defining clusters of samples with shared composition. However, this paradigm poorly accounts for samples whose composition is a mixture of cluster-characterizing ones and which therefore lie in between them in the cluster space. This paper addresses unsupervised learning of 2-way clusters. It defines a mixture model that allows 2-way cluster assignment and describes a variant of generalized k-means for learning such a model. We demonstrate applicability to microbial 16S rDNA sequencing data from the Human Vaginal Microbiome Project.


---

## Body

## 1. Introduction

Microbiome analysis [1] by sequencing of ubiquitous genes, most commonly 16S rRNA, is a standard, cost-effective way to characterize the composition of a microbial sample. Standard analysis tools facilitate quantifying the fraction of sequence reads from each bacterial species in a sample [2]. Interpretation of composition vectors across a collection of samples typically relies on dimensionality reduction followed by clustering in the lower-dimensionality space [3]. This allows identification of functionally meaningful subsets of samples with characteristic microbiota. The Human Microbiome Project [4] and its derivatives such as the Human Vaginal Microbiome Project [5] have collected and thus analyzed large numbers of samples towards elucidating the structure and composition of microbiota across physiological and pathological states.Similar to variation in microbial genomes across different human individuals, variants along the nuclear genomes have been summarized by a small number of dimensions [6]. However, in contrast to analyses of microbiome samples, analyses of inherited genetic variation standardly assume and observe samples to be spread across a continuum in the reduced space, rather than be clustered [7]. Samples in between clusters are interpreted as originating from intermediate locales along a geographic cline [8] or as representing different levels of a mixture between cluster-specific populations.In this paper, we formally tackle the problem of clustering while allowing elements to belong to two clusters. Specifically, we will describe in detail a model for clustering inℝd. We construct a model that generalizes k-means clustering by allowing data points to be assigned to a point in the space along the line between two assigned clusters [9]. Each cluster is still modeled as a Gaussian with uniform, spherical covariances; the key difference is the presence of a parameter u∈0,1 for each 2-way-assigned data point xi, which determines the proportional assignment of xi between its two cluster representatives. We first describe the 2-way model’s inputs, parameters, and outputs. We then give the objective function, an algorithmic description, and a series of performance metrics. Next, we evaluate the performance on simulated data, describing benchmarks for optimal performance. Finally, we apply the model to real data of 16S rDNA sequencing from 1500 midvaginal bacterial samples by the Vaginal Human Microbiome Project.


## 2. Methods

### 2.1. 2-Wayk-Means

The model characterizes a mixture where points are each sampled either from ak-mixture of uniform, spherical Gaussian distributions or from pairwise weighted averages of these Gaussians.Formally, we describe a generative model for a setX of data points xi1n∈ℝd. The model involves k∈ℤ+ clusters. The jth cluster is parametrized by its mean μj∈ℝd. To simulate xi, the model first chooses a pair of cluster indices (j,j′) along with a weighting ui∈0,1. xi is drawn from a Gaussian distribution whose parameters are ui-weighted averages of two representative clusters. Specifically, xi~Nxi;μ˜ijj′,Σ such that μ˜ijj′=uiμj+1−uiμj′ and Σ∈ℝd×d is the given uniform, spherical covariance matrix.The inference problem involves the inputs of dataX and number of clusters k, seeking output of the generative model parameters, that is, the vectors of assignments C=c1,…,cn and weights U=u1,…,un.


### 2.2. Generalizedk-Means

Given inputx1,…,xn∈ℝd and cardinality k∈ℕ, k-means traditionally provides us with the following objective:
(1)∑i=1nminj∈kxj−cj22,where c1,…,cj are the cluster representatives. The k-means objective can be generalized as the following:
(2)minC,Φ∑i=1nxi−Cϕi22,where Φ=ϕ1∣ϕ2∣⋯∣ϕn∈0,1k×n are the cluster assignments and C=c1∣c2∣⋯∣ck∈0,1d×k are the cluster representatives.A common generalization ofk-means is to permit each ϕi to have s nonzero entries (in our case, we set s=2). An algorithm for this generalized objective is simply to hold C fixed while performing sparse regression on Φ and then hold Φ fixed and use ordinary least squares (OLS) to find C.In our case, because we only allow pointsxi to lie uniformly between two cluster representatives, the two nonzero entries in a given ϕi are restricted to some ui∈0,1 and 1−ui∈0,1. Our problem is instead the following:
(3)minC,Φ∑i=1nxi−Cϕi22,subject to
(4)ϕi0≤2,ϕi1=1,ϕi≥0.


### 2.3. 2-Wayk-Means Algorithm

Our goal is to find a nonnegative 2-sparse solution for eachϕi. To do so, we can minimize over all k2 cluster representative possibilities. This 2-sparse solution gives us indices j,j′ which correspond with the two cluster representatives. This corresponds with the following objective:
(5)minui∈ℝ,cj′,cj∈Cxi−uicj+1−uicj′22,subject to ui∈0,1.For a givencj and cj′, minimizing with respect to uijj′ reveals a global minimum at
(6)cj′−cjTcj′−xicj′−cj22.After minimizing with respect touijj′, we project uijj′ to the region 0,1. We set uijj′=0 if the minimizer is less than 0 and set uijj′=1 if the minimizer is greater than 1. This allows us to achieve the minimum value of ui over the domain 0,1 for xi.After minimizing the assignmentΦ, we then use OLS to pick optimal C as specified before. Formally, OLS produces a vector ciT that minimizes the squared residual error between an input matrix ΦT and vector xiT.
(7)minciTxiT−ΦTciT22.Taking the gradient and setting equal to zero yields the following formula:(8)ciT=ΦΦT−1ΦxiT.Thus, we perform OLS for all vectorsciT at once with matrix multiplication:
(9)CT=ΦΦT−1ΦX.Thus, this gives us representativesc1,…,ck that minimize the residual error between the cluster representatives and data points subject to Φ. We then alternate this process for r rounds until convergence.


### 2.4. Performance Metrics

We use the 2-wayk-means objective as a performance metric in measuring the accuracy of a model in unsupervised examples.
(10)objX,k,r=minC,Φ∑i=1nxi−Cϕi22,where Φ has at most two nonzero entries with values ui∈0,1 and 1−ui∈0,1.Additionally, we also use four different error rates to measure the accuracy of 2-wayk-means on test cases. Let ci∗,μj∗, and uijj′∗ be the ground truth instance parameters, that is, respectively, true 2-way cluster assignment of xi, center of cluster j, and 2-way weighting for xi between clusters j,j′.e
r
r
f
x defines the 0-1 error rate for 2-way cluster assignment:
(11)errfx=∑i=1n1ci≠ci∗n.e
r
r
μ defines the squared deviation from optimal μ∗:
(12)errμ=∑∀jμj−μj∗2.e
r
r
u defines the squared deviation from optimal uijj′∗. WLOG, we assume uijj′=maxu,1−u, where u is the variable drawn from 0,1:
(13)erru=∑i=1nuijj′∗−uijj′2n.


## 2.1. 2-Wayk-Means

The model characterizes a mixture where points are each sampled either from ak-mixture of uniform, spherical Gaussian distributions or from pairwise weighted averages of these Gaussians.Formally, we describe a generative model for a setX of data points xi1n∈ℝd. The model involves k∈ℤ+ clusters. The jth cluster is parametrized by its mean μj∈ℝd. To simulate xi, the model first chooses a pair of cluster indices (j,j′) along with a weighting ui∈0,1. xi is drawn from a Gaussian distribution whose parameters are ui-weighted averages of two representative clusters. Specifically, xi~Nxi;μ˜ijj′,Σ such that μ˜ijj′=uiμj+1−uiμj′ and Σ∈ℝd×d is the given uniform, spherical covariance matrix.The inference problem involves the inputs of dataX and number of clusters k, seeking output of the generative model parameters, that is, the vectors of assignments C=c1,…,cn and weights U=u1,…,un.


## 2.2. Generalizedk-Means

Given inputx1,…,xn∈ℝd and cardinality k∈ℕ, k-means traditionally provides us with the following objective:
(1)∑i=1nminj∈kxj−cj22,where c1,…,cj are the cluster representatives. The k-means objective can be generalized as the following:
(2)minC,Φ∑i=1nxi−Cϕi22,where Φ=ϕ1∣ϕ2∣⋯∣ϕn∈0,1k×n are the cluster assignments and C=c1∣c2∣⋯∣ck∈0,1d×k are the cluster representatives.A common generalization ofk-means is to permit each ϕi to have s nonzero entries (in our case, we set s=2). An algorithm for this generalized objective is simply to hold C fixed while performing sparse regression on Φ and then hold Φ fixed and use ordinary least squares (OLS) to find C.In our case, because we only allow pointsxi to lie uniformly between two cluster representatives, the two nonzero entries in a given ϕi are restricted to some ui∈0,1 and 1−ui∈0,1. Our problem is instead the following:
(3)minC,Φ∑i=1nxi−Cϕi22,subject to
(4)ϕi0≤2,ϕi1=1,ϕi≥0.


## 2.3. 2-Wayk-Means Algorithm

Our goal is to find a nonnegative 2-sparse solution for eachϕi. To do so, we can minimize over all k2 cluster representative possibilities. This 2-sparse solution gives us indices j,j′ which correspond with the two cluster representatives. This corresponds with the following objective:
(5)minui∈ℝ,cj′,cj∈Cxi−uicj+1−uicj′22,subject to ui∈0,1.For a givencj and cj′, minimizing with respect to uijj′ reveals a global minimum at
(6)cj′−cjTcj′−xicj′−cj22.After minimizing with respect touijj′, we project uijj′ to the region 0,1. We set uijj′=0 if the minimizer is less than 0 and set uijj′=1 if the minimizer is greater than 1. This allows us to achieve the minimum value of ui over the domain 0,1 for xi.After minimizing the assignmentΦ, we then use OLS to pick optimal C as specified before. Formally, OLS produces a vector ciT that minimizes the squared residual error between an input matrix ΦT and vector xiT.
(7)minciTxiT−ΦTciT22.Taking the gradient and setting equal to zero yields the following formula:(8)ciT=ΦΦT−1ΦxiT.Thus, we perform OLS for all vectorsciT at once with matrix multiplication:
(9)CT=ΦΦT−1ΦX.Thus, this gives us representativesc1,…,ck that minimize the residual error between the cluster representatives and data points subject to Φ. We then alternate this process for r rounds until convergence.


## 2.4. Performance Metrics

We use the 2-wayk-means objective as a performance metric in measuring the accuracy of a model in unsupervised examples.
(10)objX,k,r=minC,Φ∑i=1nxi−Cϕi22,where Φ has at most two nonzero entries with values ui∈0,1 and 1−ui∈0,1.Additionally, we also use four different error rates to measure the accuracy of 2-wayk-means on test cases. Let ci∗,μj∗, and uijj′∗ be the ground truth instance parameters, that is, respectively, true 2-way cluster assignment of xi, center of cluster j, and 2-way weighting for xi between clusters j,j′.e
r
r
f
x defines the 0-1 error rate for 2-way cluster assignment:
(11)errfx=∑i=1n1ci≠ci∗n.e
r
r
μ defines the squared deviation from optimal μ∗:
(12)errμ=∑∀jμj−μj∗2.e
r
r
u defines the squared deviation from optimal uijj′∗. WLOG, we assume uijj′=maxu,1−u, where u is the variable drawn from 0,1:
(13)erru=∑i=1nuijj′∗−uijj′2n.


## 3. Results

### 3.1. Example Run for 2-Wayk-Means

We find it illuminating to demonstrate the performance of 2-wayk-means versus vanilla k-means on a cartoon example.In Figure1, we simulated n=500 data points in ℝ2 from three clusters, with respective means μ1=0,0,μ2=1,1,andμ3=2,0 and covariance matrices Σ=0.001I. Data points are drawn into pairwise clusters by choosing two cluster representatives without replacement from the following prior probabilities:
(14)Pc1=0.2,Pc2=0.5,Pc3=0.3.Figure 1
n
=
500 simulated data points. The white diamonds are cluster centers for three simulated clusters (red cluster, bottom left; green cluster, top; and blue cluster, bottom right). Points are colored as a linear combination of the clusters they lie between (according to u).We initialize the cluster representatives with vanillak-means. Vanilla k-means achieves the results in Figure 2. Statistics for vanilla k-means are given as follows:
(15)obj:41.1025,errfx:0.112,errμ:0.4039,erru:0.2435.Figure 2
n
=
500 simulated data points after k-means. The white diamonds are cluster centers for three simulated clusters (red cluster, bottom left; green cluster, top; and blue cluster, bottom right). White triangles are cluster centers determined by u means. Colors are u values determined by k-means.k-means predicts the 2-way cluster assignments of ≈11% points incorrectly due to skewing the cluster means toward the middle of the graph. 2-way k-means, however, significantly improves on all error rates. After 10 rounds of 2-way k-means, we achieve the results in Figure 3.
(16)obj:5.822,errfx:0.032,errμ:0.0495,erru:0.2084.Figure 3
n
=
500 simulated data points after 10 rounds of 2-way k-means. White diamonds are cluster centers determined by 2-way k-means (10 rounds). Colors are u values determined by 2-way k-means (10 rounds).For every statistic, the results are clearly an improvement on standardk-means. The ≈3% error rate on cluster assignment still exists because 2-way k-means points closest to cluster representatives may be assigned to an incorrect secondary cluster.


### 3.2. Benchmarks

#### 3.2.1. Sparsity (Avg. of 10 Trials, 10 Rounds Each)

Our sparsity test was conducted by keeping cluster prior probabilities and cluster centersμ constant while varying the number of data points (ratio of α means n=500α). From Figure 4, we see that the algorithm performs consistently well under a variety of conditions, but too few data points can hurt performance to an extent.Figure 4
Error rates whenn=α500, and cluster priors and centers are fixed.


#### 3.2.2. Cluster Separation (Avg. of 10 Trials, 10 Rounds Each)

We test the error rate as a function of the Euclidean distances ofμ (ratio of α means μ1=α0,0,μ2=α1,1,andμ3=α2,0). From the results in Figure 5, we can see that a certain threshold is required for proper performance of the algorithm. This makes sense, as when α=0.01, the clusters are almost on top of each other and difficult to distinguish. Additionally, as the cluster centers are moved farther apart, the ℓ2 norm between the cluster representative determined by the algorithm and the actual cluster representative increases (but this is to be expected).Figure 5
Error rates as a function of the Euclidean distances ofμ, where μ1=α0,0,μ2=α1,1,andμ3=α2,0.


#### 3.2.3. Variance (Avg. of 10 Trials, 10 Rounds Each)

We increase the variance of the clusters while fixing cluster prior probabilities, data points, and cluster centers (ratio ofα means Σ=α0,0.0001,0,0.001). From the results in Figure 6, we can see that large variance hurts proper performance of the algorithm. Analogous to cluster separation, as when α=100, the clusters are too close to distinguish.Figure 6
Error rates as a function of cluster varianceΣ, where Σ=α0,0.001,0,0.001.


### 3.3. Real Data

Publicly available sequence data for the Human Microbiome Project (HMP) study SRP002462, described as metagenomic sequencing of 16S rDNA from vaginal and related samples from clinical and twin subjects, was downloaded from the NCBI SRA database [10]. The downloaded sets of data correspond to two separate submissions: SRA169809 (1608/1608 samples were downloaded) and SRA273234 (34/133 samples were downloaded), for a total of 1642 samples.The SFF files were processed and cleaned using the microbial community analysis software mothur [11], based on a standard protocol developed for 454 sequence data processing and quality control [12]. The dissimilarities between the samples were calculated using the Clayton-Yue dissimilarity measure. The data was subsampled to 5000 sequences per sample (this step results in dropping out 136 samples that had less than 5000 reads in total) 500 times to produce the distance file, which was used to calculate principal coordinates. Figure 7 shows the graph of ~1500 data points after PCoA. After implementing the 2-way k-means algorithm [13], we initialized with k-means k=5 and ran 2-way k-means for 5 rounds on the data.Figure 7
1500 data points graphed after PCoA.Unfortunately, the nonlinear arches between the clusters pushed the cluster representatives slightly outside the clusters. Nonetheless, the algorithm was still an improvement overk-means. We note that after k-means, the 2-way objective had a value of 108.0 while our 2-way k-means algorithm converged on an objective of ≈51.0 after 5 rounds. Additionally, the algorithm gives us a characterization of the samples lying between two clusters. The results can be seen in Figure 8.Figure 8
~1500 data points after 5 rounds of 2-wayk-means. The dark points are closer to the cluster representatives while the lighter points are more between them. The cluster representatives are the larger blue points that are slightly outside the clusters (compensating for the nonlinear arches between clusters).


### 3.4. Discussion

We first get the most abundant operational taxonomic unit (OTU) in each sample (down to the genus level) and the closest cluster assignment for each sample. We use this to observe which OTUs are the most common to each cluster. We can find the closest sample to each data point by simply taking theargmaxu for each data point xi.From Table1, we see that four of the five clusters have a unique and most abundant OTU, while cluster c3 has a variety of abundant types. Aside from the top four OTUs, separating the data into discrete clusters obscures how the rest of the OTUs can be characterized.Table 1
The most abundant OTU per cluster. Because the 16S rDNA data maps multiple sequences to the same genus level, we use subscripts to denote different OTUs with the same genus.
Family; genus
c
1
c
2
c
3
c
4
c
5
Lactobacillaceae;Lactobacillus1
5
534
17
0
1
Lactobacillaceae;Lactobacillus2
0
0
3
269
0
Bifidobacteriaceae;Gardnerella
227
0
13
0
1
Lachnospiraceae; unclassified1
0
0
4
0
150
Lactobacillaceae;Lactobacillus3
2
0
44
2
0
Lactobacillaceae;Lactobacillus4
2
4
32
1
1
Leptotrichiaceae;Sneathia1
8
0
33
0
0
Prevotellaceae;Prevotella1
2
0
30
0
0
Prevotellaceae;Prevotella2
3
0
16
0
1
Unclassified; unclassified2
2
1
15
0
2
Prevotellaceae;Prevotella3
0
0
4
0
1
Leptotrichiaceae;Sneathia2
1
0
2
0
1
Lachnospiraceae; unclassified3
0
0
1
0
3
Streptococcaceae;Streptococcus1
0
0
18
0
0
Veillonellaceae; unclassified4
0
0
0
0
0
Streptococcaceae;Streptococcus2
0
1
15
0
0
Mycoplasmataceae;Mycoplasma
0
1
7
0
0
Bifidobacteriaceae;Bifodobacterium
0
0
9
0
0
Fusobacteriaceae;Fusobacterium
0
0
7
0
0
Enterobacteriaceae; unclassified5
0
0
8
0
0By using each data point’s cluster-pair assignment, we further separate the data intok2−k clusters. Let cjj′ designate the data points that are between clusters j and j′ but are nearer to cluster j than cluster j′. We take the most abundant OTUs in each sample and the cluster pair for each sample. We can then find the most abundant OTUs for each cluster pair.Table2 shows the structure of the most abundant OTU types for each 2-way cluster cjj′ defined before. Once again, we find that clusters c1j, c2j, c4j, and c5j are all dominated by the same single OTU from before. Yet, observing clusters c3j′ provides us with a more in-depth understanding of the diverse cluster c3.Table 2
The most abundant OTUs per cluster pair.
Genus
c
12
c
13
c
14
c
15
c
21
c
23
c
24
c
25
c
31
c
32
c
34
c
35
c
41
c
42
c
43
c
45
c
51
c
52
c
53
c
54
Lactobacillus
1
5
0
0
0
113
69
57
295
0
15
2
0
0
0
0
0
0
1
0
0
Lactobacillus
2
0
0
0
0
0
0
0
0
0
0
3
0
14
60
27
168
0
0
0
0
Gardnerella
75
32
95
25
0
0
0
0
10
0
3
0
0
0
0
0
1
0
0
0
Unclassified1
0
0
0
0
0
0
0
0
0
0
1
3
0
0
0
0
54
23
23
50
Lactobacillus
3
0
2
0
0
0
0
0
0
6
7
31
0
0
0
2
0
0
0
0
0
Lactobacillus
4
1
1
0
0
0
4
0
0
4
4
23
1
0
1
0
0
1
0
0
0
Sneathia
1
0
8
0
0
0
0
0
0
22
2
4
5
0
0
0
0
0
0
0
0
Prevotella
1
0
2
0
0
0
0
0
0
9
4
17
0
0
0
0
0
0
0
0
0
Prevotella
2
0
0
0
3
0
0
0
0
15
0
0
1
0
0
0
0
0
0
1
0
Unclassified2
2
0
0
0
0
1
0
0
8
6
0
1
0
0
0
0
1
0
1
0
Prevotella
3
0
0
0
0
0
0
0
0
3
0
0
1
0
0
0
0
1
0
0
0
Sneathia
2
0
1
0
0
0
0
0
0
1
1
0
0
0
0
0
0
0
0
1
0
Unclassified3
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
1
1
1
0
Streptococcus
2
0
0
0
0
0
0
0
0
2
2
14
0
0
0
0
0
0
0
0
0
Unclassified4
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Streptococcus
2
0
0
0
0
0
1
0
0
2
2
11
0
0
0
0
0
0
0
0
0
Mycoplasma
0
0
0
0
1
0
0
0
5
2
0
0
0
0
0
0
0
0
0
0
Bifodobacterium
0
0
0
0
0
0
0
0
1
1
7
0
0
0
0
0
0
0
0
0
Fusobacterium
0
0
0
0
0
0
0
0
1
2
4
0
0
0
0
0
0
0
0
0
Unclassified5
0
0
0
0
0
0
0
0
0
2
5
1
0
0
0
0
0
0
0
0Interestingly, we see that the makeups ofc31, c32, c34, and c35 are remarkably different. We immediately see that the top four OTUs are all predominantly contained in cluster pairs that include their single main cluster. In addition, we notice that the samples with abundant Sneathia1, Prevotella2, and unclassified types are predominantly contained in c31. c32 contains samples with a variety of abundant OTUs. Lactobacillus3, Lactobacillus4, Prevotella1, Streptococcus1, Streptococcus2, and Bifodobacterium are abundant in samples that are predominantly contained in c34. Finally, almost no samples are in the cluster pair c35, aside from a few Sneathia1 types.In this way, 2-wayk-means also opens up a wealth of information on the relationships between samples. In particular, it now makes more sense to characterize the samples as being in 6 different clusters: c1, c2, c31, c34, and c5. We also see that certain clusters have mixed relationships, while others have almost no interaction. Without 2-way k-means, this would not be immediately obvious.


## 3.1. Example Run for 2-Wayk-Means

We find it illuminating to demonstrate the performance of 2-wayk-means versus vanilla k-means on a cartoon example.In Figure1, we simulated n=500 data points in ℝ2 from three clusters, with respective means μ1=0,0,μ2=1,1,andμ3=2,0 and covariance matrices Σ=0.001I. Data points are drawn into pairwise clusters by choosing two cluster representatives without replacement from the following prior probabilities:
(14)Pc1=0.2,Pc2=0.5,Pc3=0.3.Figure 1
n
=
500 simulated data points. The white diamonds are cluster centers for three simulated clusters (red cluster, bottom left; green cluster, top; and blue cluster, bottom right). Points are colored as a linear combination of the clusters they lie between (according to u).We initialize the cluster representatives with vanillak-means. Vanilla k-means achieves the results in Figure 2. Statistics for vanilla k-means are given as follows:
(15)obj:41.1025,errfx:0.112,errμ:0.4039,erru:0.2435.Figure 2
n
=
500 simulated data points after k-means. The white diamonds are cluster centers for three simulated clusters (red cluster, bottom left; green cluster, top; and blue cluster, bottom right). White triangles are cluster centers determined by u means. Colors are u values determined by k-means.k-means predicts the 2-way cluster assignments of ≈11% points incorrectly due to skewing the cluster means toward the middle of the graph. 2-way k-means, however, significantly improves on all error rates. After 10 rounds of 2-way k-means, we achieve the results in Figure 3.
(16)obj:5.822,errfx:0.032,errμ:0.0495,erru:0.2084.Figure 3
n
=
500 simulated data points after 10 rounds of 2-way k-means. White diamonds are cluster centers determined by 2-way k-means (10 rounds). Colors are u values determined by 2-way k-means (10 rounds).For every statistic, the results are clearly an improvement on standardk-means. The ≈3% error rate on cluster assignment still exists because 2-way k-means points closest to cluster representatives may be assigned to an incorrect secondary cluster.


## 3.2. Benchmarks

### 3.2.1. Sparsity (Avg. of 10 Trials, 10 Rounds Each)

Our sparsity test was conducted by keeping cluster prior probabilities and cluster centersμ constant while varying the number of data points (ratio of α means n=500α). From Figure 4, we see that the algorithm performs consistently well under a variety of conditions, but too few data points can hurt performance to an extent.Figure 4
Error rates whenn=α500, and cluster priors and centers are fixed.


### 3.2.2. Cluster Separation (Avg. of 10 Trials, 10 Rounds Each)

We test the error rate as a function of the Euclidean distances ofμ (ratio of α means μ1=α0,0,μ2=α1,1,andμ3=α2,0). From the results in Figure 5, we can see that a certain threshold is required for proper performance of the algorithm. This makes sense, as when α=0.01, the clusters are almost on top of each other and difficult to distinguish. Additionally, as the cluster centers are moved farther apart, the ℓ2 norm between the cluster representative determined by the algorithm and the actual cluster representative increases (but this is to be expected).Figure 5
Error rates as a function of the Euclidean distances ofμ, where μ1=α0,0,μ2=α1,1,andμ3=α2,0.


### 3.2.3. Variance (Avg. of 10 Trials, 10 Rounds Each)

We increase the variance of the clusters while fixing cluster prior probabilities, data points, and cluster centers (ratio ofα means Σ=α0,0.0001,0,0.001). From the results in Figure 6, we can see that large variance hurts proper performance of the algorithm. Analogous to cluster separation, as when α=100, the clusters are too close to distinguish.Figure 6
Error rates as a function of cluster varianceΣ, where Σ=α0,0.001,0,0.001.


## 3.2.1. Sparsity (Avg. of 10 Trials, 10 Rounds Each)

Our sparsity test was conducted by keeping cluster prior probabilities and cluster centersμ constant while varying the number of data points (ratio of α means n=500α). From Figure 4, we see that the algorithm performs consistently well under a variety of conditions, but too few data points can hurt performance to an extent.Figure 4
Error rates whenn=α500, and cluster priors and centers are fixed.


## 3.2.2. Cluster Separation (Avg. of 10 Trials, 10 Rounds Each)

We test the error rate as a function of the Euclidean distances ofμ (ratio of α means μ1=α0,0,μ2=α1,1,andμ3=α2,0). From the results in Figure 5, we can see that a certain threshold is required for proper performance of the algorithm. This makes sense, as when α=0.01, the clusters are almost on top of each other and difficult to distinguish. Additionally, as the cluster centers are moved farther apart, the ℓ2 norm between the cluster representative determined by the algorithm and the actual cluster representative increases (but this is to be expected).Figure 5
Error rates as a function of the Euclidean distances ofμ, where μ1=α0,0,μ2=α1,1,andμ3=α2,0.


## 3.2.3. Variance (Avg. of 10 Trials, 10 Rounds Each)

We increase the variance of the clusters while fixing cluster prior probabilities, data points, and cluster centers (ratio ofα means Σ=α0,0.0001,0,0.001). From the results in Figure 6, we can see that large variance hurts proper performance of the algorithm. Analogous to cluster separation, as when α=100, the clusters are too close to distinguish.Figure 6
Error rates as a function of cluster varianceΣ, where Σ=α0,0.001,0,0.001.


## 3.3. Real Data

Publicly available sequence data for the Human Microbiome Project (HMP) study SRP002462, described as metagenomic sequencing of 16S rDNA from vaginal and related samples from clinical and twin subjects, was downloaded from the NCBI SRA database [10]. The downloaded sets of data correspond to two separate submissions: SRA169809 (1608/1608 samples were downloaded) and SRA273234 (34/133 samples were downloaded), for a total of 1642 samples.The SFF files were processed and cleaned using the microbial community analysis software mothur [11], based on a standard protocol developed for 454 sequence data processing and quality control [12]. The dissimilarities between the samples were calculated using the Clayton-Yue dissimilarity measure. The data was subsampled to 5000 sequences per sample (this step results in dropping out 136 samples that had less than 5000 reads in total) 500 times to produce the distance file, which was used to calculate principal coordinates. Figure 7 shows the graph of ~1500 data points after PCoA. After implementing the 2-way k-means algorithm [13], we initialized with k-means k=5 and ran 2-way k-means for 5 rounds on the data.Figure 7
1500 data points graphed after PCoA.Unfortunately, the nonlinear arches between the clusters pushed the cluster representatives slightly outside the clusters. Nonetheless, the algorithm was still an improvement overk-means. We note that after k-means, the 2-way objective had a value of 108.0 while our 2-way k-means algorithm converged on an objective of ≈51.0 after 5 rounds. Additionally, the algorithm gives us a characterization of the samples lying between two clusters. The results can be seen in Figure 8.Figure 8
~1500 data points after 5 rounds of 2-wayk-means. The dark points are closer to the cluster representatives while the lighter points are more between them. The cluster representatives are the larger blue points that are slightly outside the clusters (compensating for the nonlinear arches between clusters).


## 3.4. Discussion

We first get the most abundant operational taxonomic unit (OTU) in each sample (down to the genus level) and the closest cluster assignment for each sample. We use this to observe which OTUs are the most common to each cluster. We can find the closest sample to each data point by simply taking theargmaxu for each data point xi.From Table1, we see that four of the five clusters have a unique and most abundant OTU, while cluster c3 has a variety of abundant types. Aside from the top four OTUs, separating the data into discrete clusters obscures how the rest of the OTUs can be characterized.Table 1
The most abundant OTU per cluster. Because the 16S rDNA data maps multiple sequences to the same genus level, we use subscripts to denote different OTUs with the same genus.
Family; genus
c
1
c
2
c
3
c
4
c
5
Lactobacillaceae;Lactobacillus1
5
534
17
0
1
Lactobacillaceae;Lactobacillus2
0
0
3
269
0
Bifidobacteriaceae;Gardnerella
227
0
13
0
1
Lachnospiraceae; unclassified1
0
0
4
0
150
Lactobacillaceae;Lactobacillus3
2
0
44
2
0
Lactobacillaceae;Lactobacillus4
2
4
32
1
1
Leptotrichiaceae;Sneathia1
8
0
33
0
0
Prevotellaceae;Prevotella1
2
0
30
0
0
Prevotellaceae;Prevotella2
3
0
16
0
1
Unclassified; unclassified2
2
1
15
0
2
Prevotellaceae;Prevotella3
0
0
4
0
1
Leptotrichiaceae;Sneathia2
1
0
2
0
1
Lachnospiraceae; unclassified3
0
0
1
0
3
Streptococcaceae;Streptococcus1
0
0
18
0
0
Veillonellaceae; unclassified4
0
0
0
0
0
Streptococcaceae;Streptococcus2
0
1
15
0
0
Mycoplasmataceae;Mycoplasma
0
1
7
0
0
Bifidobacteriaceae;Bifodobacterium
0
0
9
0
0
Fusobacteriaceae;Fusobacterium
0
0
7
0
0
Enterobacteriaceae; unclassified5
0
0
8
0
0By using each data point’s cluster-pair assignment, we further separate the data intok2−k clusters. Let cjj′ designate the data points that are between clusters j and j′ but are nearer to cluster j than cluster j′. We take the most abundant OTUs in each sample and the cluster pair for each sample. We can then find the most abundant OTUs for each cluster pair.Table2 shows the structure of the most abundant OTU types for each 2-way cluster cjj′ defined before. Once again, we find that clusters c1j, c2j, c4j, and c5j are all dominated by the same single OTU from before. Yet, observing clusters c3j′ provides us with a more in-depth understanding of the diverse cluster c3.Table 2
The most abundant OTUs per cluster pair.
Genus
c
12
c
13
c
14
c
15
c
21
c
23
c
24
c
25
c
31
c
32
c
34
c
35
c
41
c
42
c
43
c
45
c
51
c
52
c
53
c
54
Lactobacillus
1
5
0
0
0
113
69
57
295
0
15
2
0
0
0
0
0
0
1
0
0
Lactobacillus
2
0
0
0
0
0
0
0
0
0
0
3
0
14
60
27
168
0
0
0
0
Gardnerella
75
32
95
25
0
0
0
0
10
0
3
0
0
0
0
0
1
0
0
0
Unclassified1
0
0
0
0
0
0
0
0
0
0
1
3
0
0
0
0
54
23
23
50
Lactobacillus
3
0
2
0
0
0
0
0
0
6
7
31
0
0
0
2
0
0
0
0
0
Lactobacillus
4
1
1
0
0
0
4
0
0
4
4
23
1
0
1
0
0
1
0
0
0
Sneathia
1
0
8
0
0
0
0
0
0
22
2
4
5
0
0
0
0
0
0
0
0
Prevotella
1
0
2
0
0
0
0
0
0
9
4
17
0
0
0
0
0
0
0
0
0
Prevotella
2
0
0
0
3
0
0
0
0
15
0
0
1
0
0
0
0
0
0
1
0
Unclassified2
2
0
0
0
0
1
0
0
8
6
0
1
0
0
0
0
1
0
1
0
Prevotella
3
0
0
0
0
0
0
0
0
3
0
0
1
0
0
0
0
1
0
0
0
Sneathia
2
0
1
0
0
0
0
0
0
1
1
0
0
0
0
0
0
0
0
1
0
Unclassified3
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
1
1
1
0
Streptococcus
2
0
0
0
0
0
0
0
0
2
2
14
0
0
0
0
0
0
0
0
0
Unclassified4
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Streptococcus
2
0
0
0
0
0
1
0
0
2
2
11
0
0
0
0
0
0
0
0
0
Mycoplasma
0
0
0
0
1
0
0
0
5
2
0
0
0
0
0
0
0
0
0
0
Bifodobacterium
0
0
0
0
0
0
0
0
1
1
7
0
0
0
0
0
0
0
0
0
Fusobacterium
0
0
0
0
0
0
0
0
1
2
4
0
0
0
0
0
0
0
0
0
Unclassified5
0
0
0
0
0
0
0
0
0
2
5
1
0
0
0
0
0
0
0
0Interestingly, we see that the makeups ofc31, c32, c34, and c35 are remarkably different. We immediately see that the top four OTUs are all predominantly contained in cluster pairs that include their single main cluster. In addition, we notice that the samples with abundant Sneathia1, Prevotella2, and unclassified types are predominantly contained in c31. c32 contains samples with a variety of abundant OTUs. Lactobacillus3, Lactobacillus4, Prevotella1, Streptococcus1, Streptococcus2, and Bifodobacterium are abundant in samples that are predominantly contained in c34. Finally, almost no samples are in the cluster pair c35, aside from a few Sneathia1 types.In this way, 2-wayk-means also opens up a wealth of information on the relationships between samples. In particular, it now makes more sense to characterize the samples as being in 6 different clusters: c1, c2, c31, c34, and c5. We also see that certain clusters have mixed relationships, while others have almost no interaction. Without 2-way k-means, this would not be immediately obvious.


## 4. Conclusion

The complexity of microbial populations is unfolding as microbiome data becomes increasingly available. Yet, standard methodologies oversimplify microbial compositions by pigeonholing them into discrete clusters. This paper further refines the models for microbial abundance across groups of samples. We allow samples to be presented as a weighted average of two clusters, rather than belonging to only one. This may be motivated biologically, as the sample often reflects a mixture of two sources of microbiota, each well represented by a cluster. An alternative explanation is that the averaged sample represents an intermediate, potentially temporary state of the microbial composition, between the more stable ones represented by the clusters themselves.Technically, we formalize this model as a generalization ofk-means. We derive a simple algorithm to infer such a structure and validate its benchmarks on simulated data.Applying our algorithm to real data from the Human Vaginal Microbiome Project provides empirical support to the 2-way model. We showed that while most of the samples lie in six clusters: four well-defined clusters and two subclusters. Furthermore, while previously, a sizable fraction of samples in between clusters was ignored, the 2-way model characterized the entire distribution. Using 2-wayk-means, we can tell that a large portion of the previously unclustered samples, which lie in between two clusters, contains shared properties. In addition, we see that certain clusters have mixed relationships, while others have almost no interaction.


## 5. Further Research

In addition, this paper leaves several open questions and opportunities for further research:(i)
How can we efficiently characterize a 2-way distribution with nonspherical covariance matrices?(ii)
How can we efficiently characterize ak-way distribution?(iii)
How can we efficiently characterize a 2-way distribution with nonlinear paths between cluster representatives?Addressing these questions will further help us understand the composition of microbial populations.



---

*Source: 5284145-2017-09-05.xml*