The first and foremost step in any time series analysis pertains to visualization of temporal trends of the constituent entities (for example taxonomic groups in a microbial ecosystem). This workflow can be used to visualize and identify high, medium, and low abundant taxa. In addition, it allows identification of 'core,' 'persistent,' and 'transient' microbial groups which serve as important characteristic constituents of the ecosystem. The core microbiome refers to those taxa which are present across all timepoints. On the other hand, the persistent microbiota refers to the ones that are present across extended time points, but not in all. In contrast, the transient group comprises of those sets of taxa which show frequent trends of appearance and disappearance. It should be noted that although the threshold parameters used for defining the 'core,' 'transient,' and 'persistent' have been taken from a previously reported study (Caporaso et al., 2011), they are prone to biases due to sequencing depth.

### Workflow-2

### **Compare temporal trends between selected taxa**

Analysis of time series data often requires trend comparison of a custom set of taxonomic groups. For example, the group may be a set of taxa previously known to show a characteristic behavior. The current workflow allows easy graphical comparison of two

or more user selected taxa over the sampled timeline. TIME also allows comparison of trends in microbial abundance using a simple 'select and plot' operation, wherein users can choose microbial taxa (at a specified taxonomic level) using a simple auto-complete search or from a dropdown selection. One or more taxa can then be appended to or removed from the existing plot, thereby providing an easy way to study a selected set of microbes. Most of the microbiome datasets consist of sets of highly abundant as well as rare taxa. This poses a major problem while plotting and visualizing multiple taxa together in a single plot, with highly abundant taxa dominating the scale, thereby making it difficult to decipher patterns for the rare groups. In case the selected taxa have different abundance scales, users can utilize the 'log scaling' option for comparing their trends. Particular time stretches of interest can also be zoomed for in depth analysis.

### Workflow-3

#### **Identify temporally stable/unstable taxa**

One of the important steps in time series analysis pertains to identification of stationary entities corresponding to the ones which have mean, standard deviation and variance constant over time. In microbial time series studies, identification of stationary taxa is especially crucial to detect inter microbial competition (David et al., 2014). As demonstrated in an earlier study (David et al., 2014), the presence of competition among the resident taxa is expected to cause sustained growth of some of them leading to their non-stationary behavior. Since in most cases the microbiota are in stable state, only a few taxonomic groups are expected to be non-stationary. A significant test of non-stationarity hence can be considered as a hint for a restoring force governing bacterial dynamics. However, fluctuations due to diet and environment may also affect stationarity of taxa and hence a cautious interpretation of results is required. Additionally, the similarities in phylogeny of non-stationary taxa may also provide clues pertaining to resource competition as genetically similar taxa are more likely to exhibit resource competition.

### Workflow-4

### **Identify variations in taxonomic groups between two time ranges**

Time series experiments involving perturbation events (like administration of antibiotics) are likely to disrupt the microbial community structure. In such analysis, it might be of interest to identify and visualize the exact temporal effect of the perturbation on the resident microbial groups. This workflow allows identification of taxa which undergo noticeable changes between two selected time ranges along with statistical inferences. Taxonomic behaviors like gradual increase or decrease in abundances can be easily inferred from the tabular summary generated using this workflow.

### Workflow-5a

#### **Cluster groups of taxa having similar behavior over time**

An important goal while analyzing microbial time series data pertains to identification of groups of taxa which show similar trends over a time stretch. Similar temporal behavior by different bacterial taxa could arise due to reasons like symbiotic relationship between two or more bacteria. On the other hand, it is also important to know which bacterial taxa behave in temporally opposite ways, since such behavior might be an indicator of some underlying interaction or competition among them. Taxonomic groups depicting similar behavior in a selected timeframe are identified using Dynamic Time Warping (DTW) algorithm (described under "Materials and Methods" section). The output can be visually explored using interactive tree and trend plots. Each branch of the tree corresponds to a set of taxa having similar time series trends. Users can select a branch (a group of taxa having similar temporal patterns) or an individual terminal node (taxa) from the tree and visually explore the time series trends using the assistive plot.

### Workflow-5b

#### **Explore pair-wise relationship among taxonomic groups**

Visualization of correlation and other similarity indices between the resident taxonomic groups often helps to gather meaningful insights. This workflow allows users to select Pearson correlation or modified DTW (referred to as TIME-DTW) index and use it to generate heatmaps. Such heatmaps are useful for visual pattern mining and the corresponding distance metric can be exported for further advanced network analysis.

### Workflow-6

#### **Explore inter taxa interactions using causality network**

The existence of a strong correlation in the abundance of two or more taxonomic groups across a time scale may not always be ascertained to causation. A recent study has utilized 'module networks' to understand causality relationships among bacteria (Lu et al., 2017). A causation event can be ascertained between two taxonomic groups when the past values of one taxon are observed to have some information about the future values of the other. This analysis is performed in 'TIME' using a Granger causality algorithm (described in details under "Materials and Methods" section). The global community behavior over the whole sampled timeline is captured using interactive causality networks and trend plots. Each node in the network can be queried for its causality using interactive operations. While right clicking on a node (corresponding to a taxon) highlights the nodes (or taxa) that are affected ('Granger caused') by it, left clicking on the same highlights the nodes (or taxa) responsible for affecting ('Granger causing') its temporal changes.

### Workflow-7

### **Cluster time points based on similar community patterns**

Many microbial time series datasets are often observed to have a typical composition of constituent entities which gives rise to seasonality or periodicity of microbial communities. These similarities and differences in the proportion of the constituent taxonomic groups give rise to 'community states' in the microbiome. Such 'community states' could be useful for obtaining insights into the microbial dynamics (Gajer et al., 2012). The interactive hybrid trend plot and heatmap generated using this workflow is useful for visualizing the temporal changes in the community structure.

### Case Studies on Publicly Available Time Series Microbiome Data

We demonstrate the applicability and utility of each workflow using three publicly available time series microbiome datasets. 'Caporaso-Dataset' corresponds to a longitudinal metagenomic time series data of gut microbiome samples from an American healthy male and female subject, collected at regular intervals spanning a long time period (Caporaso et al., 2011). A second time series metagenomic dataset ('Dethlefsen-Dataset') corresponds to a study evaluating the effects two doses of antibiotic treatments on the gut microbiome of three adult American females (Dethlefsen and Relman, 2011). The third dataset ('Gajer-Dataset') corresponds to a temporal sampling of vaginal microbiome of 32 reproductive age women over a period of 16-weeks (Gajer et al., 2012). All the above datasets are pre-loaded into the 'TIME' application for users' convenience.

### Case Study 1: Analysis of Microbial Perturbation from Microbiome Time Series Data

In order to demonstrate the applicability of 'TIME' in analysis of perturbation, 'Dethlefsen-Dataset' (antibiotic treatment) was selected and various relevant workflows were used for analysis. The 'Dethlefsen-Dataset' had an associated metadata mapping for the time points corresponding to the different states (for all three individuals – D, E, and F), namely before antibiotic treatment ('PreCp'), during the two doses ('FirstCp' and 'SecondCp'), the week immediate post the two treatments ('FirstWPC' and 'SecondWPC'), gap between the doses ('Interim') and the time points post treatment ('PostCp'). A drastic drop in diversity and richness specifically at the points of perturbation ('FirstWPC' and 'SecondWPC'), could be visually inferred using the diversity plots generated using TIME (**Figure 1A**). The Figure also shows the slow but incomplete recovery in diversity post perturbation, which is in line with the reported findings (Dethlefsen and Relman, 2011). The core taxa identified using 'Workflow-1' (at 'genus' level) also indicates an interindividual variation among the three subjects (**Figure 1B**), with only four genera to be consistently common across all (Bacteroides, Coprococcus, Roseburia, and Dorea). In order to identify the taxa which are most affected by the antibiotic treatment on 'Sample E' (as a representative example), the Workflow-4 was employed after selecting two time stretches, namely, 'before Cp1' (Period 1 ranging from time point 0–59) and 'after Cp1' (Period 2 ranging from time point 65–124). This analysis identified the affected genera sorted by the log fold change in the mean abundances between period 1 and period 2. Haemophilus, Butyrivibrio, Eubacterium, Turicibacter, and Parabacteroides were identified to be the top five affected genera upon antibiotic treatment based of log-fold abundance (**Figure 2**) but none of them were found to be statistically significant (when evaluated with Wilcoxon Rank-Sum Test using P-values corrected for multiple testing). Subsequently, to gather a deeper insight into the pattern of the affected genera during perturbation or genera similarly affected during perturbation, Workflow-5a was used (selecting sample 'E,' time point as 'FirstCp' and a rare taxa cutoff of 0.5) to generate the DTW tree (**Figure 3A**). Visual inference of the tree revealed three clear clusters (**Figure 3A**), each of which were used to generate their corresponding trend plots (**Figure 3B**). While Cluster 2 seemed to contain genera whose abundance is most strongly decreased by antibiotic treatment, Cluster 1 contained the moderately affected ones. On the other hand, Cluster 3 consisted of genera which increased post perturbation, possibly due to the reduced abundances of taxa belonging to Clusters 1 and 2.

### Case Study 2: Insights into Microbial Inter-Dependencies Using Causality Networks

To analyze the effect of stationary genera and its relation to causality, the female subject (at genus level) from 'Caporaso-Dataset' (the 6 months spanning time series sampling) was used to generate the causality network using Workflow-6 (keeping a rare taxa cutoff of 0.5). The non-stationary genera information was overlaid on the network using one of the features in TIME which highlights the corresponding names (**Figure 4**). While a majority of the genera were seen to be stationary, a few exhibited non-stationary behaviors, an observation similar to an earlier study on a different gut microbiome dataset (David et al., 2014). Further, the majority of the non-stationary genera belonged to the phylum Firmicutes, strengthening the hypothesis of phylum level (genetically similar) resource competition (David et al., 2014). However, owing to the complexity of the community interactions in a gut microbiome, further experimental validations are required to support this hypothesis. In order to infer the effect of a non-stationary genus on others, we chose two non-stationary genera nodes (Faecalibacterium and Clostridium) from the causality network. While Faecalibacterium is a well documented commensal gut bacterium, a number of species belonging to Clostridium are known to have several pathogenic effects on human. Right clicking on these nodes enables one to highlight the edges connecting the genera affected ('Granger caused') by them and correspondingly displays the trend plot of all the associated taxonomic groups. A quick look into the edge connections showed that most of the genera affected by Faecalibacterium are non-stationary as compared to the ones affected by Clostridium. This observation suggests a differential influence of one taxon over others.

### Case Study 3: Importance of Time Series Community Analysis

Microbial communities in different body sites have been reported to exhibit differences in their compositions. These compositions are also known to change over time. For example, studies on temporal variation of human gut microbiome have reported the presence of periodic as well as non-periodic diversity patterns (Caporaso et al., 2011). Such similar temporal patterns arise due to a comparable microbial community composition across these time points. Workflow-7 of 'TIME' is dedicated to identify such community clusters and visualize their variations across the timeline. 'Caporaso-Dataset' and 'Gajer-Dataset' (corresponding to gut and vaginal time series microbiome, respectively) were used as a part of this analysis pipeline (results are summarized in **Figures 5**, **6**, respectively). In order to consider the effect of only the 'non-rare taxa' (taxa which occur in at least 70% of

samples), a rare taxa cutoff of 0.7 was selected and a bi-directional clustering was done (for time points and taxa). Each of the two 'time-point clusters' (**Figure 5**) represent a group of time points having similar microbial distributions, called 'community states' (see section "Materials and Methods" for details). A comparison of the female and male gut microbiome time series ('Caporaso-Dataset') using the above workflow revealed a clear bias of one of the two 'community states' in the male (**Figure 5B**) while almost an equal distribution of the two 'community states' was found in the female (**Figure 5A**). In male, while the dominant cluster had mainly the genera Bacteroides and Parabacteroides as distinguishable marker, other genera namely Prevotella and Campylobacter were observed to be the prominent contributors of the less dominant cluster (**Figure 5B**). On the other hand, the female microbiome had one cluster prominently dominated by Akkermansia, with no single clearly dominant member in the other (**Figure 5A**). To explore community states in a different body site, vaginal microbiome from subject-1 of 'Gajer-Dataset' was considered and analyzed using Workflow-7. A clear periodic pattern in the 'community states' was observed (**Figure 6B**), probably due to the prominent changes in the menstrual cycle and related hormonal changes in reproductive age females. While one 'community state' showed a dominance of Lactobacillus iners, the other showed a dominance of Atopobium. The genera Atopobium is known to be associated with bacterial vaginosis, while lactic acid producing bacteria (like L. iners) are known to prevent pathogen colonization by creating an acidic environment (Gajer et al., 2012). The generated heatmap (**Figure 6A**) as well as trend comparison plot using Workflow-2 (**Figure 6C**) indicate an antagonistic behavior between the above two taxa.

### DISCUSSION

The various workflows in 'TIME' allow visualizing time series microbiome data as well as analyzing them to obtain meaningful biological insights. It is to be noted that a few key points need to be considered before interpreting the generated outputs and building hypotheses based on such datasets (Weiss et al., 2017). For instance, the microbial abundance files used as input for analysis represent the count of clustered sequences (OTUs) across several time points corresponding to one or more

FIGURE 2 | Demonstrating the utility of 'Workflow-4' in identifying the taxonomic groups completely eliminated and the ones mostly affected by antibiotic treatment during the first dosage period ('FirstCp' of 'Dethlefsen-Dataset' used in the case study). Two time periods ('BeforeCp' and 'AfterCp') were chosen based on the metadata.

sources (samples). Regardless of strict experimental designs, not all sources as well as time points are sampled/sequenced at similar depths due to sampling constraints as well as sequencing limitations. Hence, samples sequenced at lower depth may display biased diversity estimates and consequently affect the downstream analyses. For example, workflow 4 in TIME can predict differential abundant taxa between two time stretches with increased confidence if the sequencing depths are sufficiently high and even since samples with higher number of sequence will have better estimates of abundances. Similarly, if some time points are sampled deeper than the others, it makes interpretation of transient and rare taxa difficult (workflows 1

shown in the legend while bold labels indicate non-stationary taxa.

and 3) without normalization. In addition, presence of sparse OTUs represents uncertainty in counts owing to limitations in the sequencing detection ability (since they are below the detection threshold). A majority of microbiome studies consider either a relative normalization route (OTU counts scaled to proportions) or a rarefaction based normalization step (each sample is sub-sampled to an even depth), both of which are implemented in TIME for convenience. Use of rarefaction curves can provide guidance on choosing a suitable rarefaction depth for normalization and lower the false discovery rates (Weiss et al., 2017). However, it should be kept in mind that rarefying a data might impact a number of downstream analysis workflows due to removal of a subset of the data. Moreover, time series data involving perturbation events, if normalized using rarefaction, might subdue the effect of the perturbation itself. Relative normalization on the other hand, is also prone to create several artifacts (Stämmler et al., 2016). Both rarefied as well as relatively normalized data are compositional, therefore, fluctuations in abundance of one taxon might lead to spurious fluctuations in abundance of other taxa resulting in false correlations (Weiss et al., 2017). A lack of knowledge of absolute abundance can thus impact the interpretation of the results of the analyses. For example in workflow 3, although a taxon might change in abundance and appear to be non-stationary, it may actually be not changing but taxa around it may be changing in relative abundance. Moreover, relative abundance based approaches ignore the possibility that the altered abundance itself could be a key identifier of a disease state (Vandeputte et al., 2017). It may also be noted that both relative and absolute abundances are required for obtaining a comprehensive understanding of time series microbiome data (Props et al., 2017). Additionally, data obtained from appropriately designed experiments (e.g., using replicates for each time point) will increase confidence on the obtained results. Advanced experimental protocols have also been reported (Stämmler et al., 2016) which helps in normalizing the biases arising due to differential microbial loads across samples.

The incorporated Granger causality based interaction networks in 'TIME' provides a way to capture the overall global microbial community behavior and is ideal for datasets having evenly sampled time-points. Variations of Granger causality have been applied earlier to decipher ecological relationships

FIGURE 4 | (A,B) Represent two composite plots for Granger causality graphs (A1 and B1) and Trend plots (A2 and B2) corresponding to genera Faecalibacterium and Clostridium respectively. Granger causality graph (A1 and B1) for the constituent taxa in the female subject of 'Caporaso-Dataset' used in the case study generated using 'Workflow-6.' The trend plots (A2 and B2) for two genera namely Faecalibacterium and Clostridium along with the genera caused (or affected) by them are displayed below the corresponding circular graphs. The arrows in the graph represent the causality relationships between the source and target nodes. The genera color labels correspond to their respective phyla as shown in the legend while bold labels indicate non-stationary taxa.

(Detto et al., 2012) and in gene expression networks (Yang et al., 2017) with reasonable success. However, not all Granger causal interactions correctly predict biological causality and are merely statistical predictions. It should also be noted that such predictions do not provide explanations regarding the origin of the interactions and could be due to an indirect influence. For instance, one time series may be a strong predictor of another time series because both are shaped together by a common underlying cause. Hence, like any other statistical prediction, a cautious interpretation of each predicted interaction is required to be made before building any hypothesis. Incorporation of functional data like metabolic co-dependencies (Levy et al., 2015) might help to strengthen the basis of a predicted interaction.

### CONCLUSION

The various workflows implemented in 'TIME' can help end users not only to perform a number of analyses, but also gain meaningful insights from the interactive visualizations. Analysis on a few well known publicly available datasets illustrate the utility of the options available in 'TIME.' For example, apart from obtaining information regarding the temporal effect of antibiotic treatment on human gut microbiome, 'TIME' could identify similarly perturbed groups of microbial genera. Additionally, the inter-microbial competition among the pathogens and commensals could also be inferred from the causality networks and stationarity analysis. In another example, the periodic changes in community structure of the vaginal microbiome were illustrated using the 'community state' analysis workflow. Although the scope of the case studies presented here is limited in this communication, the workflows can be further utilized to gain additional insights. We expect 'TIME' to be a valuable contribution in the field of microbial time series data analysis and visualization.

### MATERIALS AND METHODS

'TIME' web application uses Python and JavaScript to execute the backend algorithms and for browser based data processing, respectively. We used the DyGraphs<sup>1</sup> (DyGraphs Java Script,

<sup>1</sup>http://dygraphs.com/

2017) module for rendering time series line charts since it has the ability to handle large datasets seamlessly. Other visualizations are implemented using D3.js library (Bostock, 2017) with extensive interactive operations.

### Input Format

User data (consisting of the microbial abundance table) along with the available metadata information can be incorporated in 'TIME' using a simple form. The abundance table can either be provided as a standard 'QIIME' output (Caporaso et al., 2010) or as a tab delimited file. The metadata file is required to have information related to the source of the microbiome sample, sample names (identical to the ones in the abundance file), time stamp information along with the sample condition for each time point. A detailed description of the input files is provided in the user manual (available in the website).

### Normalization, Visual Exploration and Segregation of Microbiome Time Series Data

The microbial abundances obtained for analysis represents the count of clustered sequences belonging to the constituent taxa as operational taxonomic units (OTUs). The abundances of each OTU across different time points constitute the OTU abundance matrix. Restraints in sampling at multiple time points as well as sequencing errors result in unequal sequencing depths. 'TIME' provides methods to circumvent this limitation using either a proportion based or rarefaction based normalization. Rarefaction plots serve as one of the means to identify unequally sampled data points and subsequently can be used to normalize the OTU matrices such that all time points have similar counts. Users can generate a rarefaction curve for each metagenomic source by selecting either all the time points or a set of equidistant 5 or 10 time points. The generated curve can be used as a guide to select a suitable rarefaction normalization depth. Alternatively, users may proceed with absolute count data (without any normalization) or perform relative proportion based normalization. It is advisable to choose appropriate normalization method (refer to the "Discussion" section for more details).

The visual examination of the temporal trends is an important step in any time series analysis. In 'TIME,' all the taxa abundances at any particular taxonomic level can be viewed together as interactive line plots across the sampled timeline. An important challenge for carrying out such comparative microbial data analysis pertains to the problem of taxa abundances with

case study generated using 'Workflow-2.'

different orders of magnitude (with some taxa having very high abundances and some having extremely low counts). In other words, it is difficult to visualize the trends of the lower abundant taxa owing to the dominant influence of the very high abundant ones on the plot. 'TIME' provides two ways to tackle this problem. While one uses 'quartile segregation,' the other utilizes 'log scaling.' In quartile segregation, the different taxa are grouped into four quartiles based on their abundance information which can be viewed separately. The very high abundant and the very low abundant taxonomic groups (or potential outliers) tend to occupy the top and bottom quartiles, respectively. The remaining quartile accommodates the taxonomic groups having the intermediate abundances. This makes sure that during visual exploration the temporal trends of the low abundance taxa do not get compressed (or dominated) by the trends of the very high abundant taxa. TIME also offers the option of log scaling the abundance values so that the trends of low abundant and that of high abundant taxa can be compared on the same plot. Additionally, the tool provides an option to view the core, persistent and the transient groups of bacteria which are reported to have distinct roles in microbial ecosystems (Caporaso et al., 2011). A taxon is considered to be persistent if it is observed in more than 20% of the time points, with at least 90% of these observations being consecutive. On the other hand, the transient taxa are those which are observed in at least 60% of the time points, with at most 75% of these observations being consecutive. However, TIME provides an option to modify these parameters (prevalence threshold and consecutive observations) for definition of core, persistent and transient in 'workflow 1' to accommodate differences in wide number of datasets. In addition to the above measures, the richness and diversity of the studied microbial communities are also calculated using well known indices. While richness of a microbiome denotes the unique number of constituent taxa present in each sample (at a time point), the diversity provides a measure of how evenly the taxonomic entities are distributed. Although diversity of a microbiome can be calculated using a number of ways, the widely accepted Shannon index for diversity (Shannon, 1948) has been implemented in 'TIME.'

$$\text{Shannon } Index \ = -\sum\_{i=1}^{R} p\_i \ \ln \ p\_i$$

Where, p<sup>i</sup> refers to the proportion of the abundance of the i th taxon in the population consisting of 'R' taxa.

With implementation of each of the above methods and their corresponding visualizations, an interactive operation for selecting a 'subplot window' of the plotted timescale is presented. This feature enables users to graphically choose the start and end time points using simple mouse operations and visualize the selected time range. The 'subplot window' can then be dragged along the time-scale with a zoomed view. The moving average of a time series can also be specified using a text box available at the bottom corner of the plot window. This feature smoothes the short term fluctuations in the time series data and shows the overall trends (and cycles) across a longer timescale.

### Identification of Stationary Taxonomic Groups

Stationarity of taxa in microbial time series data is important to understand inter-microbial competition (David et al., 2014). A taxon is considered stationary if its mean, variance, covariance, and autocorrelation are constant over time, due to the absence of a unit root process. A unit root process is said to be present in a time series if its autoregressive model has an estimated coefficient close to one. The presence of a unit root indicates that a perturbation in the value of the entity in the time series has a persistent impact on its future values and hence a cause of non-stationarity. The most commonly used method for calculating stationarity is the Augmented Dickey Fuller (ADF) Test. While the null hypothesis of the ADF Test is that there is a unit root process governing the dynamics of the entity (taxon), the alternate hypothesis states that there is no unit root. The ADF test statistic is a number, the more negative it is, the stronger will be the confidence with which the null hypothesis can be rejected. 'TIME' allows identification of microbial groups detected to be stationary and non-stationary and lists the same in a searchable table with options to export the results. The stationarity information corresponding to each taxon is also used to augment other plots in the tool along with the phylogeny information.

### Generation of Inter-Microbial Causality Network

One of the important objectives in time series studies pertains to identifying causal relationships among entities. Causality aims to find the direct interactions between entities such that one entity can trigger/suppress or be triggered/suppressed by the other. It should be noted that causation should not be confused with correlation. For example, two taxa ('A' and 'B') in a microbiome time series dataset may be correlated, but may not have any causal relationships. Granger Causality (Granger, 1969) is one of the most well established statistical tests for checking causality among two time series. The basic premise of this method is that, if one variable causes another, the past values of the former must have some information about the future values of the latter (which is not available otherwise). For example in a microbiome time series data, if taxon A affects taxon B, the future values of taxon B can be better predicted using the past values of both taxa A and B, rather than using the past values of taxon B alone. In order to ascertain if taxon A influences taxon B, two regressions are performed. In one, past values of taxa A and B are used to predict the present values of taxon B. In the other, only the past values of taxon B is used to predict the present values of taxon B. If a significant increase of the goodness of fit of the former regression over the latter is observed, then taxon A is said to 'Granger cause' taxon B.

Since a typical microbiome time series data has more than two entities (taxa), 'Granger Lasso Causality' (Hlavácková- ˇ Schindler and Pereverzyev, 2015) can be used to find causal relationships among all taxa. Thus, apart from the 'Pairwise Granger causality' (described above) for all possible taxa pairs, 'Granger-Lasso' method has also been implemented in 'TIME.' The LASSO (Least Absolute Shrinkage and selection operator) is one of the most well known and widely used methods for feature selection and regularization in machine learning. LASSO works by adding a regularizing penalty to the sum of squared errors. This objective function is minimized (by optimization) for estimating the values of the coefficients of regression (thus reducing the weightage/coefficients of the unimportant predictors), thereby finding the best set of predictors for every variable. Granger-Lasso utilizes the LASSO methodology for identifying causal relationships among all entities (taxa) in multivariate microbiome time series dataset (Arnold et al., 2007). Another option allowing selection of causality pairs predicted by both 'Pairwise Granger' and 'Granger-Lasso' is implemented for improved Granger Causality predictions. All these three methods are available in 'TIME' which can be finally used to generate a directed causality network.

### Identifying Taxonomic Groups Having Similar Temporal Patterns

In time series datasets, it is not only important to evaluate temporal changes of different entities and the causal relationships among them, but also to identify entities which exhibit similar temporal patterns. The Euclidean distance based clustering of entities is unsuitable for identifying similar temporal patterns since this distance measure does not take into account the distortion across time series (Keogh and Ratanamahatana, 2005). In other words, the temporal behavior of two taxa which are out of phase is assigned a high value by Euclidian measure. On the other hand, Dynamic Time Warping (DTW) gives due importance to the phase displacement and obtains the optimal alignment between the two time series (Berndt and Clifford, 1994). DTW uses a dynamic programming based approach to align and score the similarity of the temporal patterns corresponding to two entities (taxa in the case of microbial time series). Since the DTW algorithm is relatively slow with a worst case time complexity of O(n<sup>2</sup> ), a modified DTW algorithm (Sakoe and Chiba, 1990) is implemented in 'TIME'. In this algorithm, a constraint is applied in such a way that a limited number of cells are evaluated during computing the cost matrix of the alignment, thereby making the overall computation process much faster (Salvador and Chan, 2007). 'TIME' uses the calculated pair-wise DTW distances among the different taxa for hierarchical clustering. The resulting dendrograms can be viewed as trees in standard or radial layouts. One of the limitations of the DTW distance pertains to the

inability to interpret the distance score easily as it does not fall in a definite range. Therefore, it is desirable to have a modified score with a definitive range that can be universally interpreted. In order to achieve this, we introduce a new method for calculating the distance between two time series, called the 'TIME-DTW Distance,' In a microbial time series data, one taxon can have a difference of several orders of magnitude with another, but their time series may have similar overall shape. Thus, a standard normalization step is first applied to minimize such differences. Following this, the DTW distance is calculated and normalized by the average 'sum of the absolute difference' (SAD) between each time series and its 'mirror image' (Supplementary Material). The resultant value ('TIME-DTW distance') will hence always fall between a range of zero and one. 'TIME' allows an easy and interactive way to explore the results using a 'clustered heat map.' In addition to using the TIME DTW Distance as the measure of similarity/dissimilarity, the pairwise similarity between taxa can also be viewed using Pearson Correlation coefficient. The resulting heatmaps are hierarchically clustered based on their distances along the vertical axis, and taxonomic hierarchies along the horizontal axis.

### Understanding Community Structure Based on Similarities across Time Points

Apart from understanding the temporal similarities among the resident entities (taxa), clustering of time points having similar entity distribution is expected to yield valuable insights regarding the microbial community dynamics. The identified time points having similar taxonomic distributions (i.e., phylotype proportions) can be considered as a 'community state' (Gajer et al., 2012). Jenson Shannon divergence (JSD) metric has been utilized earlier to identify such 'community state' in microbiome time series data (Gajer et al., 2012). In 'TIME,' a modification of the method is implemented to make it applicable for any microbiome time series. The taxa abundances are first normalized to generate probability distributions, which are then used to calculate the JSD among the different time points. Thus, a pairwise JSD matrix is obtained for all time points. Since, the K-medoids clustering algorithm (Jin and Han, 2011) is known

### REFERENCES

to be robust to noise and outliers (as compared to K-means), it was utilized for clustering the time points using the generated JSD matrix. The number of clusters can be chosen by the user based on visual inspection. Along with clustering the different samples based on their microbial community structure, it is often useful and sometimes necessary to find the drivers of the cluster, i.e., the most dominant taxa among the clusters. Keeping this in view, 'TIME' also provides the option to view the (normalized) relative abundances of the taxa across different clusters and different time points as a heatmap, which helps in visual exploration and determination of the distinctive/driver taxa or groups of taxa.

### AUTHOR CONTRIBUTIONS

KB, BK, and SM conceived the idea and designed the overall methodology. KB and BK implemented the algorithms and developed the web server. BK performed the case studies. BK, KB, and SM analyzed the results and drafted the manuscript. All authors read and approved the final manuscript.

### ACKNOWLEDGMENTS

We acknowledge Dr. Swadha Anand and Dr. Anirban Dutta for their help during the course of this work. We also acknowledge Mr. Pramod Nanadikar and Mr. Govind Gopal for their help in deploying the application. BK is also an industry sponsored Ph.D. student at Chemical Engineering and Process Development Division, CSIR-National Chemical Laboratory, Pune 411008 (India) and would like to acknowledge the Academy of Scientific and Innovative Research (AcSIR) for its support.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.00036/full#supplementary-material

Ernst, J., and Bar-Joseph, Z. (2006). STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics 7:191. doi: 10.1186/1471-2105-7-191

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Baksi, Kuntal and Mande. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Bayesian Semiparametric Regression Model for Joint Analysis of Microbiome Data

Juhee Lee<sup>1</sup> \* and Marilou Sison-Mangus <sup>2</sup>

*<sup>1</sup> Department of Applied Mathematics and Statistics, University of California, Santa Cruz, Santa Cruz, CA, United States, <sup>2</sup> Department of Ocean Sciences, University of California, Santa Cruz, Santa Cruz, CA, United States*

The successional dynamics of microbial communities are influenced by the synergistic interactions of physical and biological factors. In our motivating data, ocean microbiome samples were collected from the Santa Cruz Municipal Wharf, Monterey Bay at multiple time points and then 16S ribosomal RNA (rRNA) sequenced. We develop a Bayesian semiparametric regression model to investigate how microbial abundance and succession change with covarying physical and biological factors including algal bloom and domoic acid concentration level using 16S rRNA sequencing data. A generalized linear regression model is built using the Laplace prior, a sparse inducing prior, to improve estimation of covariate effects on mean abundances of microbial species represented by operational taxonomic units (OTUs). A nonparametric prior model is used to facilitate borrowing strength across OTUs, across samples and across time points. It flexibly estimates baseline mean abundances of OTUs and provides the basis for improved quantification of covariate effects. The proposed method does not require prior normalization of OTU counts to adjust differences in sample total counts. Instead, the normalization and estimation of covariate effects on OTU abundance are simultaneously carried out for joint analysis of all OTUs. Using simulation studies and a real data analysis, we demonstrate improved inference compared to an existing method.

Keywords: count data, Laplace prior, metagenomics, microbiome, regularizing prior, process convolution, negative binomial model, 16S ribosomal RNA sequencing

### 1.  Their successional dynamics could further change in response to perturbations of the host or of the surrounding environments (Turnbaugh et al., 2009; Needham and Fuhrman, 2016). Understanding how abiotic and biotic factors influence the dynamics of microbial communities is of great interest in the field of microbiome studies. Recent revolutionary advances in next-generation sequencing (NGS) technologies along with rapidly decreasing costs, have facilitated the accumulation of large datasets of 16S ribosomal RNA (rRNA) amplicon sequences across various disciplines such as medicine, biology, ecology, and environmental sciences (Woo et al., 2008). Sequencing data is usually pre-treated for quality filtering, noise removal and chimera checking through bioinformatics algorithms and the filtered sequences are clustered into Operational Taxonomic Units (OTUs), which represent similar organisms (microbial species) based on sequence homology (called OTU picking). An OTU abundance table is generated,

### Edited by:

*Michele Guindani, University of California, Irvine, United States*

#### Reviewed by:

*Yanxun Xu, Johns Hopkins University, United States Michael Pester, Deutsche Sammlung von Mikroorganismen und Zellkulturen (DSMZ), Germany*

> \*Correspondence: *Juhee Lee<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Systems Microbiology, a section of the journal Frontiers in Microbiology*

Received: *15 November 2018* Accepted: *08 March 2018* Published: *26 March 2018*

#### Citation:

*Lee J and Sison-Mangus M (2018) A Bayesian Semiparametric Regression Model for Joint Analysis of Microbiome Data. Front. Microbiol. 9:522. doi: 10.3389/fmicb.2018.00522* recording counts for OTUs in samples. Further statistical data analyses are then performed using the OTU table to answer biological and ecological questions.

Analysis of huge NGS data is computationally expensive and challenging. One of the key challenges is the normalization of counts across samples. Total counts (often called library size or sequencing depth) may vastly vary across different samples due to technical reasons. Thus, observed counts are not directly comparable across samples and cannot be used as a measure of the abundance of an OTU. Normalized counts through rarefaction or relative frequencies are commonly used for easy comparison of OTU abundance across samples. However, such ad hoc normalization procedures have been criticized from a statistical perspective since using pre-normalized quantities may undermine the performance of downstream analysis (McMurdie and Holmes, 2014). Another challenge is high dimensionality and sparsity in OTU count data. A dataset typically includes hundreds or thousands of OTUs and a majority of them has zero or very low frequencies in most of samples. For example, **Figure 1A** illustrates a heatmap of OTU counts in our motivating dataset described in section 2.3. It shows that a majority of OTUs has very low counts (gray) in a sample, and the set of OTUs having large counts (blue) vary across samples. Due to such sparsity in data, borrowing strength across OTUs through joint analysis of all OTUs is crucial for improved inference. Recently, various statistical methods including Romero et al. (2014), Chen and Li (2016), Gibbons et al. (2017), and Zhang et al. (2017) have been developed for microbiome studies using NGS data. For example, Zhang et al. (2017) used a negative binomial mixed regression model to study interactions between the microbiome and host environmental/clinical factors. Random effects are used to induce correlation among samples from a group. Common to most of recent methods including Zhang et al. (2017) is separately analyzing each OTU at a time.

We develop a Bayesian semiparametric generalized linear regression model to study the effects of physical and biological factors on abundance of microbes. The proposed method performs mode-based normalization through a hierarchical model, which enables direct modeling of OTU counts. Furthermore, the hierarchical model facilitates borrowing strength between OTUs, between samples, and between time points through joint analysis and improves inference on the effects of covariates **X** on OTU abundance which are the parameters of our primary interest. Specifically, a negative binomial (NB) distribution parameterized by a mean parameter µ and an overdispersion parameter s is assumed for OTU counts. The NB distribution flexibly accommodates overdispersion often seen in NGS data and is commonly used as a robust alternative to a Poisson distribution (Anders and Huber, 2010). The expected count µ of an OTU is decomposed as a product of factors, a baseline mean count g and a nonnegative function η(**X**) of covariates that describes their effects on the mean count. We use the log link function for η(**X**) and assume that change in a covariate has a multiplicative effect on mean count, where the associated coefficient quantifies the size and direction of the effect. We consider a Laplace prior for the coefficients, a shrinkage prior that is essential in a high dimensional regression setting. Shrinkage priors in regression yield sparse point estimates of the coefficients, where many of the coefficients have values close to zero and few have large values. The sparse estimates improve out-of-sample prediction and produce more interpretable models (Park and Casella, 2008). In addition, shrinkage priors such as a Laplace prior in a regression problem mitigate potential problems by multicollinearity and yield improved coefficient estimates when covariates are high-dimensional and potentially highly correlated (Polson and Scott, 2012). For baseline mean counts, we develop a nonparametric model to combine all OTUs for joint analysis. Baseline mean counts may vary across samples and OTUs. Also, as in our motivating data for which samples were taken over time, there may be temporal dependence in baseline mean counts. To tackle the problem, we further decompose the baseline count g into sample size factor (r), OTU size factor (α0), and OTU and time factor (αt), that is, <sup>g</sup> <sup>=</sup> <sup>r</sup> <sup>×</sup> <sup>α</sup><sup>0</sup> <sup>×</sup> <sup>α</sup><sup>t</sup> . Due to the overparametrization of the baseline mean abundance, individual factors are not identifiable. To avoid identifiability issues, we place the regularizing priors with mean constraints (Li et al., 2017) for sample size factor r and OTU size factor α0. In addition, we model a temporal dependence structure between the baseline expected counts for an OTU through a convolutional Gaussian process (Higdon, 1998). The process convolution approach is often used as an alternative approach of the Gaussian process to construct a dependent process due to its efficient computation (Lee et al., 2005; Liang and Lee, 2014). Through simulation studies, we show that estimates of individual parameters r, α0, and α<sup>t</sup> are not fully interpretable under the proposed model, but baseline mean counts g are identifiable. The model also provides a posterior distribution of g for uncertainty quantification.

The rest of the paper is organized as follows. In section 2 we describe the proposed model and discuss the prior formulations and the resulting posterior inference. We perform simulation studies to assess the proposed model and perform comparison with an existing method that analyzes one OTU at a time. We then apply the proposed model to an ocean microbiome dataset. Section 3 presents the performance of the proposed model from the simulation experiment and the ocean microbime data. Section 4 concludes the paper with a discussion on limitations and possible extensions.

### 2. MATERIALS AND METHODS

### 2.1. Bayesian Semiparametric Regression Model

Suppose that samples are taken at n different time points, 0 ≤ <sup>t</sup><sup>i</sup> <sup>≤</sup> <sup>T</sup>, <sup>i</sup> <sup>=</sup> 1, ... , <sup>n</sup>, and with <sup>K</sup><sup>i</sup> replicates at time point <sup>t</sup><sup>i</sup> . We consider count yt<sup>i</sup> ,k,<sup>j</sup> of OTU j in replicate k taken at time ti , where <sup>i</sup> <sup>=</sup> 1, ... , <sup>n</sup>, <sup>k</sup> <sup>=</sup> 1, ... , <sup>K</sup><sup>i</sup> , and <sup>j</sup> <sup>=</sup> 1, ... , <sup>J</sup>. A sample is thus indexed by t<sup>i</sup> and k. We let the total number of samples N = P<sup>n</sup> <sup>i</sup>=<sup>1</sup> <sup>K</sup><sup>i</sup> . Let **<sup>Y</sup>** <sup>=</sup> [yt<sup>i</sup> ,k,j ] denote the N × J matrix of counts, where yt<sup>i</sup> ,k,j is integer-valued and nonnegative. Also, suppose that covariates **<sup>X</sup>**t<sup>i</sup> <sup>=</sup> (Xt<sup>i</sup> ,1, ... , Xt<sup>i</sup> ,P) ′ are recorded at

time t<sup>i</sup> . For example, covariates are physical and biological factors in our motivating data.

2.1.1. Sampling Model

Count data by NGS methods is often modeled through a Poisson distribution. The assumption under the Poisson distribution that the variance is equal to the mean is often too restrictive to accommodate overdispersion that variation in data exceeds the mean. The negative binomial (NB) distribution is a popular and convenient alternative to address the overdispersion problem and is widely recognized as a model that provides improved inference to NGS count data (for example, see Robinson and Smyth, 2007; Anders and Huber, 2010). A NB distribution can be characterized by mean and overdispersion parameters. We suppress index i for simpler notation and assume a NB model for count yt,k,<sup>j</sup> of OTU

 $j$  in predicate  $k$  at time  $t$ ,

$$\mathcal{Y}\_{t,k,j} \stackrel{indep}{\sim} \text{NB}(\mu\_{t,k,j}, s\_j),\tag{1}$$

where mean count µt,k,<sup>j</sup> > 0 and overdispersion parameter s<sup>j</sup> > 0. The model in Equation (1) implies that count of OTU j in replicate k at time t has mean E(yt,k,<sup>j</sup> <sup>|</sup> <sup>µ</sup>t,k,<sup>j</sup> ) <sup>=</sup> <sup>µ</sup>t,k,<sup>j</sup> and variance Var(yt,k,<sup>j</sup> <sup>|</sup> <sup>µ</sup>t,k,<sup>j</sup> ,sj) <sup>=</sup> <sup>µ</sup>t,k,<sup>j</sup> <sup>+</sup> <sup>µ</sup> 2 t,k,j sj . The model allows different dispersion levels across OTUs through OTU-specific overdispersion parameters s<sup>j</sup> . In the limit as s<sup>j</sup> → 0, the model in Equation (1) yields the Poisson distribution with mean µt,k,<sup>j</sup> . We assume a gamma distribution for a prior distribution of s<sup>j</sup> , sj iid ∼ Ga(a<sup>s</sup> , <sup>b</sup>s), <sup>j</sup> <sup>=</sup> 1, ... , <sup>J</sup>, with fixed <sup>a</sup><sup>s</sup> and <sup>b</sup><sup>s</sup> .

Frontiers in Microbiology | www.frontiersin.org

### 2.1.2. Model for Regression

We next model the mean count µt,k,<sup>j</sup> of yt,k,<sup>j</sup> . We decompose the mean count into factors, a baseline mean count and a function of covariates, <sup>µ</sup>t,k,<sup>j</sup> <sup>=</sup> <sup>g</sup>t,k,jηj(**X**t). Here parameter <sup>g</sup>t,k,<sup>j</sup> denotes the baseline mean abundance of OTU j in sample (t, k) and ηj(**X**t) is a function of covariates **X**<sup>t</sup> for OTU j to model the covariate effects. We construct a generalized regression model by letting log(ηj(**X**t)) <sup>=</sup> **<sup>X</sup>** ′ tβj , where <sup>β</sup><sup>j</sup> <sup>=</sup> (βj1, ... , <sup>β</sup>jP) ′ is a P-dimensional vector of regression coefficients of OTU j (Lawless, 1987; McCullagh and Nelder, 1989). The coefficient βj,<sup>p</sup> quantifies the effect of covariate p X<sup>p</sup> on the mean abundance of OTU j. A vector β<sup>j</sup> close to the zero vector produces a value of ηj(**X**t) close to 1, and the mean count remains similar to the baseline mean count gt,k,<sup>j</sup> , implying insignificant covariate effects. A negative (positive) of βj,<sup>p</sup> implies a negative (positive) association between mean counts and the p-th covariate, and a larger value of Xj,<sup>p</sup> decreases (increases) the mean count, while holding the other covariates constant. We consider a Laplace prior on β<sup>j</sup> . Specifically, we express the Laplace distribution as a scale mixture of normals and assume for <sup>j</sup> <sup>=</sup> 1, ... , <sup>J</sup> and <sup>p</sup> <sup>=</sup> 1, ... , <sup>P</sup>,

$$\begin{aligned} \alpha\_{j,p} \mid \sigma\_j^2, \phi\_{j,p} \stackrel{indep}{\sim} \text{N}(0, \sigma\_j^2 \phi\_{j,p}), \quad \phi\_{j,p} \stackrel{indep}{\sim} \text{Exp}(\frac{\lambda\_j^2}{2}),\\ \lambda\_j^2 \stackrel{iid}{\sim} \text{Ga}(a\_\lambda, b\_\lambda), \quad \sigma\_j^2 \stackrel{iid}{\sim} \text{IG}(a\_\sigma, b\_\sigma), \end{aligned} \tag{2}$$

where aλ, bλ, a<sup>σ</sup> , and b<sup>σ</sup> are fixed. σ 2 j and φj,<sup>p</sup> denote the global and local shrinkage parameters, respectively, for OTU j. After integrating φj,<sup>p</sup> out, the prior distribution of βj,<sup>p</sup> is the Laplace distribution with location parameter 0 and scale parameter q σ 2 j /λ<sup>j</sup> , that is, <sup>p</sup>(βj,<sup>p</sup> <sup>|</sup> <sup>λ</sup> 2 j , σ 2 j ) <sup>∝</sup> exp <sup>−</sup>λ<sup>j</sup> <sup>|</sup>βj,p|/ q σ 2 j . Compared to a normal distribution that is a common choice for to experimental artifacts. For example, counts of an OTU even in the replicates taken at a time point may vastly differ. Sample specific size factors rt,<sup>k</sup> account for different total counts in different samples and expected counts normalized by rt,<sup>k</sup> are comparable across samples. Factor α0,<sup>j</sup> explains variabilities in baseline mean abundances of OTUs and αt,<sup>j</sup> models temporal dependence of the mean counts for an OTU, respectively. Factors α0,<sup>j</sup> and αt,<sup>j</sup> are not indexed by replicate k and account for stochastic change over time in normalized baseline expected counts of OTU j. Collecting all, we write the mean count as

$$
\mu\_{t,k,j} = \lg\_{t,k,j} \eta\_j(\mathbf{X}\_t) = r\_{t,k} \alpha\_{0,j} \alpha\_{t,j} \eta\_j(\mathbf{X}\_t), \tag{3}
$$

The model for gt,k,<sup>j</sup> in Equation (3) is overparameterized and the individual parameters are not identifiable. To avoid potential identifiability issues, many of NB models rely on some form of approximation for the baseline mean counts. For example, one may find the maximum likelihood estimates (MLEs) of baseline mean abundance under some constraints and plug in those estimates to infer the mean abundance levels µt<sup>i</sup> ,<sup>j</sup> of OTUs (Witten, 2011). Plugging in MLEs is simple but may not be robust. In particular, the inference is greatly affected by a small change in a few OTUs that have large counts. Moreover, the errors introduced in the baseline mean count estimation will not be reflected in the inference. Several approaches to robustify the estimates are proposed (for example, see Anders and Huber, 2010; Witten, 2011). To circumvent the identifiability issue and provide uncertainty quantification for estimation of gt,k,<sup>j</sup> , we take an alternative in Li et al. (2017) by imposing regularizing priors with mean constraints for rt,<sup>k</sup> and α0,<sup>j</sup> . We let the logarithm of the factors r˜t,<sup>k</sup> = log(rt,<sup>k</sup> ) and <sup>α</sup>˜ 0,<sup>j</sup> <sup>=</sup> log(α0,j), and assume the regularizing prior distribution with mean constraints,

$$\begin{split} \tilde{r}\_{l,k} \mid \psi^r, \eta^r, w^r, \nu\_r^2, c\_r \stackrel{\text{iid}}{\sim} & \quad \sum\_{\ell=1}^{L'} \psi\_\ell^r \left\{ w\_\ell^r \phi(\eta\_\ell^r, \nu\_r^2) + (1 - w\_\ell^r) \phi\left(\frac{c\_r - w\_\ell^r \eta\_\ell^r}{1 - w\_\ell^r}, \nu\_r^2\right) \right\}, \\ \tilde{\alpha}\_{0,j} \mid \psi^a, \eta^a, w^a, \nu\_a^2, c\_a \stackrel{\text{iid}}{\sim} & \quad \sum\_{\ell=1}^{L''} \psi\_\ell^a \left\{ w\_\ell^a \phi(\eta\_\ell^a, \nu\_a^2) + (1 - w\_\ell^a) \phi\left(\frac{c\_a - w\_\ell^a \eta\_\ell^a}{1 - w\_\ell^r}, \nu\_a^2\right) \right\}, \end{split} \tag{4}$$

the prior of βj,p, the Laplace distribution has more concentration around zero but allows heavier tails. The regularized regression through the Laplace prior more shrinks the coefficients of insignificantly related covariates into zero and less pulls the coefficients of important covariates toward zero. Shrinkage of β estimates through the model in Equation (2) mitigates possible issues due to multicollinearity and efficiently improves estimation of β in a high dimensional setting (Polson and Scott, 2012).

#### 2.1.3. Model for Baseline Mean Count

We next build a prior probability model for the baseline mean count gt,k,<sup>j</sup> of OTU j in sample (t, k). We assume <sup>g</sup>t,k,<sup>j</sup> <sup>=</sup> <sup>r</sup>t,kα0,jαt,<sup>j</sup> to separate sample (rt,<sup>k</sup> ), OTU (α0,j), and OTU-time (αt,j) factors. Sample total counts <sup>y</sup>t,k,· <sup>=</sup> <sup>P</sup><sup>J</sup> j=1 yt,k,<sup>j</sup> may greatly differ for different samples possibly due where φ(η, v 2 ) is the probability density function of the normal distribution with mean η and variance v 2 , constraints for the mixture weights P<sup>L</sup> r <sup>ℓ</sup>=<sup>1</sup> <sup>ψ</sup> r ℓ = P<sup>L</sup> α <sup>ℓ</sup>=<sup>1</sup> <sup>ψ</sup> α <sup>ℓ</sup> <sup>=</sup> 1 with 0 < ψ<sup>r</sup> <sup>ℓ</sup> < 1 and 0 < ψ<sup>α</sup> <sup>ℓ</sup> < 1 , 0 < w r <sup>ℓ</sup> < 1, and 0 < w α <sup>ℓ</sup> < 1 for all ℓ. Mixture models as in Equation (4) are often used as a basis to approximate any distribution. Each component in Equation (4) is further composed of a mixture of two normals, N(ηℓ, v 2 ) and N (c−wℓηℓ) (1−wℓ) , v 2 with weights w<sup>ℓ</sup> and 1 − wℓ, respectively, and the mean of the component isc. In consequence, the prior and posterior of <sup>r</sup>˜ and <sup>α</sup>˜ under the model in Equation (4) satisfy their prespecified mean constraints c<sup>r</sup> and cα, respectively. Li et al. (2017) showed that the model in Equation (4) flexibly accommodates various features in a distribution such as skeweness or multi-modality while satisfying the constraints. Furthermore, the model based normalization through Equation (4) enables joint analysis of all OTUs and can further improve estimation of the covariate effects. With the regularizing priors, baseline mean counts gt,k,<sup>j</sup> are identifiable, while rt,<sup>k</sup> , α0,<sup>j</sup> , and αt,<sup>j</sup> are not directly interpretable. More importantly, the parameters of primary interest ηj(**X**t) can be uniquely estimated and βj,p's keep their interpretation as parameters that quantify the effects of covariates on mean counts. We used an empirical approach to fix the mean constraints c<sup>r</sup> and cα. Sensitivity analyses were conducted to assess the robustness to the specification of c<sup>r</sup> and cα and show that the model provides reasonable estimates of gt,k,<sup>j</sup> and moderate changes in the values of c<sup>r</sup> and c<sup>α</sup> minimally change the estimates. More details of the specification of c<sup>r</sup> and cα are discussed in section 3.1. We fix the numbers of mixture components, L r and L α and variances v 2 r and v 2 α . We let η r ℓ iid ∼ N(c<sup>r</sup> , ω 2 r ) and η α ℓ iid <sup>∼</sup> N(cα, <sup>ω</sup> 2 α ), where ω 2 r and ω 2 α are fixed. We assume ψ <sup>r</sup> <sup>=</sup> (<sup>ψ</sup> r 1 , ... , ψ r L <sup>r</sup> ) ∼ Dir(d<sup>r</sup> , ... , dr) and ψ <sup>α</sup> <sup>=</sup> (ψ α 1 , ... , ψ α L <sup>α</sup> ) <sup>∼</sup> Dir(dα, ... , <sup>d</sup>α), with fixed <sup>d</sup><sup>r</sup> and <sup>d</sup>α. We let w r ℓ iid ∼ Be(a<sup>r</sup> , <sup>b</sup>r), <sup>ℓ</sup> <sup>=</sup> 1, ... , <sup>L</sup> r and w α ℓ iid <sup>∼</sup> Be(aα, <sup>b</sup>α), <sup>ℓ</sup> <sup>=</sup> 1, ... , L <sup>α</sup> with fixed a<sup>r</sup> , b<sup>r</sup> , aα, and bα.

Recall that samples are collected over time points t1, ... , t<sup>n</sup> in [0, T] and αt,<sup>j</sup> accounts for temporal dependence in the baseline mean count for an OTU. We let <sup>α</sup>˜t,<sup>j</sup> <sup>=</sup> log(αt,j) a function in time t and use a stochastic process to model temporal dependence among µt,k,<sup>j</sup> . The Gaussian process (GP) is one of the most popular stochastic models for the underlying process in spatial and spatio-temporal data (for example, see Cressie, 1992; Banerjee et al., 2014 among many others). The GP effectively represents the underlying phenomenon in a variety of applications, but it has some drawbacks such as a complex computation that requires a matrix decomposition and problematic estimation of the parameters in its covariance function, potentially leading to difficulties in exploring the posterior distribution (Lee et al., 2005; Liang and Lee, 2014). To alleviate such difficulties of GP models while still maintaining their flexibility and adaptiveness, we use a convolution approach with a kernel function developed in Higdon (1998, 2002). For each OTU, we specify the latent process θj(t) to be nonzero only at the time points u1, ... , u<sup>M</sup> in [0, T]. Specifically, we consider the GP convolution model,

$$\tilde{\alpha}\_{t\bar{\jmath}} = \sum\_{m=1}^{M} Z(t - \mu\_m) \theta\_{m\bar{\jmath}\*}$$

where {u1, ... , <sup>u</sup>M} a set of basis points in [−<sup>t</sup> ′ 1 , T + t ′ 2 ] with t ′ 1 , t ′ <sup>2</sup> <sup>&</sup>gt; 0, and <sup>Z</sup>(<sup>t</sup> <sup>−</sup> <sup>u</sup>m) a Gaussian kernel centered at <sup>u</sup>m, Z(t − um) = √ 1 2πγ <sup>2</sup> exp{ −(t−um) 2 2γ 2 }. The number of basis points M, their locations u<sup>m</sup> and the range parameter γ can be treated as random variables by placing prior distributions, e.g., consider a gamma prior for γ . For simplicity, we fix them as follows. We first choose a value for M and let u<sup>m</sup> evenly spaced over time [−t ′ 1 , T + t ′ 2 ]. Following Xiao (2015), we let the range parameter γ <sup>2</sup> <sup>=</sup> ((2<sup>T</sup> <sup>+</sup> <sup>t</sup> ′ <sup>1</sup> + t ′ 2 )/M) 2 , that is, the range parameter depends on the value of M. Through simulations, we studied the impact of different values of M on the posterior inference of gt,k,<sup>j</sup> . A discussion is included in section 3.1. Given the number of basis points M, we assume θm,<sup>j</sup> | τ 2 j indep <sup>∼</sup> N(0, <sup>τ</sup> 2 j ) and τ 2 j iid ∼ IG(a<sup>τ</sup> , b<sup>τ</sup> ), <sup>m</sup> <sup>=</sup> 1, ... , <sup>M</sup> and <sup>j</sup> <sup>=</sup> 1, ... , <sup>J</sup>.

We implement posterior inference on the parameters ˜<sup>θ</sup> <sup>=</sup> (βj , σ 2 j , λ 2 j , <sup>φ</sup>j,p,r˜t,<sup>k</sup> , ψ r ,w r ℓ , η r ℓ , <sup>α</sup>˜ 0,<sup>j</sup> , <sup>α</sup>˜t,<sup>j</sup> , ψ α ,w α ℓ , η α ℓ , θ<sup>j</sup> , τ 2 j ,sj) via a Markov chain Monte Carlo (MCMC) method based on Metropolis-Hastings algorithm and Gibbs sampling. Each of the parameters is iteratively updated conditional on the currently computed values of all other parameters to simulate a sample from the posterior distribution. The parameters <sup>r</sup>˜ and <sup>α</sup>˜ <sup>0</sup> jointly determine baseline mean counts and joint updating of <sup>r</sup>˜ and <sup>α</sup>˜ <sup>0</sup> may greatly improve the mixing. In our ocean microbiome data, some discretized covariates are missing. We treat them as random variables by assuming a uniform distribution over possible categories, and impute their values in MCMC simulation. Full details of our MCMC algorithm are given in Supplementary section 1. We diagnose convergence and mixing of the described posterior MCMC simulation using trace plots and autocorrelation plots of imputed parameters. For the upcoming simulation examples and the data analysis, we found no evidence of practical convergence problems. An R package of the code used for simulations and the analysis of the ocean microbiome dataset in the following sections is available from the authors website https://users.soe.ucsc.edu/~juheelee/.

### 2.2. Simulation Experiment: Data Generation and Comparative Study

We conducted simulation studies to assess the performance of our model. We compared the model to an alternative model, the negative binomial mixed model (NBMM) in Zhang et al. (2017). We assumed a sample of J = 200 OTUs. We used the same time points (ti) and numbers of replicates (Ki) of our ocean microbiome data as shown in **Figure 1B**. We let β TR <sup>j</sup>,<sup>p</sup> = 0 with probability 0.85. For β TR j,p 6= 0 we simulated <sup>β</sup> TR j,p from either of N(−1.5, 0.05<sup>2</sup> ) or N(1.5, 0.05<sup>2</sup> ) with equal probability, where N(a, b 2 ) denotes the normal distribution with mean a and variance b 2 . It implies that a covariate has no effect on OTU abundance with probability 0.85 or may significantly affect mean abundance with the remaining probability 0.15. To specify r TR t,k and α TR 0,j , we did not assume any distribution and used their classical estimates from our ocean microbiome data; following Witten (2011), we first computed estimates of sample size factors r ′ ti ,k and OTU size factors α ′ 0,j using the ocean microbiome data, r ′ ti ,<sup>k</sup> = yt<sup>i</sup> ,k,·/y··· and α ′ 0,j = 1 N P<sup>n</sup> i=1 PK<sup>i</sup> k=1 yti ,k,j/r ′ ti ,k where yti ,k,· = P<sup>J</sup> j=1 yti ,k,<sup>j</sup> and y··· = P<sup>J</sup> j=1 y··,<sup>j</sup> . We then randomly sampled from the pool of r ′ t,k and α ′ 0,j to specify the true values. To simulate temporal dependence in OTU abundance, we let α˜ TR ti ,<sup>j</sup> = at<sup>i</sup> ,<sup>j</sup> cos(2π(t˜ <sup>i</sup> − bt<sup>i</sup> ,j)) + ct<sup>i</sup> ,j(t˜ <sup>i</sup> − t˜ ⋆ ) 2 . Here t˜ <sup>i</sup> denotes time ti in year and t˜ ⋆ the median of t˜ i . We let at,<sup>j</sup> iid <sup>∼</sup> N(0.15, 0.1<sup>2</sup> ), bt,<sup>j</sup> iid <sup>∼</sup> N(0, 0.5<sup>2</sup> ), and ct,<sup>j</sup> iid <sup>∼</sup> N(0.1, 0.1<sup>2</sup> ) to have different patterns for OTUs. For some OTUs, <sup>α</sup>˜ TR ti ,j are illustrated in red squares in **Figures 4E–G**. We generated s TR j iid ∼ Ga(1, 10). We used the covariate matrix of the ocean microbiome data illustrated in **Figure 2** for the simulation study. For the missing covariates in the data, we generated a value of possible categories with equal probability. We finally simulated OTU counts yt<sup>i</sup> ,k,j from the

negative binomial distribution yt<sup>i</sup> ,k,j indep <sup>∼</sup> NB(<sup>µ</sup> TR ti ,k,j ,s TR j ), where µ TR ti ,k,<sup>j</sup> = r TR ti ,k α TR 0,j exp(α˜ TR ti ,<sup>j</sup> <sup>+</sup> **<sup>X</sup>** TR <sup>t</sup> β TR j ).

For comparison, we used the negative binomial mixed model (NBMM) in Zhang et al. (2017). Similar to the proposed model, the NBMM uses a negative binomial distribution with mean µ NBMM and shape parameter θ NBMM to model OTU counts and assumes log(µ NBMM t,k,j ) = log(yt,k,· )+<sup>β</sup> NBMM 0,<sup>j</sup> <sup>+</sup>**X**t<sup>β</sup> NBMM <sup>j</sup> <sup>+</sup>**Z**t,k**<sup>b</sup>** NBMM j where **X**<sup>t</sup> and **Z**t,<sup>k</sup> are the covariate matrices for fixed effects and random effects, respectively. It assumes random effects **b** NBMM j ∼ N(**0**, 9). By letting the replicates at a time point share the same random effect, OTU abundances in the replicates at a time point are correlated. The NBMM normalizes OTU counts by sample total counts. That is, sample total counts yt,k,· are used as an offset to adjust for the variability in total counts across samples. Similar to other existing methods, the NBMM performs separate analyses of OTUs. An iterative weighted least squares algorithm is developed to produce the MLEs under the NBMM and implemented in a R function glmm in R package BhGLM.

### 2.3. Ocean Microbiome Data: Data Description and Preprocessing

We applied the proposed statistical method to ocean microbiome data. Seawater samples were collected weekly at the end of Santa Cruz Municipal Wharf (SCW), Monterey Bay (36.958 <sup>o</sup>N, <sup>−</sup>122.017 <sup>o</sup>W), with an approximate depth of 10 m. SCW is one of the ocean observing sites in Central and Northern California (CenCOOS), where harmful algal bloom species [HAB species: Alexandrium (Ax), Dinophysis (Dp), Pseudo-nitzschia (Pn) etc.] are monitored weekly along with nutrient measurements [ammonia (NH4), silicate (Si), nitrate (N), phosphate (P)], temperature (T), domoic acid (DA) concentration, and chlorophyll (Chl). Details of phytoplankton net tow sampling of measuring phytoplankton abundance, measurement of physical (nutrients and temperature) and biological parameters (chlorophyll α and DA concentration) are described in Sison-Mangus et al. (2016). Pseudo-nitzschia, Dinophysis, and Alexandrium cells were counted with a Sedgewick rafter counter under the microscope. Data for physical and biological factors are available from the website link http://www.sccoos.org/query/? project=Harmful%20Algal%20Blooms&study[]=Santa%20Cruz %20Wharf. Among the 10 variables, the concentration levels of Alexandrium, Dinophysis, Pseudo-nitzschia, and domoic acid have highly right-skewed distributions and are discretized into categories based on their biological properties for our analysis. The ranges of the concentration levels for the discretization are in Supplementary Table 1 and **Figures 2A–J** illustrates all covariates included for analysis. The values of −1, 0, 1, 2, 3, and 4 represent missing values and the categories of None, Low, Medium, High, and Very High for the discretized variables, respectively. Due to high right skeweness, categories corresponding to high concentration levels have low frequencies. Values of the Dinophysis concentration level are missing at 20 time points among 55 points used for analysis. Sample correlations between the factors are relatively strong. **Figures 2K,L** shows scatterplots for some selected pairs of the factors.

For bacterial RNA samples, three depth-integrated (0, 5, and 10 ft) water samples were collected at a total of 55 time points between April 2014 and November 2015. Two or three samples are sequenced at each time point. The numbers of replicates at the time points are illustrated in **Figure 1B**. Microbial RNA in the samples was extracted for 16S rRNA sequencing. Postprocessing of sequences was performed using the Quantitative Insights Into Microbial Ecology (QIIME 1.9.1) pipeline. A total of nearly 39,823 OTUs were obtained in data after removing singletons. We restricted our attention to OTUs that have greater than or equal to five counts on average. The rule leaves in the end J = 263 OTUs for the 150 samples for the analysis. A heatmap of the counts in the filtered data is shown in **Figure 1**. The primary goal of the study is to investigate the effects of physical and biological factors on abundance levels of OTUs, while accounting for baseline abundance levels of OTUs in samples.

### 3. RESULTS

### 3.1. Simulation Experiment: Model Fitting and Comparison

To fit the proposed model for the simulated data designed in section 2.2, we specified hyperparameter values of the model as follows; for the Laplace prior of <sup>β</sup>j,p, we let <sup>a</sup><sup>λ</sup> <sup>=</sup> <sup>b</sup><sup>λ</sup> <sup>=</sup> 0.5 for a gamma prior of λ 2 j (with mean aλ/b<sup>λ</sup> and variance aλ/b 2 λ ) and a<sup>σ</sup> = b<sup>σ</sup> = 0.3 for an inverse gamma prior for common variances σ 2 j . For the regularizing priors of r˜t<sup>i</sup> ,<sup>k</sup> and <sup>α</sup>˜ 0,<sup>j</sup> , we fixed <sup>d</sup><sup>α</sup> <sup>=</sup> <sup>d</sup><sup>r</sup> <sup>=</sup> 10, <sup>a</sup><sup>r</sup> <sup>=</sup> <sup>b</sup><sup>r</sup> <sup>=</sup> <sup>a</sup><sup>α</sup> <sup>=</sup> <sup>b</sup><sup>α</sup> <sup>=</sup> 1, <sup>ω</sup> 2 <sup>r</sup> <sup>=</sup> <sup>ω</sup> 2 <sup>α</sup> = 1.0, v 2 <sup>r</sup> = 1, and v 2 <sup>α</sup> = 2.0. We also fixed the number of mixture components for the regularizing priors L<sup>r</sup> = 30 and L<sup>α</sup> = 50. To specify values of the mean constraints c<sup>r</sup> and cα, we took an empirical approach. We used the simulated yt<sup>i</sup> ,k,j , computed estimates of rt<sup>i</sup> ,k,<sup>j</sup> and αj,0 as described in section 2.2 and fixed the mean constraints at the means of the logarithm of the estimates, respectively. Note that the specified values of c<sup>r</sup> and c<sup>α</sup> were very different from the means of their true values. For the process convolution prior of OTU-time factor <sup>α</sup>˜t<sup>i</sup> ,j , we chose a value of M such that the kernel function at a basis point is not entirely located in a place where no sample is obtained. We let the number of basis <sup>M</sup> <sup>=</sup> 13 and basis <sup>u</sup>m, <sup>m</sup> <sup>=</sup> 1, ... , <sup>M</sup> evenly spaced between −10 and T<sup>i</sup> + 10. For overdispersion parameter s<sup>j</sup> we let a<sup>s</sup> = 1 and b<sup>s</sup> = 2. To run MCMC simulation, we initialized the parameters by simulating with their prior distributions. We then implemented posterior inference using MCMC simulation over 25,000 iterations, discarding the first 10,000 iterations as burn-in and choosing every other sample as thinning.

**Figure 3** illustrates the comparison of posterior estimates of βj,<sup>p</sup> to their true values β TR j,p for some selected covariates. In the figure, dots and blue dashed lines represent posterior means βˆ j,p of βj,<sup>p</sup> and their 95% credible intervals, respectively. βˆ <sup>j</sup>,ps are around the 45 degree line (red dotted line) for most of (j, p) and most of the interval estimates captures the true values. It implies that the proposed model reasonably recovers β TR j,p . For categories 3 and 4 of X<sup>4</sup> in **Figures 3I,J**, the credible intervals tend to be wider due to their low frequencies in the data as shown

FIGURE 2 | Ocean microbiome data. Bar plots of discretized covariates, concentration levels of *Alexandrium* (Ax, *X*1) and *Dinophysis* (Dp, *X*2), Pseudo-nitzchia (Pn, *<sup>X</sup>*3), domoic acid (DA, *<sup>X</sup>*4) in (A–D). The values of <sup>−</sup>1, 0, 1, 2, 3, and 4 represent a missing value, none and low, medium, high, and highest concentration levels, respectively. Histograms of continuous covariates, concentration levels of ammonia (NH4, *X*5), nitrate (N, *X*6), phosphate (P, *X*7), silicate (Si, *X*8), water temperature (T, *X*9), concentration level of chlorophyll (Chl, *X*10) in (E–J). The variables are standardized to have mean 0 and variance 1 prior to analysis. A scatterplot of the concentrations of ammonia and silicate and a scatterplot of the concentrations of phosphate and silicate are shown in (K,L), respectively.

in **Figure 2D**. The insert plot in each panel illustrates a scatter plot of βˆ <sup>j</sup>,<sup>p</sup> for (j, p) with β TR <sup>j</sup>,<sup>p</sup> = 0 . It shows that the proposed regression model with the Laplace prior effectively shrinks βj,<sup>p</sup> with β TR <sup>j</sup>,<sup>p</sup> = 0 to zero, as is desired in our simulation setup. Supplementary Figure 1 has plots for all covariates.

**Figures 4A–C** illustrate plots of g TR t,k,j vs estimates of gt,k,<sup>j</sup> with their means (black dots) and 95% credible intervals (blue vertical lines) for some selected OTUs, j = 8, 34, and 48. Recall that we do not attempt to recover the true values of individual rt<sup>i</sup> ,k , α0,<sup>j</sup> , and αt,<sup>j</sup> , but we rather aim to reasonably recover the true baseline mean counts, g TR ti ,k,<sup>j</sup> = r TR ti ,k α TR 0,j exp(α˜ TR ti ,j ). In the figure the estimates are tightly around the 45 degree line, providing evidence that reasonable estimates of baseline mean counts are obtained under the proposed model. **Figure 4D** has a histogram of averaged differences between baseline mean count estimates and their true values, D<sup>j</sup> = P<sup>n</sup> i=1 PK<sup>i</sup> k=1 (gˆti ,k,j − g TR ti ,k,j )/N. The averaged differences are around zero, implying that the proposed model provides reasonable estimates of baseline mean counts for most of OTUs. We further examined individual parameters. **Figures 4E–G** shows the comparison of estimates of <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> to their true values over time for the same OTUs in **Figures 4A–C**. Black dots and blue vertical lines represent estimates of posterior means of <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> and their 95% credible intervals, respectively. Red squares represent their true values. From the figure, the estimates of <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> are consistently greater than their true values at all time points, but capture their overall temporal trend. **Figure 4H** illustrates a scatterplot of r˜ TR t,k and their posterior estimates of r˜t,<sup>k</sup> , where dots and blue vertical

FIGURE 3 | Simulation 1—proposed model. Comparison of the true values β TR *j*,*p* and the posterior estimates of β*j*,*p* under the proposed model for some selected covariates. Dots and blue dashed lines represent estimates of posterior means βˆ *j*,*p* of β*j*,*p* and 95% credible intervals (CIs) of β*j*,*p*, respectively. The insert plot in each panel is a scatter plot of βˆ *j*,*p* and β TR *j*,*p* for (*j*, *p*) with β TR *j*,*p* <sup>=</sup> 0.

intervals denote estimates of posterior means and 95% credible intervals, respectively, and the gray horizontal line is atc<sup>r</sup> used for analysis. Different from the estimates of <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> , the estimates of r˜t,<sup>k</sup> fall below the 45 degree line approximately by the same distance for all OTUs. It shows that estimates of <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> and r˜t,<sup>k</sup> have discrepancies from their true values but in the opposite direction and the model can produce reasonable estimates of gt,k,<sup>j</sup> as seen in **Figures 4A–D**. The true overdispersion parameters s TR j are reasonably well estimated as shown in **Figure 4I**. We check the posterior predictive distribution of Yt,k,<sup>j</sup> . The posterior predicted values of Yt<sup>i</sup> ,k,<sup>j</sup> with their 95% predictive intervals for OTUs j = 8, 34, and 48 are compared to their observed values in Supplementary Figure 2. The figure indicates a reasonable model fit.

In addition, we conducted a sensitivity analysis to the specification of mean constraints c<sup>r</sup> and c<sup>α</sup> for the priors of r˜ and <sup>α</sup>˜ <sup>0</sup>. We used different values for <sup>c</sup><sup>r</sup> and <sup>c</sup><sup>α</sup> and compared the estimates of gt,k,<sup>j</sup> to their truth. Supplementary Figures 3a– c has histograms of averaged differences D<sup>j</sup> between gˆt<sup>i</sup> ,k,<sup>j</sup> and g TR ti ,k,j for different specification of c<sup>r</sup> and cα. The histograms show

FIGURE 4 | Simulation 1—proposed model. Panels (A–C) illustrate plots of the true baseline mean counts *g* TR *t*,*k*,*j* vs their estimates *g*ˆ *t*,*k*,*j* for some selected OTUs *<sup>j</sup>* <sup>=</sup> 8, 34, 48. Panel (D) shows a histogram of averaged differences between *<sup>g</sup>* TR *t*,*k*,*j* and *g*ˆ *t*,*k*,*j* for each OTU. Panels (E–G) show plots of estimates of <sup>α</sup>˜0,*<sup>j</sup>* + ˜α*t*,*<sup>j</sup>* over time for OTUs *<sup>j</sup>* <sup>=</sup> 8, 34, 48. Panel (H) has a scatterplot of <sup>ˆ</sup> *r*˜ *t*,*k* vs. *r*˜ TR. Panel (I) has a scatterplot of *<sup>s</sup>*<sup>ˆ</sup> vs. *<sup>s</sup>* TR. Dots represent posterior mean estimates and blue vertical dotted lines 95% credible intervals. Red squares represent the true values.

minor change in estimates of gt<sup>i</sup> ,k,<sup>j</sup> under different specifications of c<sup>r</sup> and cα. An sensitivity analysis to the specification of the number <sup>M</sup> of basis points in the GP convolution model for <sup>α</sup>˜t,<sup>j</sup> was also performed. We used M = 8, 13, and 18 and examined estimates of the baseline mean counts, gt,k,<sup>j</sup> . Supplementary Figures 3a,d,e has histograms of averaged differences D<sup>j</sup> for each of M. The results indicate that the baseline mean counts are reasonably estimated for a range of values of M in the simulation study.

For comparison, we used the NBMM to the simulated data. Since the NBMM does not accommodate missing covariates, we used **X** TR to fit the NBMM. **Figure 5** compares the MLEs βˆNBMM j,p of βj,<sup>p</sup> to the true values for the same covariates used in **Figure 3**. Dots and blue vertical lines represent the MLEs under the NBMM and their 95% confidence intervals, respectively. Comparing **Figure 5** to **Figure 3**, the NBMM produces poor estimates. The MLEs are biased for some covariates (e.g., **Figure 5A**). Also, confidence intervals under the NBMM often

FIGURE 5 | Simulation 1—NBMM. Comparison of the true values β TR *j*,*p* and maximum likelihood estimates <sup>β</sup>ˆNBMM *j*,*p* of β*j*,*p* under the negative binomial mixed model (NBMM) for some selected covariates. Dots and blue dashed lines represent <sup>β</sup>ˆNBMM and their 95% confidence intervals, respectively. The insert plot in each panel is a scatter plot of <sup>β</sup>ˆNBMM *j*,*p* and β TR *j*,*p* for (*j*, *p*) with β TR *j*,*p* <sup>=</sup> 0.

fail to capture the true values and their interval estimates under the NBMM tend to be much wider than those under the proposed model. Normalization through observed sample total counts and inducing correlation in replicates through iid (independent and identically distributed) random effects under the NBMM may lead to poor estimation of the baseline mean abundance for the simulated data, resulting in deterioration of coefficient estimation. In addition, separate analyses of OTUs under the NBMM do not allow to strengthen estimates through combining information across OTUs. Comparing the insert plots in **Figure 5** to those in **Figure 3**, βˆNBMM <sup>j</sup>,<sup>p</sup> with β TR <sup>j</sup>,<sup>p</sup> = 0 tends to more widely spread out from zero and often their confidence intervals fail to capture zero. Supplementary Figure 4 has plots of βj,<sup>p</sup> for all covariates. Supplementary Figures 4, 5 shows the comparison of the estimates θˆNBMM of overdispersion parameters under the NBMM to their true values. Note that θ NBMM is the inverse of s in our model. The NBMM underestimates s<sup>j</sup> for many OTUs, and yields poor predicted values, implying the lack of a fit.

We further examined the performance of the proposed model through additional simulation studies, Simulations 2 and 3 in Lee and Sison-Mangus Bayesian Model for Microbime Data

Supplementary section 2. In these simulations, we studied the robustness of the model when different simulation setups are used to simulate data. In Simulation 2, we assumed no temporal dependence among OTU abundance and generated independent samples from a normal distribution for <sup>α</sup>˜t,<sup>j</sup> . We assumed that all β TR j,p has nonzero effects for all OTUs and simulated βj,<sup>p</sup> from a mixture of normals. The performance of our model is almost the same as in Simulation 1 (see Supplementary Figures 6–8). Interestingly, the NBMM that assumes iid random effects performs poorly for β estimation. In Simulation 3, we simulated α˜ TR t,j from a discontinuous function. The results show that when the temporal dependence pattern is not smooth as assumed for the GP, estimates of the baseline mean counts under the proposed model are slightly deteriorated but the model produces reasonable inference on βj,<sup>p</sup> (see Supplementary Figures 9–11). A more detailed summary of the additional simulations is given in Supplementary section 2.

### 3.2. Ocean Microbiome Data: Model Fitting and Comparison

We specified hyperparameters similar to those in the simulations and analyzed the microbiome data in section 2.3. The MCMC simulation was run over 25,000 iterations. The first 15,000 iterations were discarded as burn-in and every other sample was kept as thinning and used for inference. **Figure 6** illustrates inference on covariate effects for some selected OTUs, j = 16, 34, and 49, taxonomically belonging to Alteromonadales, Halomonas sp., and Alteromonadales in the Gamma-proteobacteria phyla, respectively. Dots and vertical solid lines represent the posterior mean βˆ <sup>j</sup>,<sup>p</sup> and 95% credible interval estimates, respectively. Similar to the results of the simulation study, the credible intervals for high and highest levels of the discretized covariates tend to be wider due to their low frequencies in the data. From **Figure 6A**, on average the medium concentration level of domoic acid (DA, X4) and the concentration level of nitrate (N, X6) significantly decrease the mean abundance of OTU 16 by a multiplicative factor of exp(−0.572) = 0.564 and exp(−0.260) = 0.771, respectively. One may infer that the medium concentration level of domoic acid is significantly associated with lower expected counts for the OTU compared to those with category none of the domoic acid concentration level. A similar argument can be applied to the inference on the nitrate concentration level. Interestingly, we observed statistically significant reduction in abundance from many OTUs belonging to Gammaproteobacteria including those OTUs for increasing domoic acid concentration (not shown). The resulting inference was further validated through a lab experiment. Most notably, four bacterial cultured isolates belonging to Gamma-proteobacteria (three among them are Alteromonadales) were observed to be severely retarded in growth after 2 days of exposure to increasing domoic acid of 0 to 150 µg/ml in the experiment (Sison-Mangus, unpublished data). This demonstrates that the proposed model successfully identifies important OTUs in ocean bacterial community dynamics for further investigation. More results are presented in Supplementary section 3. Supplementary Figures 12a–c illustrates the posterior estimates of baseline expected counts <sup>α</sup>˜ 0,<sup>j</sup> + ˜αt,<sup>j</sup> normalized by sample size factors for the OTUs. From the figure, the baseline expected counts vary over time for those OTUs and the temporal pattern of OTU j = 34 is different from those of OTUs j = 16 and 49. Histograms of the posterior mean estimates βˆ <sup>j</sup>,<sup>p</sup> of βj,p, are illustrated in Supplementary Figure 13. The figure does not show clear overall tendency in the direction of association between covariates and OTU counts. Posterior inference on sample size factors rt<sup>i</sup> ,k and OTU specific overdispersion parameters s<sup>j</sup> is illustrated in Supplementary Figures 12d,e.

For comparison, we fitted the NBMM to the data. Since the NBMM does not account for missing values, we use the maximum a posteriori estimates of the missing values under the proposed for the NBMM. We used the R function glmm and the algorithm produced warning messages on convergence for 32 OTUs. **Figure 7** illustrates βˆNBMM j,p (dots) with their 95% confidence intervals (solid vertical lines) for OTUs j = 16, 34, and 49. Inference on the covariate effects is different from that under the proposed. For example, domoic acid (DA) levels do not have significant effects on the mean counts for OTU j = 16 and 49 from **Figures 7A,C**. Comparing **Figures 7A,C** to **Figure 7**, the NBMM produces wider interval estimates for βj,p. Histograms of the MLEs of βj,p, βˆNBMM j,p under the NBMM are shown in Supplementary Figure 14. The histograms are much dispersed than those under the proposed model shown in Supplementary Figure 13. Estimates βˆ <sup>j</sup>,<sup>p</sup> and βˆNBMM j,p for all covariates are also compared in Supplementary Figure 15. From the figure, the NBMM yields extremely large or small values for βˆ <sup>j</sup>,<sup>p</sup> for some OTUs, possibly due to the convergence problem. The insert plots show that for regions of small values of βˆ <sup>j</sup>,p, the estimates under the proposed are more shrunken toward zero than those under the NBMM, similar to the results in section 3.1. The overdispersion parameter estimates under the NBMM tend to be smaller than those under the proposed (shown in Supplementary Figure 12f), which may lead to different predictive distributions of OTU counts.

### 4. DISCUSSION AND CONCLUSIONS

In this paper, we developed a Bayesian semiparametric regression model for joint analysis of microbiome data. We formulated the mean counts of OTUs as a product of factors and built models for the factors. We utilized the regularizing priors with mean constraints to avoid possible idenfiability issues, and the process convolution model to capture the temporal dependence structure in the baseline mean abundance of an OTU. The flexible model developed for baseline abundance enables joint analysis of all OTUs in the data and allows borrowing information across OTUs, across samples, and across time points. The model produces accurate estimates of the baseline mean counts and yields improved estimates of the effects of the covariates. We incorporated the Laplace distribution, a sparsity inducing shrinkage prior for the coefficients and the proposed model produces sparse estimates that is more desirable when the problem is high-dimensional and covariates are highly

correlated. We compared the proposed model to a comparable frequentist model that does separate analyses for individual OTUs. The comparisons through the simulation study and real data analysis show better performance of the proposed model.

Although we focused on the analysis of NGS count data, the proposed model is quite general and can be applied for analyses of any count data. Future work will explore alternative approaches to model the effects of covariates on the mean counts. For example, one may consider a nonparametric model using linear combinations of basis functions (Kohn et al., 2001) to flexibly capture shape in the response function. In such a case, an elaborate development of the prior model may be needed to produce a robust inference since both the baseline mean counts and the covariate effects are nonparametrically modeled. Other possible extensions are to include a variable selection method such as a stochastic search variable selection (George and McCulloch, 1993) if it is reasonable to assume that some covariate effects are exactly zero, and to let coefficients vary over time if covariate effects evolve with time. For time varying coefficients, we may use the random walk process in Leybourne (1993) to

induce relationship between βj,p,t−<sup>1</sup> and βj,p,<sup>t</sup> . Considering the high dimensionality in OTU data, posterior computation may need to be carefully handled. Also, prior information may be needed to produce sensible inference due to sparsity in data.

led the collaboration with MS-M for statistical analysis. MS-M provided the ocean microbiome data, participated the statistical model development, provided biological interpretation of the resulting inference and edited the manuscript.

### AUTHOR CONTRIBUTIONS

JL developed the statistical model and conducted simulation studies and data analysis. She also prepared the first draft and

### FUNDING

This work was supported by NSF grant DMS-1662427 (JL) and NOAA-ECOHAB PROGRAM (Grant No. NA17NOS4780183, ECOHAB #905) (MS-M and JL). Collection of environmental data from the Santa Cruz Municipal Wharf was supported by Cal-PReEMPT with funding from the NOAA-MERHAB (#206), ECOHAB and IOOS programs.

### ACKNOWLEDGMENTS

We gratefully acknowledge Raphael Kudela who provided the environmental data in this study, Michael

### REFERENCES

Kempnich and Sanjin Mehic for doing the water sampling and processing the water samples for DNA extraction.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.00522/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Lee and Sison-Mangus. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# SplinectomeR Enables Group Comparisons in Longitudinal Microbiome Studies

#### Robin R. Shields-Cutler<sup>1</sup> , Gabe A. Al-Ghalith<sup>2</sup> , Moran Yassour3,4 and Dan Knights1,5 \*

<sup>1</sup> BioTechnology Institute, College of Biological Sciences, University of Minnesota, Minneapolis, MN, United States, <sup>2</sup> Bioinformatics and Computational Biology, University of Minnesota, Minneapolis, MN, United States, <sup>3</sup> Broad Institute of Massachusetts Institute of Technology, Harvard University, Cambridge, MA, United States, <sup>4</sup> Center for Computational and Integrative Biology, Massachusetts General Hospital and Harvard Medical School, Boston, MA, United States, <sup>5</sup> Department of Computer Science and Engineering, University of Minnesota, Minneapolis, MN, United States

Longitudinal, prospective studies often rely on multi-omics approaches, wherein various specimens are analyzed for genomic, metabolomic, and/or transcriptomic profiles. In practice, longitudinal studies in humans and other animals routinely suffer from subject dropout, irregular sampling, and biological variation that may not be normally distributed. As a result, testing hypotheses about observations over time can be statistically challenging without performing transformations and dramatic simplifications to the dataset, causing a loss of longitudinal power in the process. Here, we introduce splinectomeR, an R package that uses smoothing splines to summarize data for straightforward hypothesis testing in longitudinal studies. The package is open-source, and can be used interactively within R or run from the command line as a standalone tool. We present a novel in-depth analysis of a published large-scale microbiome study as an example of its utility in straightforward testing of key hypotheses. We expect that splinectomeR will be a useful tool for hypothesis testing in longitudinal microbiome studies.

Keywords: bioinformatics, microbiome analysis, R packages, computational biology methods, permutation tests, longitudinal data analysis

## INTRODUCTION

Biological studies in humans are subject to significant variability and noise, often great enough to obscure all but the most dramatic differences. Longitudinal studies are powerful in these cases, allowing researchers to observe (and account for) both inter- and intra-individual variability, or measure changes in response to an intervention in real time (Gonzalez et al., 2012; Gerber, 2014). As the costs of DNA sequencing have decreased, microbiome researchers have a greater opportunity to perform such longitudinal studies. While longitudinal data with multiple timepoints always provide more information than single-timepoint data, the computational tools to analyze longitudinal microbiome studies with multiple timepoints per subject lag behind. A number of practical concerns often complicate analysis of longitudinal microbiome data: time points are usually not in sync or differ in number between subjects, longitudinal variation may not follow a normal distribution, and timeseries data may follow arbitrary curves, for example during the maturation of the infant microbiome. To overcome these challenges, researchers in many studies have collapsed samples across time points to average individuals' signals or they have

#### Edited by:

Jessica Galloway-Pena, University of Texas MD Anderson Cancer Center, United States

#### Reviewed by:

Ekaterina Smirnova, University of Montana, United States Richard Allen White III, RAW Molecular Systems (RMS) LLC, United States

> \*Correspondence: Dan Knights

#<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Systems Microbiology, a section of the journal Frontiers in Microbiology

Received: 15 November 2017 Accepted: 06 April 2018 Published: 23 April 2018

#### Citation:

Shields-Cutler RR, Al -Ghalith GA, Yassour M and Knights D (2018) SplinectomeR Enables Group Comparisons in Longitudinal Microbiome Studies. Front. Microbiol. 9:785. doi: 10.3389/fmicb.2018.00785

summarized first with multivariate approaches that condense the initial observations (e.g., David et al., 2014; Zhou et al., 2015; Yassour et al., 2016). These approaches have been sufficient to make important discoveries in published studies, but there may still be opportunities to gain statistical power by using additional information content and directionality of the temporal axis.

To address this gap in data analysis, we introduce splinectomeR for direct hypothesis testing of categorical variables in longitudinal studies. SplinectomeR's implementation is straightforward and complements recently developed mixed-effects models that are used for discovering differentiating taxa (Chen and Li, 2016). At the core of the tests is the loess spline that uses weighted local polynomials to model data that may not follow any classical model or shape (as is common in real biological data) (Cleveland, 1979; Cleveland and Devlin, 1988). Null distributions are generated by permutation of the data, similar to methods implemented in multivariate tests such as PERMANOVA (Anderson, 2001). Lastly, its implementation as an R package makes it practical and easy to adopt by microbiome researchers.

splinectomeR contains three key functions: permuspliner, which tests for an overall significant difference between two groups' summary splines over the longitudinal time course; sliding\_spliner, which interpolates across the group spline and tests for significance between two groups at each interval to illuminate regions of time where significant differences exist; and trendyspliner, which tests for a significant non-zero overall trend in a single population over time. To demonstrate the utility of these tests, we performed additional testing of key hypotheses in a large published cohort of 37 infant microbiomes sampled over the first 3 years of life (Yassour et al., 2016). SplinectomeR is open-source and freely available for download and installation on GitHub at https://github.com/RRShieldsCutler/splinectomeR.

### METHODS AND IMPLEMENTATION

splinectomeR contains three primary functions that test specific hypotheses about longitudinal trends (see **Figure 1** for schematic diagram). Each function uses loess splines to smooth longitudinal data before performing the specific statistical test. The input is a properly formatted data frame: each quantitative measurement or metadata category is its own column, and each row is a separate observation. This is a common structure for bioinformatics metadata, including in microbiome analyses, and therefore tests may be performed with little or no reformatting required. The standalone command-line version of the scripts requires a tab-delimited file of the same structure. We refer the reader to the online package vignettes (included as Supplemental Data Files 1 and 2) for examples of proper input data and for manual data wrangling techniques.

### Permuspliner

The objective of the permuspliner function is to test whether two groups of individuals follow more different trajectories over time (or along any continuous axis) than would be expected by random chance. It compares two groups over time without collapsing the timepoints to a single averaged point. Since differences between groups may not be consistent over time, and responses may even invert relative to one another over time, collapsing could mask or nullify important differences between groups. This problem can be avoided by considering the entire time series in a statistical comparison. Permuspliner fits a spline to the average time series of each group of individuals, and then measures the absolute area between the splines to determine the observed group difference.

Specifically, the input data is first split into the two groups to be compared, A and B. In some studies, participants may only appear once or twice in the dataset; to filter out samples from low-prevalence participants, the user may set a threshold with "cut\_low = n." A loess spline is fit to each group's total time series, if it also meets the minimum data sparseness threshold set by the "cut\_sparse = n" parameter. Lastly, the direction of the null hypothesis can be set with "test\_direction = [more/less]"; this enables the user to test for not only group differences that are greater than expected by chance (test\_direction = "more," the default) but also differences smaller than expected by chance (test\_direction = "less"). To calculate the observed group distance, points are interpolated along each spline and the sum of the areas across these points yields the group distance. Splinectomer uses 1000 points by default but this can be edited in the function arguments (with "ints = n") to account for a longer time series or greater/lower resolution. This area determination is then recalculated after permuting the group assignments: to which group each participant belongs (A vs. B) is shuffled randomly without replacement. Permuting the group labels, as opposed to the participant data points, preserves the underlying distributions and patterns of the individuals' timeseries curves. The permutation is repeated (default setting is 999 permutations) to generate a null distribution over the random between-group distances, from which an empirical p-value is calculated and reported by comparison to the observed distance. Because the random distribution is built from the observed values, it inherently reflects the noise and variability of the dataset, and is therefore tailored to each study's unique character. SplinectomeR also includes plotting functions for visualizing the permuted distribution.

### Sliding Spliner

As a complement to the permuspliner test, the sliding spliner function allows the user to test the data series at defined intervals and ask whether the two groups of interest are significantly different at a any point in time during the time series. This is often challenging in clinical datasets where sampling is not coordinated between individuals, so analysis requires artificial binning or averaging across time in order to compare enough data at any one time point. Here, each individual—as opposed to the whole group, as in the permuspliner test above—is summarized by a spline, thus filling gaps across their time series. Every interpolated interval can then be tested as a complete distribution of all participants without binning or unnecessarily dropping samples. Low-prevalence participants are again filtered out with the "cut\_low" parameter. Since the splines are not extrapolated beyond the start and end of the individual's time

series, the start and end regions may become less dense as fewer subjects have early and/or late samples. To account for this, the "test\_density = n" argument sets the minimum number of participants required in each group to perform and report a p-value (default = 3). The output from this function is a table containing a p-value for each time interval, summarized by the companion plotting function where intervals at which the groups significantly diverge are visualized, and the points are scaled according to the data density—larger circles on the plot mean more data were used to calculate the p-value at that interval than for smaller circles.

### Trendyspliner

While testing for association between a categorical variable and a longitudinal variable is straightforward using regressions and correlations, methods for quickly testing whether a response is increasing or decreasing over time are less established. The trendyspliner function tests whether a set of responses in one group changes in a non-zero direction over the time series (or other continuous independent axis). A spline is fit to the data and interpolated across the number of set intervals. Non-zero change is measured as the area between the group spline and a linear baseline that is established from the start point of the group spline. Thus, if the pattern of observations does not meaningfully increase or decrease over time, the spline will not diverge from the baseline and the areas will remain small. To generate the null distribution, the time series within each individual is then permuted, and the spline is recalculated along with the area to the permuted baseline. The permutation is repeatedly executed to generate a random distribution of areas from which the two-sided p-value is calculated by comparison to the observed value. In some biological measurements, individuals' initial values may be variable (e.g., body weight, height). To normalize these differences and improve the ability to detect an increase or decrease in these situations, the user can elect to "mean\_center" the observations before calculations, which shifts each individuals' mean over time to the group mean from all individuals. As in the other modes, a plotting function allows the user to visualize the permuted splines in the context of the real data.

### User Customization

In each of the above functions, the user can alter several additional parameters for specific applications. These include the spline span parameter, a standard spline parameter that determines how large the local smoothing neighborhood is, and

therefore the degree of smoothing; the number of permutations, which influences the sensitivity of the test and p-value, as more permutations will allow a lower p-value to be detected but will also increase run time (default is 999 permutations); and the number of intervals over which to divide the data (larger values also increase run time and memory but provide finer resolution and more precise comparisons). These arguments allow flexibility for more advanced users with particular needs and unique data shapes.

### Package Implementation

The splinectomeR code was written in R version 3.4.0, and the package was built in RStudio using devtools and roxygen2 to generate and populate the package documentation (Wickham and Chang, 2017; Wickham et al., 2017). SplinectomeR is opensource and freely available on GitHub at https://github.com/ RRShieldsCutler/splinectomeR. The figures generated by the secondary plotting functions are ggplot2 objects and can be saved in most image formats at a size and resolution specified by the user (e.g., through base R drivers or the "ggsave" function) (Wickham, 2009).

### EXAMPLE ANALYSIS

To demonstrate how splinectomeR detects group changes over a time series, we first generated a simulated data set. The response variable was perturbed at one or two regions of the time series, and at three magnitudes of change (**Figures 1B,C**). A linear model is a poor fit to these data shapes, and testing for absolute change from beginning to end is not sensitive to the internal dynamics. However, splinectomeR uses the entire dataset to compare between the baseline and perturbed data. As the magnitude of change increases, the permuspliner test yields decreasing p-values as expected in both the single region (p = 0.3, p = 0.005, p = 0.001 for 1x, 2x, and 4x shifts, respectively; **Figure 1B**) and double region perturbations (p = 0.79, p = 0.36, p = 0.002 for 1x, 2x, and 4x shifts, respectively; **Figure 1C**). Non-linear changes of this sort during a time series may be of great biological interest, although their statistical significance is difficult to test with existing tools.

Freely available datasets were used to test and further demonstrate the splinectomeR functions, as documented in the package vignettes that are available to view online in HTML format at https://rrshieldscutler.github.io/splinectomeR/. A proof-of-concept analysis was performed on the ChickWeights dataset in the R "datasets" package, and is available as a vignette on the website. To demonstrate splinectomeR's utility on a more complex dataset, we used the publicly available OTU tables and associated metadata from a published longitudinal study of infant microbiomes by Yassour et al. (2016).

### Analysis of Longitudinal Microbiome Data

We tested splinectomeR's utility on the taxonomic and metadata profiles from Yassour et al. (2016), to evaluate whether we could statistically support patterns described by the authors and investigate novel hypotheses. A more extensive analysis including all code used for data formatting is available as an online vignette and as Supplemental Data File 2.

In several figures in the original publication, the authors draw visual comparisons between taxon abundances in antibiotic exposed versus non-exposed infants. These stream plots are powerful in displaying the inter-individual diversity and inspired us to use splinectomeR to perform statistical hypothesis testing on the time series data for differences between the infant groups. SplinectomeR's permuspliner function can test whether a taxon's abundance pattern over time is significantly different between antibiotic-exposed and non-exposed infants. We used splinectomeR to calculate that the difference in Bacteroidaceae abundance is not statistically significant across the overall time series (p = 0.28). The permutated differences support this conclusion, as the output shows in **Figure 2A**. However, the results indicate that the groups may be diverging toward the end of the time series, as the observed distance is higher relative to most of the permuted values. We tested this with the sliding spliner function, which generates a series of p-values across the longitudinal scale. The results, as the function's output plot shows in **Figure 2B**, indeed show that there is a temporal pocket of significance surrounding the 30-month time point. We were also able to confirm the finding that the genus Bacteroides is significantly different between vaginal and cesarean born infants (p = 0.04), and that this difference is most pronounced in the first year of life (see Supplemental Data File 2).

Because these tests are implemented as R functions, they can be used programmatically for multiple hypothesis testing, such as testing all dominant bacterial families for significant differences between antibiotic exposure status. We performed this test on the present dataset (see vignette for full analysis), revealing that Porphyromonadaceae is the most discriminatory family (p = 0.05). The built-in plotting function (permuspliner.plot.permsplines) shows that antibiotic-exposed infants do indeed develop a notably higher abundance of Porphyromonadaceae (**Figure 2C**). The permuted splines lie predominantly between the two observed curves, supporting the conclusion that this difference is larger than expected by chance, and this observation becomes greater over time. The distance plot further supports this conclusion (**Figure 2D**), showing that initial distances are not greater than chance but become significant after approximately 12 months of age.

In their published analysis, the authors found a significant difference in abundance of the butyrate-producing Clostridium groups IV and XIV by antibiotic status at the final time point (36 months). Given butyrate's important roles in gut homeostasis (e.g., Pryde et al., 2002; Ridaura et al., 2013; Zhang et al., 2016), it is worth investigating whether this difference exists over the infants' first 3 years of development, or is just established at 36 months. To test this, we used the permuspliner function on a summarized table containing the following putatively butyrate-producing genera from Clostridium groups IV and XIV present in the OTU table: Clostridium, Coprococcus, Dorea, Lachnospira, Roseburia, Ruminococcus, and Faecalibacterium (Louis and Flint, 2009; Lopetuso et al., 2013; Van den Abbeele et al., 2013). When all of the longitudinal

FIGURE 2 | Permuted spline tests for statistical significance in longitudinal microbiome data. (A) The permuspliner output plot shows that the difference between Bacteroidaceae abundance between antibiotic exposed and non-exposed infants (distance between group splines shown as red line) is not significantly greater than 95% of the permuted values shown in translucent black. (B) The plot output of the sliding spliner function shows the p-value at each specified interval (shown with default 100 intervals) derived from the distribution of points from individuals' smoothed splines. Dotted line indicates p = 0.05. The number of infants with data at a given interval is used to scale the data point size, as some entered and exited the study later or earlier, respectively. (C) Porphyromonadaceae abundance over time distinguishes antibiotic exposed (group spline in blue) and non-exposed infants (group spline in red; 999 permutations, p = 0.053). (D) Group distance plot, as in (A), for the Porphyromonadaceae comparison, showing that permuted splines support a trend toward a greater true statistical difference after 12 months of life.

data are included in this comparison, we find that antibiotic exposed infants do indeed have significantly lower Clostridium group IV/XIV abundance compared to non-exposed infants over time (**Figure 3A**). Notably, our test using the entire 36-month time series yields a lower p-value than that reported using just the 36 month data: p = 0.003 vs. p = 0.037, respectively, and the resulting plot suggests that the most divergent time point is near 1 year of age (**Figures 3A,B**). This demonstrates the utility and simplicity of the splinectomeR tests; we are able to directly include 3 years of observations, strengthen support for this finding, and suggest directions for new hypotheses.

To test and demonstrate the third function in the splinectomeR package, we hypothesized that the infants' alpha diversity would significantly increase over time. We can use the permuspliner test to show that the alpha diversity is highly similar between the antibiotic (p = 0.92) and birth mode (p = 0.98) groups, but to statistically test whether it changes over time, we used the trendyspliner function. As shown in **Figure 4**, the randomly permuted splines generate a linear and zero-slope distribution, while the observed spline increases steadily over the time series, confirming our hypothesis that infant microbiome alpha diversity increases over time (p = 0.01). In this case, the trend is evident and expected; many biological datasets, however, have slight trends that are hard to interpret, in which case the trendyspliner test provides a straightforward permutation-based statistical test to confirm whether the deviation is greater than expected from chance.

### Summary

Yassour et al. (2016) present a prodigious dataset with dense longitudinal data that details the human gut microbiome's complex dynamics over the first 3 years of life. The approaches presented above provide statistical support for observations and conclusions the authors reported, and allow us to test and develop additional hypotheses from the dataset. Researchers analyzing new longitudinal microbiome data with multiple samples per subject may benefit from these analytical capabilities provided by splinectomeR.

### DISCUSSION

Longitudinal studies hold great promise for understanding the effects of interventions and environmental stimuli in the context of a naturally variable population. Analysis of these complex data has been impeded by a lack of clear, simple methods for directly comparing observations across multiple individuals and many time points without averaging or summarizing across time points. As we have demonstrated here using a largescale longitudinal microbiome study, the splinectomeR package performs straightforward tests that are easy to interpret and will allow researchers to test important hypotheses from within R or a command line interface.

The approaches here are not without limitations; reliance on the loess spline means that the tests may be impacted by outliers, particularly in sparse datasets. User-definable arguments for sparseness cut offs and spline resolution (the smoothing parameter or spar) can help minimize these effects. Additionally, the size of longitudinal studies mean that the tests with many permutations can be slow to complete, though still easily performed on a standard workstation running R. From a standard metadata table, splinectomeR tests are run with a single command in R or on the command line. We provide the user with interpretable results in the form of pre-formatted plots that can be saved at publication quality and re-generated from the results object stored by the function.

### CONCLUSION

In summary, we have shown how a new open-source R package, splinectomeR, can quickly assess statistical significance

in large longitudinal microbiome studies by summarizing longitudinal group data with splines and using randomly permuted distributions to evaluate the probability that the observed magnitude of differences between groups, or of trends over time, is due to chance. By providing three distinct types of hypothesis tests, we can explore overall changes in abundances or other biological measures, and compare longitudinal trends between groups of interest. Altogether, we are able to perform these tests in a way that takes full advantage of longitudinal data and maintains individual observations, thus leveraging all possible data points. Longitudinal studies generating "big data" with multi-omics approaches are now commonplace, but we lack appropriate tools to interpret these data. We offer splinectomeR as an open-source solution to testing key hypotheses in complex longitudinal data. SplinectomeR may also simplify analysis for longitudinal studies beyond the microbiome research field.

### AUTHOR CONTRIBUTIONS

RS-C, GA-G, and DK conceived and designed the spline tests. RS-C and GA-G wrote the R code, and RS-C built the R package and conceived, and performed the example experiments. RS-C

### REFERENCES

and MY analyzed the experiments. RS-C and DK wrote the manuscript.

### FUNDING

This research was supported, in part, by NIH Grant R01AI121383 (to DK) and Minnesota Partnership for Biotechnology and Genomics Grant #13.28 (to DK).

### ACKNOWLEDGMENTS

We are grateful to all members of the Knights lab for testing preliminary code and providing feedback on various aspects of the project.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.00785/full#supplementary-material

**Conflict of Interest Statement:** DK serves as CEO and holds equity in CoreBiome, a company involved in the commercialization of microbiome analysis. The University of Minnesota also has financial interests in CoreBiome under the terms of a license agreement with CoreBiome. These interests have been reviewed and managed by the University of Minnesota in accordance with its Conflict-of-Interest policies.

The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Shields-Cutler, Al-Ghalith, Yassour and Knights. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# On the Use of Diversity Measures in Longitudinal Sequencing Studies of Microbial Communities

Brandie D. Wagner 1,2 \*, Gary K. Grunwald<sup>1</sup> , Gary O. Zerbe<sup>1</sup> , Susan K. Mikulich-Gilbertson<sup>3</sup> , Charles E. Robertson<sup>4</sup> , Edith T. Zemanick <sup>2</sup> and J. Kirk Harris <sup>2</sup>

*<sup>1</sup> Department of Biostatistics and Informatics, Colorado School of Public Health, University of Colorado, Anschutz Medical Campus, Aurora, CO, United States, <sup>2</sup> Department of Pediatrics, School of Medicine, University of Colorado, Anschutz Medical Campus, Aurora, CO, United States, <sup>3</sup> Department of Psychiatry, School of Medicine, University of Colorado, Anschutz Medical Campus, Aurora, CO, United States, <sup>4</sup> Department of Molecular, Cellular and Developmental Biology, School of Medicine, University of Colorado, Anschutz Medical Campus, Aurora, CO, United States*

#### Edited by:

*Jessica Galloway-Pena, University of Texas MD Anderson Cancer Center, United States*

#### Reviewed by:

*Bradley Stevenson, University of Oklahoma, United States Richard R. Rodrigues, Oregon State University, United States*

> \*Correspondence: *Brandie D. Wagner<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Systems Microbiology, a section of the journal Frontiers in Microbiology*

Received: *23 January 2018* Accepted: *01 May 2018* Published: *22 May 2018*

#### Citation:

*Wagner BD, Grunwald GK, Zerbe GO, Mikulich-Gilbertson SK, Robertson CE, Zemanick ET and Harris JK (2018) On the Use of Diversity Measures in Longitudinal Sequencing Studies of Microbial Communities. Front. Microbiol. 9:1037. doi: 10.3389/fmicb.2018.01037* Identification of the majority of organisms present in human-associated microbial communities is feasible with the advent of high throughput sequencing technology. As substantial variability in microbiota communities is seen across subjects, the use of longitudinal study designs is important to better understand variation of the microbiome within individual subjects. Complex study designs with longitudinal sample collection require analytic approaches to account for this additional source of variability. A common approach to assessing community changes is to evaluate the change in alpha diversity (the variety and abundance of organisms in a community) over time. However, there are several commonly used alpha diversity measures and the use of different measures can result in different estimates of magnitude of change and different inferences. It has recently been proposed that diversity profile curves are useful for clarifying these differences, and may provide a more complete picture of the community structure. However, it is unclear how to utilize these curves when interest is in evaluating changes in community structure over time. We propose the use of a bi-exponential function in a longitudinal model that accounts for repeated measures on each subject to compare diversity profiles over time. Furthermore, it is possible that no change in alpha diversity (single community/sample) may be observed despite the presence of a highly divergent community composition. Thus, it is also important to use a beta diversity measure (similarity between multiple communities/samples) that captures changes in community composition. Ecological methods developed to evaluate temporal turnover have currently only been applied to investigate changes of a single community over time. We illustrate the extension of this approach to multiple communities of interest (i.e., subjects) by modeling the beta diversity measure over time. With this approach, a rate of change in community composition is estimated. There is a need for the extension and development of analytic methods for longitudinal microbiota studies. In this paper, we discuss different approaches to model alpha and beta diversity indices in longitudinal microbiota studies and provide both a review of current approaches and a proposal for new methods.

Keywords: microbiome, Hill's numbers, repeated measures, alpha diversity, beta diversity, Shannon index, mixed model

## INTRODUCTION

Identification of the majority of organisms present in humanassociated microbial communities is now feasible with the advent of high throughput sequencing technology. Several studies have shown large subject-to-subject variability (Flores et al., 2014) as well as many different factors that might contribute to variability in microbiome studies, i.e., diet, region, exposure, genetics, etc. Given the highly personalized microbiome, valuable information is likely to come from studies following subjects over time. The use of longitudinal study designs is important to better understand the contribution of the microbiome to human health (Flores et al., 2014). Complex study designs with longitudinal sample collection require analytic approaches to account for this additional source of variability, and to allow examination of changes within subjects.

Extension and development of analytic methods are needed for longitudinal microbiota studies (Gerber, 2014). Current approaches include extending models applied to individual taxa to address repeated measures over time (Chen and Li, 2016; Fang et al., 2016; Wagner et al., 2017) but not much attention has been given to discussion and extension of ecological community indices, which are useful for describing the community biodiversity. The majority of analytic methods for these measures were developed for studying one community over time in the field of ecology. This paper, therefore, focuses on the application and development of methods for diversity indices in order to model multiple communities (i.e., subjects) over time.