 Methodology article
 Open Access
 Published:
Inferring microbial interaction network from microbiome data using RMN algorithm
BMC Systems Biologyvolume 9, Article number: 54 (2015)
Abstract
Background
Microbial interactions are ubiquitous in nature. Recently, many similaritybased approaches have been developed to study the interaction in microbial ecosystems. These approaches can only explain the nondirectional interactions yet a more complete view on how microbes regulate each other remains elusive. In addition, the strength of microbial interactions is difficult to be quantified by only using correlation analysis.
Results
In this study, a rulebased microbial network (RMN) algorithm, which integrates regulatory OTUtriplet model with parametric weighting function, is being developed to construct microbial regulatory networks. The RMN algorithm not only can extrapolate the cooperative and competitive relationships between microbes, but also can infer the direction of such interactions. In addition, RMN algorithm can theoretically characterize the regulatory relationship composed of microbial pairs with low correlation coefficient in microbial networks. Our results suggested that Bifidobacterium, Streptococcus, Clostridium XI, and Bacteroides are essential for causing abundance changes of Veillonella in gut microbiome. Furthermore, we inferred some possible microbial interactions, including the competitive relationship between Veillonella and Bacteroides, and the cooperative relationship between Veillonella and Clostridium XI.
Conclusions
The RMN algorithm provides the reconstruction of gut microbe networks, and can shed light on the dynamical interactions of microbes in the infant intestinal tract.
Background
Microbes are the most abundant and diverse organisms on earth and their interactions are crucial in understanding both the ecology and the evolution of microorganisms. Microbial interactions, including mutualism, competition, parasitism and commensalism, are difficult to quantify as the underlying processes usually cannot be observed directly and are often too complex for laboratory experiments [1]. However, recent advances in highthroughput sequencing technology have made large scale surveys of microbial communities feasible. Metagenomic studies and networkbased approaches have yield detailed information on the composition of microbial communities, which in turn pave the way to study the structure of microbial ecosystems and their dynamics [2–4].
Elucidating competitive and cooperative relationship is a challenge in generating a microbial interaction network because of the direction of such interactions [5]. Competition and cooperation are the two most studied microbial interactions in the recent times with the former dominating the latter in various microbial communities [6–9]. Recent studies have also shown that competitive interactions can drive the evolution of cooperation in microbial ecosystems [9]. Therefore, identifying competitive and cooperative relationships between microbes is profound importance; however, the directional nature of such interactions also poses as a difficult challenge in network construction.
There are several approaches in constructing a microbial network. One commonly used method is the similaritybased network construction where the cooccurrences of two species over multiple timeseries samples are measured to infer their interaction [2, 3, 10–12]. In such networks, nodes correspond to organisms and an edge between two nodes represents the significant relationship of two taxa across a set of time series samples. Although other derived approaches [13–15] can identify the pairwise relationships between microbes by using correlation estimation, they do not identify the nature and the strength of these relationships [1]. In general, these methods cannot capture the complexity of microbial interactions and cannot elucidate how microbes regulate each other.
To model complex relationships (one species influencing multiple others), methods employing mathematical and statistical models have been developed recently. For example, the generalized Lotka–Volterra framework was used to model the dynamics of microbial populations in order to estimate parameters governing speciesspecies interactions [16–18]. Moreover, other studies have used nonparametric regression models to infer the dynamic relationships between three microbial populations [19, 20]. These regression models use smooth functions to describe microbial relationships and estimate their relative effects. Because these regression analyses do not use a priori knowledge to build dynamic microbial interactions, the prediction from those models often lack biological explanations, especially when the microbial community under examination is a complex one [21]. Therefore, how to construct a meaningful microbial interaction network remains a challenge in microbial ecology.
Due to the importance of competitive and cooperative relationships between microbes in complex ecological communities [8, 22, 23], an inference model for constructing microbial regulatory network should consider these interactions simultaneously. However, all of the studies mentioned above do not incorporate these features into their construction of regulatory networks. Association rule mining can be a promising way to infer complex relationships. A rulebased method, named smooth response surface (SRS) algorithm, was developed to construct complex regulatory networks [24]. Basing on a piecewise linearquadratic polynomial, the SRS algorithm uses a threedimensional model surface (3DMS) to successfully infer the relationships between the target, activator, and repressor nodes.
In this paper, we adopt a similar rationale as the SRS algorithm to detect the cooperative and competitive relationships between microbes for constructing microbial regulatory networks. Our new methodology, named the rulebased microbial network (RMN) algorithm, is an advanced version of the rulebased approach. We used Tanh functions for modeling nonlinear responses as regulatory effects of microbes. To our understanding, Tanh functions have been used to model nonlinear responses in some microbial studies, such as dose influences of pathogens [25], regulatory effects of microbes [26], influences of microbial compositions (phospholipids fatty acids) [27]. The RMN algorithm uses a nonlinear regulatory OTUtriplet (NRO) model to capture cooperative and competitive relationships and consequently infers the directional interaction between microbes. We applied our new approach to simulated data and showed that the RMN algorithm was capable to reconstruct microbial regulatory networks. Furthermore, we used Koenig’ data (2011) to construct the interaction network for the microbial community in infant gut [28]; and our results suggest that Bifidobacterium, Streptococcus, Clostridium XI, and Bacteroides play essential roles in regulating the abundance of Veillonella in an infant gut microbiome.
Methods
Sample collection
To construct microbial regulatory networks in infant intestinal tracts, we used the pyrosequencing reads of 16S rRNA genes from 61 infant gut samples, which were treated with four diet states at four time periods [28]. These infant gut samples were collected over 61 time points. The average length of 16S rRNA gene sequences were 237 bps [28], we used only 16S rRNA gene sequences with length longer than 200 bps in this study.
Taxonomic classification and OTU assignment
The Ribosomal Database Project (RDP) classifier [29] was used to classify 16S rRNA pyrosequencing reads and to perform an operational taxonomic unit (OTU) assignment. After accomplishing the OTU assignment with a threshold of 80 % confidence, we carried out a further filtering step by removing operational taxonomic units (OTUs) that appeared in fewer than half of the total time points.
Relative abundance analysis
To estimate the relative difference between individual OTUs at each time point, we calculated the relative abundance of individual OTUs. First, we defined the pyrosequencing read number of the i ^{th} OTU at j ^{th} time point as Y _{ ij }. Second, we defined the sum of the pyrosequencing read numbers at j ^{th} time point as Y _{ j }. Finally, we divided Y _{ ij } by Y _{ j } to obtain the relative proportion X _{ ij } of the i ^{th} OTU at j ^{th} time point_{.} After all pyrosequencing read numbers were transformed into relative proportions, we removed the l ^{th} OTU if both the maximum relative proportion of all X _{ lj } was less than 0.1 and the coefficient of variation of all X _{ lj } was larger than 3.4. In addition, we filtered out the data at the k ^{th} time point if both the maximum of all X _{ ik } was less than 0.1 and the total number of OTUs was less than 60 % at the k ^{th} time point. Finally, the missing values were filled through the Bayesian principal component analysis (BPCA) with K = 5 [30].
Data standardization
The relative proportion of each OTU at each time point was standardized in order to estimate the difference between OTUs at each time point. The relative proportion X _{ ij } of the i ^{th} OTU at the j ^{th} time point was transformed such that it fell within the interval [0,1]. We defined P _{ ij } as (X _{ ij } − X _{min}) /(X _{max} − X _{min}), where P _{ ij } is the standardized relative proportion (SRP) of the i ^{th} OTU at the j ^{th} time point, while X _{min} and X _{max} are the minimum and maximum relative proportions respectively.
Network construction
We developed a rulebased microbial network (RMN) algorithm to construct microbial regulatory networks. The RMN algorithm is based on a nonlinear regulatory OTUtriplet model, named the NRO model for short. Figure 1 shows the flowchart of the RMN algorithm, and the details are as follows:
(i) To investigate how each OTU was involved in microbial regulatory network, we constructed a nonlinear regulatory OTUtriplet model by using triplets of OTUs across all time points (or samples). The nonlinear regulatory OTUtriplet model has three assumptions: First, a cooperation–competition system is assumed to exist in the microbial regulatory network. An OTUtriplet is defined as {O _{ m }, O _{ c }, O _{ t }}, where O _{ m } , O _{ c }, and O _{ t } represent the m ^{th} , c ^{th} and t ^{th} OTUs respectively. Here, O _{ m } and O _{ c } are assumed to be the cooperator and the competitor of O _{ t } respectively. In other words, the relationship between the m ^{th} and the t ^{th} OTUs is assumed to be cooperative, while that between the c ^{th} and the t ^{th} OTUs is assumed to be competitive. Note that O _{ m }, O _{ c } and O _{ t } represent different microbes in genus level. Second, it is assumed that the standardized relative proportion (SRP) of O _{ t } is high when the SRP of O _{ m } is high and the SRP of O _{ c } is low. Conversely, the SRP of O _{ t } is assumed to be low when the SRP of the O _{ m } is low and the SRP of the O _{ c } is high. Third, the relationship between O _{ m }, O _{ c } and O _{ t } can be modeled by using a threedimensional model surface (3DMS) (Fig. 2) which is based on a piecewise nonlinearquadratic polynomial:
In (1), P _{ m } and P _{ c } are the standardized relative proportions (SRPs) of O _{ m } and O _{ c } respectively while P(P _{ m, } P _{ c }) is the inferred SRP of O _{ t }.
(ii) To find suitable OTUtriplets close to the 3DMS for the OTUtriplet regulatory model, we used a lackoffit function to analyze all possible OTUtriplets. The lackoffit function is defined as:
where P ^{^} _{ tn } is the inferred SRP of O _{ t } at the n ^{th} time point (e.g. P(P _{ m, } P _{ c }) at the n ^{th} time point ); P _{ tn } is the observed SRP of O _{ t } at the n ^{th} time point; and P ^{−} _{ tn } represents the average of P _{ tn } across all k time points. A small L(O _{ m }, O _{ c }, O _{ t }) value indicates a good fit between an OTUtriplet {O _{ m }, O _{ c }, O _{ t }} and the 3DMS. In other words, a smaller L(O _{ m }, O _{ c }, O _{ t }) value reveals a stronger relationship among O _{ m }, O _{ c }, and O _{ t }. Similar rationale has been successfully applied for constructing gene regulatory networks [24].
(iii) To check the reliability of OTUtriplets, we used an adjustment function to fix imperfect measurements at one or more time points (or samples). The adjustment function is defined as:
where L _{ (n) }(O _{ m }, O _{ c }, O _{ t }) is the value of lackoffit function when the information from the n ^{th} time point is removed. A smaller D(O _{ m }, O _{ c }, O _{ t }) value indicates the information of the OTUtriplet relationship is reliable. Finally, an integrated function, I(O _{ m }, O _{ c }, O _{ t }) is defined as the value of L(O _{ m }, O _{ c }, O _{ t }) multiplied by the value of (1+ D(O _{ m }, O _{ c }, O _{ t })). A smaller I(O _{ m }, O _{ c }, O _{ t }) value indicates a stronger OTUtriplet relationship among O _{ m }, O _{ c }, and O _{ t. }
(iv) Finally, we defined two criteria, L _{ d }(O _{ m }, O _{ c }, O _{ t }) = L(O _{ m }, O _{ c }, O _{ t })/(1 L(O _{ m }, O _{ c }, O _{ t })) and D(O _{ m }, O _{ c }, O _{ t }), to filter out unreliable OTUtriplets, when 0 ≦ L _{ d }(O _{ m }, O _{ c }, O _{ t }) ≦ L and D(O _{ m }, O _{ c }, O _{ t }) < D. Note the L and D are constants. We then ranked reliable OTUtriplets basing on the sorted I(O _{ m }, O _{ c }, O _{ t }) values in the ascending order, and selected OTUtriplets with lower I(O _{ m }, O _{ c }, O _{ t }) values as candidates for network construction.
Performance analysis
To evaluate the performance of the RMN algorithm, we simulated a microbial regulatory network consisting of ten gut bacteria, three latent factors, and ten noise factors as follows:
We used Tanh functions to create simulated a microbial regulatory network because Tanh functions have been used to model regulatory responses of pathogens and microbes [25–27]. In addition, the similar idea has been used to create simulated gene networks for evaluating the performance of GASA method [31]. In detail, we defined this synthetic regulatory model in Fig. 3, where O _{ i }(s) is the abundance proportion of the i ^{th} OTU with i =1,2,…,10; F _{ j }(s) is the abundance proportion of j ^{th} latent factor with j = 1, 2, 3; and N _{ k }(s) is the abundance proportion of the k ^{th} noise factor with k = 1,2,…,10. Note the value of O _{ i }(s) is between 0 and 1, while F _{ j }(s) is a random variable between 0 and 1. N _{ k }(s) is a random variable representing the noise level in the simulated data (when N _{ k }(s) = 0, no noise exists in the simulated data). Three noise levels are assumed in the simulated data: low noise if the value of N _{ k }(s) is equal to the variance of O _{ i }(s) divided by R _{ k }(s); medium noise if N _{ k }(s) is equal to five times the low noise level; and high noise if N _{ k }(s) is equal to ten times the low noise level. Here, R _{ k }(s) represents the signalnoise ratio and its range was between 4 and 12. Note that R _{ k }(s) is an integral and is generated by using timeseries data of gut bacteria [19]. In detail, R _{ k }(s) is the average proportion of gut bacteria divided by its standard deviation.
We evaluated the performance of the RMN algorithm by calculating the following measures: true positive rate (TPR), true negative rate (TNR), Fmeasure, and Accuracy. TPR is defined as TP/(TP + FN), where TP and FN are the numbers of true positive interactions (links) and false negative interactions (links) respectively. TNR is defined as TN/(TN + FP), where TN and FP are the numbers of true negative interactions (links) and false positive interactions (links) respectively. We define Fmeasure as 2*TPR*TNR/ (TPR + TNR), and accuracy as (TP + TN)/(TP + FN + TN + FP).
Results and discussion
Performance of RMN algorithm using simulation data
We generated simulated data under different levels of noise and evaluated the performance of our RMN algorithm by using the above mentioned four evaluation measures. Then, we analyzed the occurrence rate (%) of 50000 simulated microbial regulatory networks when any three of four evaluation measures showed the efficiency superior or equal to the defined value. (see Fig. 4). All four evaluation measures suggest high performance from our RMN algorithm. Comparing different noise levels, RMN algorithm has shown the best efficiency when there is no noise in the simulated data. Even at medium noise level, over 82 % of the results show efficiency values higher than 0.64. In addition, the assessment from those four evaluation measures shows 0.64 efficiency value on more than 74 % simulated data. There are no significant differences between the performance results when the noise is varied from zero to medium level (see Fig. 5(a)). Although the performance of our RMN algorithm at the high noise level is worse than those from lower noise levels, its efficiency value is still over 0.68. More interestingly, according to TNR, there were fewer false positive interactions (links) when inferred from the high noise data than that from lower noise data (see Fig. 5(a)). In addition, best efficiency under different noise interferences exceeds 0.82 (see Fig. 5(b)). During developing RMN algorithm, we used two conditions, L _{d}(O _{m}, O _{c}, O _{t}) and D(O _{m}, O _{c}, O _{t}), to filter out unreliable OTUtriplets, and we tried to select the suitable thresholds for L _{d}(O _{m}, O _{c}, O _{t}) and D(O _{m}, O _{c}, O _{t}) by analyzing the simulation data (Table 1). The best efficiency for all evaluations reached over 0.82 when the threshold at L _{ d }(O _{ m }, O _{ c }, O _{ t }) and D(O _{ m }, O _{ c }, O _{ t }) were used (see Table 1). Moreover, the efficiency for all evaluations reached over 0.71 when the threshold at L _{ d }(O _{ m }, O _{ c }, O _{ t }) and D(O _{ m }, O _{ c }, O _{ t }). All in all, our efficiency evaluations demonstrate that the performance of our RMN algorithm is robust under three latent factors and different noise levels.
Inferring the regulatory relationship of microbial pairs with low correlation coefficient
Recent study shows that correlations between the abundances of microbes do not imply microbial interactions [18]. In order to investigate whether OTUtriplets obtained by RMN algorithm are significantly correlated based on their correlation coefficient and whether RMN algorithm can find the regulator relationship behind standard correlationbased approaches, we analyzed the optimal resulting triplets with true positive links (filtered by RMN algorithm). Here we used Spearman's rank correlation coefficient for the estimation of nonlinear correlation coefficient. The Pr value represents a probability of RMN algorithm prediction for OTU pairs whose nonlinear correlation coefficient is less than 0.5. In detail, Pr value (%) was defined as (O _{ m } O _{ t }_count + O _{ c } O _{ t } _count) *100 %/Total_count, where O _{ m } O _{ t }_count is the number of resulting triplets with a true positive link between cooperator OTU and target OTU pair whose nonlinear correlation coefficient is less than 0.5; O _{ c } O _{ t } _count is the number of resulting triplets with a true positive link between competitor OTU and target OTU pair whose nonlinear correlation coefficient is greater than −0.5; Total count is the number of resulting triplets with a true positive link. As Table 2 shown, under different levels of noise inference, Pr values for most analyzed results were more than 15 %. That may imply the resulting triplets contain OTU pairs with low nonlinear correlation coefficient. For example, more than 15 % OTU pairs with low nonlinear correlation coefficient (less than 0.5) could be found when the threshold was set at 0≦L _{ d }(O _{ m }, O _{ c }, O _{ t }) ≦0.8 and D(O _{ m }, O _{ c }, O _{ t }) < 0.8. Based on the threshold, the performances of RMN algorithm were more than 0.82 (see Table 1). In addition, more than 20 % OTU pairs with low nonlinear correlation coefficient (less than 0.5) could be found when the threshold was set at 0≦L _{ d }(O _{ m }, O _{ c }, O _{ t }) ≦3.8 and D(O _{ m }, O _{ c }, O _{ t }) < 0.8. Based on the threshold, the performances of RMN algorithm were more than 0.71 (see Table 1). In other words, the RMN algorithm may characterize the regulatory relationship composed of OTU pairs with low correlation coefficient, and has a choice to find the regulator relationship behind standard correlationbased approaches.
The construction of microbial regulatory network from infant gut data
To further demonstrate the capability of our RMN algorithm, we constructed microbial interaction networks from a set of real microbial community time series data. The data used consist of 61 infant gut samples treated with four diet states at four time periods [28].
After performing taxonomic classification and OTU assignment, 101 OTUs were obtained (Additional file 1). However, after the filtering process removing OTUs that appeared in less than half of the total time points, only 14 OTUs met the requirement (Additional file 2). In addition, after performing relative abundance analysis and filling missing values, 10 OTUs with 43 time points were obtained. Finally, according to Fig. 3c from Koenig, et al. [28], only 10 OTUs with 36 time points between four steps were selected for analyzing using our RMN algorithm (Additional file 3). Consistently, among 10 OTUs, 4 OTUs (OTU15, OTU16, OTU34, and OTU94) were found to be involved in competitive and cooperative interactions in Freilich’s study [32]. Since microbiome data are usually enormous and inherently noisy, the threshold was set at 0≦L _{ d(} O _{ m }, O _{ c }, O _{ t }) ≦3.8 and D(O _{ m }, O _{ c }, O _{ t }) < 0.8 for better prediction sensitivity (see Table 3). Before filtering by this threshold, top 100 possible OTUtriplets from 10 OTUs analysis were listed in Additional file 4. After filtering by this threshold, only 3 reliable OTUtriplets from 5 OTUs (Additional file 5) were used to construct microbial regulatory network. The resulting microbial interaction network among five OTUs is displayed in Fig. 6. Gramnegative anaerobic cocci, Veillonella (OTU100), is found to interact cooperatively with three OTUs, namely Clostridium XI (OTU30), Bifidobacterium (OTU16) and Streptococcus (OTU94), but competitively with Bacteroides (OTU15). Among them, Bifidobacterium, Streptococcus, and Bacteroides have been identified in competitive and cooperative metabolic interactions [32]. In particular, we observed a positive correlation between Veillonella and Bifidobacterium , both residing in human intestines and oral mucosa as anaerobic commensal organisms [33]. This directional interaction also suggests the commensal relationship between Veillonella and Bifidobacterium. In addition, Veillonella was also observed to have a positive correlation with Streptococcus. They have been found to interact cooperatively during the production and degradation of Lactic acid [34–40]. Recently, positive correlations were also found between the abundances of Veillonella and Streptococcus [41]. In addition, the constructed network also suggests that competitive relationship may exist between Veillonella and Bacteroides. Veillonella and Bacteroides are gutassociated obligate anaerobic genera found in maternal faeces, breast milk and neonatal faeces [42]; and the difference in their efficiency on oligosaccharide consumptions may suggest they occupy different metabolic niches [43]. In addition, our network suggests there might be a cooperative relationship between Veillonella and Clostridium XI Clostridium glycolicum; and this cooperative relationship might suggest their simultaneous involvement in the ear infections in infants. In fact, those two gut bacteria have been found to be associated with bacterial infections, such as bacteremia [44, 45]. Finally, microbial interaction network from Fig. 6 were analyzed to check if the RMN algorithm can predict microbial interactions from low correlation data. As shown, it confirmed that OTU pairs with low correlations can be predicted by RMN algorithm (Additional file 6). In general, our constructed network can extrapolate cooperative and competitive relationships between microbes in infant guts.
Conclusions
In this paper we presented RMN algorithm, a rulebased algorithm using the OTUtriplet model with parametric weighting function, to construct microbial regulatory networks. RMN algorithm can theoretically extract directional interactions to delineate the cooperative and competitive relationships between microbes from highthroughput sequencing time series data. The efficiency estimation shows the high performance of RMN algorithm on simulated data with three latent factors and different noise levels. We also applied our framework to identify new relationships in infant gut samples. RMN algorithm is both computationally feasible and capable of detecting biologically significant progresses embedded in a microbial community.
Although simple pairwise relationship has been used to describe complex forms of microbial interactions, advanced techniques are still required to infer complex microbial networks. Recent studies show that the cooperative and competitive relationships between microbes coexist in complicated ecological communities [8, 22, 23]. In addition, triplet relation, like activatorrepressortarget pattern, has been adopted to construct gene regulatory network form time course gene expression data [24]. Basing on such a rationale, we developed RMN algorithm to describe the relationship among every triplet of microbe as the cooperatorcompetitortarget. RMN algorithm is used to search for triplet relations of OTUs where the abundances of the target should be high while that for the cooperator is high and that for the competitor is low (conversely, the abundances of the target should be low while that for the cooperator is low and that for the competitor is high). Our evaluation results suggest RMN algorithm has high performance when subjected to simulated data with different levels of noise (see Figs. 4 and 5). As presented in the results, our algorithm has characterized wellknown and possible microbial interactions from time series microbiome data of infant gut samples (see Fig. 6).
In summary, networks reconstructed by using simple similaritybased method often fail to capture the complex topology of a real system, and such networks may be full of false positive links. Furthermore, biologically significant relationships are often difficult to be identified from such networks. Our RMN algorithm can produce low false positive links even for highly noisy data (see Fig. 5(a)). Moreover, RMN algorithm is designed to explore complex networks consisted of nonlinear, decentralized and directed interactions among microbes. However, having said all these, the drawback of RMN algorithm is that it may not be able to detect simple linear correlation between two OTUs. Another potential concern is that we tested our method using a simulated data with the Tanh functions, the circularity in the model construction and validation can cause bias on our results. Overall, RMN algorithm can theoretically predict microbial interactions from low correlation data (see Tables 1 and 2). Our method should be a promising starting step to identify novel microbial interactions that normally cannot be found by using similaritybased methods. Therefore, by integrating similaritybased methods and our RMN algorithm, one can potentially gain a more accurate picture of microbial interactions towards a better understanding of microbial dynamics.
Abbreviations
 RMN:

Rulebased microbial network
 SRS:

Smooth response surface
 3DMS:

Threedimensional model surface
 NRO:

Nonlinear regulatory OTUtriplet
 RDP:

Ribosomal Database Project
 OTU:

Operational taxonomic unit
 OTUs:

Operational taxonomic units
 BPCA:

Bayesian principal component analysis
 SRP:

Standardized relative proportion
 SRPs:

Standardized relative proportions
 TPR:

True positive rate
 TNR:

True negative rate
References
 1.
Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol. 2012;10(8):538–50. doi:10.1038/nrmicro2832.
 2.
Zhou J, Deng Y, Luo F, He Z, Yang Y. Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2. mBio. 2011;2(4). doi:10.1128/mBio.0012211.
 3.
Deng Y, Jiang YH, Yang Y, He Z, Luo F, Zhou J. Molecular ecological network analyses. BMC Bioinf. 2012;13:113. doi:10.1186/1471210513113.
 4.
Sun Y, Cai Y, Mai V, Farmerie W, Yu F, Li J, et al. Advanced computational algorithms for microbial community analysis using massive 16S rRNA sequence data. Nucleic Acids Res. 2010;38(22), e205. doi:10.1093/nar/gkq872.
 5.
Freilich S, Kreimer A, Meilijson I, Gophna U, Sharan R, Ruppin E. The largescale organization of the bacterial network of ecological cooccurrence interactions. Nucleic Acids Res. 2010;38(12):3857–68. doi:10.1093/nar/gkq118.
 6.
Lo Giudice A, Brilli M, Bruni V, De Domenico M, Fani R, Michaud L. Bacteriumbacterium inhibitory interactions among psychrotrophic bacteria isolated from Antarctic seawater (Terra Nova Bay, Ross Sea). FEMS Microbiol Ecol. 2007;60(3):383–96. doi:10.1111/j.15746941.2007.00300.x.
 7.
PerezGutierrez RA, LopezRamirez V, Islas A, Alcaraz LD, HernandezGonzalez I, Olivera BC, et al. Antagonism influences assembly of a Bacillus guild in a local community and is depicted as a foodchain network. ISME J. 2013;7(3):487–97. doi:10.1038/ismej.2012.119.
 8.
Foster KR, Bell T. Competition, not cooperation, dominates interactions among culturable microbial species. Curr Biol. 2012;22(19):1845–50. doi:10.1016/j.cub.2012.08.005.
 9.
Celiker H, Gore J. Competition between species can stabilize publicgoods cooperation within a species. Mol Syst Biol. 2012;8:621. doi:10.1038/msb.2012.54.
 10.
DuranPinedo AE, Paster B, Teles R, FriasLopez J. Correlation network analysis applied to complex biofilm communities. PLoS One. 2011;6(12), e28438. doi:10.1371/journal.pone.0028438.
 11.
Zhou J, Deng Y, Luo F, He Z, Tu Q, Zhi X. Functional molecular ecological networks. mBio. 2010;1(4). doi:10.1128/mBio.0016910.
 12.
Xia LC, Steele JA, Cram JA, Cardon ZG, Simmons SL, Vallino JJ et al. Extended local similarity analysis (eLSA) of microbial community and other time series data with replicates. BMC systems biology. 2011;5 Suppl 2:S15. doi:10.1186/175205095S2S15.
 13.
Fuentes S, van Nood E, Tims S, Heikampde Jong I, ter Braak CJ, Keller JJ, et al. Reset of a critically disturbed microbial ecosystem: faecal transplant in recurrent Clostridium difficile infection. ISME J. 2014;8(8):1621–33. doi:10.1038/ismej.2014.13.
 14.
Xu Z, Hansen MA, Hansen LH, Jacquiod S, Sorensen SJ. Bioinformatic approaches reveal metagenomic characterization of soil microbial community. PLoS One. 2014;9(4), e93445. doi:10.1371/journal.pone.0093445.
 15.
Zhang Z, Geng J, Tang X, Fan H, Xu J, Wen X, et al. Spatial heterogeneity and cooccurrence patterns of human mucosalassociated intestinal microbiota. ISME J. 2014;8(4):881–93. doi:10.1038/ismej.2013.185.
 16.
Stein RR, Bucci V, Toussaint NC, Buffie CG, Ratsch G, Pamer EG, et al. Ecological modeling from timeseries inference: insight into dynamics and stability of intestinal microbiota. PLoS Comput Biol. 2013;9(12), e1003388. doi:10.1371/journal.pcbi.1003388.
 17.
Marino S, Baxter NT, Huffnagle GB, Petrosino JF, Schloss PD. Mathematical modeling of primary succession of murine intestinal microbiota. Proc Natl Acad Sci U S A. 2014;111(1):439–44. doi:10.1073/pnas.1311322111.
 18.
Fisher CK, Mehta P. Identifying keystone species in the human gut microbiome from metagenomic timeseries using sparse linear regression. PLoS One. 2014;9(7), e102451. doi:10.1371/journal.pone.0102451.
 19.
Trosvik P, Rudi K, Naes T, Kohler A, Chan KS, Jakobsen KS, et al. Characterizing mixed microbial population dynamics using timeseries analysis. ISME J. 2008;2(7):707–15. doi:10.1038/ismej.2008.36.
 20.
Trosvik P, Stenseth NC, Rudi K. Convergent temporal dynamics of the human infant gut microbiota. ISME J. 2010;4(2):151–8. doi:10.1038/ismej.2009.96.
 21.
Trosvik P, Skanseng B, Jakobsen KS, Stenseth NC, Naes T, Rudi K. Multivariate analysis of complex DNA sequence electropherograms for highthroughput quantitative analysis of mixed microbial populations. Appl Environ Microbiol. 2007;73(15):4975–83. doi:10.1128/AEM.0012807.
 22.
Korb J, Foster KR. Ecological competition favours cooperation in termite societies. Ecol Lett. 2010;13(6):754–60. doi:10.1111/j.14610248.2010.01471.x.
 23.
Mitri S, Xavier JB, Foster KR. Social evolution in multispecies biofilms. Proc Natl Acad Sci U S A. 2011;108 Suppl 2:10839–46. doi:10.1073/pnas.1100292108.
 24.
Xu H, Wu P, Wu CF, Tidwell C, Wang Y. A smooth response surface algorithm for constructing a gene regulatory network. Physiol Genomics. 2002;11(1):11–20. doi:10.1152/physiolgenomics.00060.2001.
 25.
Yang SX. A neural network model for dose–response of foodborne pathogens. Appl Soft Comput. 2003;3(2):85–96.
 26.
Mozhayskiy V, Tagkopoulos I. Horizontal gene transfer dynamics and distribution of fitness effects during microbial in silico evolution. BMC Bioinf. 2012;13 Suppl 10:S13.
 27.
Schryver JC, Brandt CC, Pfiffner SM, Palumbo AV, Peacock AD, White DC, et al. Application of nonlinear analysis methods for identifying relationships between microbial community structure and groundwater geochemistry. Microb Ecol. 2006;51(2):177–88.
 28.
Koenig JE, Spor A, Scalfone N, Fricker AD, Stombaugh J, Knight R, et al. Succession of microbial consortia in the developing infant gut microbiome. Proc Natl Acad Sci U S A. 2011;108 Suppl 1:4578–85. doi:10.1073/pnas.1000081107.
 29.
Wang Q, Garrity GM, Tiedje JM, Cole JR. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol. 2007;73(16):5261–7. doi:10.1128/AEM.0006207.
 30.
Oba S, Sato MA, Takemasa I, Monden M, Matsubara K, Ishii S. A Bayesian missing value estimation method for gene expression profile data. Bioinformatics. 2003;19(16):2088–96.
 31.
Chen CM, Lee C, Chuang CL, Wang CC, Shieh GS. Inferring genetic interactions via a nonlinear model and an optimization algorithm. BMC Syst Biol. 2010;4:16.
 32.
Freilich S, Zarecki R, Eilam O, Segal ES, Henry CS, Kupiec M, et al. Competitive and cooperative metabolic interactions in bacterial communities. Nat Communicants. 2011;2:589. doi:10.1038/ncomms1597.
 33.
Verma R, Dhamija R, Ross SC, Batts DH, Loehrke ME. Symbiotic bacteria induced necrotizing pancreatitis. JOP Journal of the pancreas. 2010;11(5):474–6.
 34.
Chalmers NI, Palmer Jr RJ, Cisar JO, Kolenbrander PE. Characterization of a Streptococcus sp.Veillonella sp. community micromanipulated from dental plaque. J Bacteriol. 2008;190(24):8145–54. doi:10.1128/JB.0098308.
 35.
Palmer RJ, Diaz PI, Kolenbrander PE. Rapid succession within the Veillonella population of a developing human oral biofilm in situ. J Bacteriol. 2006;188(11):4117–24.
 36.
Kara D, Luppens SB, Cate JM. Differences between single and dualspecies biofilms of Streptococcus mutans and Veillonella parvula in growth, acidogenicity and susceptibility to chlorhexidine. Eur J Oral Sci. 2006;114(1):58–63. doi:10.1111/j.16000722.2006.00262.x.
 37.
Mikx FHM, van der Hoeven JS. Symbiosis of Streptococcus mutans and Veillonella alcalescens in mixed continuous cultures. Archives of Oral Biology. 1975;20(7):407–10. doi:http://dx.doi.org/10.1016/00039969(75)902241.
 38.
Mikx FH, van der Hoeven JS, Konig KG, Plasschaert AJ, Guggenheim B. Establishment of defined microbial ecosystems in germfree rats. I. The effect of the interactions of streptococcus mutans or Streptococcus sanguis with Veillonella alcalescens on plaque formation and caries activity. Caries Res. 1972;6(3):211–23.
 39.
McBride BC, Van der Hoeven JS. Role of interbacterial adherence in colonization of the oral cavities of gnotobiotic rats infected with Streptococcus mutans and Veillonella alcalescens. Infect Immun. 1981;33(2):467–72.
 40.
van der Hoeven JS, Toorop AI, Mikx RH. Symbiotic relationship of Veillonella alcalescens and Streptococcus mutans in dental plaque in gnotobiotic rats. Caries Res. 1978;12(3):142–7.
 41.
Rigsbee L, Agans R, Shankar V, Kenche H, Khamis HJ, Michail S, et al. Quantitative profiling of gut microbiota of children with diarrheapredominant irritable bowel syndrome. Am J Gastroenterol. 2012;107(11):1740–51. doi:10.1038/ajg.2012.287.
 42.
Jost T, Lacroix C, Braegger CP, Rochat F, Chassard C. Vertical motherneonate transfer of maternal gut bacteria via breastfeeding. Environ Microbiol. 2014;16(9):2891–904. doi:10.1111/14622920.12238.
 43.
Marcobal A, Barboza M, Froehlich JW, Block DE, German JB, Lebrilla CB, et al. Consumption of human milk oligosaccharides by gutrelated microbes. J Agric Food Chem. 2010;58(9):5334–40. doi:10.1021/jf9044205.
 44.
Fisher RG, Denison MR. Veillonella parvula bacteremia without an underlying source. J Clin Microbiol. 1996;34(12):3235–6.
 45.
Elsayed S, Zhang K. Clostridium glycolicum bacteremia in a bone marrow transplant patient. J Clin Microbiol. 2007;45(5):1652–4. doi:10.1128/JCM.0184506.
Acknowledgements
We thank for the comments from the anonymous reviewers. The experimental data provided by Jeremy E. Koenig et al. is also appreciated. This work was supported by the National Science Council of Taiwan (Grant No: MOST 1032313B001004).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
KNT and DW designed the analyses. KNT and SHL collected the data. KNT and WCL performed the analyses. KNT and DW wrote the paper. KNT and DW were the principal investigators and conceived the analyses. All authors read and approved the final manuscript.
Additional files
Additional file 1:
The table of OTUs for all the samples at the genus level with OTU assignment before filtering process. (XLS 81 kb)
Additional file 2:
The table of OTUs for all the samples at the genus level with OTU assignment after filtering process. (XLS 34 kb)
Additional file 3:
The table of OTUs for all the samples at the genus level with data standardization. (XLS 26 kb)
Additional file 4:
List of top 100 possible OTUtriplets ranking by I ( O _{ m } , O _{ c } , O _{ t } ) values in the ascending order. These OTUtriplets were selected without using the threshold set at 0≦L _{d}(O _{m}, O _{c}, O _{t}) ≦3.8 and D(O _{m}, O _{c}, O _{t}) < 0.8. (XLS 28 kb)
Additional file 5:
List of the reliable OTUtriplets ranking by I ( O _{ m } , O _{ c } , O _{ t } ) values in the ascending order. These OTUtriplets were selected when the threshold was set at 0≦L _{d}(O _{m}, O _{c}, O _{t}) ≦3.8 and D(O _{m}, O _{c}, O _{t}) < 0.8. (XLS 19 kb)
Additional file 6:
The nonlinear correlation coefficient of OTU pairs predicted by RMN algorithm. (XLS 20 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Microbiome
 Nextgeneration sequencing
 Microbial regulatory network
 Cooperative and competitive relationship
 OTUtriplet model
 RMN algorithm