Suppose that we have abundance measurements from n independent microbiome samples and p OTUs/ASVs, denoted by X = (X ₁, X ₂, … , X _i , … , X _n )^T (1 ≤ i ≤ n), where X _i = (x _i1, x _i2, … , x _ij , … , x _ip )^T (1 ≤ j ≤ p) and x _ij is the (normalized) abundance of the j ^th OTU/ASV in the i ^th sample. Let Y = (y ₁, y ₂, … , y _i , … , y _n )^T (1 ≤ i ≤ n) denote the vector for the outcome variable, such as clinical outcomes and environmental gradients. Additionally, we may have q covariates, such as age and biological sex, which are denoted by Z _n×q = (Z ₁, Z ₂, … , Z _i , … , Z _n )^T, where Z _i = (z _i1, z _i2, … , z _ik , … , z _iq )^T (1 ≤ k ≤ q) are the measurement of the q covariates in the ith sample. Moreover, we may have a rooted phylogenetic tree G capturing the phylogenetic relatedness of the OTUs/ASVs. G has p leaves (terminal vertices with a degree of 1) and one node (an internal vertex with a degree greater than 1) called root. The p leaves correspond to the p OTUs/ASVs while the root is theoretically assumed to be the last common ancestor of all vertices in the phylogenetic tree. In a path connecting a leaf and the root, the vertices closer to the root are regarded as “ancestors” of vertices that are farther from; thus, this ancestral relationship describes the relative closeness of vertices to the root of G. The aim for RFtest is to test the association between Y _n×1 and X _n×p while adjusting Z _n×q .

The tree of life underpins our understanding towards microorganisms (Washburne et al., 2018). Closely related microorganisms share similar biological traits and association signals tend to be clustered with respect to their phylogenetic relationship (Xiao et al., 2017; Xiao et al., 2018a; Xiao et al., 2018b). We therefore aim to utilize the phylogenetic information in the random forest test to improve its power. We incorporate such phylogenetic information by augmenting the OTU/ASV-level abundance data with the abundances of the internal nodes of the phylogenetic tree G. This is achieved by creating an n-by-m feature matrix W _n×m = (w _il ) _n×m for the m internal nodes in G, where the features accumulate the abundance of OTUs/ASVs belonging to the same ancestor in G. As each leaf corresponds to one OTU/ASV in microbiome and there exists exactly one path between each leaf and the root, the total abundance of all OTU/ASV leaves that shares a specific common ancestor or internal node l is well-defined. Thus, we have w i l = ∑ j ∈ A x i j (1) where w _il is the collective abundance of the l ^th internal node of the i ^th sample and A is the set of OTUs/ASVs whose ancestor is the l ^th internal node.

The RFtest uses the generalization error rate estimate (Breiman, 2001) of random forest as a test statistic, and uses permutation to calculate p-values. Specifically, random forest is firstly grown using the “ranger” package (Wright and Ziegler, 2017) in the R platform (Team, 2020) using Y _n×1 as outcome variable and X _n×p and W _n×m as input features, and the observed out-of-bag (OOB) error rate T _obs is used as the test statistic. The OOB error is the average error for each observation calculated using predictions from the trees that do not contain in their respective bootstrap sample. Here, we use the probabilistic prediction for classification and the OOB error is essentially a Brier’s score (Malley et al., 2012). Regression and classification trees are used for continuous and binary Ys, respectively. When there are no covariates, it permutes the outcome Y _n×1 B times and calculates the OOB error rate T ˜ b   ( b = 1 ,   … ,   B ) based on the permuted Y _n×1. The p-value is calculated using: p ‐ value = [ # ( T ˜ b ≤ T o b s ) + 1 ] / ( B + 1 ) (2) where #( T ˜ b ≤ T _obs) is the number of permuted datasets satisfying T ˜ b ≤ T _obs.

When covariates are present, RFtest accommodates covariates using the following steps. Firstly, Y _n×1 is regressed on covariate Z _k (1 ≤ k ≤ q) using linear model if Y is continuous: y i = β ^ 0 + β ^ 1 z i 1 + … + β ^ q z i q + e i = β ^ 0 + ∑ k = 1 q β ^ k z i k + e i (3) and using logistic regression model if Y is binary: logit ( P ( y i = 1 ) ) = β ^ 0 + ∑ k = 1 q β ^ k z i k (4) where β ^ 0 and β ^ k (1 ≤ k ≤ q) are the estimated coefficients, and e _i are regression residuals. Next, for a continuous Y, we generate Y ˜ b using residual permutation. The observed error rate T _obs is calculated based on the input features X _n×p and W _n×m and the adjusted outcome Y _adj = (e ₁, e ₂, … , e _i , … , e _n )^T (1 ≤ i ≤ n). Thereafter, the permutated Y ˜ b =   ( y ˜ 1 b , y ˜ 2 b ,   … ,   y ˜ i b ,   … ,   y ˜ n b ) T is generated by y ˜ i b = e ˜ i b (5) where e ˜ i b is the permutated regression residuals for the ith sample. For a binary covariate Y, Y ˜ b is generated using a (0, 1) random number generator according to adjusted probabilities of logit ( P ( y ˜ i b = 1   |   ∑ i y ˜ i b =   ∑ i y i ) ) = β ^ 0 + ∑ k = 1 q β ^ k z i k (6) where we conditioned on the number of observed cases. Finally, we calculate the error rate T ˜ b under permutation based on Y ˜ b similarly. Consequently, p-value can be obtained using (Eq. 2).

We implemented the random forest test in the package “RFtest” on the R platform, which is available on GitHub (https://github.com/Lujun995/Random-forest-test-RFtest).

Simulations were conducted under various scenarios to study whether RFtest would control type I error rates at desired levels and whether it would be a powerful testing approach compared with competing methods. Instead of using a parametrical model such as the Dirichlet-multinomial model (Chen and Li, 2013), the microbiome data were directly resampled from a large gut microbiome study by Hale et al. (2017). Briefly, the study compared the fecal microbiome profiles of patients with adenomas versus healthy controls. 16s rRNA sequences were analyzed using IM-TORNADO pipeline (Jeraldo et al., 2014), OTUs were clustered at 97% identity, and singletons were removed (Hale et al., 2017). After rarefaction to 20,000 counts per sample, the adenoma dataset contained 439 samples and 2,100 OTUs, where we resampled n = 50 samples, i.e., X _50×p , without replacement for each simulated dataset. We then constructed the outcome variable Y _50×1 under six scenarios, following the strategy by Zhao et al. (2015). Let S denote the set that comprises OTUs associated with Y. We generated the continuous and binary outcome Y = (y ₁, y ₂, … , y _i , … , y ₅₀)^T (1 ≤ i ≤ 50) based on y i = β 0 + z i + β     s c a l e [ ∑ j ∈ S ( x i j ) ] + ε i , (7)

and log it ( P ( y i = 1 ) ) = β 0 + z i + β     s c a l e [ ∑ j ∈ S ( x i j ) ] , (8) where β ₀ is a constant, β is an adjustable effect size, ε _i ∼ N (0, σ ²), and the “scale” function standardizes the data to have mean 0 and standard deviation 1. We used β ₀ = 10 for a continuous Y and β ₀ = 0 for a binary Y, ε _i ∼ N (0, 1).

The first scenario (S0) was used to study the type I error rate of RFtest by setting the effect size β = 0 under three cases, including no covariates [z _i = 0 and ε _i ∼ N (0, 1)], one covariate independent of X (z _i ∼ N (0, 1) and ε _i ∼ N (0, 9)), and one covariate associated with X {z _i = scale[ Σ j ∈ S ( x i j ) ] + N (0, 1) and ε _i ∼ N (0, 9)}, respectively. In the third case, S consisted of OTUs from an abundant lineage S _A , which contributed to 15% of the total OTU number and 21% of the total abundance.

The other five scenarios (S1–S5) were used to evaluate the power of RFtest. No covariates were included (z _i = 0) in these scenarios. In S1, we investigated different signal types (phylogenetically clustered vs. non-phylogenetically clustered) and different signal densities (5% vs. 15%). For phylogenetically clustered signals, the signal OTUs for 5% and 15% densities were from two abundant lineages S _B and S _A , respectively, where S _B was contained in S _A described above and contributed to 5% of the total OTU number and 11% of the total abundance. For non-clustered signals, the signal OTUs were randomly selected and OTUs for 5% density were also contained in those for 15%. We further substituted the term Σ j ∈ S ( x i j ) in Eq. 7 and Eq. 8 with Σ j ∈ S ( x i j / x ⋅ j ¯ ) , where x ⋅ j ¯ = 1 n ∑ i = 1 n ( x i j ) , to avoid several OTUs dominating the overall signal strength.

The scenario S2 was designed to further validate the results of clustered signals in S1 using different lineages. We studied seven disjoint major lineages ( S _I : I = A, C, D, E, F, G, H) in the dataset of Hale et al. (2017), which spanned 80% of the total OTU number and more than 80% of the total abundance. Each lineage possessed 5%–20% of the total OTU number and 1%–40% of the overall abundance. The simulations in this scenario were conducted under β = 2.25 for a binary Y and β = 0.75 for a continuous Y.

The scenario S3 was to evaluate the power of the RFtest when the outcome variable was non-linearly associated with the signal OTUs. We applied a non-linear link function f _link to x _ij . Specifically, y i = β 0 + β     s c a l e [ ∑ j ∈ S f l i n k ( x i j ) ] + ε i (9) for a continuous Y, and log it ( P ( y i = 1 ) ) = β 0 + β     s c a l e [ ∑ j ∈ S f l i n k ( x i j ) ] (10) for a binary Y, where f _link(x _ij ) = log₂(x _ij + 1) (x _ij ≥ 0).

The scenario S4 studied a complex association between Y and X where there was interaction between two sets of signal OTUs. Particularly, for a continuous Y, it was generated via y i = β 0 + β     s c a l e [ ∑ j ∈ S ( x i j ) ] · s c a l e [ ∑ j ' ∈ S ' ( x i j ' ) ] + ε i (11)

and for a binary Y, it was generated using log it ( P ( y i = 1 ) ) = β 0 + β   s c a l e [ ∑ j ∈ S ( x i j ) ] · s c a l e [ ∑ j ' ∈ S ' ( x i j ' ) ] , (12) where β was fixed at 1.33 and 5 for a continuous and binary Y, respectively, and S and S ' were two disjoint sets comprising 15% and 13% of total OTUs, respectively. For phylogenetic signals, we let S = S _A and S ’ = S _C , where S _A had been characterized in S0 and S _C was another major lineage accounting for 12% of the total abundance. For non-phylogenetic signal, the terms ∑ j ∈ S (x _ij ) and ∑ j ' ∈ S '   (x _ij’) in Eq. 11 and Eq. 12 were normalized using ∑ j ∈ S   or   S ( x i j / x ⋅ j ¯ ) , where x ⋅ j ¯ = 1 n ∑ i = 1 n ( x i j ) . The sample size used in this scenario ranged from 50 to 250 as detection of an interaction generally requires a relatively large sample size.

The last scenario (S5) was used to assess the robustness of RFtest to outliers. Firstly, the outcome variable Y was generated according to the procedure in S1, using clustered or non-clustered signals with a density of 15%. Subsequently, the order of OTUs in 0, 1, or 3 samples was randomly shuffled, yielding 0, 1, or 3 outliers; therefore, these outliers would possess distinct microbiome profiles.

The source code of this section is available at GitHub (https://github.com/Lujun995/RFtest-Simulations).

The competing methods include the optimal microbiome regression-based kernel association test (optimal MiRKAT) (version 1.1.1, https://cran.r-project.org/package=MiRKAT) (Zhao et al., 2015), the adaptive microbiome-based sum of powered score test (aMiSPU) (version 1.0, https://cran.r-project.org/package=MiSPU) (Wu et al., 2016) and optimal microbiome-based association test (OMiAT) (version 6.0, https://github.com/hk1785/OMiAT) (Koh et al., 2017). While multiple distance or dissimilarity functions could be used in MiRKAT, we followed the example in the “MiRKAT” package (Zhao et al., 2015) and selected weighted and unweighted UniFrac distance (Lozupone and Knight, 2005; Lozupone et al., 2007) and Bray–Curtis dissimilarities (Bray and Curtis, 1957), which have been widely used in microbiome studies. All the results were averaged over 1,000 simulation runs.

In this section, we aimed to compare the results of RFtest, optimal MiRKAT, aMiSPU, and OMiAT in real-world examples. We re-analyzed the relationship between outcome variables and microbiome profiles in two published datasets. The first example was taken from a study on the throat microbiome (Charlson et al., 2010). That study investigated the effect of smoking on human microbiota in the upper respiratory tract. While detailed information of sample collecting and data processing procedures can be accessed from Charlson et al. (2010), a summary is provided here. Nylon-flocked swabs were taken from the nasopharynx and oropharynx of 62 healthy subjects, including 33 non-smokers and 29 smokers. From each swab, DNA was extracted using the QIAamp DNA Stool Minikit (Qiagen) and the V1–V2 region of the 16S rRNA was amplified. Thereafter, this 16S rRNA was sequenced using a 454 Life Sciences Genome Sequencer FLX instrument (Roche). The sequence reads were denoised (Quince et al., 2009), analyzed using the QIIME pipeline (Caporaso et al., 2010), and clustered into OTUs at 97% similarity using UCLUST (Edgar, 2010).

In the original study (Charlson et al., 2010), the association between smoking and the respiratory tract microbiome was tested by Permutational Multivariate Analysis of Variance (Anderson, 2001). based on weighted and unweighted UniFrac distances (Lozupone and Knight, 2005; Lozupone et al., 2007). A difference in microbial community structure was reported between smokers and non-smokers (p < 0.05). In the present study, we re-analyzed the microbiome data and found consistent results with previous studies (Charlson et al., 2010; Zhao et al., 2015; Wu et al., 2016). When no covariate was considered, the p-value estimated by the RFtest was 0.001 while those of the optimal MiRKAT, the OMiAT, and the aMiSPU were 0.006, 0.008, and 0.009, respectively. When biological sex was included as a confounder, the estimated p-values became 0.002, 0.009, 0.010, and 0.005 for the RFtest, the optimal MiRKAT, the OMiAT, and the aMiSPU, respectively. The RFtest provided more significant p-values in general, while all competing methods rejected the null hypotheses at a significance level of 0.01.

Another relevant example was taken from a study of the distance–decay relationship in microbial ecology (Xue et al., 2021). This relationship can be portrayed as relatedness of microbial communities decreases as their spatial distance increases (Hanson et al., 2012). In brief, surface soil was collected intact from a paddy field in Wenling, Zhejiang Province, China (28°21′ N, 121°15′ E) in November 2017. From the sample, a soil cube (2.0 cm × 2.0 cm × 2.0 cm) was selected and further divided into 4 × 4 × 4 cubes of which each had sides 0.5 cm in length. DNA samples were extracted from these sub-cubes, and the V4–V5 region of the 16S rDNA genes was amplified and subsequently sequenced using an Illumina HiSeq platform. After removal of adaptors and quality control, 16s rDNA sequences were aligned using USEARCH11 (https://www.drive5.com/usearch/) and OTUs were clustered at 97% identity using UPARSE (Edgar, 2013). Finally, the microbial communities were rarefied to 41,752 sequences per sample.

As one of the original findings (Xue et al., 2021), a decreased community similarity, measured by 1 − Bray–Curtis dissimilarity (Bray and Curtis, 1957) between microbial communities, was observed as the spatial distance increased in the 64 sub-cubes (Mantel test, p = 0.001). Herein, we re-examined this distance–decay association using the RFtest via an assessment of microbial changes along each spatial axis of the xyz-coordinate defined in the study of Xue et al. (2021). We found a similar result that the microbiome was associated with the x- and y-axes, and p-values by the RFtest were 0.001, 0.001, and 0.310 for the x-, y-, and z-axes, respectively. Those of the optimal MiRKAT were 0.011, 0.001, and 0.618, respectively; those of the OMiAT were 0.001, 0.001, and 0.265; and those of aMiSPU were 0.006, 0.001, and 0.135. While all methods discovered a statistically significant association between microbial changes and the x-axis, the RFtest reported a more significant p-value than the optimal MiRKAT and the aMiSPU, rejecting the null hypotheses at a significance level of 0.01.