TECHNOLOGY AND CODE published: 28 January 2021 doi: 10.3389/fmicb.2021.609048 kernInt: A Kernel Framework for Integrating Supervised and Unsupervised Analyses in Spatio-Temporal Metagenomic Datasets Elies Ramon 1* , Lluís Belanche-Muñoz 2 , Francesc Molist 3 , Raquel Quintanilla 4 , Miguel Perez-Enciso 1,5 and Yuliaxis Ramayo-Caldas 4 1 Plant and Animal Genomics, Statistical and Population Genomics Group, CSIC-IRTA-UAB-UB Consortium, Centre for Research in Agricultural Genomics (CRAG), Bellaterra, Spain, 2 Department of Computer Science, Polytechnic University Edited by: of Catalonia, Barcelona, Spain, 3 Schothorst Feed Research B.V., Lelystad, Netherlands, 4 Animal Breeding and Genetics Isabel Moreno Indias, Program, IRTA, Caldes de Montbui, Spain, 5 Catalan Institution for Research and Advanced Studies (ICREA), Barcelona, University of Málaga, Spain Spain Reviewed by: Xiaoquan Su, Qingdao University, China The advent of next-generation sequencing technologies allowed relative quantification Jun Chen, of microbiome communities and their spatial and temporal variation. In recent years, Mayo Clinic, United States Yang Dai, supervised learning (i.e., prediction of a phenotype of interest) from taxonomic University of Illinois at Chicago, abundances has become increasingly common in the microbiome field. However, a gap United States exists between supervised and classical unsupervised analyses, based on computing *Correspondence: ecological dissimilarities for visualization or clustering. Despite this, both approaches Elies Ramon eramongurrea@gmail.com face common challenges, like the compositional nature of next-generation sequencing data or the integration of the spatial and temporal dimensions. Here we propose a kernel Specialty section: This article was submitted to framework to place on a common ground the unsupervised and supervised microbiome Evolutionary and Genomic analyses, including the retrieval of microbial signatures (taxa importances). We define Microbiology, two compositional kernels (Aitchison-RBF and compositional linear) and discuss how a section of the journal Frontiers in Microbiology to transform non-compositional beta-dissimilarity measures into kernels. Spatial data is Received: 22 September 2020 integrated with multiple kernel learning, while longitudinal data is evaluated by specific Accepted: 07 January 2021 kernels. We illustrate our framework through a single point soil dataset, a human dataset Published: 28 January 2021 with a spatial component, and a previously unpublished longitudinal dataset concerning Citation: Ramon E, Belanche-Muñoz L, pig production. The proposed framework and the case studies are freely available in the Molist F, Quintanilla R, Perez-Enciso M kernInt package at https://github.com/elies-ramon/kernInt. and Ramayo-Caldas Y (2021) kernInt: A Kernel Framework for Integrating Keywords: microbiome, metagenomics, kernel, supervised, unsupervised, spatio-temporal, SVM, kPCA Supervised and Unsupervised Analyses in Spatio-Temporal Abbreviations: ANN, artificial neural network; ASV, amplicon sequence variant; JSK, Jensen-Shannon Kernel; kPCA, Metagenomic Datasets. kernel principal components analysis; MDS, multidimensional scaling; MKL, multiple kernel learning; NGS, next-generation Front. Microbiol. 12:609048. sequencing; NMSE, normalized mean squared error; OTU, operational taxonomic unit; PCA, principal components analysis; doi: 10.3389/fmicb.2021.609048 PCoA, principal coordinates analysis; RBF, radial basis function; RF, random forests; SVM, support vector machines Frontiers in Microbiology | www.frontiersin.org 1 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis INTRODUCTION supervised methods is the identification of microbial signatures (i.e., taxa that are predictive of a certain phenotype), which The microbiome is defined as the ensemble of microorganisms may enable a biological interpretation of the results. RF are and their genomes in a given environment. Microorganisms endowed with variable importance measures that can be used to are present in ecological niches as diverse as soil, oceans, this effect, while there is not such straightforward heuristic for freshwater, plants, and animals, but a large fraction of these ANN, although several possible strategies exist (Ibrahim, 2013). taxa cannot be cultivated with culture-dependent methods. The Another supervised method, selbal (Rivera-Pinto et al., 2018), is advent of next-generation sequencing (NGS) revolutionized this focused on the identification of microbial signatures based on field by allowing the massive sequencing and quantification of balances (i.e., the geometric means of data from two groups of microbial habitats. taxa), and has the particularity of being purely compositional. Proper analysis of microbiome data is challenging for a variety As microbial communities are highly dynamic systems, it of reasons. Abundance data obtained with NGS is multivariate, is important to address their spatial and/or temporal variation sparse and compositional in nature (Gloor et al., 2017). Also, (Berg et al., 2020). In spatial-structured studies, repeated samples microbial communities are very dynamic biological systems, of different sites (e.g., body sites, depth layers) are obtained from thus justifying spatial or time-course studies (Bodein et al., the same individuals or entities, thus raising the question of how 2019; Berg et al., 2020). The first approach on the field used to integrate them. A more general challenge is the integration statistical tools from standard ecological studies (Gloor et al., of datasets coming from different sources (e.g., “omics”), which 2017). For example, one of the first steps in nearly all microbiome may have different data types. Several statistical methods have studies consists in computing alpha and beta-diversities. Beta- been proposed to solve this question in the microbiome field. diversity measures, e.g., Bray-Curtis or Unifrac, quantify the Some examples are Link-HD (Zingaretti et al., 2020), mixKernel difference in diversity between samples from different habitats. (Mariette et al., 2018), and MOFA (Argelaguet et al., 2018), all They are used for clustering analysis or, more commonly, focused in the unsupervised learning setting. In most supervised for visualization techniques like principal coordinates analysis methods, this integration is usually performed at the input (PCoA) or multidimensional scaling (MDS). However, this data level (early integration), for example by concatenating the approach has been challenged, as the abundance data obtained by datasets; or after the model is built (late integration), combining NGS has a particular nature. The total number of reads delivered their scores as in ensemble methods. However, early integration is bounded by an uninformative sum: the library size (i.e., the may be not possible if data nature differs across sources number of total reads per sample). Library size is uninformative (Schölkopf et al., 2004). The case of the longitudinal studies because it does not contain information about the population. (which follow the evolution over time of microbial communities) Instead, it is arbitrarily fixed by the sequencing process and is more complex. Typically, longitudinal data is modeled by may vary by orders of magnitude across samples (McMurdie fitting a function (e.g., polynomial interpolation, splines) to the and Holmes, 2014). This kind of data is called compositional data points over time. To date, there exist few analytical tools for and deserves a specific mathematical treatment (Gloor et al., this kind of data in the microbiome field. Two examples can be 2017). In the case of metagenomics, extensive research is being found at Bodein et al. (2019) and Coenen et al. (2020), but they done to translate current statistical techniques to this paradigm are restricted to unsupervised analysis. (Gloor et al., 2017; Silverman et al., 2017; Rivera-Pinto et al., Difficulties like the compositionality of data or how to 2018). One example is the proposal of using the compositional accommodate the spatial and temporal dimensions affect Aitchison distance instead of the classic beta-diversity measures supervised and unsupervised methods alike. However, there is a (Quinn et al., 2018). gap between the most widely used supervised learning methods In machine learning, the aforementioned clustering, and the unsupervised analyses typical of the microbiome field ordination and visualization techniques belong to the so-called (Figure 1A). Libraries like QIIME2 juxtapose traditional analyses unsupervised learning. Supervised learning, which is focused (e.g., PCoA) with many different and powerful prediction on prediction, is not so widespread in microbiome analysis yet, algorithms, but both branches remain independent from a but the number of studies using this kind of approach is rapidly mathematical point of view. It is true that some beta-diversity growing in the last years (Zhou and Gallins, 2019). Due to this dissimilarity-based engines can be used as classifiers (Su et al., rise in popularity, widely used libraries for microbiome analysis 2020; Shenhav et al., 2019). However, as these tools are strongly like QIIME2 (Bolyen et al., 2019) now include plugins for focused on distinguishing among a limited number of bacteria- supervised learning in their toolbox. Typical available methods related conditions, they are not aimed at regression problems, include random forests (RF), artificial neural networks (ANN), nor do they give any information about the microbial signatures. support vector machines (SVM), and ridge regression (Qu et al., We consider that carrying out all aforementioned analyses in 2019; Zhou and Gallins, 2019; Namkung, 2020). Among the a common mathematical framework would provide a new, aforementioned, RF are popular in the microbiome context and holistic view to microbiome studies. With all this in mind we tend to outperform other methods (Zhou and Gallins, 2019; propose a generic and flexible kernel framework (Figure 1B) Namkung, 2020). ANN have shown excellent performance in as a way to handle unsupervised and supervised microbiome some cases but are susceptible to overfitting, especially if sample analyses (including the retrieval of microbial signatures), while size is greatly exceeded by the number of taxa, as is often the case paying special attention to data compositionality and spatial in metagenomics and metataxonomics. A desirable feature for and temporal integration. Kernel methods are a family within Frontiers in Microbiology | www.frontiersin.org 2 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis FIGURE 1 | Metagenomic analysis workflow. (A) Current state-of-the art: supervised and unsupervised learning are completely independent. (B) Kernel framework: the pivotal position of the kernel matrix is clearly observed. In gray, several tasks not performed during the present work but that merit future research. Frontiers in Microbiology | www.frontiersin.org 3 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis machine learning methods that share the use of kernel functions where xi and xj represent the taxonomic abundances in two or, simply, kernels. Some of these methods have been already different individuals, D is the number of different taxa, G(.) is applied to some specific problems or areas within microbiome the geometric mean, and γ > 0 is a hyperparameter that has analysis (Zhan et al., 2017; Mariette et al., 2018; Zhou and Gallins, to be tuned. This non-linear kernel derives from the Aitchison 2019) but their potential has not been fully exploited. In this distance, which is Euclidean and therefore (Eq. 1) is a valid work, we propose two new compositional kernels and discuss kernel. The logarithm term can be identified as the compositional how to translate non-compositional, but nonetheless widespread, clr-transformation (Gloor et al., 2017) over the original data. beta-diversity matrices to the kernel framework. We perform Analogously, we define the compositional linear kernel as: supervised and unsupervised analyses from the same kernel D matrix, and show how to extract microbial signatures. Spatial and xjk     X xik longitudinal data are also treated with specific kernel tools. This cLin (xi , xj ) = log log (2) G(xi ) G(xj ) kernel framework is illustrated with three case studies: a single k=1 point soil metagenomic dataset, a human dataset with a spatial Although cRBF is related to Aitchison distance and has the component, and a previously unpublished longitudinal dataset advantage of non-linearity, cLin is easier to interpret and allows concerning pig gut microbiota. An R package implementing the the retrieval of the microbial signatures. proposed methods, along with the analyzed datasets, is freely available at https://github.com/elies-ramon/kernInt. Non-compositional Kernels The most widely beta-diversity measures are Bray-Curtis, Unifrac and Jensen-Shannon (Gloor et al., 2017). Bray-Curtis and Jensen- Shannon are computed from taxonomic tables, while Unifrac MATERIALS AND METHODS additionally needs a phylogenetic tree. The Jensen-Shannon is metric and has a kernel counterpart that is already described in Kernels for Microbiome Data Bai and Hancock (2011) as the Jensen-Shannon Kernel (JSK): A real symmetric two-place function is a kernel iff, for every finite set of objects x1 , . . . , xN , it generates a positive semi-definite  ! 1 XD 2 xik JSK xi , xj = 1 −  matrix of dimension N × N: the kernel matrix (Schölkopf et al., xik ln 2 k=1 xik + xjk 2004; Shawe-Taylor and Cristianini, 2004). Probably the most ! widely known and used kernel functions are the linear and radial XD 2 xjk basis function (RBF) kernels, both defined for real vectors. + xjk ln (3) k=1 xik + xjk Intuitively, a kernel can be understood as a measure of the similarity between xi and xj . As objects x1 , . . . , xN are never provided that xi and xj contain relative frequencies. The Bray- represented explicitly, kernels can be designed for non-standard Curtis dissimilarity is semimetric, and so we propose using data types if a notion of what is considered “similar” in that Jaccard, a similar distance (Gardener, 2014), instead. The Jaccard given context exists (Schölkopf et al., 2004). Each kernel provides distance is paired with a well-known kernel (Bouchard et al., a different grasp of the dataset. Furthermore, as similarity 2013) and has a variant suitable for quantitative data. The measures, kernels are related (but opposite) to the beta-diversities quantitative Jaccard (also known as Ružička) kernel is defined in widely used in microbiome analyses. However, although every Gardener (2014) as: beta-diversity distance or dissimilarity is paired with a similarity XD min xik , xjk  counterpart, not all of them fulfill the aforementioned conditions and are, therefore, kernels. qJac(xi , xj ) = (4) k=1 max xik , xjk  We now present two compositional and two non- compositional kernels, all of them available in kernInt. In All aforementioned kernels have an asymptotic computational addition, users have the option of entering any kernel matrix, complexity of O(N 2 D). pre-computed with a kernel of their choice. In this work we are restricted only to kernels that can be obtained from Kernel Methods and Framework taxonomic abundance tables, but further insights can be found Kernel methods share the use of symmetric and positive semi- in the Discussion. definite matrices (i.e., kernel matrices), and not the original data, as input. That limits the potential similarity measures Compositional Kernels that one can use to only valid kernels, but also guarantees Here we define two kernels analogous to the linear and RBF that every matrix generated can be processed by the kernel kernels, but specific for compositional data. We introduce the method. Furthermore, using kernels places all different analyses Aitchison-RBF kernel as: in a common mathematical ground (see Figure 1B), which we refer as the kernel framework. For phenotype prediction, we use SVM, a classical method that can perform regression and cRBF(xi , xj )  classification (both binary and multi-class). For the unsupervised     x 2  analyses we use kernel principal components analysis (kPCA), xik = exp −γ D jk P k=1 log G(xi ) − log G(xj ) (1) a kernelized version of the standard algorithm. In both Frontiers in Microbiology | www.frontiersin.org 4 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis cases, kernInt allows the user to choose the values of the number of time points and 1t the time increment, then: hyperparameters and (in the case of SVM) to perform a T complete cross-validation and performance evaluation using an X fRBF (f , g) = exp(−γ (f (ti ) − g(ti ))2 ) (8) independent test set. i=1 Spatial Data T The kernel framework is particularly well suited for the X integration of spatial or heterogeneous data types (Schölkopf fLin (f , g) = 1t f (ti ) g (ti ) (9) i=1 et al., 2004; Mariette et al., 2018). This is because the integration can be done directly at the kernel matrices level. Let K1 , . . . , KM The discrete approach is sound in cases with few data points, be the kernel matrices computed from M different sources of when the modeling is less reliable. However, contrarily to (Eqs 6, data coming from the same individuals. Then, we can obtain a 7), these expressions cannot deal with irregular sampling times consensus kernel matrix K∗ : or missing data. M In multivariate scenarios, for instance microbiome data, many features are simultaneously sampled over time. Let fk and gk X K∗ = βz Kz (5) z=1 model taxon k in two individuals, being D the total number of taxa. The aforementioned kernels can be combined as in: with the restriction βz ≥ 0. The optimal βz values can be D obtained through an optimization process, which is known as Y fRBF 0 (f , g) = fRBF(fk , gk ) (10) multiple kernel learning (MKL) (Schölkopf et al., 2004). In k=1 unsupervised scenarios, a consensus matrix K∗ can be obtained by choosing the β coefficients that maximize average similarity of D X K∗ with all Kz matrices (Mariette et al., 2018). fLin0 (f , g) = fLin(fk , gk ) (11) k=1 Temporal Data With a computational complexity of O(N 2 TD) if A time series is an ordered set of repeated samples indexed (Eqs 8, 9) are used. by time, in the form {xi , ti }. The natural way to summarize It should be noted that the kernel approach allows the this type of data is through a function, which can be obtained integration of data that is both spatial and temporal-structured. using polynomial interpolation or splines. When data contains kernInt first handles the temporal dimension using a kernel for the time series of several individuals, it is commonly referred as longitudinal data (fLin or fRBF) over each space point, and then longitudinal data. integrates the spatial dimension by performing MKL over the The functional RBF kernel (Chen et al., 2013) translates the fLin or fRBF kernel matrices coming from the same individual. RBF kernel to accept real functions as input. Therefore, evolution over time among individuals is compared and used afterward for phenotype prediction or unsupervised tasks. Let f (t) and g(t) Microbial Signature be univariate functions, so that they represent the variation of In a broad sense, the “microbial signature” is the collection of a single feature in two different individuals between the time taxa associated with a trait of interest that has a high predictive interval [ta , tb ]. Then, the kernel definition is: value in the context of a given model (Rivera-Pinto et al., 2018). It can be retrieved from a linear SVM using the orientation of  Ztb  the separating hyperplane (Guyon et al., 2002): if the plane is 2 orthogonal to a particular feature dimension, then that feature fRBF f , g = exp −γ f (t) − g (t) dt   (6) is maximally informative. This method takes into account the ta correlation between taxa. As cLin is a translation of the linear In an analogous way, the functional linear kernel is defined as: kernel for compositional data, using (Eq. 2) we can retrieve the microbial signatures, which should be understood as the taxa Ztb importances after the clr-transformation. The same occurs when fLin (f , g) = f (t) g (t) dt (7) assessing the variable influence on the principal components ta in kPCA. A general permutation technique is proposed in Mariette et al. (2018), but using cLin permits obtaining the These kernels allow irregular sampling intervals and missing taxa influence in the same straightforward way than standard time points, but suffer of the cost of computing numerically principal components analysis (PCA). the integral (e.g., if an algebraic solution is not possible). The linearity also permits extending the microbial signature Computations can be simplified if fLin and fRBF are defined for retrieval, when using SVM, to the longitudinal and spatial cases. discrete functions, so the modeling of time series as continuous When performing MKL, as long as the cLin kernel is strictly functions is skipped. In this case, f (t) and g(t) may directly applied to all sampled sites, the global importance of a given denote the original objects {xi , ti }, so each time value directly taxon among all sites can be computed as the weighted sum maps to a certain value of the feature variable x. If T is the total (using the optimal β coefficients) of its partial importance in Frontiers in Microbiology | www.frontiersin.org 5 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis each site. In the longitudinal case, the global importance of each out at the ASV (3,577 ASVs were obtained) and at the Genera taxon k can be obtained from (Eq. 9) by addition of the partial taxonomic levels. importances over all T time points. Experimental Set-Up Case Studies and Data Pre-processing Analyses across the three datasets included a comparison with We illustrate our framework with three case studies: a single point the original reports (for Soil and Smokers datasets), as well dataset, a dataset with a spatial component, and a longitudinal as contrast with results from RF. The cLin and cRBF kernels dataset. The latter is previously unpublished while the rest of were applied directly to the raw counts, as they handle data the data is public. in an inherently compositional manner. Before computing both kernels, a number under the detection limit was added to Soil Dataset all dataset entries to handle zeroes (Quinn et al., 2018). An Bacterial composition of soil varies significantly at a alternative normalization of data, the cumulative sum scaling biogeographical scale, and is related to chemical and (Paulson et al., 2013) was performed prior to applying the environmental factors. Here we reanalyzed a single point non-compositional Jensen-Shannon and Jaccard kernels. That dataset by Lauber et al. (2009), who used 16S small-subunit way the compositional and non-compositional kernels could be ribosomal (16S rRNA) gene pyrosequencing to profile the compared. In the rest of cases (RF and longitudinal) we used the bacterial communities of different soils across North and compositional clr-transform over data. RF were obtained with South America. Authors reported that soil pH was significantly the R package randomForest (Liaw and Wiener, 2002), while the correlated with beta-diversity distances between samples. They kernel approach was carried out using kernInt (which relies on also found correlation with alpha diversity, which was highest the kernlab package for computing kPCA and SVM). A step-by- in soils with near-neutral pHs. To perform our analysis, we step guide with examples can be found at the kernInt package retrieved the taxonomic abundances as well as the associated vignette: https://elies-ramon.github.io/kernInt/. metadata from Qiita https://qiita.ucsd.edu/ (ID: 103). The Unsupervised analyses were carried out using the whole number of operational taxonomic unit (OTUs) was 7,396, datasets. Instead, for the supervised analyses, each dataset was while the number of soil samples was 89. As a part of the split at random into the training set (80% of data) and the test pre-processing, we excluded sample number 89, with only 1 read, set (20%). Optimal hyperparameters’ values (number of trees in which was also not included in the original paper. RF, cost in SVM, and γ for RBF-like kernels) and β coefficients for MKL were obtained by 5 × 5 cross-validation on the training Smokers Dataset set. Hyperparameters’ ranges are in Supplementary Table 1. Charlson et al. (2010) analyzed the impact of cigarette smoking Once the best values were found, the final model was built on the global airway microbial population. Bacterial communities using the whole training set. We repeated the whole process 40 were profiled using 454 pyrosequencing of the 16S rRNA gene times, each time with different 80/20 randomly split training/test in four airway sites: the left and right sides of nasopharynx and partitions, to obtain an error distribution. Performance over the oropharynx. Authors reported that composition was primarily test set was computed using normalized mean squared error determined by airway site, with individuals exhibiting minimal (NMSE) for regression and Accuracy for classification. We lateral or temporal variation. They used RF to predict the measured with the microbenchmark package the running time smoking status from the taxonomic abundances. We retrieved of computing the SVM models on a 64-bit Ubuntu 20.04 LTS the dataset (metadata and taxonomic abundances) from Qiita workstation with Intel(R) Core(TM) i5-6300U CPU at 2.40 GHz (ID: 524) to perform our analysis. Of the original 70 individuals, and 12 GiB of RAM (see Supplementary Figure 1). For the sake we discarded those that reported airway illness or antibiotic of comparison, the running time of several RF implementations usage in the 3 months prior to sampling. Thus, we analyzed the (including the randomForest package) can be found at same 62 individuals of the original work (29 smokers and 33 Wright and Ziegler (2017). non-smokers). Number of different OTUs was 2,817. For the Smokers and Pig case studies, additional conside- rations had to be taken into account. In the Smokers dataset, in Pig Dataset addition to the kPCA analysis, we computed the similarity among Here we present a previously unpublished dataset, which kernel matrices of different body sites with the mixKernel package evaluates the relationship of pre-weaning diarrhea with the (Mariette et al., 2018). We compared the performance of data early gut microbiota colonization in piglets. Gut microbiota was integration via MKL (the kernel approach) with that of RF when profiled in 153 piglets during their first week of life. Between using early and late integration approaches. In the former case, days 8 and 21 (weaning day), 79 out of the 153 piglets had the input of the RF was the concatenated data of the four sites. diarrhea and were treated with antibiotics. Swab sampling was Instead, in the latter case we used the forests created for each site done within 5 min after farrowing (day 0) and at days 3 and separately to vote for the final decision (Li et al., 2018). 7 post-farrowing. DNA was extracted from fecal samples and In the Pig dataset, to make sure that the training and test sets profiled using Illumina sequencing of 16S rRNA gene in each were completely independent, piglets from the same litter (full of the three time points. The cleaned sequences were processed sibs) were always placed either in one or other set. Performance of into amplicon sequence variants (ASVs). Further details are fLin and fRBF was contrasted to those of RF and their analogous described in Supplementary Method 0. Analyses were carried non-longitudinal kernels (cLin and cRBF) using all available days Frontiers in Microbiology | www.frontiersin.org 6 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis at once. For the non-longitudinal methods, 80% of the piglets basic side, orders Solibacterales and RB41 on the acid side). were used to train the model, using their three time points data in We used MiRKAT to test the significance of the association of separate rows, with time included as an additional variable. The the pH with both the top ten and 5% most important taxa, remaining piglets were reserved to test the model, but using only according to the cLin kernel. In both cases, we obtained very low one of their time points (either day 0, 3, or 7) chosen at random p-values (<10−8 ). and discarding the rest. This way, both longitudinal and non- longitudinal approaches had the same test set size. Longitudinal Smokers Data kernels fLin and fRBF were computed using (Eq. 9) and (Eq. 8), We predicted smoking status from the taxonomic abundances. as only three time points were available and we preferred not to At first models were built using the four sites separately, as in the interpolate the day’s in-between. Also, using the expression for original study. Authors used RF and reported a median accuracy discrete functions we could obtain the microbial signatures. The of 64% on the right and 65% on the left oropharynx (i.e., throat), information of all taxa was combined as in (Eqs 11, 10) and the and 71% on the right and 68% on the left nasopharynx. We training/test partitions were carried out as in the normal case. In re-computed the RF accuracies with our data pre-processing, a second step, the dataset was decomposed by sampling times and and obtained very similar results (Figure 3A), with the only the analysis was carried out for days 0, 3, and 7 separately using exception of the right nasopharynx (new median accuracy: 66%). RF, cLin and cRBF in the usual way. Regarding the kernels, the worst one was cLin (Supplementary Microbial signatures from SVMs were obtained from the Figure 5), which nonetheless gave similar accuracies to RF. The hyperplane normal vector w. The importance of taxon k is best kernel was the Jaccard kernel (Figure 3A), which improved computed by kernInt as (wk )2 (Guyon et al., 2002). When using substantially the RF accuracies, especially in the throat. Then, we RF, we used the mean decrease in node impurity (for regression combined the spatial-structured samples of the same individuals tasks) and mean decrease in Gini index (for classification). Both to test if accuracy increased when using an integrative approach RF and SVM give absolute values of taxa importance, so they were (Figure 3A). For the kernels, we first used MKL to combine converted to relative values. We used the R package MiRKAT the kernel matrices at the airway level (nasopharynx on one (Zhan et al., 2017) to test if the association of the target phenotype hand and oropharynx on the other) and, finally, we integrated with the signatures we obtained was statistically significant. all sites. This decreased the error substantially and delivered the best classification result, with a median accuracy of 92%. As for the RF, we tested both the early integration approach and RESULTS the late integration approach, and found that the latter granted better predictions. At best, integration of the four sites delivered Soil Data a median accuracy of 83%. The results for the rest of kernels can The cLin kPCA over the bacterial abundances is shown in be found in Supplementary Figure 5. In all cases, integration of Figure 2A. The remaining kPCAs, which gave a similar profile, the four datasets using our MKL proposal increased the accuracy can be found at Supplementary Figure 2. Soil samples are in comparison to the individual models, and doing so gave better clearly separated by their pH, in agreement with the original or equivalent results that those of RF integration approaches. results. The U-shaped projection is typical of data structured by a The only exceptions to this trend were the nasopharynx and gradual transition with few overlapping OTUs at the endpoints oropharynx models delivered by cLin (but not the model with the (Supplementary Figure 3). The peak diversity in near-neutral four sites combined). soils in contrast with extreme pHs may also have some effect Next, we recovered the overall microbial signature (i.e., across (Supplementary Figure 4). In addition, we used SVMs with the four sampling sites). The importance distribution is not as the four kernels described above to predict the pH of each skewed as in the Soil dataset: here the top 5% taxa accounted soil site from the bacterial abundances. This was not done in for the 62% of overall importance. The association of this subset the original work and so we used RF, a non-kernel, alternative of taxa with the target phenotype was highly significant (p- method, as benchmark. Results are shown in Figure 2B. The best value < 10−8 ). Top ten taxa are shown in Figure 3B. Neisseria sp. compositional kernel was cLin, having a median error of ∼0.09; large impact in discriminating smokers from non-smokers was and the best non-compositional one was JSK, with a median error already reported in the original work, especially in oropharynx of ∼0.10. In comparison, RF had a higher median error, almost models. The rest of highlighted taxa in Figure 3B were also noted the double of cLin, around 0.17. to have a role, either in models from nasopharynx alone or from To go further in the interpretation of the results, we both airways sites (Charlson et al., 2010). This mostly agrees analyzed the microbial signature retrieved from RF and cLin- with our results when the sampling sites are analyzed separately SVM. The distribution of the importances was highly skewed. (Supplementary Figure 6). For subsequent analyses we kept only 5% of the taxa, which Following the original work, differences in bacterial accounted for around the 90% (RF) and 95% (SVM) of total communities among the body sites were also analyzed. We importance, with the two methods having 42% of OTUs in present results for the Jaccard kernel in Figures 3C,D, while common. Top ten relevant taxa are shown in Figure 2C (RF) the rest are in Supplementary Figure 7. Figure 3C shows and Figure 2D (SVM). In agreement with the kPCA results, the similarity across kernel matrices derived from left and prediction is primarily driven by few OTUs of extreme pH right nasopharynx and oropharynx. The highest similarity ecosystems (e.g., genera Rubrobacter and Balneimonas on the was achieved within matrices of the same airway site but Frontiers in Microbiology | www.frontiersin.org 7 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis FIGURE 2 | (A) Compositional linear kPCA over the 88 soils. Color represents pH, while point size stands for the number of different observed taxa. (B) pH prediction error distribution over the 40 replicates. (C,D) Top relevant taxa for pH prediction according to RF and cLin-SVM. Standard error across the 40 replicates is marked with error bars. different laterality. As in the original paper (Supplementary worked best at the ASV level (with a median accuracy around Figure 8), using a kPCA (Figure 3D) we could discriminate 76%) than in Genera data (median accuracy ∼70%). Although between nasopharynx and oropharynx sites (first PC) but not aggregating taxa to the genus level is a relatively common between left and right. practice –see e.g., Rivera-Pinto et al. (2018)–, in our case using a coarser taxonomic resolution decreased the accuracy. Within the Pig Data non-longitudinal approach, we obtained similar accuracies using Evolution of gut microbiota from 153 healthy piglets over their RF and kernels, and both were close to the median accuracy of first week of life was used to predict the occurrence of pre- the random model (50.1%). To further understand the results, weaning diarrhea. In Figure 4A we compared the performance the analysis was carried out in days 0, 3, and 7 separately using of the longitudinal kernels (fLin and fRBF) vs. their analogous RF, cLin and cRBF kernels. Figure 4B reveals that all models non-longitudinal kernels (cLin and cRBF) plus RF when using from days 0 and 3 had no predictive power. Accuracy increased all available days at once. The longitudinal approach clearly dramatically after day 7 to a maximum of 73% for cRBF (ASV outperformed the non-longitudinal approach at both Genera level), only slightly worse than its analogous longitudinal kernel and ASVs levels. fRBF had a better performance than fLin, and fRBF. We used the kernel machine test of MiRKAT to further Frontiers in Microbiology | www.frontiersin.org 8 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis FIGURE 3 | (A) Non-smoker/smoker accuracies from taxonomic data: RF (green) vs. the quantitative Jaccard kernel (red). NoseL, NoseR, OroL, and OroR models are obtained from single datasets, the Nasopharynx and Oropharynx panels contain the information from the right and left sides, and All the combination of the four datasets. Early and Late refer to the early integration and late integration approaches for RF. (B) Top ten global cLin-SVM importances across all body sites. (C) Similarity across the Jaccard kernel matrices of the four sites (Nasopharynx Right and Left and Oropharynx Right and Left). (D) Jaccard kPCA of the taxonomic abundances. Color code represent airway site, whereas shape indicates the laterality of the samples. confirm that days 0 and 3 were not significantly associated for a 46% (fLin) and 58% (cLin) of the total importance, with an with phenotype, while day 7 was. As expected, only the kernel overlap between RF and cLin in day 7 of 2/3 of the ASVs. The matrices of day 7 delivered significant p-values (Genera: cLin association of the 5% most important taxa (global and day 7) with p-value < 10−6 , RBF p-value < 10−7 ; ASV: cLin and cRBF the phenotype was statistically significant according to MiRKAT p-values < 10−8 ) after Bonferroni correction. (p-values < 10−8 ). In a second step we analyzed the kPCA and microbial signatures, after discarding all models without predictive power. Figures 5A,B show the fLin and cLin (day 7) kPCA, while DISCUSSION fRBF and cRBF are in Supplementary Figure 9. In all cases a partial separation between healthy and sick piglets, with a large The kernel framework allows performing a great diversity of area of overlap, is observed. Genera relevance on prediction of analyses in a common ground, while allowing a great flexibility on pre-weaning diarrhea is shown in Figures 5C–E. We discuss how data is approached. However, within the microbiome field, the microbial signature at the Genera level, as around 2/3 of previous application of kernel methods has been mostly restricted the ASVs lack species assignation (Supplementary Figure 10). to specific areas. Zhan et al. (2017) proposed a kernel-based According to fLin, beneficial genera like Lactobacillus and semi-parametric regression method for testing the association of Bacteroides had the higher overall importance during the first the human microbiota communities with multiple phenotypes. week. In day 7, it was striking the great importance given by Their method was implemented in the R package MiRKAT. cLin to Desulfovibrio, and secondarily to Streptococcus. RF also In turn, Mariette et al. (2018) combined metagenomic data highlighted the butyrate-producing genus Dorea. Distribution of and environmental measures of the TARA ocean expedition the microbial signature at the ASV level was skewed, but again, using unsupervised MKL with the mixKernel package. In some much less than in the Soil case study. The top 5% ASVs accounted reports that compare the performance of different supervised Frontiers in Microbiology | www.frontiersin.org 9 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis FIGURE 4 | (A) Accuracy for RF, non-longitudinal kernels (cRBF, cLin) and Longitudinal kernels (fLin, fRBF) in the prediction of neonatal diarrhea from all available days. (B) Accuracy for RF and non-longitudinal cLin and cRBF kernels from metagenomic data of days 0, 3, and 7 post-birth separately. In both panels, the red dashed line marks the accuracy of the random model. methods in microbiome data, SVM often appear along RF or current state-of-the-art in the following points: (i) Proposal ANN (Qu et al., 2019; Zhou and Gallins, 2019; Namkung, 2020). and implementation of specific kernels for microbiota, while Thus, kernel methods were mostly used in an isolated way, QIIME2 currently provides default kernels for real vectors (the without exploiting the kernel framework ability to integrate a linear, RBF, polynomial and sigmoid kernels). (ii) As far as we great range of analyses while giving a unitary view. Another know, SVM is available in QIIME2 but kPCA is not; therefore, advantage of this framework is that it can handle virtually any it is not possible performing both supervised and unsupervised data type. However, to our best knowledge, it has not been analysis under the same mathematical point of view (Figure 1). previously applied to longitudinal microbiome studies. Finally, in (iii) Integration of spatial and temporal samples: QIIME2 previous works there was a lack of kernels that took into account does not have a specific handling of spatial (and, potentially, the compositional nature of metagenomic datasets. Here we multi-omic) data, while kernInt allows performing unsupervised, addressed all these questions, while also providing some examples supervised and retrieval of microbial signatures in this kind of of how previous kernel-based tools like MiRKAT and mixKernel datasets. On the other hand, the QIIME2 “longitudinal” plugin can fit into our framework. implements several tools for longitudinal data, but the option of When comparing kernInt to a popular package for performing supervised learning from the variation of microbiota microbiome analysis like QIIME2, it becomes apparent that over time is absent. the former is more specific in its scope. kernInt is not concerned Throughout this work, we summarized the microbiome with sequence alignment, taxonomic assignation and quality analyses in three branches: unsupervised learning (represented control as QIIME2 is, but with the analysis once the abundance by kPCA), supervised learning (SVM) and identification of table is obtained. Both packages are aimed at community ecology phenotype-associated microbial signatures. The Soil case study analysis (in QIIME2: alpha and beta diversities, PCoA, etc.) clearly illustrated how all three types are intertwined and and supervised learning areas. While kernInt does not have the complementary. In agreement with the original publication, both great range of methods available in QIIME2, it improves the SVM and kPCA results showed that taxonomic abundances and Frontiers in Microbiology | www.frontiersin.org 10 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis FIGURE 5 | Above: kPCA of fLin in panel (A) and of cLin (day 7) in panel (B). Below: Microbial signatures at the Genera level. Global importance for the first week is in panel (C). Importances for the day 7 according to the cLin kernel and RF are in panels (D,E). Standard error across the 40 replicates is marked with error bars. pH are strongly related. This granted a quite low prediction error the kernel methods tend to have an advantage over variable- (up to a median NMSE of 0.09) but, by itself, does not explain the oriented methods (e.g., in supervised learning: ridge regression, underlying mechanism connecting microbial abundance and pH. decision trees, RF, etc.) and over ANN (for the reasons stated in Microbial signature revealed that the SVM learning is driven by the “Introduction” section) when faced with N < < < D data. few taxa of opposite pH ecosystems. For instance, RB41 belong This is a common scenario in metagenomics when working at to the phylum Acidobacteria. The Rubrobacter genus contains the OTU or ASV level, but not necessarily in coarser taxonomic well known extremophiles and, like the Balneimonas (renamed resolutions. This is illustrated with the different behavior of Microvirga) genus, has preference for clearly alkaline soils (Dahal kernels with respect to RF in Figure 4B (see ASV vs. Genera and Kim, 2017; Chen et al., 2018). Furthermore, the arch in the results). In the other cases, SVM were consistently better (or at kPCA projection indicated that communities from acid and basic least equivalent) to RF in all the case studies that we analyzed. habitats did not overlap (Morton et al., 2017). Taken together, This disagrees with some previous reports in the microbiome these complementary views point that soil microbial structure is area, e.g., Zhou and Gallins (2019). However, it should be noted shaped by a gradual niche differentiation strongly modulated by that SVM performance depends on the kernel used, and these the pH. This agrees with previous findings on this dataset (Lauber reports used generic linear and RBF kernels. Even when using et al., 2009; Morton et al., 2017) but appears in a more concise and kernels specific for metagenomic data, we observed differences unified way using the kernel framework. among their mean NMSE or accuracies as large as fifteen In comparison to other methods, the kernel framework did percentage points. At the same time, our results suggest that not only allow a holistic view of data, but also gave good results in there is not a single kernel that systematically achieves the best each learning area. Concerning supervised learning, in general, performance in every problem. We found that cLin was the best Frontiers in Microbiology | www.frontiersin.org 11 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis one in the first case study, quantitative Jaccard in the second and application to time-structured samples. Pre-weaning diarrhea is fRBF in the third. In this scenario, we consider that the linear- an important issue in pig production, as the antibiotic treatment like kernels like cLin are a safe starting point. They allow for increases both the emergence of resistances and the economic the retrieval of the microbial signatures, are faster to compute costs. It is already known that gut colonization starts immediately and easier to interpret than non-linear kernels, and with high- after birth, and it evolves from a highly variable to a more stable dimensional data (>103 –104 taxa) they tend to match the RBF and homogeneous ecosystem over the first weeks. However, most kernel (usually considered the gold standard) in performance of the current studies in pig production ignore early dynamics (Hsu et al., 2003; Keerthi and Lin, 2003). RBF may be useful if in gut microbiota (Mach et al., 2015; Han et al., 2018; Massacci the number of different taxa is low, or when a strong non-linear et al., 2020). We wanted to test if pre-weaning diarrhea could relationship is suspected. The weakness of the compositional be anticipated as soon as the first week of life. In this sense, kernels that we proposed is that they cannot handle zeroes our results suggest that the first stages of intestinal microbiota without pre-processing; instead, zeroes pose no problems to the convey some valuable information indeed. kPCAs showed a quantitative Jaccard kernel. How to deal with zeroes is, currently, partial separation between piglets affected of diarrhea vs. healthy an open topic of research in compositional analysis (Weiss piglets, and by using longitudinal kernels we achieved a moderate et al., 2017; Quinn et al., 2018). If there is not enough a priori accuracy of 76%. However, it was unclear if this accuracy was to information that permits selecting the kernel beforehand, visual be attributed to a different taxa evolution in the two groups over assessment of the candidate kernel matrices via kPCA could be the first week, or to a single time point with a great predictive of some help. A more rigorous approach is to perform nested value. The day-by-day prediction clarified this issue, and showed cross-validation (Cawley and Talbot, 2010) to avoid overfitting that day 7 achieved a median accuracy of 73% while the rest when selecting both the candidate hyperparameters’ values and of points lacked predictive power. Even so, longitudinal kernels the best kernel for a given problem. Finally, phylogenetic kernels were able to slightly improve prediction (76% vs. 73% at the ASV were beyond the scope of this work, nor the available datasets had level, and 69% vs. 64% using Genera), so global taxa evolution the phylogenetic trees needed to compute them. However, they may also have a small role. may be derived from (Eqs 1, 2) by replacing the clr term with This is also seen in the underlying microbial signatures of other transformations, e.g., the PhILR transformation (Silverman the global first week (fLin) vs. day 7 (cLin). To be noted, et al., 2017). A phylogeny-based kernel was also proposed in Xiao in day 7 the most important genus was sulfate-reducing and Chen (2017). bacteria Desulfovibrio, which is known to have a relevant role Concerning our unsupervised analyses, we observed that the during pig gut colonization (Mach et al., 2015). Instead, the main structure revealed by the original MDS/PCoA (ordination global (longitudinal) model was mainly led by Lactobacillus by pH in the Soil dataset, and by body site in the Smoker dataset) and Bacteroides. Relationship of both genera to pre-weaning was conserved in our kPCAs. On the other hand, microbial diarrhea is well sustained in literature. Lactobacillus spp. are well signatures obtained with SVM had a biological interpretation. In known probiotic bacteria, while members of Bacteroides genus general, the most important taxa retrieved from SVM coincided are associated with increased infants gut microbial diversity with those of RF (40–65% of overlap depending on the dataset), (Stewart et al., 2018). Furthermore, both play an important role and could be recovered too when dealing with spatial and on mammals’ gut microbial colonization (Sawicki et al., 2017; temporal-structured datasets. However, we acknowledge that a Wexler and Goodman, 2017) and are dominant in healthy pigs drawback of these signatures (though they handle well the cases compared with diarrhea-affected piglets (Song et al., 2017), which of multicollinearity) is that they are based on linear kernels. In gives confidence in the reliability of our findings. turn, RF can take into account both non-linearity and complex In summary, our kernel framework successfully places the interactions among taxa. In any case, the informativeness most important analyses in the microbiome field on a common of a microbial signature can be assessed by the prediction ground, takes into account the compositionality of data, and performance of the SVM model that generated it. is flexible enough to integrate spatial and temporal dimensions Apart of the aforementioned advantages of the kernel of the datasets. framework, we also showed how it can accommodate datasets with spatial and/or temporal components. We illustrated the integration of spatial-structured samples with the Smokers DATA AVAILABILITY STATEMENT dataset. The analysis in the original work was carried out in each The data analyzed in this study is subject to the following sampling site independently, with a maximum median accuracy licenses/restrictions: This manuscript utilizes proprietary data. of 71%. Here we showed how combining the body sites using Requests to access these datasets should be directed to YR-C/ MKL increased the median accuracy to 92%. Therefore, our IRTA/yuliaxis.ramayo@irta.cat. results remark the relevance of using an integrative approach to improve the accuracy of phenotype prediction when spatial- structured samples of the same individuals are available. ETHICS STATEMENT In addition to the package and framework proposal, we analyzed a previously unpublished dataset profiling the The animal study was reviewed and approved by the Central microbiota evolution and pre-weaning diarrhea incidence in 153 Authority for Scientific Procedures on Animals of Netherlands– piglets. Through this dataset we illustrated the kernel framework Centrale Commissie Dierproeven (CCD). Frontiers in Microbiology | www.frontiersin.org 12 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis AUTHOR CONTRIBUTIONS from the Spanish Ministry of Science and Innovation. We acknowledge further financial support from the Spanish Ministry YR-C, MP-E, RQ, and ER contributed to conception and design of Economy and Competitiveness through the “Severo Ochoa of the study. FM was in charge of the pig data sampling. YR-C and Programme for Centres of Excellence in R&D” 2016–2019 (SEV- MP-E supervised the overall research, while LB-M supervised the 2015-0533),” and from the EU through the BFU2016–77236-P machine learning part. ER performed the all analysis and wrote (MINECO/AEI/FEDER, EU). the first draft of the manuscript. YR-C, MP-E, and LB-M revised and wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version. ACKNOWLEDGMENTS FUNDING The authors warmly thank all technical staff from Schothorst Feed Research. This work was funded by projects PID2019-108829RB-I00, AGL2016–78709-R, and AGL2017–88849-R awarded by the Spanish Ministry of Economy and Competitiveness. ER has SUPPLEMENTARY MATERIAL funding from a FI-AGAUR Ph.D. studentship grant, with the support of the Secretaria d’Universitats i Recerca de la Generalitat The Supplementary Material for this article can be found de Catalunya and the European Social Fund. YR-C is recipient of online at: https://www.frontiersin.org/articles/10.3389/fmicb. a Ramon y Cajal post-doctoral fellowship (RYC2019-027244-I) 2021.609048/full#supplementary-material REFERENCES Gardener, M. (2014). Community ecology: analytical methods using R and Excel. United Kingdom: Pelagic Publishing Ltd. Argelaguet, R., Velten, B., Arnol, D., Dietrich, S., Zenz, T., Marioni, J. C., et al. Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V., and Egozcue, J. J. (2017). (2018). Multi-Omics Factor Analysis–a framework for unsupervised integration Microbiome Datasets Are Compositional: And This Is Not Optional. Front. of multi-omics data sets. Mole. Syst. Biol. 14:e8124. Microbiol. 8:2224. doi: 10.3389/fmicb.2017.02224 Bai, L., and Hancock, E. R. (2011). Graph clustering using the jensen-shannon kernel. Guyon, I., Weston, J., Barnhill, S., and Vapnik, V. (2002). Gene selection for cancer In International Conference on Computer Analysis of Images and Patterns, classification using support vector machines. Mach. Learn. 46, 389–423. Berlin: Springer, 2011, 394–401. Han, G. G., Lee, J. Y., Jin, G. D., Park, J., Choi, Y. H., Kang, S. K., et al. (2018). Berg, G., Rybakova, D., Fischer, D., Cernava, T., Vergčs, M. C. C., Charles, T., Tracing of the fecal microbiota of commercial pigs at five growth stages from et al. (2020). Microbiome definition re-visited: old concepts and new challenges. birth to shipment. Scient. Rep. 8, 1–9. Microbiome 8, 1–23. Hsu, C. W., Chang, C. C., and Lin, C. J. (2003). A practical guide to support vector Bodein, A., Chapleur, O., Droit, A., and Lê Cao, K. A. (2019). A generic multivariate classification. Available online at: https://www.csie.ntu.edu.tw/~{}cjlin/papers/ framework for the integration of microbiome longitudinal studies with other guide/guide.pdf [accessed on June 12, 2020] data types. Front. Genet. 10:963. doi: 10.3389/fgene.2019.00963 Ibrahim, O. M. (2013). A comparison of methods for assessing the relative Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, importance of input variables in artificial neural networks. J. Appl. Sci. Res. 9, et al. (2019). Reproducible, interactive, scalable and extensible microbiome data 5692–5700. science using QIIME 2. Nat. Biotechnol. 37, 852–857. Keerthi, S. S., and Lin, C. J. (2003). Asymptotic behaviors of support vector Bouchard, M., Jousselme, A. L., and Doré, P. E. (2013). A proof for the positive machines with Gaussian kernel. Neural. Comput. 15, 1667–1689. doi: 10.1162/ definiteness of the Jaccard index matrix. Int. J. Approx. Reas. 54, 615–626. 089976603321891855 doi: 10.1016/j.ijar.2013.01.006 Lauber, C. L., Hamady, M., Knight, R., and Fierer, N. (2009). Pyrosequencing-based Cawley, G. C., and Talbot, N. L. (2010). On over-fitting in model selection and assessment of soil pH as a predictor of soil bacterial community structure at the subsequent selection bias in performance evaluation. J. Mach. Learn. Res. 11, continental scale. Appl. Environ. Microbiol. 75, 5111–5120. doi: 10.1128/AEM. 2079–2107. 00335-09 Charlson, E. S., Chen, J., Custers-Allen, R., Bittinger, K., Li, H., Sinha, R., et al. Li, Y., Wu, F. X., and Ngom, A. (2018). A review on machine learning principles (2010). Disordered microbial communities in the upper respiratory tract of for multi-view biological data integration. Brief. Bioinform. 19, 325–340. cigarette smokers. PLoS one 5:e15216. doi: 10.1371/journal.pone.0015216 Liaw, A., and Wiener, M. (2002). randomForest: Breiman and Cutler’s random Chen, R. W., Wang, K. X., Wang, F. Z., He, Y. Q., Long, L. J., and Tian, X. P. (2018). forests for classification and regression. R Package Vers. 4, 6–10. Rubrobacter indicoceani sp. nov., a new marine actinobacterium isolated from Mach, N., Berri, M., Estellé, J., Levenez, F., Lemonnier, G., Denis, C., et al. Indian Ocean sediment. Int. J. Systemat. Evolut. Microbiol. 68, 3487–3493. doi: (2015). Early-life establishment of the swine gut microbiome and impact on 10.1099/ijsem.0.003018 host phenotypes. Environ. Microbiol. Rep. 7, 554–569. doi: 10.1111/1758-2229. Chen, H., Tang, F., Tino, P., and Yao, X. (2013). “Model-based kernel for efficient 12285 time series analysis,” in Proceedings of the 19th ACM SIGKDD international Mariette, J., Villa-Vialaneix, N. (2018). Unsupervised multiple kernel learning for conference on Knowledge discovery and data mining, eds J. He, Y. Koren, R. heterogeneous data integration. Bioinformatics 34, 1009–1015. doi: 10.1093/ Ghani, R. Parekh, I. Dhillon, and T. E. Senator (United States: Association for bioinformatics/btx682 Computing Machinery), 392–400. Massacci, F. R., Berri, M., Lemonnier, G., Guettier, E., Blanc, F., Jardet, D., et al. Coenen, A. R., Hu, S. K., Luo, E., Muratore, D., and Weitz, J. S. (2020). A Primer (2020). Late weaning is associated with increased microbial diversity and for Microbiome Time-Series Analysis. Front. Genet. 11:310. doi: 10.3389/fgene. Faecalibacterium prausnitzii abundance in the fecal microbiota of piglets. Anim. 2020.00310 Microb. 2, 1–13. Dahal, R. H., and Kim, J. (2017). Microvirga soli sp. nov., an alphaproteobacterium McMurdie, P. J., and Holmes, S. (2014). Waste not, want not: why rarefying isolated from soil. Int. J. Syst. Evolut. Microbiol. 67, 127–133. doi: 10.1099/ijsem. microbiome data is inadmissible. PLoS Comput. Biol. 10:e1003531. doi: 10.1371/ 0.001582 journal.pcbi.1003531 Frontiers in Microbiology | www.frontiersin.org 13 January 2021 | Volume 12 | Article 609048 Ramon et al. A Kernel Framework for Microbiome Analysis Morton, J. T., Toran, L., Edlund, A., Metcalf, J. L., Lauber, C., and Knight, R. (2017). childhood from the TEDDY study. Nature 562, 583–588. doi: 10.1038/s41586- Uncovering the horseshoe effect in microbial analyses. Msystems 2, 166–e116. 018-0617-x Namkung, J. (2020). Machine learning methods for microbiome studies. Su, X., Jing, G., Sun, Z., Liu, L., Xu, Z., McDonald, D., et al. (2020). Multiple-Disease J. Microbiol. 58, 206–216. doi: 10.1007/s12275-020-0066-8 Detection and Classification across Cohorts via Microbiome Search. Msystems Paulson, J. N., Stine, O. C., Bravo, H. C., and Pop, M. (2013). Differential abundance 5, 150–e120. analysis for microbial marker-gene surveys. Nat. Methods 10, 1200–1203. doi: Weiss, S., Xu, Z. Z., Peddada, S., Amir, A., Bittinger, K., Gonzalez, A., et al. (2017). 10.1038/nmeth.2658 Normalization and microbial differential abundance strategies depend upon Qu, K., Guo, F., Liu, X., Lin, Y., and Zou, Q. (2019). Application of machine data characteristics. Microbiome 5:27. learning in microbiology. Front. Microbiol. 10:827. doi: 10.3389/fmicb.2019. Wexler, A. G., and Goodman, A. L. (2017). An insider’s perspective: Bacteroides as 00827 a window into the microbiome. Nat. Microbiol. 2, 1–11. Quinn, T. P., Erb, I., Richardson, M. F., and Crowley, T. M. (2018). Understanding Wright, M. N., and Ziegler, A. (2017). ranger: A Fast Implementation of Random sequencing data as compositions: an outlook and review. Bioinformatics 34, Forests for High Dimensional Data in C++ and R. J. Statist. Softw. 77, 1–17. 2870–2878. doi: 10.1093/bioinformatics/bty175 Xiao, J., and Chen, J. (2017). Phylogeny-based kernels with application to Rivera-Pinto, J., Egozcue, J. J., Pawlowsky-Glahn, V., Paredes, R., Noguera-Julian, microbiome association studies. In New Advances in Statistics and Data Science. M., and Calle, M. L. (2018). Balances: a New Perspective for Microbiome Cham: Springer, 217–237. Analysis. mSystems 3, 53–e18. Zhan, X., Tong, X., Zhao, N., Maity, A., Wu, M. C., and Chen, J. (2017). A small- Sawicki, C. M., Livingston, K. A., Obin, M., Roberts, S. B., Chung, M., and sample multivariate kernel machine test for microbiome association studies. McKeown, N. M. (2017). Dietary fiber and the human gut microbiota: Genet. Epidemiol. 41, 210–220. doi: 10.1002/gepi.22030 application of evidence mapping methodology. Nutrients 9:125. doi: 10.3390/ Zhou, Y. H., and Gallins, P. (2019). A review and tutorial of machine learning nu9020125 methods for microbiome host trait prediction. Front. Genet 10:579. doi: 10. Schölkopf, B., Tsuda, K., and Vert, J. P. (2004). Kernel methods in computational 3389/fgene.2019.00579 biology. New York: MIT press. Zingaretti, L. M., Renand, G., Morgavi, D. P., and Ramayo-Caldas, Y. (2020). Link- Shawe-Taylor, J., and Cristianini, N. (2004). Kernel methods for pattern analysis. HD: a versatile framework to explore and integrate heterogeneous microbial Cambridge: Cambridge university press. communities. Bioinformatics 36, 2298–2299. doi: 10.1093/bioinformatics/ Shenhav, L., Thompson, M., Joseph, T. A., Briscoe, L., Furman, O., btz862 Bogumil, D., et al. (2019). FEAST: fast expectation-maximization for microbial source tracking. Nat. Methods 16:627. doi: 10.1038/s41592-019- Conflict of Interest: The authors declare that the research was conducted in the 0431-x absence of any commercial or financial relationships that could be construed as a Silverman, J. D., Washburne, A. D., Mukherjee, S., and David, L. A. (2017). A potential conflict of interest. phylogenetic transform enhances analysis of compositional microbiota data. Elife 6:e21887. Copyright © 2021 Ramon, Belanche-Muñoz, Molist, Quintanilla, Perez-Enciso Song, D., Peng, Q., Chen, Y., Zhou, X., Zhang, F., Li, A., et al. (2017). and Ramayo-Caldas. This is an open-access article distributed under the terms Altered gut microbiota profiles in sows and neonatal piglets associated of the Creative Commons Attribution License (CC BY). The use, distribution or with porcine epidemic diarrhea virus infection. Scient. Rep. 7, reproduction in other forums is permitted, provided the original author(s) and the 1–10. copyright owner(s) are credited and that the original publication in this journal Stewart, C. J., Ajami, N. J., O’Brien, J. L., Hutchinson, D. S., Smith, D. P., Wong, is cited, in accordance with accepted academic practice. No use, distribution or M. C., et al. (2018). Temporal development of the gut microbiome in early reproduction is permitted which does not comply with these terms. Frontiers in Microbiology | www.frontiersin.org 14 January 2021 | Volume 12 | Article 609048