Introducing METATEA. METATEA aims to be a toolbox for easy analyses of the microbial community using amplicon sequence data (and whole genome sequence data in the future) at the sequence level and for the detection of antagonistic relationships. Since intraspecies variation in even one gene can change the virulence of a pathogen12, we developed pipelines for analysis at the sequence level (Fig. 2). Analyzing 428 IBD samples as a case, we had challenges from sequencing errors. Sequencing errors may increase the number of unique sequences, which are weakly proportional to the total numbers of reads (Supplementary Table S1); we assumed such erroneous reads to generally occur once or twice in a sequence file. Therefore, we developed filtering step and normalization step for the removal of rare sequences and the adjustment. Also, various quality control steps were developed to accurately handle reads with Ns and short/long reads at the sequence level. For example, most ambiguous bases (Ns) from the pyrosequencing platform can be corrected13 using METATEA. For the detection of microbial antagonism, METATEA implements the calculation of ECV (see Methods). Being more user-friendly, METATEA supports matrix output format. The format was designed to easily perform further analyses using other statistical tools like R. METATEA is written in Perl, and we support both Windows-compatible and Linux-compatible executables. METATEA Linux version was compiled from C code for relatively fast analyses of large data sets. METATEA is free, and its detailed manual is available from the project website (https://sourceforge.net/projects/metatea/).
Comparison between taxon-based and sequence-based analyses. We determined the proportional threshold (≥ 0.1%) and extracted 3,202 unique sequences and00 calculated their proportions. To validate our method, the microbial community structure was compared to that in the previous study11, from which our data set was derived. Although the previous study had used different methods and samples relative to our current one, the regression line in Supplementary Fig. S1 showed our present results to be in accordance with the previous one’s with respect to community composition (r2 = 0.989 in HC and 0.973 in CD).
Generally, analysis results at the sequence level agreed with, and even outperformed, the results at family/genus/OTU levels. First, increase/decrease of bacteria, according to the disease subtypes, at the sequence level generally matched the results of a previous study11 at the genus or family level (Supplementary Fig. S2). For example, five bacterial sequences in the Pasteurellaceae family increased in CD, as reported in a previous study11. However, Bacteroidales is comprised of various bacteria, whose increase was offset by its decrease. Second, the risk of CD was predicted using random forests at both sequence and OTU (1% and 3% dissimilarity) levels. Prediction at the sequence level outperformed that at the OTU level in terms of the area under the receiver operating characteristic (ROC) curve (AUC) shown in Fig. 3. Third, the co-occurrence of CD-associated organisms in the previous study generally coincided with the correlation network in this study (Supplementary Fig. S3), despite the use of different samples and methods. For example, both sequence identifiers S0276 (Fusobacterium periodonticum) and S0345 (Veillonella dispar) were increased in CD and showed possible agreement (Supplementary Table S2). However, we could also discover a few exceptions, like S0066 (Clostridium symbiosum; Lachnospiraceae) and S0346, which were increased in CD and showed possible agreement, although Lachnospiraceae was co-excluded in CD in the previous study11.
Analyses of the microbial community composition. The most frequent sequence identifiers in our IBD data set were S0078 (Bacteroides vulgatus), S0052 (Bacteroides fragilis), S0029 (E. coli), S0051 (Faecalibacterium prausnitzii-like), and S0152 (Bacteroides dorei) in descending order (Supplementary Table S2). The microbial sequences were the main coefficients in our PCA analysis (Supplementary Fig. S4a). Generally, HC samples clustered in the center, slightly toward sequence S0078, in the PCA analysis. PCA results coincided with the heat map results (Supplementary Fig. S4b). Single bacterial sequences did not show high proportions in HC samples (≥ 50%), whereas certain CD and UC samples did show high proportions (≥ 60%) of sequence S0029 or S0046. Bacterial community structures in CD (pMANOVA: P < 0.001) and UC (pMANOVA: P < 0.001) were significantly different from that in HC. Several intraspecies variations in disease subtypes were identified, and their average proportions changed differently according to disease subtypes (Supplementary Table S2). Particularly, intraspecies variation of B. vulgatus, E. coli, F. prausnitzii, Parabacteroides distasonis, and Kineothrix alysoides in average proportion, according to disease subtypes, were prevalent, whereas variation of B. dorei, Haemophilus parainfluenzae, Prevotella copri, and Fusobacterium nucleatum were less obvious (Supplementary Fig. S5).
Detection of disease-associated community members. Multiple Mann-Whitney U tests were performed using METATEA to identify sequences associated with disease subtypes (Supplementary Table S3). Generally, H. parainfluenzae, Dialister pneumosintes, F. periodonticum, and many other sequences were significantly increased in CD. H. parainfluenzae-like, F. prausnitzii-like, Sutterella massiliensis-like, and many others were significantly increased in UC. We newly discovered Veillonella atypica, Pseudomonas lini, and Hyphomicrobium zavarzinii to possibly be related to CD, IC, and UC, respectively, probably due to high-resolution analyses at the sequence level. Lachnospiraceae showed variations in disease subtypes (Supplementary Fig. S6). Intraspecies variations in Z scores were apparent, particularly among sequences of low proportions (Supplementary Table S2). For example, although F. prausnitzii (S0376) and some F. prausnitzii-like sequences (S0051) were increased in HC, other F. prausnitzii-like sequences, such as S0397 and S2816, were increased in CD and UC, respectively (Supplementary Table S3). Generally, sequences with high Z scores (≥ 2.3264) in CD and UC compared to those in HC might be those of bacteria associated with various sites of the body or opportunistic environmental bacteria, rather than members of the normal microbiota in the gut (Supplementary Table S4). Some of these aberrant bacteria of putatively external origins are reportedly related to various diseases like endocarditis14, carcinoma15, and abscess16. For more accurate prediction of CD risk, using more biomarkers at the identical sequence level may be an appropriate approach rather than using a small number of powerful biomarkers. For example, sequences like S0362 (Roseburia inulinivorans) can be the most powerful biomarkers to predict CD risk (Fig. 3b). However, even the highest proportion of S0362 cannot confirm HC status (Supplementary Fig. S7).
Comparison between ECV, PCC, and SCC. ECV was designed to selectively identify highly exclusive relationships among all microbial pairs (Fig. 4). This ability is very important in terms of microbial data because most microbial pairs can show widely scattered inversely proportional correlations as shown in Fig. 5 around the median values of ECV, PCC, and SCC. ECV successfully identified strong antagonistic (almost non-coexisting) relationships between K. pneumoniae (S0390) and F. prausnitzii-like (S0764, S1650, and S25900) sequences around the maximum values (Fig. 5a). Although the negative relationships between K. pneumoniae and F. prausnitzii-like sequences are to be experimentally verified, K. pneumoniae increased in disease samples, whereas F. prausnitzii increased in control samples in previous studies17,18.
We then compared ECV, PCC, and SCC for their capabilities to identify antagonism. Widely used PCC is one of the robust means to detect positive and negative linear correlations. Although PCC finds well co-occurring microbial members having positive correlations at the maximum values (~ 1), strong negative relationships (~-1) were undetectable using PCC (Fig. 5b). Among the microbial pairs (B. vulgatus and B. dorei as well as closely related sequences) with extremely low correlation coefficients, none had the coefficients lower than − 0.3 (Supplementary Table S5). SCC assesses monotonic relationships using a rank-based measure (whether linear or not), while PCC assesses linear correlations. SCC identified microbial pairs with inversely proportional correlations better around the minimum values than PCC (Fig. 5c). However, strong antagonistic relationships could not be identified even with SCC. For example, S0535 (Haemophilus pittmaniae) and S1732 (E. coli-like) showed strong antagonistic relationships, and their SCC was 0.039. Pairs (Lachnospiraceae and Oscillibacter ruminantium-like; Blautia wexlerae and Escherichia coli-like; Bacteroides uniformis and Escherichia coli-like) with correlation coefficients from the minimum value were − 0.367, -0.346, and − 0.336.
Relationships between community members. Relationships across different sequences of the same species may be associated with disease subtypes and be useful for predicting CD risk. For example, IBD samples were categorized into two groups, as shown in Supplementary Fig. 8a, and high ratio of S0103 (B. vulgatus-like) to S0092 (B. vulgatus-like) might be related to CD. We successfully detected mutually exclusive relationships using ECV at three different OTU levels (0, 1, and 3% dissimilarity); S0390 (K. pneumoniae), 1%OTU281, and 3%OTU010 generally had significantly exclusive relationships with a variety of bacteria (Supplementary Table S6-S8; Supplementary Fig. S9) and increased in CD than in HC (Fig. 6). Particularly, high proportions of S0390 might be related to reduced bacterial diversity (Supplementary Fig. S10). Different microbes showed different ECV values with certain disease-associated microbes. We regarded 15 sequences (average proportion ≥ 0.1%) as putative probiotics that occurred at least three times in HC than in CD, and 8 sequences were removed by the ECV formula (see Methods), The 7 putative probiotic sequences, S0342 (F. prausnitzii), S0362 (R. inulinivorans), S0463 (Blautia luti), S0464 (Clostridium amygdalinum), S0661 (Lachnospiracea_incertae_sedis), S0706 (Lachnospiracea_incertae_sedis), and S0707 (unclassified_Lachnospiraceae) in general showed exclusive relationships with disease-associated sequences (Supplementary Table S9). Particularly S0342 showed strong mutually exclusive relationships with E. coli whereas S0464 did with F. periodonticum (Supplementary Table S10). Additionally, ECV was designed to be less affected by sample sizes, and pairs of S0390 and other sequences showed higher ECV values in HC than in CD (Supplementary Table S11). For example, ECV value of S0089 (Bacteroides koreensis-like) and S0390 in HC was higher than in CD, suggesting putative synergistic effects of S0390 whether it is direct or not (Supplementary Fig. S8b). Additional studies using more clinical data would validate our observation.