https://orcid.org/0000-0001-7501-0519
ABSTRACT
A large amount of sequencing data is produced in microbial community ecology studies using the high-throughput sequencing technique, especially amplicon-sequencing-based community data. After conducting the initial bioinformatic analysis of amplicon sequencing data, performing the subsequent statistics and data mining based on the operational taxonomic unit and taxonomic assignment tables is still complicated and time-consuming. To address this problem, we present an integrated R package-‘microeco’ as an analysis pipeline for treating microbial community and environmental data. This package was developed based on the R6 class system and combines a series of commonly used and advanced approaches in microbial community ecology research. The package includes classes for data preprocessing, taxa abundance plotting, venn diagram, alpha diversity analysis, beta diversity analysis, differential abundance test and indicator taxon analysis, environmental data analysis, null model analysis, network analysis and functional analysis. Each class is designed to provide a set of approaches that can be easily accessible to users. Compared with other R packages in the microbial ecology field, the microeco package is fast, flexible and modularized to use and provides powerful and convenient tools for researchers. The microeco package can be installed from CRAN (The Comprehensive R Archive Network) or github (https://github.com/ChiLiubio/microeco).
INTRODUCTION
Applications of the high-throughput sequencing technique in microbial ecology produce large amounts of sequencing data. Nowadays, sequencing data analysis in microbial community ecology can be classified arbitrarily into the bioinformatic analysis and the following statistical analysis. Bioinformatic analysis is a typically computation-intensive work, and the unix/linux server is suitable for this operation. After obtaining the necessary operational taxonomic unit (OTU) table, taxonomic assignments and phylogenetic tree based on the amplicon sequencing data, the downstream data analysis tends to be diverse, complicated and time-consuming. Although there are some statistical and plotting approaches in the bioinformatic platform, such as QIIME (Caporaso et al. 2010), mothur (Schloss et al. 2009), RDP (http://rdp.cme.msu.edu/) and SiLVAngs (https://ngs.arb-silva.de/silvangs/), more powerful, flexible and comprehensive tools are necessary to perform the data mining effectively.
R programming language and its package systems for statistical computing stand out for their power and flexibility in scientific applications. Besides the classical ecological packages, e.g. vegan (Oksanen et al. 2019), several R packages have been designed to perform the complicated analysis of sequencing-based microbial community data, such as phyloseq (Mcmurdie and Holmes 2013), microbiome (https://github.com/microbiome/microbiome), microbiomeSeq (http://www.github.com/umerijaz/microbiomeSeq), ampvis2 (https://madsalbertsen.github.io/ampvis2/reference/index.html), MicrobiomeR (https://github.com/vallenderlab/MicrobiomeR) and Rhea (Lagkouvardos et al. 2017). In terms of the software design, the phyloseq package is well developed based on the S4 class. The microbiome, microbiomeSeq, MicrobiomeR and ampvis2 packages all depend on the phyloseq class of phyloseq package and cover a range of functions in microbial community ecology research (Table S1, Supporting Information). However, they all lack the encapsulations relevant to analysis approaches and lack some significant and cutting-edge approaches (Table S1, Supporting Information). Thus, those package frameworks are insufficient for the current data mining pipeline of microbial community research. In the context of data analysis with multiple statistical and plotting approaches, a modularized software is urgently needed for users to reduce time cost on the learning and usage. The Rhea package has a brief and clear framework for the workflow, but it lacks many important methods and is not a traditional function-based package.
Here we created a package named microeco based on the R6 class system (Chang 2019), which provides the encapsulated object-oriented programming paradigm. Classes were created for each part of analysis approaches to make the package framework modularized, clear and brief (Fig. 1). Each class encapsulates a series of functions and diverse algorithms. We integrated some frequently used approaches to make the package supports a wide range of microbial community analysis, such as LEfSe (Segata et al. 2011), redundancy analysis (RDA), co-occurrence network analysis, functional prediction and null model analysis (Stegen et al. 2013). We also developed some unique and advanced functions, such as the taxa composition analysis in venn diagram and functional redundancy calculation for both prokaryotes and fungi.
Flowchart of the processing performed in microeco. The workflow starts with the OTU table (or other traditional community-species abundance table), taxonomic table, sample information table and phylogenetic tree (necessary for phylogenetic analysis). The designed classes include the basic class, i.e. microtable class, which is used for data preprocessing, and the extending classes, and each class represents one type of analysis. All the extending classes were named with the prefix ‘trans_’, representing the transformation of various data. The forward arrows represent the analysis orders. The reverse arrows mean that the operation may generate a new class or change the object in the corresponding class.
METHODS
The core processing workflow of the package is presented in Fig. 1. Here we introduced the main functions and characteristics of each class and the uniqueness and usefulness of some important functions in the classes.
microtable class
The microtable class is the basic class that is designed for creating the object storing the input files and for the subsequent data preprocessing. Compared with the phyloseq class in the phyloseq package, the microtable class is brief. For example, the file format required in the object creation is the user-friendly data.frame class for all the OTU table, taxonomy table and sample information table (see supplemental material Example 1). The built-in functions for data preprocessing includes tidy_dataset(), filter_pollution(), rarefy_samples(), merge_samples(), merge_taxa(). By applying function tidy_dataset(), all the basic files in the object can be trimmed and have consistent information. If the numbers of sequences across samples vary greatly, rarefy_samples() is recommended to make the sequence number equal for each sample to reduce the effects of sequencing depth on the alpha- and beta-diversity calculation. Running cal_abund() can automatically calculate the taxon abundance at each taxonomic rank and return a list (taxa_abund) containing all the tables. The function cal_alphadiv() is used to calculate the alpha diversity, including Chao1, Shannon–Wiener index, Simpson index and phylogenetic alpha diversity (Faith 1992). Furthermore, the function cal_betadiv() is available for obtaining the distance matrix of the beta diversity, including Bray–Curtis, Jaccard and UniFrac (Lozupone and Knight 2005; Jun Chen 2018).
trans_abund and trans_venn classes
To visualize taxa abundance changes across samples or groups, we created the trans_abund class to implement the plotting approaches for barplot, boxplot, heatmap and pie chart. The trans_venn class is designed for venn diagram to decipher the shared or particular taxa among groups, including the analysis of groups more than five. Compared with other venn analysis packages, such as VennDiagram (Hanbo Chen 2018), we also added the percentage of OTUs’ number or sequences’ number for each particular or shared fraction. Generally, researchers are interested in the taxa compositions of particular and common fractions of the venn analysis results. Thus, we developed the unique function trans_venn_com() to analyze the taxa compositions of each fraction (Mendes et al. 2011; Andrew et al. 2012). This function can transform the composition of OTUs in each fraction into the community format table and return a new microtable object for plotting the taxa composition of each fraction rapidly (see supplemental material Example 2).
trans_alpha class
The trans_alpha class is designed to transform the alpha diversity table of microtable object into other formats for statistics and plotting. The function cal_diff() can be performed for the differential test among groups by the Kruskal–Wallis rank sum test or ANOVA (analysis of variance) with multiple comparisons for all the alpha diversity indexes. The following function, plot_alpha(), is used to present the alpha diversity data and add the significance by either multiple comparisons or paired comparisons directly.
trans_beta class
A crucial part in the ecological analysis of microbial community data is the beta diversity, which can be carried out with trans_beta class. The trans_beta class currently implements several commonly used unconstrained ordination approaches, such as principal component analysis, principal coordinates analysis (PCoA) and non-metric multidimensional scaling. Besides, the distances among samples from the same group or different groups can indicate the community dissimilarity within or between groups (Lim et al. 2015). The function cal_group_distance() and plot_group_distance() can be performed for group distance transformations and plotting, respectively. The clustering analysis and plotting was also implemented in function plot_clustering(). Additionally, to simplify the use of permutational multivariate analysis of variance (perMANOVA), we developed the function cal_manova() for the overall comparisons, paired group comparisons and multiple factor comparisons of perMANOVA.
trans_diff class
The identification of the indicator taxon is significant for explaining the biological mechanisms underlying the difference in community structure among different groups. With the increase of the complexity of sequencing data in microbial community ecology, evaluating the taxa with significant differences across groups is a challenge. So some tools, e.g. LEfSe (Segata et al. 2011), combine the advantages of supervised machine learning and differential abundance test to identify important indicator taxa that distinguish the communities of one group from others. The trans_diff class currently implements three well-known approaches: LEfSe (Segata et al. 2011), random forest (An et al. 2019) and metastat (White, Nagarajan and Pop 2009). LEfSe combines linear discriminant analysis (LDA) and the differential abundance test to find indicator taxa. Compared with the Python version, the LEfSe approach in trans_diff class is reimplemented using R code to reduce the software dependence. Meanwhile, plotting of the results (barplot of LDA score and cladogram) are also easy to adjust using the package. We integrated the random forest analysis and differential abundance test to increase the convenience and power of the random forest approach compared with the randomForest package. The approach metastat is particularly useful in differential abundance test between two groups (White, Nagarajan and Pop 2009). We implemented the metastat method to reduce the workload of the file format transformation and added the corresponding plotting function.
trans_env class
Evaluating the effects of environmental factors on microbial community structure is critical in inferring the underlying mechanisms dominating community assembly. In the class trans_env, the redundancy analysis (RDA) and distance-based RDA (db-RDA) are implemented in cal_rda() based on the vegan function rda(). In order to make the visualization of RDA result easy, we reinforced the transformation and plotting approaches in the trans_rda() and plot_rda(), respectively. The Mantel test between environmental data and the beta-diversity distance matrix can be conveniently calculated using cal_mantel() for all environmental variables and all distance matrix. Correlation heatmap is an important method for presenting the correlations between taxa abundance and environmental factors. This can be done with two steps using trans_env class. First, the correlations between taxa and environmental factors can be calculated with cal_cor() for all samples or performed within different groups. Then, plot_corr() can be performed for the correlation heatmap using ggplot2 package style or the pheatmap package style with clustering plot, even though the correlation analysis is employed for different groups.
trans_nullmodel class
In recent decades, the integration of phylogenetic analysis and null model promotes the inference of niche and neutral influences on community assembly more powerfully by adding a phylogeny dimension (Stegen et al. 2013). The trans_nullmodel class provides an encapsulation, including the calculation of the phylogenetic signal, beta mean pairwise phylogenetic distance (betaMPD), beta mean nearest taxon distance (betaMNTD), beta nearest taxon index (betaNTI), beta net relatedness index (betaNRI) and Bray–Curtis-based Raup–Crick (RCbray). The approach for phylogenetic signal analysis is mantel correlogram (Liu et al. 2017). The algorithms of betaMNTD and betaMPD have been optimized to be faster than those in the picante package (Kembel et al. 2010). The combinations between RCbray and betaNTI (or betaNRI) can be used to infer the strength of each ecological process dominating the community assembly under the specific hypothesis (Stegen et al. 2013; Liu et al. 2017). This can be achievable by the function cal_process() to parse the percentage of each inferred process.
trans_network class
The analysis of co-occurrence patterns in microbial ecology is a hot topic and is usually analyzed by applying the network analysis. In trans_network class, three network construction approaches are provided: correlation network (Deng et al. 2012), SPIEC-EASI (SParse InversE Covariance Estimation for Ecological Association Inference) network (Kurtz et al. 2015) and probabilistic graphical models (PGM) based network (Tackmann, Matias Rodrigues and Mering 2019). For the correlation-based network (step 1 in Fig. 2), three correlation calculation methods are optional, including cor.test function, WGCNA package (Langfelder and Horvath 2008) and SparCC (Friedman and Alm 2012). The optimization of the correlation coefficient is optional by using the random matrix theory-based method (Deng et al. 2012). Another network construction approach is SPIEC-EASI in SpiecEasi R package (Kurtz et al. 2015). The third network approach we implemented is the PGM network, which can automatically invoke the Julia language in the system and the FlashWeave package (Tackmann, Matias Rodrigues and Mering 2019). Saving a network with the function save_network() and the default setting can generate the network.gexf file, which can be opened with Gephi software (https://gephi.org) for network visualization. The module classification, network topological attributes and node topological attributes calculations were all developed based on the igraph package (Csardi and Nepusz 2006). The fast greedy modularity optimization algorithm was implemented to classify modules with the cluster_fast_greedy() function in igraph. The first principal component of OTU abundance matrix from each module was extracted as the eigengenes of modules to reveal higher order organization information based on PCA analysis (Deng et al. 2012). The eigengene analysis of modules was provided in the function cal_eigen(). Furthermore, we developed the function subset_network() for extracting a part of the network and function cal_sum_links() for summing the linkages among taxa at any taxonomic rank.
Flowchart of processing in the trans_network class. The workflow starts with the trans_network object creation using the microtable object. The numbers represent the main operations in the class. (1): correlation and significance calculations; (2): construct correlation-based network; (3): construct SPIEC-EASI (SParse InversE Covariance Estimation for Ecological Association Inference) network with SpiecEasi R package; (4): construct probabilistic graphical model based network; (5): calculate network topological properties; (6): save network as gexf format file; (7): calculate node topological properties; (8): plot the results from (7); (9): sum the linkages among different taxa; (10): plot the results from (9); (11): extract a sub-network from the network according to the nodes; (12): calculate sub-network topological properties and (13): eigengene analysis of modules. The detailed code of each step was provided in supplemental material Example 3.
trans_func class
Functional analysis is an attractive part in the data analysis of microbial community, mainly coming from its interpretability on the biological questions. In microbial community ecology, exploring functional redundancy is a challenge because of the difficulty in identifying the functional traits of OTUs or species. To make this analysis feasible in the trans_func class, we developed the function cal_spe_func() to automatically match the taxonomic assignments of OTUs against the taxonomic information of the FAPROTAX database (Louca, Parfrey and Doebeli 2016) for prokaryotes or against the taxonomic information of the FUNGuild database (Nguyen et al. 2016) for fungi (step 1 in Fig. 3). FAPROTAX database was created based on the information from the Bergey's Manual of Systematic Bacteriology and related literatures, so the result of prokaryotic trait identification was reliable and conservative. Subsequently, the function cal_spe_func_perc() can be performed to calculate the percentages of taxa (abundance unweighted) or individuals (abundance weighted) associated with each trait in the communities or network modules (step 2 or 3 in Fig. 3). We considered higher percentage of a functional trait indicates higher functional redundancy of this trait in the community or network module.
Flowchart of processing in the trans_func class. The workflow starts with the trans_func object creation using the microtable object. The numbers represent the main operations in the class. (1): identify the traits of OTUs for prokaryotes or fungi; (2): calculate functional redundancy for communities; (3): calculate functional redundancy for modules; (4): plot the results from (2) or (3); (5): perform correlations between functional redundancy and environmental factors; (6): plot the correlation heatmap from (5); (7): predict metagenomic pathway abundance of prokaryotic community with Tax4Fun package; (8): create a new microtable object according to the pathway abundance table; (9): calculate pathway abundance at each KEGG (Kyoto Encyclopedia of Genes and Genomes) level; (10): perform the LEfSe analysis based on metabolic pathway abundance. The detailed code of each step was provided in supplemental material Example 4.
To predict and analyze the functional potential of communities easily, we integrated the Tax4Fun R package (Aßhauer et al. 2015) for the metagenomic metabolic pathway prediction (step 6 in Fig. 3) and FAPROTAX software (http://www.loucalab.com/archive/FAPROTAX) to predict the functional potential of communities on the biogeochemical cycles. After obtaining the functional potential of communities, it is easy to conduct the following statistical analysis by using microeco package (steps 7–9 in Fig. 3).
RESULTS AND DISCUSSION
The case studies in microeco package involved a 16S rRNA gene amplicon sequencing dataset from Chinese wetland soils (An et al. 2019) and an ITS amplicon sequencing dataset from the rhizosphere microbiome research (Gao et al. 2019). The detailed sample information was described in the sample information tables (sample_info_16S and sample_info_ITS). The 16S rRNA gene amplicon sequencing dataset was employed to demonstrate the examples in this article. Samples in the dataset can be classified into Chinese inland wetlands (IW), coastal wetlands (CW) and Tibet plateau wetlands (TW). A detailed online tutorial was provided to show usage of the microeco package (https://chiliubio.github.io/microeco/). Here, we focused on several important examples and the corresponding codes (see supplemental material), which can be applied to achieve the results in Figs 4 –6. The subset of this example data analysed here involved 90 samples, 30 in each group (CW, IW and TW). We first created the microtable object and performed the basic preprocessing operations as the codes indicate in supplemental material Example 1. In total, 13 296 OTUs remained after removing the taxa not assigned in the Kingdom ‘k__Archaea’ or ‘k__Bacteria’ or taxa with the assignment information ‘Mitochondria’ or ‘Chloroplast’. The OTU table was rarefied to make the sequence number same across samples (10 000 sequences in each sample). The taxa abundance, alpha diversity and distance matrix of beta diversity were all calculated for the downstream analysis (supplemental material Example 1).
Venn diagram (A), barplot of OTU compositions in venn analysis (B), alpha diversity (C), PCoA based on Bray–Curtis distance (D), betaNRI (E), LEfSe (F), abundance of the indicator taxa from LEfSe (G), correlation heatmap between indicator taxa and environmental factors (H). *: P < 0.05, **: P < 0.01; ***: P < 0.001. (a): the percentages represent the sequence number percentages of OTUs in this fraction. (b): the frequency in y axis represents the percentage of the OTUs number in each genus. (d): the ellipse represents the 95% confidence interval. (e): betaNRI, beta net relatedness index. Horizontal dashed lines denote values of +2 or −2, and the values beyond which an individual β-deviation value was considered statistically signifcant. (f): LDA score is the linear discriminant analysis score in LEfSe. (h): MAP, mean annual precipitation; MAT, mean annual temperature; TOC, total organic carbon; TN, total nitrogen.
Visualization of network (A), chord plot of summed linkages among phyla (B), roles of nodes in terms of within-module and among-module connectivity (C) and correlation heatmap between environmental factors and module eigengenes (D). *: P < 0.05, **: P < 0.01; ***: P < 0.001. (a): M represents the module, and M1–M5 represents module 1–5. (d): MAP, mean annual precipitation; MAT, mean annual temperature; TOC, total organic carbon; TN, total nitrogen.
The functional redundancy of network modules (A), correlation heatmap between environmental factors and functional redundancy of prokaryotic communities (B) and barplot of LDA score for metabolic pathways (C). *: P < 0.05, **: P < 0.01; ***: P < 0.001. The functional redundancy defined here is the percentage of OTUs (abundance unweighted) with a trait in the community or network module. (a): M represents the module, and M1–M5 represents module 1–5. (b): MAP, mean annual precipitation; MAT, mean annual temperature; TOC, total organic carbon; TN, total nitrogen. (c): LDA score is the linear discriminant analysis score in LEfSe.
First, several operations were shown in Fig. 4 and in supplemental material Example 2. We performed the venn diagram and analyzed the OTU compositions for each shared or particular fractions (Fig. 4A). The sequence number percentages of OTUs in venn diagram were added to illustrate the abundance information of OTUs in each fraction (Fig. 4A). The results indicated that the shared OTUs among CW, IW and TW and the abundance of those shared OTUs account for the most number and abundance of total OTUs. To explore the taxa composition in each fraction, we transformed the results and obtained the OTU compositions at genus level for each fraction (Fig. 4B). This approach can help explain whether some taxa (here genera) are comprised of particular or shared OTUs. For example, in the shared fraction CW-IW-TW, the Planctomyces genus has a relatively high percentage than other groups. The next example was the paired comparisons of alpha diversity between groups, which indicates that the Chao1 in IW is significantly higher than that in CW (Fig. 4C). To further show whether there is large differences on the beta diversity among the groups, the PCoA based on the Bray–Curtis distance (Fig. 4D) was conducted as shown in supplemental material Example 2. We also used trans_nullmodel to calculate the betaNRI to evaluate whether there is different phylogenetic turnover patterns within different groups. The result indicated that the strength of basic phylogenetic turnover in TW is significantly higher than those in CW and IW (Fig. 4E). LEfSe (Fig. 4F and G) were used to identify the indicator taxa, demonstrating that the indicator taxa were distinct across groups. For instance, the Proteobacteria and Firmicutes phyla were significantly enriched in CW (Fig. 4F and G). To further explore how environmental factors influence those indicator taxa, we employed trans_env class to correlate the abiotic factors to the abundance of these taxa and performed the correlation heatmap (Fig. 4H). The clustering of taxa and environmental factors were also evaluated to reveal the extent to which some taxa have more similar niche preferences than others. The results clearly illustrated that the correlation heatmap can offer useful aids to visualize how different indicator taxa may be affected by environmental factors (Fig. 4H).
To further study the species co-occurrence patterns, we applied trans_network class to the correlation network analysis (supplemental material Example 3). The correlation network was constructed with the specific threshold (spearman correlation: R = 0.7, P = 0.01) and visualized with gephi (Fig. 5A). We performed chord plot to show the linkage number within and between phyla (Fig. 5B). To explore the node importance in network, the node roles were assigned according to the within-module and among-module connectivity (Fig. 5C). This approach has been widely applied in the network analysis and was validated to show a comprehensive evaluation on the node roles in network (Deng et al. 2012; Shi et al. 2016). Then the module eigengenes were calculated and correlated to the abiotic factors (Fig. 5D), revealing that OTU compositions of modules may be dominated by different environmental factors. It is noted that, under the condition of strong phylogenetic clustering imposed by abiotic filtering, the correlation network may be comprised of a large number of edges unassociated with the biotic interactions. So, it is better to use correlation network to explain the species niche similarities and the underlying mechanisms instead of the biotic interactions (Röttjers and Faust 2018).
Functional analysis is of great significance in unraveling the mechanisms structuring species co-occurrence patterns and community assembly. Previous studies scarcely ever decipher the functional profiles of microbial co-occurrence network modules, mainly resulting from the methodological difficulty. We used FAPROTAX database to identify species traits and analyzed the functional redundancy in the network modules (Fig. 6A). The results showed that M2 (the second largest module) has a high percentage of populations (OTUs) with aerobic chemoheterotrophy function. There was no N-cycle and S-cycle related populations found in M4. Then we calculated functional redundancy in communities and correlated the functional redundancy to the abiotic factors to explore how abiotic factors may affect the functional redundancy (Fig. 6B). The results indicated that conductivity and MAT (mean annual temperature) both had extremely significant positive correlations with the functional redundancy associated with sulfite respiration and photoautotrophy. Finally, we used the integrated function to predict the abundance of metabolic pathways and performed LEfSe among CW, IW and TW (Fig. 6C). The highest LDA score in CW, IW and TW are genetic information processing, metabolism and environmental information processing, respectively, meaning that the enriched and indicative pathways in each region reflected the unique biological and gene patterns related with the community structures. This predictive and meaningful result may point toward more promising and deeper consideration on the design and data mining of metagenomic sequencing research.
It is necessary to note that statistical and plotting approaches are extremely diverse in microbial ecology with the increase of complexity and depth of research, especially in the recent decade. Even though we have packaged many frequently-used and unique approaches in microeco package, it is impossible to consider all approaches and meet all the requirements in this field. Our primary philosophy in the design of the microeco package is simplicity without loss of power. So the key point is how to make researchers to perform the data mining in the field fast and conveniently. We selected the advanced and representative approaches to increase the power and classified those approaches into several classes approximately to make the framework clear.
One of the goals creating the microeco project is to help researchers close the gap that is ever-present between the difficulty in learning of analysis methods and the demand on the fast data mining with the combination of many approaches, especially for the novice programmers. Currently, it is achievable to use microeco package to explore data in many hot and significant scientific questions, e.g. species co-occurrence patterns (Fig. 2) and functional profile (Fig. 3). Another advantage using microeco package is that it is highly reproducible for the investigators and convenient to describe the complex approach in publications.
CONCLUSIONS
The microeco package features (i) high modularity in that each class covers each part of the analysis, and the structure and functions of the package are easy to understand, remember, search and use; (ii) high flexibility in that diverse algorithms or approaches are implemented in each part to make the analysis flexible; (iii) fast acting in that some algorithms have been optimized; and (iv) power in that several cutting-edge approaches are implemented, such as the LEfSe, network analysis, null models and functional redundancy calculation.
AUTHORS’ CONTRIBUTIONS
CL, YMC, XZL and MJY conceived the ideas and designed the methodology; CL and MJY created the package; CL, YMC, XZL and MJY wrote the manuscript.
ACKNOWLEDGEMENTS
This work was supported by the National Natural Science Foundation of China (31670503), the 13th Five-year Informatization Plan of Chinese Academy of Sciences (XXH13503-03-106) and China Biodiversity Observation Networks (Sino BON).
Conflicts of Interest
None declared.
