Development of the test statistic

Suppose two samples R and C, each with size n_R and n_C, were drawn independently from population $$ \mathfrak{P}$$ of human microbiomes. In HMAS, these samples may correspond to the sets of microbiomes of volunteers in the reference group ($$ {\left\{{\overrightarrow{y}}_i^R\right\}}_{i=1}^{n_R}$$) and case group ($$ {\left\{{\overrightarrow{y}}_i^C\right\}}_{i=1}^{n_C}$$), respectively; each microbiome is a vector of which the elements represent relative frequencies $$ y_j^S$$ of j = 1,2,⋯,t, where t is the number of operational taxonomic units (OTUs) at different taxonomic levels from phylum to strain and S∈{R, C}. The summary statistics are vectors $$ \overrightarrow{r}$$ and $$ \overrightarrow{c}$$, of which the elements are the mean frequencies r_j and c_j of the OTUs in the pooled data obtained from R and C, respectively. Now consider the microbiome of an individual q, $$ {\overrightarrow{y}}^q$$, of which the elements are simple binary measures of presence/absence of OUT j (i.e., $$ y_j^q\in \left\{0,1\right\}$$), and suppose we want to determine whether $$ {\overrightarrow{y}}^q$$ is a member of R or C using $$ \overrightarrow{r}$$ and $$ \overrightarrow{c}$$. First, calculate a distance d for OUT j using the absolute difference between $$ y_j^q$$ and c_j and the absolute difference between $$ y_j^q$$ and r_j as follows:

Assuming that OTUs are independent and invoking the central limit theorem for the large number of OTUs (t>50) examined in the HMAS, z-score of d_j across all OTUs will follow the standard normal distribution, N(0, 1).

Since the t = max(t_R, t_C, t_y) is large, the variance of d can be estimated reliably by the sample variance s². The test statistic Z was inspired by Homer et al. [14] and Braun et al. [16]. The authors used a similar test statistic calculated from the single nucleotide polymorphism (SNP) genotyping data to identify an individual’s genotype in GWAS samples. I modified the original test statistic in order to use the binary taxonomic data in HMAS. Because an individual microbiome randomly drawn from population $$ \mathfrak{P}$$ should be equally distant from R and C, μ₀ was presumably expected to be zero, i.e., N(0, 1) was a putative null distribution. Thus, under the null hypothesis $$ H_0:Z=0({\overrightarrow{y}}^q$$ is a random draw from $$ \mathfrak{P}$$), the alternative hypothesis $$ H_R:Z<0({\overrightarrow{y}}^q$$ is a member of R) or $$ H_C:Z>0({\overrightarrow{y}}^q$$ is a member of C) can be tested with an appropriate significance level α. For example, Z>1.65 rejects H₀ in favor of H_C at α = 0.05 (one-tailed test).

Distribution simulations

I examined the feasibility of the test statistic Z in identifying the presence of an individual’s microbiome in an HMAS sample with simulated datasets. Suppose that the OTUs correspond to species-level affiliations of microorganisms. Then, the presence of species j in $$ \mathfrak{P}$$ is a Bernoulli random variable with parameter p_j. To model the probability distribution of p_j $$ \left(\overrightarrow{p}\right)$$, I used four different Beta(π₁, π₂) distributions. For $$ \mathfrak{P}$$ with a small number of high-frequency species, π₁ = 0.1 and π₂ = 1.0 were assumed, and for $$ \mathfrak{P}$$ with a large number of high-frequency species, π₁ = 1.0 and π₂ = 0.1 were assumed. For $$ \mathfrak{P}$$ with a high number of high-frequency species as well as a high number of low-frequency species, π₁ = 0.1 and π₂ = 0.1 were assumed, which might be more realistic than the above two distributions in that the high-frequency species may correspond to constitutional or autochthonous prokaryotic populations and the low-frequency species may correspond to opportunistic or heterochthonous prokaryotic populations across individual microbiomes. In addition, a Beta(1, 1) (uniform) distribution was also used, which represents our ignorance of the distribution of p_j according to the principle of indifference.

Because the individual microbiome in R or C is a set of t random draws of y_j~Bernoulli(p_j), i.e., $$ \overrightarrow{y}\sim Bernoulli\left(\overrightarrow{p}\right)$$, the summary statistics for samples R and C were simulated using $$ \overrightarrow{r}\sim Binomial(n_R,\overrightarrow{p})/n_R$$ and $$ \overrightarrow{c}\sim Binomial(n_C,\overrightarrow{p})/n_C$$, respectively. To construct density curves of true positives for samples R and C, random draws of $$ {\overrightarrow{y}}^{R+}$$ and $$ y_j^{C+}$$ from samples R and C were used to calculate the test statistics Z^R+ and Z^C+ $$ (n_Z^{R+}=n_Z^{C+}=100)$$, respectively. Density curves were estimated using the Gaussian kernel density estimator. To construct density curves of the true null distribution, random draws of $$ {\overrightarrow{y}}^P$$ from $$ \mathfrak{P}$$ were used to calculate test statistics Z^P $$ (n_Z^P=100)$$. Python code for the simulation is available at https://colab.research.google.com/drive/1dOZi8OSo5qHmF7JPyGAP1I_BQdBVSSiC?usp=sharing.

Overview of simulation results

A simulation study was started with the number of OTUs t =2,000, which roughly reflects the number of species (including uncultivated candidate species) in the human gut microbiome [17], and with n_R = n_C = 10, which could correspond to small-scale HWAS, under the assumption that p_j follows uniform distribution (Fig 1). The null distribution estimated using Z^P was very close to N(0, 1) and crisply separated from the distributions of true positives (Z^R+ and Z^C+), indicating the feasibility of the microbiome-based identification. The distributions of true positives moved toward the null distribution with increasing n or with decreasing t but moved away from the null distribution with decreasing n or with increasing t. Similar distribution patterns were observed for all the underlying distributions of p_j assumed (S1–S3 Figs).

Fig 1

Distributions of the test statistic Z under the assumption that the population OTU frequencies follow a uniform (Beta(1, 1)) distribution.

Density curves for true positives of samples R (Z^R+) and C (Z^C+) are denoted by green and red lines, respectively. Density curves of simulated null distribution and standard normal distribution are denoted by black and gray lines, respectively. Single and double asterisks represent type II error probabilities β<0.05 and β<0.01, respectively.

For numerical interpretation of the simulation results, I focused on the probability of type II error at α = 0.05 (β^R = P(Z^R+>z_α|H_R) or (β^C = P(Z^C+<z_α|H_C); the probability that the test statistic is not in the H₀ rejection range, given that the alternative hypothesis is true). The power of the test (1−β) is important to a privacy attacker because it would be especially difficult for the attacker with a high β to determine whether the individual under investigation belongs to a particular HMAS sample due to the high false negative rate. α level critical values were obtained from percentiles of Z^P distribution because the true null distribution might diverge from N(0, 1). As expected in the density curves, β was far less than 0.01 for the experiments with a large t or small n (S1–S4 Tables).

To formulate the guidelines for human microbiome privacy, the initial simulation study was expanded for virtual HMAS samples with n_R = n_C = 10, 20,…,100,…,1000 and with t = 20, 30,…,100,…,1000,…,20000. The method was in general slightly more powerful for p_j under the uniform distribution (Fig 2) than for p_j under other beta distributions (S4 and S5 Figs). In the resulting contour plots, the yellowish area represents higher β (i.e., lower test power); thus, the HMAS samples located in the dark blue area are considered to be vulnerable to privacy risk. The power of the test decreased notably with increasing HMAS sample size (n) or with a decreasing number of OTUs (t) from which the HMAS summary statistics are calculated, and it is possible to notice a borderline where β decreases considerably, which could help in assessing the privacy risk of the HMAS data.

Fig 2

Contour plot representations of the type II error probabilities (β) for true positives of samples R and C under the assumption that the population OTU frequencies follow a uniform (Beta(1, 1)) distribution.

Sample size and the number of OTUs are log-scaled. Dotted line denotes suggested minimal guidelines for HMAS privacy.

Effect of sample size

Samples of a large size approximate the population $$ \mathfrak{P}$$ because $$ {\lim }_{n_{R\to \infty }}{r}_j=p_j$$ or $$ {\lim }_{n_{C\to \infty }}{c}_j=p_j$$ (hence, $$ {\lim }_{{n_R,n}_{C,\to \infty }}{d}_j=0$$), and the distribution of Z^R+ or Z^C+ would overlap with the null distribution as shown in the selected density curves, indicating that a large sample size will make the classifying method ineffective. Contrarily, for the samples of a small size, $$ \overline{d}$$ significantly deviates from μ₀ even if the difference between r_j and p_j or between c_j and p_j is very small. Note that the sample size does not increase with the number of technical replicates used in HMAS data, since the sample points in technical replicates are not independent.

For an unequal sample size (e.g., n_R = 1,000 and n_C = 10), the samples with a larger size would approximate the null distribution, but the test power for identifying a microbiome in a sample with a smaller size would not be affected (S6 Fig). Considering that taxonomic profiles of individual microbiomes in the case sample could be more homogeneous than those in the control sample, the effective size of the case sample might be smaller than the census size of the case sample. This would result in much higher statistical power for testing whether an individual is a member of the case sample, which is a primary interest of the attackers.

Effect of the number of OTUs

With the sample size fixed, an increased number of OTUs would increase the test power. This situation can happen when the HMAS data contains taxonomic profiles with a resolution finer than species level. Under the arbitrary speculation that each species is comprised of ca. 10 strains (t = 20,000), the method was very powerful (β≪0.01) even for a moderate sample size. This is because the denominator of the test statistic shrinks by increasing t compared to the sample mean $$ \overline{d}$$, which results in a large Z^R+ or Z^C+. In the same vein, a reduced number (e.g., t = 20) of OTUs (roughly corresponding to near phylum-level) would decrease the power of the test.

Effect of correlation among OTUs

The occurrence of many OTUs in the human microbiome could be correlated, since the microbiome itself is an ecological community [18]. If the correlation among OTUs is significant, the violation of the assumption that OTUs are independent could change the test power. I analytically evaluated the effect of the OUT correlation on test power. If OTUs are not independent, the variance of $$ \overline{d}$$ includes covariance (Cov) terms as follows:

Note that t_e = t if $$ \tilde{\rho }=0$$, 1≤t_e<t if $$ 0<\tilde{\rho }\le 1$$, and t<t_e if $$ -1/(t-1)\le \tilde{\rho }<0$$ (t_e→∞ as $$ \tilde{\rho }\to -1/(t-1)$$). For example, if three of 10 OTUs show a strong positive correlation, two of these OTUs cannot contribute to the number of OTUs as equally as other independent OTUs. Thus, the effect of a positive $$ \tilde{\rho }$$ could be similar to the effect of a decreased number of OTUs, resulting in decreased test power as in the case of linkage disequilibrium among SNPs [16]. However, $$ \tilde{\rho }$$ can also be negative if negative interactions among OTUs are more prevalent than positive interactions among OTUs, which could subsequently increase test power. Even very small negative correlations (e.g., $$ \tilde{\rho }=-0.001$$) among a moderate number of OTUs (t = 1000) would increase t_e to 10⁶, while very small positive correlations (e.g., $$ \tilde{\rho }=+0.001$$) decrease t_e to 500. The negative correlation overwhelmingly affects t_e because t_e is a reciprocal function of $$ \tilde{\rho }$$ when t is fixed (S7 Fig). Nonetheless, because $$ \tilde{\rho }$$ cannot be measured from summary statistics unless provided or assumed, the effect of $$ \tilde{\rho }$$ should be investigated more comprehensively in the future.

Additional considerations

I used the binary vector as a query microbiome because it was considered that the binary data is less prone to experimental variations. If a frequency vector is used, the $$ \overline{d}$$ could increase because small differences between $$ |y_j^q-r_j|$$ and |$$ y_j^q-c_j$$| can accumulate across a large number of OTUs. However, it would not always result in a larger Z^R+ or Z^C+ because the increased sample variance could counteract the increase in $$ \overline{d}$$. Thus, the use of a frequency vector might not guarantee more test power, rather it could make $$ \overline{d}$$ more prone to errors in OTU frequencies.

It was assumed that R and C were samples drawn independently from the population of human microbiomes ($$ \mathfrak{P}$$). Violation of this assumption could shift the location of the null distribution (μ₀). Suppose that populations underlying samples R and C (denoted by $$ \mathcal{R}$$ and $$ \mathcal{C}$$, respectively) are all different from the population underlying $$ {\overrightarrow{y}}^q\left(\mathfrak{P}\right)$$. If the systemic difference between $$ \mathcal{R}$$ and $$ \mathcal{C}$$ is significant, the difference between summary statistics can deviate from zero, i.e., E(Δ) = E(r−c)≠0. Let $$ d_j^{*}$$ be the distance metric that includes Δ term as follows:

Then, according Braun et al. [16] and using $$ P(y_j^q=0|p_j)=1-p_j$$ and $$ P(y_j^q=1|p_j)=p_j$$, the expected value of the distance metric under the null hypothesis is as follows:

Because −1≤E(Δ)≤1 and $$ 0<E\left(p\right)\le 1,{\mu }_0^{*}$$ ranges from -1 to 1 which is quite small compared to the distance between the distributions of Z^P and Z^R+ or between the distributions of Z^P and Z^C+ in the cases where the test power is very high (e.g., t≥2000 and n = 10 or t≥20000 and n = 100) (Fig 1). Moreover, as E(p) deviates from its two extreme values (0 or 1) which are highly unrealistic in microbiome composition, the 2E(Δ)E(p) term counteracts the deviation of $$ {\mu }_0^{*}$$ from zero (i.e., $$ {\mu }_0^{*}\to 0$$ as E(p)→0.5). Thus, I considered that the violation of the assumption that underlying populations are all different would not alter the main results of this study. In fact, if R (reference/control group) and C (case/treatment group) are assumed to be samples drawn from $$ \mathfrak{P}$$ and $$ \mathcal{C}$$, respectively, the test statistic would become very similar to that used for a two-tailed z-test (or t-test) which can be used to test the null hypothesis that $$ y_j^q$$ is not a member of C.

Although it might be virtually impossible for the HMAS community to develop an almighty shield that protects data against all possible privacy breaching methods, we do not need to overreact to HMAS privacy concerns since the taxonomic profile of the individual’s microbiome could not be a permanent identifier, while a subset of an individual’s microbiome may endure for most of the lifetime. Nonetheless, the HMAS community should have guidelines for reducing the privacy risk, which could be kept minimal in order not to obstruct our understanding of the human microbiome.

Concluding remarks

This study showed the possibilities of privacy breaches using the summary statistics drawn from HMAS data. The key finding was that the publication of species composition data obtained from small-scale HMAS (e.g., t≈1,000 and n_R = n_C = 10) easily exposes the privacy of the victims. The HMAS samples with moderate size (n_R = n_C≈100) and t≈1000 were on the vicinity of the borderline. I suggest these figures as a basis for HMAS data release policy while acknowledging a sophisticated test statistic that employs better distance metric along with the violation of underlying assumptions could improve the power of the privacy breaching methods. I propose minimal guidelines for HMAS data release as follows: i) increase the sample size as much as manageable, ii) do not publish species composition data even in the form of summary statistics if the sample size is less than 10, iii) avoid publishing subspecies- or strain-level data unless the sample size is far larger than 100. I also suggest that HMAS researchers evaluate the privacy risk at the time of study design using the simulation method presented in this study. Because the test statistic used was relatively robust against the variations in the models for the distribution of p_j, HMAS researchers can simply use uniform distribution and their intended t, n_R, and n_C to estimate the ‘minimal’ power of the privacy-breaching method.

This study focused only on attribute disclosure attacks under the summary statistic scenario [15]. However, HMAS data privacy concerns are not limited to summary statistics; privacy breaches can occur in various manners at any stage of HMAS data management as described by Wagner et al. [12]. Thus, it is timely that the HMAS community begins comprehensive discussions on HMAS data privacy risks and developing privacy-preserving algorithms for data storage and release.