Whole Genome Shotgun (WGS) metagenomics is increasingly used to study the structure and functions of complex microbial ecosystems, both from the taxonomic and functional point of view. Gene inventories of otherwise uncultured microbial communities make the direct functional profiling of microbial communities possible. The concept of community aggregated trait has been adapted from environmental and plant functional ecology to the framework of microbial ecology. Community aggregated traits are quantified from WGS data by computing the abundance of relevant marker genes. They can be used to study key processes at the ecosystem level and correlate environmental factors and ecosystem functions. In this paper we propose a novel model based approach to infer combinations of aggregated traits characterizing specific ecosystemic metabolic processes. We formulate a model of these Combined Aggregated Functional Traits (CAFTs) accounting for a hierarchical structure of genes, which are associated on microbial genomes, further linked at the ecosystem level by complex co-occurrences or interactions. The model is completed with constraints specifically designed to exploit available genomic information, in order to favor biologically relevant CAFTs. The CAFTs structure, as well as their intensity in the ecosystem, is obtained by solving a constrained Non-negative Matrix Factorization (NMF) problem. We developed a multicriteria selection procedure for the number of CAFTs. We illustrated our method on the modelling of ecosystemic functional traits of fiber degradation by the human gut microbiota. We used 1408 samples of gene abundances from several high-throughput sequencing projects and found that four CAFTs only were needed to represent the fiber degradation potential. This data reduction highlighted biologically consistent functional patterns while providing a high quality preservation of the original data. Our method is generic and can be applied to other metabolic processes in the gut or in other ecosystems.
Whole Genome Shotgun (WGS) metagenomics is increasingly used to study the structure and functions of complex microbial ecosystems. Thanks to a huge research effort, these high-throughput approaches are constantly improving and have already provided valuable information. Recent achievements [1–3] for analyzing metagenomic reads have focused on genome or species reconstruction. They constitute an alternative to amplicon sequencing approaches such as 16S rDNA and provide potentially more robust tools for community structure assessment. As a complement to these taxonomy oriented approaches, WGS metagenomics allows the direct functional profiling of microbial communities through gene inventories of otherwise uncultured microbial communities. The functional annotation of these genes, as well as the evaluation of their abundances requires the mapping of metagenomic reads on one or several functional databases which may include annotated, de novo assembled gene catalogs. It is a challenging task for which bioinformatic tools are still being developed (see [4] for a very recent review). As early as 2005, Tringe et al. [5] demonstrated that functional binning could discriminate environments and determine the functional potential of microbial communities. Since then, several other groups used functional data to shed light on the processes taking place in the ecosystem, even in the absence of additional expression data (see e. g. [6,7]), and developed frameworks to functionally characterize and compare ecosystems [4,8]. Investigating complex microbial communities from a purely functional point of view provides an interesting insight because their taxonomic diversity may be very high and varying in time while, owing to environmental selection pressure, their ecosystemic functions are often ubiquitous and much more stable in time. The functional profiling of microbial communities has motivated the adaptation of concepts borrowed from environmental and plant functional ecology to the framework of microbial ecology, such as community aggregated traits [9,10], which can be quantified through abundance computation of specific genes or predefined pathways from WGS data. These taxon-free approaches allow the study of key processes at the ecosystem level and correlate environmental factors and ecosystem functions [11]. Choosing relevant descriptors for ecosystem functions raises several issues. Using arbitrarily predefined pathways helps simplifying the systems description by creating combined descriptors, however it may not reflect all the information in the data. On the other hand, unsupervised statistical analyses such as e. g. principal component analysis (PCA) are quite often used to extract low-dimensional description of these high-dimensional datasets. However, they do not exploit all the richness of available information about biological processes at the individual or collective level, and may lead to irrelevant or hardly interpretable results. To overcome these issues, we propose a novel model based approach, mixing satistical data analysis and system biology. We focus on specific ecosystemic metabolic processes, described by a possibly high number of genes or functional markers abundances. Our aim is to infer a limited number of combined functional traits characterizing these processes, as well as their intensity in the ecosystem. Our approach relies on Non-negative Matrix Factorization (NMF), a popular machine learning technique in data and image analysis. NMF, together with PCA belongs to a wide family of data analysis methods designed to solve blind source separation problems. These questions are widely encountered and cover all situations where a mixture of several signals originating from different unknown sources is observed (here gene abundances in samples) and one wishes to separate them and identify the sources (the combined functional traits) and the mixture coefficients. NMF was previously used for genomic data mining in the context of microarray data analysis [12,13]. More recently, it was introduced in the context of metagenomics for reads binning [14], or analyzing datasets from various ecosystems. In particular Jiang and co-authors used this approach as a “soft” clustering tool in different ecosystems. They studied pathway abundances in diverse environmental ecosystems [15], compared habitats in marine ecosystems using protein families profiles [16] and human body sites using phylogenetic and functional data [17]. In [18] the idea of coupling source separation analyses with prior knowledge was developed in the context of bipartite network reconstruction, using prior knowledge on the network structure, and applied to the reconstruction of regulatory signals from microarray data. The originality of our work lies in the design of a constrained NMF approach where the constraints aim at selecting biologically relevant combined functional traits (the sources) to describe processes at the ecosystem level. The constraints are derived from available prior knowledge in a Bayesian perspective. Moreover, we propose a careful multicriteria selection procedure to select the relevant number of combined functional traits. As a proof of concept, we applied our approach to the modelling of ecosystemic functional traits of fiber degradation by the human gut microbiota based on 1408 samples of gene abundances.
Quantitative metagenomics allows the study of metabolic processes at the ecosystem level, by producing metagenes abundances with functional annotation. Taking advantage of this information, we consider a metabolic process occurring in an ecosystem. We assume that this process is described by a list of biochemical reactions, each associated to what we call a functional marker. A functional marker is a group of genes, or modular elements in genes, able to control the production of enzymes involved in a reaction. Typical functional markers are Kegg Orthologies (KO). The abundance of a functional marker is defined as the sum of the abundances of all the metagenes identified as belonging to the group. Our first modelling hypothesis is the existence of patterns underlying functional markers abundances in the ecosystem. Indeed, within a microbial ecosystem, genes can be looked at in a hierarchical fashion. They are first associated on microbial genomes, which are further grouped at the ecosystem level into subcommunities of microorganisms. These subcommunities involve hundreds of different bacterial species linked by complex co-occurrences or interactions. Moreover they are influenced by environmental factors (e. g. nutrients or temperature). The set of functional marker abundances within each of these subcommunities form a characteristic pattern in the ecosystem. This hypothesis is illustrated in the simplistic example of metabolic process in Fig 1. Here, the observed functional marker abundances originate from three subcommunities, resulting in three patterns. The ecological interpretation of the first pattern (green box) could be a trophic chain in which metabolite X is produced by bacterial species of type 1 and used by bacterial species of type 3 as well as an association between bacterial species of type 1 and 2 for other ecological reasons which cannot be elucidated with the observation of this particular metabolic process only. Therefore functional marker abundances in samples should reflect this structure. Our second modelling hypothesis is motivated by the fact that even if in many microbial communities microbial species composition varies over time and samples, metabolic processes at the community level are ubiquitous and much more stable in time. Consequently, we assume that the patterns mentioned above are shared by all samples of the ecosystem, in variable proportions because they were selected by specific environmental constraints. For instance, for the gut ecosystem these include anaerobic conditions, temperature or digestion residue composition coming from the host diet. As illustrated on Fig 2, when considering the abundances of functional markers in several samples, each of these patterns can be characterized by a line vector of functional marker frequencies, which we define to be a Combined Aggregated Functional Trait (CAFT), in the meaning that it represents a combination of quantifiable characteristics (functional markers) of aggregated microbial communities. Following our hypotheses, functional marker abundances in each sample are modelled as mixtures of underlying CAFTs. For each sample i = 1, …, n and marker j = 1, …, r, the observed frequency aij of marker j in sample i is the sum of the contributions of all CAFTs. This leads to a i j ≃ ∑ l = 1 k w i l h l j (1) where hl = (hl1, …, hlr) is the characterization in terms of r markers frequencies of the lth CAFT and wi = (wi1, …, wik) is the vector of the intensities or weights of the k CAFTs in sample i. Note that each hl is defined up to a multiplication by a positive constant, provided that its weights in the samples are scaled accordingly. The ℓ1-scaled version hl/ (∑j hlj) can be interpreted as the vector of relative frequencies of functional markers within the subcommunity characterized by trait l. In the rest of the paper, we denote the (n × r) data matrix by A, the (k × r) matrix whose lines are the hl by H and the (n × k) matrix whose lines are the wi by W. H will be called the trait matrix and W the weight matrix. By definition the entries of the weight and trait matrices W and H are assumed to be non-negative, we denote this condition by W ≥ 0, H ≥ 0. (2) We now show how prior knowledge on microbial genomics and metabolism can be exploited in order to select biologically relevant CAFTs. A first important consequence of our modeling framework is that genomic information can be translated into constraints on CAFTs structure. Indeed, by construction, the CAFT hl characterizes a gene pattern associated in each sample i to a subcommunity Cil. Therefore, in each of these subcommunities, the total abundance of any functional marker j is proportional to hlj. Thus, denoting nxj the number of functional markers j in the genome of a microorganism x in the ecosystem, it follows that ∑ x ∈ C i l n x j = w i l h l j. (3) From Eq (3) we observe that any constraint on the structure of the genomes (the nxj) satisfied by all x in the ecosystem induces a constraint on the CAFTs (the hlj). As an example, let us assume that two markers j and j′ are such that for each microorganism in the ecosystem, if j′ is not present in its genome, then neither is j. This implies that for each microorganism x, there is a non-negative constant Rx, j, j′ satisfying n x j = R x, j, j ′ n x j ′, (4) where we arbitrarily set Rx, j, j′ = 0 when nxj = nxj′ = 0. Then as a consequence of Eq (3) we have for each CAFT l h l j ≤ δ h l j ′, (5) where δ = supx ∈ Ecosystem Rx, j, j′. In the absence of any other information, exploiting genomic information only could lead to a large number of constraints. For instance, constraining the ratio of abundances in pairs of markers as above requires assessing the existence of constants such that Eq (4) holds for r (r − 1) ordered pairs, for each genome x. This is why we propose to combine genomic information with metabolic information, when possible, in order to select a limited but relevant set of constraints on the CAFTs. In order to do this, we formulate two additional assumptions. First, we assume that a subset of reactions in the metabolic process under consideration are well known, so that a list of metabolites that are consumed or produced by each reaction can be extracted. Second, we assume that metabolites are partitioned in two classes. The first one corresponds to metabolites that are known to be exported out of the cell, from experimental evidence or published data. The second one gathers metabolites that are actually known to stay within microbial cells, as well as metabolites for which this is strongly suspected. For convenience, these two classes will be called the extracellular and intracellular metabolites, although this may be misleading in some cases. Following similar steps as in eqs (3–5), a subset M1 of the class of intracellular metabolites can be derived from genomic information on the microorganisms present in the ecosystem, as well as a positive constant δ1, m for each m in M1, such that in each CAFT hl the marker frequencies associated to the production or consumption of m should satisfy ∑ j producing m h l j ≤ δ 1, m ∑ j' consuming m h l j ′ (6) In the same way, there is a subset M2 of intracellular metabolites such that for each m in M2, each CAFT hl should satisfy ∑ j' consuming m h l j ′ ≤ δ 2, m ∑ j producing m h l j (7) The intuitive idea behind such constraints is the following. Assuming that functional markers are relevant proxies of metabolic functions, it could be expected that if a genuinely intracellular metabolite is produced in the subcommunity characterized by a given CAFT it has to be used. This inspires constraint (6), since it imposes that the number of markers associated to reactions that produce a given intracellular metabolite is zero if the number of markers associated to reactions that consume the same intracellular metabolite is zero. Conversely, it could be expected that a metabolite cannot be used if it is not produced, which is enforced by constraint (7). This justifies the absence of constraints imposed on functional markers associated to reactions consuming or producing for the class of extracellular metabolites as the latter can actually be imported or exported outside bacterial cells in the ecosystem. However, the set of metabolites and reactions under consideration is defined in the context of a specific metabolic process, given the available knowledge. Thus it may be that in this process, some intracellular metabolites accumulate in the cells. It may also be that the selected reactions ignore metabolic routes, because they are unknown or correspond to other metabolic processes, that could produce or use the metabolites, and finally it may also happen that we were wrong in assuming that some poorly studied metabolites are intracellular. Consequently, depending on metabolites, either none, one or both constraints (6) and (7) may be satisfied. In the case where both constraints apply, the frequencies of markers associated to production and consumption are both either zero or positive with consumption over production marker frequency ratio lying in [1/δ1, m, δ2, m]. The theoretical construction of the two sets M1, M2 and the positive constants δ1, m and δ2, m is detailed in the Materials and Methods section, as well as their practical approximation from known bacterial genomes. The set of constraints (6) and (7) can be expressed in matrix notation as F Δ H T ≤ 0 (8) where Δ is the set of parameters Δ = {δ1, m, m ∈ M1} ∪ {δ2, m, m ∈ M2} and FΔ is a coefficient matrix depending on Δ. Note that, although they include metabolic information and take the familiar form of linear inequality constraints (6) and (7) are constraints on gene frequencies in subcommunities, they are not constraints on metabolic fluxes. The constrained NMF algorithm as well as an unconstrained version were implemented on a toy model in order to illustrate the improvement in terms of CAFT analysis brought by the constraints. We considered a microbial ecosystem in which a metabolic process characterized by a simple 5 reactions graph takes place, each reaction being associated to a functional marker. The metabolic process was assumed to be structured in two CAFTs. The first CAFT characterizes a subcommunity in the ecosystem in which all the functional markers are present in equal proportion, while the second CAFT characterizes a subcommunity in which only three functional markers are present, in equal proportion. The metabolic graph and the true CAFTs are displayed in Fig 3A. Thus, the true CAFT composition matrix was taken equal to H ¯ = 1 1 1 1 1 1 0 0 1 1 We built a matrix A for 80 samples assumed to be mixtures of CAFT 1 and CAFT 2. We first generated a matrix W ‾ of abundances of CAFT 1 and CAFT 2 for the 80 samples, with entries randomly drawn from a uniform distribution, and computed the corresponding marker count matrix B = W ‾ H ‾. Noisy counts were generated from B according to the formula a i j = b i j (1 + 0. 2 ϵ i j) where the ϵij are independent, centered Gaussian variables with unit variance. In order to exemplify the construction of the constraints, we assumed that a set of three highly prevalent bacteria in this ecosystem was known. Bacteria are represented by the vectors x1, x2 and x3 of functional markers counts in their genomes x 1 = (1,1, 1,0, 0) x 2 = (1,0, 0,1, 0) x 3 = (0,0, 0,0, 1) We assumed that based on published data on the physiology of these bacteria and possibly other bacteria in the ecosystem that are less prevalent, we knew that metabolites U, Z and T were extracellular. For instance, these data may come from in vitro culture experiments. We also assumed that the metabolites V, X and Y were not detected in the extracellular medium in these type of experiments, and considered them as intracellular by default. We used these representative genomes, together with the metabolic graph of Fig 3A and applied the rules described in the Materials and Methods section to build a set of constraints on the CAFTs. Considering metabolite V, it is known to be produced by the reaction associated to marker a and consumed by the reactions associated to b and d. The inspection of the three genomes led to n 1 a = 1 = 1 × (n 1 b + n 1 d) in x 1, since n 1 b + n 1 d = 1 n 2 a = 1 = 1 × (n 2 b + n 2 d) in x 2, since n 2 b + n 2 d = 1 n 3 a = 0 = 0 × (n 3 b + n 3 d) in x 3, since n 3 b + n 3 d = 0 Therefore Eq (11) is satisfied on the three reference genomes. Thus, V belongs to M ˜ 1 with δ1V = max (1,1, 0) = 1. As described in the Materials and Methods section, we compensated for calibrating δ1V on a limited number of genomes instead of all the genomes in the ecosystem by applying a multiplicative factor of 2, which seems sensible here since the highly prevalent genome have at most one copy of each functional marker in their genome. We obtained the following constraints h l 1 ≤ 2 (h l 2 + h l 4) Similarly, X belongs to M ˜ 1 with hl2 ≤ 2hl3. However, Y is not in M ˜ 1 as Eq (11) is not satisfied for genome x2, here 33% of the reference genomes. In the same way, V and X belong to M ˜ 2 with hl2+hl4 ≤ 2hl1 and hl3 ≤ 2hl2, and Y is discarded. We compared the CAFTs reconstructed using the constrained (Fig 3B) and unconstrained (Fig 3C) NMF inference. The NMF inference with constraints recovers the true CAFTs, while the unconstrained inference fails to provide relevant biological results. Besides, both constrained and unconstrained NMF provide a similar reconstruction of the data, as shown in Fig 3D. Indeed, from a mathematical point of view, the reconstructed CAFTs with and without constraints span the same linear space. Therefore, while both constrained and unconstrained NMF provide usable results for biological sample analysis, only the constrained version is relevant for CAFT inference.
We presented a new NMF based approach to investigate metabolic processes in microbial ecosystems using quantitative metagenomic data. As mentioned in the introduction, NMF techniques were previously used as a “soft” clustering tool to compare samples from various ecosystems, ranging from marine environment to human body ([15,16], [17]). However, our work focused on different aspects of NMF. Rather than clustering biological samples, we aimed at extracting biologically relevant features associated to a metabolic process, based on ecological modelling. Therefore our main interest, from the mathematical point of view, was the inference of a biologically interpretable trait matrix H rather than sample analysis through the weight matrix W. We proposed to model the functional diversity of an habitat-specific microbial ecosystem as a mixture of contributions from subcommunities, each characterized by a profile of functional marker frequencies. These profiles were named Combined Aggregated Functional Traits in reference to [10], since they are measured at the community level from a random sample of microorganisms, regardless of their taxonomic identities. More precisely they result from both a combination of functional markers to form coherent and operative metabolic pathways, and an aggregation of individuals to form a subcommunity characterized by this trait. Our main contribution was the design of a novel constrained NMF model for CAFT inference. Standard NMF approaches were popularised for biological data analysis mainly because usual dimension reduction techniques such as truncated Singular Value Decomposition (SVD), or the closely related Principal Componant Analysis (PCA) often lead to negative coefficients, which do not offer straightforward interpretations. However, while SVD or PCA provide a unique and well defined reduced dataset, optimal in the ℓ2-error sense, NMF is an intrinsically ill-posed problem with multiple solutions. Imposing constraints reflecting biological knowledge greatly helps to select a relevant solution. This was illustrated through the toy example in Fig 3, where both constrained and unconstrained NMF provided solutions leading to the same reconstructed data. However, only the constraints allow the solution that makes biological sense to be picked. Our approach is in essence Bayesian, and aims at taking advantage of prior information, when available, to explore the data. We proposed a method that exploits knowledge about the genomic structure and metabolism to build relevant constraints. Genes frequencies are not randomly distributed in the metagenome and result from metabolic association encoded in microbial genomes. Thus, constraints resulting from this metabolic structuration scale up at the community level and should be accounted for in the determination of CAFT. We designed local constraints, only involving adjacent functional markers in the graph representation of the metabolic process. Hence the constraints are not strongly restrictive and still leave a considerable freedom in assembling reactions to build the CAFT matrix, without ever forcing hard reconstruction of pathways. Our method is flexible since it allows imposing constraints only on parts of the metabolic process where biological knowledge is available and considered reliable. Moreover, it could be easily adapted to include other types of prior informations than the one we considered. For instance it should be possible, in contexts where more prior information is available, to design constraints on the weight matrix W, as proposed in [18] for microarray data analysis. We considered a new approach for selecting the number of CAFTs, based on biological reproducibility instead of the numerical stability usually considered in literature. We proposed a criterion which evaluates the concordance between CAFTs computed on independent data sets, and combined it with classical procedures. We implemented our approach on metagenes abundances from 1408 samples of the human gut microbiota, in order to determine CAFTs associated to fiber degradation in the distal gut. Fiber degradation was chosen because it is a major function of the gut microbial ecosystem. It involves the anaerobic fermentation of simple sugar which has been largely documented in the literature, allowing the design of constraints related to this step. The constraints were built by selecting 190 reference microbial genomes from prevalent species of the human gut microbiome, for which marker genes of anaerobic fermentation were annotated. This involved manual curation and relied on both state of the art knowledge of metabolism and accuracy of available annotations. We proposed an approach to account for missing information. In particular, genes were missing in pathways where they were expected to be present, which could be explained by annotation errors, or by the presence of alternative pathways undocumented in the literature. It could also come from unmodelled pathways, not related to fermentation, but sharing some common steps. Then, constraining missing genes markers would be an error. Our approach carefully took into account annotation, modelling and knowledge based errors by adjusting the constraints. We found that our ecological model of CAFTs was consistent with the data, and that the metagenomic potential for fiber degradation by the human gut microbiota could be interpreted as a mixture of 4 CAFTs, shared by all samples in variable proportions. The CAFT detailed in the Results section is the one including the genes involved in methane production, which can be detected in 50% of a healthy population through a breath test. As a matter of fact, none of the three other CAFTs showed a high frequency of this set of genes responsible for the conversion of H2 and CO2 to methane. Two of them are prevalent among individuals, characterized by their richness in terms of GH and fermentation pathways. They differ slightly in GH repertoire and downstream fermentation genes, which may be linked to dietary habits. Finally a fourth profile, less prevalent among samples, might be more represented in diseased individuals. This data reduction highlighted biologically consistent functional patterns while providing a high quality preservation of the original data. Provided that prior knowledge is available, the generic framework we designed can be applied to other metabolic processes in the gut or in other ecosystems. By acknowledging fundamental data structuration our approach enables the inference of meaningful functional traits.
The columns of matrices A and W were ℓ1-normalized; the normalized matrices correspond to the repartition of the functional markers (resp. the CAFTs) abundancies among the samples. The distance matrix dA (resp. dW) of A (resp. W) defined as the 1408 × 1408 matrix whose element (i, j) is the euclidean distance between rows i and j of A (resp. W) was computed. For each biological sample i, the Pearson correlation between the distance of i to the other samples in A and in W was evaluated: C i = cor d i, j A j ≠ i, d i, j W j ≠ i and the histogram of (Ci) i = 1, …, 1408 was displayed.
Microbial communities are highly complex ecosystems, harboring a high taxonomic diversity, rapidly varying in time. Focusing on the functional traits of microorganisms instead of species offers an interesting and complementary insight since quite often, owing to environmental selection pressure, traits are less variable in time and shared by many species. By providing gene inventories of otherwise uncultured microbial communities, Whole Genome Shotgun (WGS) metagenomics made the functional profiling of microbial communities possible, since community aggregated traits can be quantified by the abundance of relevant marker genes. Considering a global metabolic process in an ecosystem, we propose a method that exploits gene abundance data to determine combinations of functional traits characterizing this process. When applied to fiber degradation in the human gut microbiota, the method shows that only four complex functional traits are needed to characterize this process. The approach is generic and could be applied to other processes or ecosystems. It allows summarizing complex datasets by a limited number of biologically interpretable functional patterns.
biochemistry
ecology and environmental sciences
metagenomics
functional genomics
metabolites
microbial ecosystems
ecology
ecosystems
metabolic processes
genetics
biology and life sciences
genomics
ecosystem functioning
metabolism
fermentation
Inferring Aggregated Functional Traits from Metagenomic Data Using Constrained Non-negative Matrix Factorization: Application to Fiber Degradation in the Human Gut Microbiota