Multi-omic approaches for host-microbiome data integration

ABSTRACT

The gut microbiome interacts with the host through complex networks that affect physiology and health outcomes. It is becoming clear that these interactions can be measured across many different omics layers, including the genome, transcriptome, epigenome, metabolome, and proteome, among others. Multi-omic studies of the microbiome can provide insight into the mechanisms underlying host-microbe interactions. As more omics layers are considered, increasingly sophisticated statistical methods are required to integrate them. In this review, we provide an overview of approaches currently used to characterize multi-omic interactions between host and microbiome data. While a large number of studies have generated a deeper understanding of host-microbiome interactions, there is still a need for standardization across approaches. Furthermore, microbiome studies would also benefit from the collection and curation of large, publicly available multi-omics datasets.

KEYWORDS:

• Multiomics
• microbiome
• network
• disease
• analysis
• inference
• host-microbiome interactions

Introduction

The human gut microbiome comprises myriad bacteria, viruses, fungi, and protozoa that interact with each other as well as their host, ultimately affecting host physiology and disease. The microbiota contribute to various functions important for health, including digestion, immune regulation, maintenance of the intestinal mucosal barrier, and protection from pathogens.¹ Alterations in the microbiota have been associated with health conditions including colorectal cancer,² inflammatory bowel disease,³ obesity,⁴ and depression.⁵

Despite its role in health, it is only with recent advances in sequencing technology that researchers have been able to gain a holistic view of the composition of the gut microbiome. Sequencing of the gut metagenome has allowed for the discovery of species that were previously difficult to culture in vitro.⁶ Although progress has been made in culturing bacteria previously thought to be “unculturable”,⁷ a large number of taxa still have not been cultured,⁸ and next-generation sequencing remains the standard for microbiome profiling.

Analysis of the microbiome has revealed a complex web of interactions between microorganisms and their host. For example, microbes interact with one another in networks that vary spatially and with disease states.⁹ Furthermore, analyses using metabolomics and metagenomics have demonstrated the ability of microbes to metabolize drugs,¹⁰ and joint analyses of the microbiome and host transcriptome have also shown that the microbiota can regulate host gene expression.¹¹ Additional studies have identified associations between host gene expression and gut microbiome composition across disease states, underscoring the relevance of multi-omics analyses of the microbiome.^12–14

Current approaches to microbiome data analysis

Current approaches to microbiome community profiling usually involve either shotgun metagenomic or 16S rRNA amplicon sequencing.¹⁵ In shotgun sequencing, DNA is first extracted from all of the cells in a sample, cleaved into small fragments, and sequenced. Computational methods are then used to align the reads against reference genomes or marker genes to infer the abundances of taxa present in the sample.¹⁵ Shotgun sequencing may also be performed with higher read depths in order to assemble bacterial genomes de novo. In 16S rRNA amplicon sequencing, only a segment of the 16S rRNA gene of bacterial genomes in a sample is amplified and sequenced. The 16S rRNA gene contains both conserved and variable regions. The conserved regions are used as a target for PCR primers, while the variable regions are used to determine the identities and abundances of taxa present in the sample.

While microbiome research has largely focused on bacteria, which comprise more than 99% of the genetic material in the microbiome, viral, fungal, and archaeal microbial components are also of interest.¹⁶ The virome, for example, includes viruses that infect host cells and bacteriophages that target bacteria. Alterations in the virome have been associated with human diseases such as SARS-CoV-2, inflammatory bowel disease, diabetes, and colorectal cancer.¹⁷ Because viruses lack universally conserved sequences such as the 16S region, viral sequencing is usually performed using shotgun or long-read sequencing. Viral genetic material is usually present in lower concentrations in samples, which may motivate the use of higher read depths when sequencing.¹⁷ The fungal microbiota, as with viruses, are affected by host factors such as diet and antibiotic use.¹⁶ Fungi have been shown to exhibit complex interactions with the host and with other members of the microbiome.¹⁸ For example, there is evidence that fungi and bacteria together contribute to early-life immune education, and Saccharomyces boulardii has been used as a probiotic to treat bacterial infections.¹⁸ Fungi also play a role in host health and immunity. Candida albicans, for example, has been shown to modulate Th17 helper T cells.¹⁶ Candida species have also been associated with health conditions such as inflammatory bowel disease and liver disease.¹⁶ Similar to bacteria, fungal sequencing can also be performed either with shotgun sequencing or through amplicon sequencing of the 18S, ITS1 and ITS2 regions.¹⁶ However, reference ITS sequences for many fungal genomes are currently lacking, limiting the ability of sequencing projects to study the fungal microbiota, and motivating the development of alternate classification strategies.¹⁹ In contrast with other components of the microbiome, archaea represent members of the microbial community that have thus far received little attention, partially due to the lack of common archaeal pathogens. While there is evidence that archaea interact with other microbial components and form stable enterotypes, there is currently little evidence of associations between archaea and human health.²⁰ Like bacteria, archaea are quantified with either shotgun or 16S amplicon sequencing.²¹

After sequencing is performed and data is processed, microbial abundances can be represented as a two-dimensional matrix of counts where each cell represents the estimated abundance of a taxon present in a particular sample. Computational methods are then used to gain biological insights from the data. A common analysis pipeline involves identifying differentially abundant taxa between different treatment groups (such as in microbiome case-control studies), which may be performed using software such as EdgeR.²² Many bioinformatics software packages can be downloaded with Bioconductor²³ in R or with Anaconda in Python, and Knight et al. provide a comprehensive overview of current best practices in microbiome analysis.²⁴

Omics layers analyzed with microbiome data

A multitude of different data layers can be analyzed along with the microbiome, such as genomics, transcriptomics, epigenomics, radiomics, proteomics, metabolomics, diet, and clinical outcomes. Each layer must be approached with its own considerations, which we review in this section.

Host transcriptomics

Host transcriptomics can provide insights into the functional interactions between host genes and microbiome by allowing researchers to quantify gene expression activity across different treatments or disease states. Transcriptomics may be performed with short- or long-read sequencing.²⁵ Transcriptomics protocols usually filter out highly-abundant ribosomal rRNA prior to sequencing by either depleting rRNA or performing poly(A) selection.²⁵ After sequencing is performed, reads are filtered based on their sequencing quality, aligned to reference genomes, and normalized for biases due to transcript length and compositionality.²⁵

Metatranscriptomics

Metatranscriptomic analysis techniques allow researchers to quantify the abundances of microbial gene transcripts in a sample, which can provide insight into the functional characteristics of the microbiota.²⁶ Metatranscriptomics experimental protocols differ based on the organisms being studied. For example, prokaryotic mRNA lacks a poly-adenine tail, which precludes the use of poly(A) selection.²⁶ After next-generation sequencing is performed, transcripts are aligned to metatranscriptomic reference genomes and quantified.²⁶

Host genetics

Studies combining host genetics with microbiome measurements have demonstrated that microbiome phenotypes are heritable.²⁷ As a result, studies may be interested in unraveling the host genetic determinants of microbiome composition. A variety of technologies exist for performing genotyping, including next-generation whole genome sequencing, whole exome sequencing, and genotyping arrays. Kockum et al. provide an in-depth overview of current genotyping strategies,²⁸ and a comprehensive discussion on the influence of host genetics and the gut microbiome is offered by Goodrich et al.²⁹

Metabolomics

Metabolomics allows researchers a better understanding of the biochemical and metabolic processes that involve the host and microbiota. Metabolomics data is usually generated using mass spectrometry, where molecules are identified by cross-referencing their mass-to-charge ratios against reference databases. The direct measurement of small molecules provides information about cellular processes and can point to mechanistic interactions between host and microbiota. Chong et al. provide a review of current techniques used to perform integrative analyses with metabolomic and microbiome data.³⁰

Metaproteomics

Metaproteomic analyses quantify the proteins produced by host and microbiome, which can provide insight into the functional role of microorganisms in host health. Complementary to metagenomics and metatranscriptomics, metaproteomic measurements reflect the activity of cellular translational and post-translational processes.³¹ As with metabolomics, metaproteomics is generally performed using mass spectrometry. Proteomics database lookup is computationally expensive, and researchers should be mindful that the results of proteomic analyses may be sensitive to the particular choice of mass spectra database used.³¹ Peters et al. provide an in-depth review on the role of proteomics analysis in microbiome studies.³¹

Diet

Diet is a major determinant of microbiome composition. Unlike other -omics layers, diet can be directly modified, allowing researchers to study the effects of dietary intervention on host physiology. Diet data collected from human studies is usually self-reported, taking the form of food-frequency questionnaires or single-day food records.³² More precise dietary measurements can be made, such as through interviews with dietary professionals, but usually at higher cost. After data collection, diet data is harmonized into a standardized form and may be converted into nutrient values prior to downstream analysis.³³

Clinical outcomes

With the increased use of electronic health record systems over the past several decades, clinical data has also been incorporated into multi-omic analyses. Clinical data includes medical diagnoses, lab results, doctors’ notes, and sensor data such as from glucose monitors, among other health-related information.³⁴ When combined with omics data, analysis of health record data can provide insight into the clinical consequences of -omics measurements, such as whether certain genotypes are associated with disease states.³⁴ Health record data is often incomplete and noisy, and therefore it must be harmonized to a standardized form prior to analysis.³⁴ Tong et al. include a comprehensive discussion on incorporating genomics with clinical information in studies of human health.³⁴

Challenges in microbiome multi-omics integration

The simultaneous analysis of many different data layers can provide greater insights into biological systems than the analysis of single layers separately.³⁵ Different datasets offer complementary views into the same biological system, elucidating relationships between different cellular processes, for example between mRNA transcript and protein levels or between particular protein structural motifs and translational pausing sites.³⁶ Integrative approaches that study the relationships between different types of biomolecules have had applications such as predicting risk of relapse in prostate cancer³⁷ and interrogating disease mechanisms in chronic kidney disease.³⁸ An extensive literature details multi-omics techniques and analysis strategies in a variety of application domains.^39,40 In multi-omics analyses, each -omics layer poses its own unique set of challenges. Metagenomic analyses in particular must overcome ambiguity in taxon assignments, compositionality and sparsity of data, variability of the microbiome over time, sensitivity of results to the analysis pipelines used, and a paucity of publicly available multi-omic microbiome datasets. We review these challenges below.

Ambiguity in taxon assignments

Microbial taxonomic labels are generally assigned imprecisely. In shotgun sequencing, abundances are inferred based on counts of short reads in sequencing experiments (usually <300bp in length) that are aligned against multiple reference genomes to determine their taxonomic origin. Due to the immense genetic variation and diversity in the human microbiome (Tierney et al. observed that 50% of genes observed in the oral and gut metagenomes were “singletons” present in one sample but not in any others), sequences may either not match any reference genome or may match to multiple reference genomes.⁴¹ Various approaches have been taken to assign sequences to taxonomies despite this ambiguity. Kraken, for example, is a popular alignment program that assigns shotgun sequencing reads to taxonomies using a maximum-score strategy.⁴² In 16S amplicon sequencing, taxonomic classification has generally been performed by clustering sequences at an arbitrary pre-defined threshold, such as 97% or 99% sequence similarity, with the resulting taxonomic assignments known as operational taxonomic units (OTUs).⁴³ As improvements in sequencing technology have led to a decrease in sequencing error rates, taxonomies can now be assigned more precisely using methods such as amplicon sequence variants (ASVs) or zero-radius OTUs (zOTUs), which offer single-nucleotide resolution in resolving amplicons, while also allowing for the use of sophisticated error models to correct sequences that might contain errors.^43–45 Nonetheless, 16S sequencing approaches have been shown to exhibit lower resolution compared to shotgun sequencing.^46,47 Classification specificity also varies based on the particular 16S subunit that is being sequenced.⁴⁸ Studies should therefore be careful to ensure that their results do not reflect sequencing or other technical artifacts instead of true biological signal.

Compositionality

Microbiome data is compositional, in that read counts for a particular sample can only provide information about the relative abundances of one taxon relative to others instead of absolute counts.⁴⁹ Values in a compositional dataset represent fractions of a whole and are therefore constrained to sum to 1, lying on an Atchinson simplex rather than on the full Euclidean space.⁵⁰ To address this, methods such as the additive or centered log ratio transform have been developed to transform compositional data into a form in which it can be readily analyzed with non-compositional analysis techniques such as linear models. Studies of the microbiome should be careful to use statistical models that account for compositionality and to avoid erroneous conclusions about the absolute abundances of taxa.

Sparsity

Microbiome data is oftentimes sparse, with a high proportion of zero counts for taxa in many samples. Zero-counts may not reflect true biological signal, but may instead be an artifact of technical factors such as metagenomic inference pipeline parameters or sequencing depth.⁵¹ Sparsity also limits the applicability of existing statistical models which do not explicitly account for the inflation of zero-counts. For example, dimensionality reduction methods such as PCA and PCoA may generate “horseshoe” patterns when applied on sparse data where samples do not share any taxa with each other.⁵² The Tweedie distribution, instead of the more common Poisson or negative binomial distribution, has also been used to model zero-inflated microbial abundances.⁵³

Variability

The gut microbiome is highly variable over time and is strongly influenced by factors such as antibiotic use and diet.⁵⁴ A study by Vandeputte et al. found greater day-to-day variation in microbiome composition within individuals than between individuals.⁵⁵ Johnson et al. also showed that microbiome profiles correlate with diet, and that these correlations vary at the individual level.⁵⁶ Furthermore, in animal models, the microbiome is also affected by caging, bedding, food, environment, and psychological stressors.⁵⁷ Studies should be careful not to draw conclusions from transient measurements of the microbiome, such as by increasing sample sizes, controlling for factors such as diet, or performing repeated sampling.³²

Small sample sizes

While traditional genome-wide association studies have benefited in recent years by the collection of sample sizes in the millions of participants over thousands of traits,⁵⁸ multi-omics studies currently have much fewer samples, usually in the tens or hundreds of participants,⁵⁹ with only a small subset also including microbiome data. Several repositories have been developed to address the lack of large microbiome datasets. For example, the Human Microbiome Project represents a recent effort to collect a large number of microbiome samples to facilitate large-scale computational analyses. As a part of the Human Microbiome Project, the Inflammatory Bowel Disease Multi’omics Database (IMDMDB) is a multi-omics repository specifically targeting inflammatory bowel disease.⁶⁰ The related Integrative Human Microbiome Project (iHMP) has also generated comprehensive integrated microbiome data in the context of preterm birth, inflammatory bowel disease, and type 2 diabetes.⁶¹ In addition, the TAILORED-Treatment consortium, completed in 2018, generated a multi-omics dataset of 1,200 clinical samples containing genotype, metabolomics, and metagenomics data. However, the HoPOIT database, in which these results are stored, is only available to specific research partners and not to the public.⁶² Despite the existence of these projects, few others also include microbiome multi-omics data. Due to the lack of integrated microbiome data in large genomic projects, researchers have turned to creative solutions. For example, Poore et al. have recently leveraged sequencing data from The Cancer Genome Atlas to infer per-sample microbial abundances from existing human WGS and transcriptomics reads.⁶³ Similarly, Dohlman et al. also used inferred microbial abundances using sequencing data associated with TCGA samples to construct the Cancer Microbiome Atlas, in a pipeline that emphasizes the decontamination of microbial reads.⁶⁴

Sensitivity to analysis pipelines

Bioinformatics analyses can be sensitive to pipeline parameters, providing different results when different differential abundance software, alignment pipelines, or reference databases are used.^31,65,66 When multiple -omics layers are processed together, this effect may be compounded further. Studies should therefore be mindful of this potential variability in results by validating their conclusions with multiple different approaches. Furthermore, when data is combined from multiple sources, such as when performing a meta-analysis or when combining data from different cohorts of individuals, researchers should be careful to reduce pipeline-related batch effects by ensuring that samples are processed uniformly.

Reliance on databases

In addition to these challenges, multi-omic bioinformatic analyses are dependent on curated databases that support the analysis of each individual data type. For example, 16S amplicon sequencing requires microbial sequences to be aligned to known rRNA sequences, which are stored in sequence databases such as SILVA.⁶⁷ Similarly, human metabolomics studies may rely on databases such as HMDB (Human Metabolome Database)⁶⁸ that provide detailed information for a multitude of metabolites; and transcriptomics studies may be reliant on pathway databases such as KEGG (Kyoto Encyclopedia of Genes and Genomes)⁶⁹ that contextualize the abundances of various transcripts. The inclusion of additional layers into analysis pipelines may require the development of additional databases, which can be time-consuming to create.

Due to the difficulty in collecting data from several -omics layers over a large number of samples, studies may choose to combine data from many different projects into a single meta-analysis. Such studies must be careful to account for batch effects. Microbiome analyses are sensitive to the taxon inference and differential abundance pipelines used.⁶⁵ Furthermore, in low-biomass samples, technical confounders such as sample procurement strategy may also introduce effects that are greater in magnitude than the treatments themselves.⁷⁰

Current approaches to microbiome multi-omics integration

While it is clear that integrating multiple omics layers provides benefits when studying the microbiome, there is no consensus on the best way of doing so. Integration may occur at various points in the analysis pipeline, and strategies vary with the study design and the questions being asked (Figure 1). In some studies, the inference of microbial counts is itself multi-omic, with abundances being estimated from metagenomic, metatranscriptomic and metaproteomic data at the outset⁷¹ (Figure 1a). In others, data preprocessing steps such as batch correction are performed using information from multiple -omics layers⁷² (Figure 1b). Others make use of computational methods that integrate data from several -omics layers simultaneously (Figure 1c). Another approach is to perform multiple separate sequencing experiments and analyze each -omics layer individually, such as in a study by Forslund et al. that investigated the effects of medications on both the gut microbiome and blood metabolome⁷³ (Figure 1d). Below we provide an overview of methods that have been used to perform microbiome analysis on multiple -omics layers, categorized based on their underlying statistical framework.

Figure 1. The integration of – omics layers can occur at various stages in the analysis pipeline. (a) information from multiple –omics layers can be combined to inform prediction of taxon or protein abundances. (b) multi-omics information may be combined after sequence alignment to inform per-sample measurements, such as by performing integrative batch correction. (c) per-sample measurements can be modeled together, such as in a linear or graphical model. (d) Analyses can also be performed on –omics layers individually and conclusions can be drawn from the combined analyses at the end of the study.

Dimensionality reduction and clustering methods

Dimensionality reduction is often the first step in any -omics analysis because it provides a quick way to visualize the overall structure of a dataset. The most common type of dimensionality reduction, Principal Component Analysis (PCA), linearly projects a high dimensional dataset onto the data axes with the highest variance so that it can be plotted in two dimensions. Principal Coordinates Analysis (PCoA, also known as Multidimensional Scaling) is another method which produces a two-dimensional embedding of the data such that distances between data points are preserved as closely as possible. PCA or PCoA plots are usually constructed separately for each data layer prior to integrative analysis in order to identify clear patterns in the data, such as whether data points are separated between disease and control groups.^74,75

Several other dimensionality reduction strategies have also been proposed, including Isomap, t-SNE, and UMAP, each performing a different transformation to embed the data onto two-dimensional space.⁷⁶ Dimensionality reduction methods have also been developed that operate on multiple omics datasets at the same time. Multi-omics Factor Analysis (MOFA), for example, is a linear multi-table integration method that finds a small set of numeric factors that best describe samples in a dataset.⁷⁷ MOFA takes as input a series of M data matrices, one for each modality over a common set of N samples, and a pre-specified number of factors F. It then decomposes each input matrix as the product of a common N×F factor matrix and a modality-specific F×D_m weight matrix, where D_m is the number of covariates measured in modality m. MOFA imposes regularization on the weight matrices, resulting in sparse sets of features that are associated with each factor. It was used by Garcia-Etxebarria et al. to identify genotypic, microbiome, and metabolomic factors associated with adenoma and colorectal cancer risk in a cohort of 120 individuals.⁷⁸ Meng et al.⁷⁹ include a detailed discussion about other integrative dimensionality reduction methods. Another common strategy when working with multiple data modalities is to summarize an entire data layer as a single numerical value that is then used in further analysis. Wang et al., for example, summarized multi-omic microbiome data into a single risk score measurement that was then used to predict disease susceptibility.⁸⁰ In the context of image analysis, Zhao et al. also summarized MRI scan data as a single number representing gray matter volume before integrating it with other data types.⁷⁴

Clustering algorithms can also be used to identify overall patterns in a dataset. Clustering can be performed either on samples, such as when clustering patients by their metabolic markers in order to infer disease subtypes,⁸¹ or on covariates, such as when identifying clusters of co-expressed genes.⁸² In clustering, elements (either samples or covariates) are partitioned into groups, called clusters, such that elements within a cluster are more similar to each other than to elements outside of their cluster. Measures of similarity between elements may be defined in a variety of ways. Common approaches include Euclidean distance, Manhattan distance, or Bray-Curtis dissimilarity.

Clustering analyses have the potential to identify patterns of microbiome composition that have clinical significance. Clustering has been used to assign humans as having different “enterotypes” based on their intestinal microbiota composition,⁸³ which may describe alterations of the microbiome that are associated with different disease states.⁸⁴ Additional studies have identified clusters in the nasal microbiome associated with disease states. For example, Lehtinen et al. identified nasal microbiome clusters that predicted viral load and symptom severity in rhinovirus infection.⁸⁵ Similarly, Abdel-Aziz et al. identified two distinct sputum microbiome clusters in asthma patients that were predictive of gene expression and protein abundances.⁸⁶

Several approaches also support the integrative clustering of multi-omics datasets,⁸⁷ which allows cluster assignments to capture more complex relationships between -omics layers. This may be helpful when finding modules of co-regulated or co-occurring biomolecules or when identifying disease subtypes based on molecular signatures. iCluster, for example, is a popular matrix factorization algorithm that uses multiple -omics layers simultaneously to infer sample cluster assignments, and has been used to identify cancer subtypes based on copy number and gene expression data.⁸⁸ Given a series of n input data matrices and a pre-specified number of clusters K, iCluster decomposes each data matrix X_i as the product of a common cluster assignment matrix Z and n layer-specific weight matrices W_i, such that the variance explained by the cluster assignments is maximized. iCluster also imposes an L1 penalty on the weight matrices, so that each cluster is explained by a sparse set of covariates.

Several studies have used multi-omics clustering approaches to identify cancer subtypes. For example, Yuan et al. used multi-omics clustering to identify subtypes of prostate and breast cancer that were associated with patient clinical outcomes and survival.⁸⁹ They proposed Patient-specific Data Fusion (PSDF), a Bayesian approach that both estimates the number of disease subtypes as well as performs cluster assignment. Similarly, Chaudary et al. used a deep-learning clustering approach based on autoencoders to identify cancer subtypes that were predictive of survival outcomes.⁹⁰ Another approach, LRACluster, uses a maximum likelihood approach combined with K-means clustering to infer cluster assignments in continuous, binary, and count data.⁹¹

Spectral clustering is another clustering technique that has seen applications in genomic datasets. Given a dataset of size n and a pre-specified number of clusters k, spectral clustering first computes a similarity matrix that describes pairwise relationships between data points. Then, an eigendecomposition is performed on the graph Laplacian of this matrix and the first k eigenvectors are taken, which represent inter-node relationships in a lower-dimensional subspace. Finally, k-means clustering is performed on the n × k matrix described by the eigenvectors, which results in the final cluster assignments. Given an appropriate similarity function, spectral clustering is able to identify clusters of arbitrary shapes in high-dimensional datasets. Von Luxburg provides an in-depth explanation on spectral clustering and describes how it is equivalent to an n-cut problem on graphs.⁹² Extensions of spectral clustering have been developed for multi-omics datasets, such as in one approach taken by Zhang et al. that finds clusters of cells in multi-modal single-cell sequencing data.⁹³

Correlation-based methods

Approaches to identifying relationships between quantitative covariates typically involve the computation of some pairwise measure of similarity, such as Pearson’s correlation coefficient. For example, one might compute pairwise correlation coefficients between the abundances of several microbes and the expression levels of several genes. Pearson’s correlation is the most common measure of similarity, but it cannot identify non-linear associations and tends to find spurious associations in compositional datasets.⁹⁴ Spearman correlation is an alternative method where covariates are first rank-transformed, so that each measurement is replaced with its integer-valued rank between 1 and N, where N is the number of samples, and the Pearson correlation is then computed on the ranks. Spearman correlations can find monotonic non-linear associations but still are prone to finding spurious associations in compositional data. Kendall’s tau is another nonparametric measure of association between two sets of quantitative measurements that operates on ranked values from each dataset. Liu et al. include an explanation of Kendall’s tau with a comparison to Spearman’s correlation.⁹⁵ Kendall’s tau has been used by Nayfach et al. to find associations between protein family abundances and clinical variables in patients with inflammatory bowel disease.⁹⁶

Correlation-based methods applied to the microbiome must be aware of the risks of compositionality. Spearman and Pearson correlation coefficients cannot be computed directly on compositional datasets without introducing spurious associations, and therefore several methods exist to transform compositional data from simplicial space into Euclidean space so that it can be analyzed with existing methods such as linear models.⁵⁰ One such commonly-used transformation, the centered-log-ratio (CLR) transform, replaces each measurement with the logarithm of the ratio of the measurement against the sample’s geometric mean, m‘_ij = log(m_ij/g(m_i)), where m_ij is the jth measurement in the ith sample and g(m_i) is the geometric mean of all measurements in the sample. It is then possible to find Pearson or Spearman’s correlations between pairs of CLR-transformed covariates.⁴⁹ Friedman et al. have also proposed SparCC, an alternative approach for finding correlations in compositional datasets, that performs better than CLR-transformed Spearman correlations on data with very few taxa but assumes sparsity of the correlation matrix, an assumption that is often useful in practice but might not always be appropriate.⁹⁷ Studies that make use of correlation-based methods should also be careful to not draw conclusions about causal relationships, which in general cannot be inferred from observational data.

In addition to correlations, a number of different measures of similarity have been proposed, including REBACCA,⁹⁸ mutual information,⁹⁹ and cosine similarity.¹⁰⁰ One study by Faust et al. used multiple measures of similarity at once, using an ensemble of Bray-Curtis dissimilarity, Kullback – Leibler divergence, Pearson correlation, and Spearman correlation to construct a microbial co-occurrence network.¹⁰¹ You et al. compared six different measures of metabolite-microbe similarity and concluded that each measure of similarity is best suited for a different scenario, with Spearman correlation being the most performant overall.¹⁰²

As an example of a correlation-based analysis, Wang et al. computed Pearson correlations between microbial, metabolite and host gene transcript abundances in colonic tumor and tumor-adjacent tissue and concluded that the genus Fusobacterium may interfere with butyrate synthesis, possibly resulting in tumorigenesis.¹⁰³ Spearman correlations were also used by Lloyd-Price et al. to construct a network of interactions between microbiome, host transcriptome, proteome, metabolome, and virome that revealed disease-associated features of interest for follow-up analysis, including a highly connected but yet-unclassified microbe with genus Subdoligranulum as well as several acylcarnitine metabolites.⁶⁰ Several additional studies have also performed multi-omics analysis by finding feature-by-feature correlations.^104,105

Canonical correlation analysis (CCA) is a method that finds sets of highly correlated covariates across paired datasets.¹⁰⁶ Given two data matrices representing data for the same samples, such as a table of microbe abundances and a matched table of transcript abundances, CCA finds subsets of features in each dataset that have the maximum correlation with each other (Figure 2a). The features in each component are also assigned weights so that they each represent a linear combination of the covariates in the data matrices. The components are then subtracted from the original datasets and the CCA procedure is applied again to generate several additional CCA components, each of which is orthogonal to all others. The results of CCA can sometimes be difficult to interpret, as the generated components typically include nonzero coefficients for all of the covariates. To improve the interpretability of the components generated, a variant called sparse CCA (sCCA) imposes an L1 penalty on the CCA procedure to limit the number of nonzero coefficients.¹⁰⁷ An example component from sparse CCA with an appropriate penalty parameter may include just a few microbes and their associated genes. Priya et al. used sparse CCA to characterize associations between the gut microbiota and host transcriptome in patients with colorectal cancer, inflammatory bowel disease, and irritable bowel syndrome, identifying common and disease-specific microbe-transcript interactions across all three disease types.¹³ Hyuk Park et al. similarly used CCA to relate microbiome abundances with gene expression in gastric cancer and identified associations between Helicobacteraceae and inflammation-related genes as well as between Pasteurellaceae and Lachnospiraceae and cancer-related genes.¹⁰⁸ In addition, multivariate generalizations of CCA such as sparse multiple CCA have been developed that integrate three or more data tables simultaneously.¹⁰⁹ As an example, Galié et al. used sparse generalized CCA (SGCCA) in a study of metabolic syndrome to identify four major clusters exhibiting crosstalk between gut microbiota, fecal metabolites, and host plasma metabolites.¹¹⁰

Figure 2. Comparison of linear modeling strategies. (a) Canonical Correlation Analysis (CCA) finds linear combination of covariates in each –omics layer that have a maximal correlation. (b) Linear regression identifies linear relationships between a response variable, such as gene expression, and many explanatory variables, such as taxon abundances. (c) Linear discriminant analysis (LDA) finds a best separating hyperplane between two sets of data points.

Procrustes analysis is another common technique to compare microbiome data and another set of measurements over a common set of samples. Procrustes analysis measures the extent to which one set of measurements can be translated, rotated, and scaled to best fit the other.¹¹¹ It has been frequently applied in the context of gut microbiome research, for example to find associations between the gut microbiome and diet,⁵⁶ and between the gut microbiome and host gene expression.¹³

Regression and classification methods

Regression and classification methods attempt to predict one set of variables from another, such as disease subtype from gene transcript abundances. The variables to be predicted, such as disease subtype, are known as the response or dependent variables. The variables from which a prediction is made, such as gene expression, are known as the covariates or independent variables. Regression and classification methods are useful for identifying which features are predictive of an outcome. For example, they may be used to identify the microbes which are most predictive of disease severity in colon cancer. Linear models represent the most common form of regression for quantitative covariates, where the response is predicted as a linear function of the covariates and is modeled to have Gaussian-distributed noise (Figure 2b).

Generalized linear models extend linear models by allowing the response to depend on the data based on a possibly nonlinear link function and to follow a non-Gaussian noise distribution. Generalized linear models were used by Tipton et al. to find associations between the oral microbiota and lab cytokine measurements.¹¹² Logistic regression is a type of generalized linear model that uses a logit link function and has been applied to find associations between gut microbiome composition and the onset of dementia.¹¹³ Polynomial regression can also be used to capture nonlinear dependencies in metagenomic and multi-omics datasets.¹¹⁴ In polynomial regression, the data table for a set of covariates is augmented with additional covariates that take the values of various powers of the original data, such as the squares and cubes of normalized microbial abundances. Linear regression is then performed on the augmented data table.

In settings where the number of predictors is greater than the sample size, as is often the case with multi-omics datasets, it is possible to perform linear regression using penalized approaches such as Lasso,¹¹⁵ Ridge,¹¹⁶ or Elastic net regularization.¹¹⁷ Regularized regression approaches simultaneously perform regression and variable selection and can limit the number of covariates in the generated model to a smaller, more interpretable subset. The HOMINID framework, for example, uses Lasso regression and stability selection to find associations between host SNPs and microbial taxa. Given datasets of host genetic variation and microbiome composition, LASSO regression is performed with host genetic variants as the model outcome and microbial taxonomic relative abundances as the covariates. Subsequently, stability selection is used to identify robust microbiome features that are associated with each host variant. It was applied on samples from the Human Microbiome Project across 15 body sites.¹¹⁸

A nested elastic net design was used by Ghaemi et al. to predict gestational age of pregnancy from cell-free transcriptomics, plasma and serum cytokine concentrations, metagenomics, blood mass cytometry, untargeted metabolomics, and targeted plasma proteomics.¹¹⁹ Elastic net models were first trained to predict gestational age individually from each of the seven data modalities. These models were then combined by training a second elastic net model that predicted gestational age from the outputs of these single-modality elastic net models. This nested model identified a possible regulatory interaction between the steroid hormone pregnenolone sulfate and myeloid dendritic cells and regulatory T cells, which are known to play a critical role in the maintenance of pregnancy.

Mills et al. used linear regression to analyze metabolomic, metagenomic, metapeptidomic, metaproteomic, and host proteomic data in a cohort of 250 patients with inflammatory bowel disease.¹²⁰ They performed univariate linear regression to predict disease severity separately from each covariate and found that Bacteroides proteins, particularly proteases originating from B. vulgatus, were strongly associated with ulcerative colitis. They then performed a functional analysis that revealed that proteases were most correlated with disease severity. Finally, they demonstrated in vivo that protease-inhibition reduced inflammation in mice inoculated with B. vulgatus, demonstrating a possible therapeutic target for ulcerative colitis.

Linear discriminant analysis (LDA) is a classification method that finds a best-fit hyperplane that separates labeled samples (Figure 2c). It can be used to determine the axes of variation that differ between two classes of individuals, such as those with and without a particular phenotype. As an example, Gomez-Llorente et al. performed an analysis using LDA to study the relationship between weight and intestinal microbiota composition in patients with asthma. They used sparse linear discriminant analysis (sPLS-DA) to construct three classifiers: obese vs. non-obese, overweight vs. non-overweight and normal-weight vs. non-normal weight.¹²¹ They found that leptin, acetate, and bacteria from the order Clostridiales were predictive of normal vs. non-normal weight. In the same study, they also identified 12 molecular features that differentiated persistent asthma from non-persistent asthma, such as creatinine and citrate concentrations.

Random forests are a commonly used classification method in which many decision trees, each individually weak and trained on a random subset of the features, are used to construct a stronger classifier that can discriminate between outcomes of interest such as disease state. Silveira et al., in a study including metagenomics, metabolomics and fluorescence microscopy, used random forests to study the role of microbes in the progression of cystic fibrosis. They identified facultative anaerobes such as Streptococcus as being keystone bacteria in the cystic fibrosis microbiome.⁷⁵ Zeybel et al. also used random forests to identify genomic features associated with hepatic steatosis in a cohort of 78 patients.¹²² Zhao et al. similarly used random forests and linear discriminant analysis in a cohort of 50 patients to study the gut microbiome in major depressive disorders using metagenome, metabolome, inflammatory factor and MRI data.⁷⁴

Network methods

An ultimate goal of multi-omics studies is to generate a comprehensive map of interactions between the various molecules and organisms within a system.¹²³ Networks, also known as graphs, provide an intuitive approach for visualizing and studying these interactions. Unlike regression-based methods that find associations between two groups of variables at a time, network-based methods aim to model all of the interactions between entities simultaneously. Such methods have been utilized to study co-occurrence relationships between microbial taxa,¹²⁴ co-abundance relationships between metabolites,¹²⁵ and co-expression relationships between genes.¹²⁶

Graphs are specified by a set of nodes, which usually refer to covariates such as gene expression values, microbial abundances, or disease states; and a set of edges, numerical values that connect pairs of nodes and which are usually a measure of similarity such as Pearson correlation or mutual information. A graph may be generated from a dataset by first finding correlations between all pairs of covariates, next identifying which correlations are significant, and then finally connecting pairs of nodes that are significantly correlated with each other.⁸³ After a graph is constructed, various analytics can be performed on it. For example, a “clique” analysis can identify highly interconnected nodes, while a “connectivity” analysis can identify distinct groups of nodes.¹²⁷ Layeghifard et al. provide a detailed review of general network approaches applied to microbiome datasets.¹²⁸ In some applications, known pathway information can also be integrated with graph inference.¹²⁹ Liu et al. provide a comprehensive treatment of network methods for the analysis of microbiome data.¹³⁰

Cantoni et al. used correlation-based networks to study the gut microbiome, blood immune cell concentrations, and circulating metabolites in patients with multiple sclerosis.¹³¹ They constructed a pairwise association network separately for MS patients and controls, connecting covariates if they had a pairwise correlation coefficient greater than 0.7. They found that memory Th1 cells were correlated with many metabolites in controls but not in MS patients; they also identified a significant pathway linking B. thetaiotaomicron, Th17 cells, S-adenosylmethionine (associated with Th17 cell activation), and meat consumption.

Pfalzer et al. used a network approach to study the interactions between transcriptome, metabolome, and microbiome in obese mouse models of colon cancer.¹³² Three groups of mice were either fed a low-fat (LF) diet, made obese with a high-fat (HF) diet, or made obese with a leptin receptor mutation (DbDb). Differential gene expression analysis was performed twice, between HF and LF mice and DbDb and LF mice, and a co-expression network was constructed individually for both sets of differentially expressed genes. Gene modules (sets of correlated genes) were identified from both networks and correlations were then computed between these gene modules and microbial and metabolic abundances using the WCGNA R library.¹³³ Overall, the study implicated Akt-signaling genes, microbial genera including Clostridium and Sarcina, and adenosine concentration as being associated with tumor burden.

Bayesian networks are a promising approach to studying interactions between covariates and diseases.¹³⁴ Bayesian networks model the entire joint distribution of covariates as a directed graph, specifying the dependence relationships of variables to each other and allowing for the inference of the value of one variable given several others. While Bayesian network inference is computationally intractable for large numbers of covariates, approximation-based techniques have had success in inferring Bayesian network structure from measured data. Hill-climbing, for example, is a greedy approach where edges are sequentially added to the Bayesian network so as to maximize a measure of fit, such as the Bayesian Information Criterion, at each step. Su et al. present a treatment of Bayesian networks in the context of identifying disease-related genes.¹³⁵ Krishnan et al. used Bayesian networks to study multi-omics data in nonalcoholic fatty liver disease (NAFLD).¹³⁶ They constructed a Bayesian network describing gene-gene interactions and their interactions with pathophysiological conditions and identified potential regulatory genes associated with NAFLD.

Spatial and temporal multi-omics

Advances in spatial analysis have enabled the integration of multi-omics data along with spatial information about the microbiome. For example, Shi et al. have demonstrated a novel technique, HiPR-FISH, which is able to produce spatial maps of microbial taxa within the gut and may be of interest in animal studies of the microbiome.¹³⁷ In another application of spatial multi-omics, Garg et al. developed a method to integrate metabolomic, 16S metagenomic, and spatial information to study the organization of the lung microbiome in cystic fibrosis patients.¹³⁸ A lung was procured from a recently deceased CF patient and sectioned into 86 segments from which metabolomic and metagenomic data was collected. Spatial maps of the multi-omics data were generated which demonstrated that microbial composition varies within the lung and that this variation may be related to the spatial variance of disease severity.

When multi-omics measurements are taken from samples at many different time points, it becomes possible to learn about the temporal shifts in microbial interactions with the various biomolecules in a system. As gut microbial composition is highly variable in the short term and fluctuates with factors such as diet and flare status in patients with IBS or IBD, temporal analyses allow for the explicit modeling of microbiome variability over time.^139,140 One study by Mihindukulasuriya et al. studied the temporal dynamics in the gut virome in the context of diet, bacterial microbiome, metabolome, and host genome and transcriptome. For each of 50 study participants who were either healthy or diagnosed with IBS-C or IBS-D, two consecutive stool samples were sequenced for viral DNA. For a subset of 28 participants, colonic transcriptomics data was also collected. Their results indicated that the viral microbiome was stable in stool samples among both IBS and healthy individuals, despite the high variability of bacterial composition. They also identified immune-related genes associated with the virome and disease-specific associations with specific phage populations, indicating that the virome might play a role in regulating bacterial composition.¹⁴¹ In another time-series analysis, Ruiz-Perez et al. fit a dynamic Bayesian network to model the time evolution of multi-omics microbiome data in inflammatory bowel disease, incorporating prior knowledge about multi-omics interactions to constrain the model to biologically plausible interactions. The Bayesian network successfully predicted taxon abundances at future time points. It also inferred several metabolite-to-microbiome relationships which were validated experimentally.¹⁴²

Temporal microbiome analyses have also been conducted in wastewater systems to obtain epidemiological insights into the prevalence of pathogenic agents in a community and to study the efficacy of sanitation practices.¹⁴³ Herold et al. employed a multi-omics time-series analysis to investigate the effect of disturbances on microbial communities in wastewater sludge. In order to identify the niches of microbial populations in wastewater, they collected metagenomic, metatranscriptomic, metaproteomic and metametabolomic data weekly over the course of 14 months. Multidimensional Scaling (MDS) was used to identify bacterial niches. For each of several functional categories (for example, “nucleic acid metabolism”), metatranscriptomic abundances were correlated with metaproteomics levels. Four functionally distinct niches were identified, each with different responses to environmental changes, suggesting that future research directions may involve untangling these distinct subpopulations.¹⁴⁴

Multi-omic taxon and biomolecule inference

The previous methods have used multi-omics data to study the interactions between different biomolecules. However, multiple -omics datasets can also be used at an earlier stage in the analysis pipeline in order to improve inference of the biomolecules themselves. gNOMO provides one such technique.¹⁴⁵ Protein identification is commonly achieved through mass spectrometry, a technique that typically allows the quantification of only a subset of the entire set of proteins in a sample.¹⁴⁶ gNOMO uses metatranscriptomics and metagenomics information within the sample to predict the set of proteins that may be encoded by bacteria and the host using a program called Prodigal. This set of likely proteins is used as the mass spectrum lookup database for mass spectrometry.

Other approaches have combined metagenomic and metatranscriptomic information to improve the inference of microbial abundances. For example, Heintz-Buschart et al. observed that a significant proportion of reads sequenced with shotgun metagenomic sequencing did not map to a fully-sequenced reference genome.⁷¹ In order to improve taxon inference, they performed co-assembly with metagenomic and metatranscriptomics reads, which resulted in longer contig lengths and higher read usage. Furthermore, they also inferred metaproteomic abundances using a protein search database informed by the observed metagenomic and metatranscriptomic sequences.

Metabolic modeling

Microbiome metabolic modeling is an approach that constructs a mathematical representation of the entire set of biochemical reactions that occur within a microbial ecosystem. Metabolic models allow researchers to reason about inter-microbe and microbe-environment interactions. Metabolic models can be used to predict steady-state microbial abundances and the rates of metabolite consumption and secretion by bacteria in a system.¹⁴⁷ Magnúsdóttir et al. developed the AGORA (Assembly of Gut Organisms through Reconstruction and Analysis) resource which contains genome-scale metabolic constructions for 773 human gut bacteria.¹⁴⁸ AGORA was initially applied to identify a growth medium for Bacteroides caccae. In another application, a metabolic model using AGORA was used to identify differences in microbial metabolism between patients with IBD and healthy controls.¹⁴⁹

Metabolic models support the integration of multi-omics data by imposing constraints based on metagenomic, metabolomic, metaproteomic, or diet information.¹⁴⁹ Yizhak et al., for example, developed IOMA (integrative omics-metabolic analysis), a method that integrates metabolic reconstructions with matched proteomics and metabolomics measurements in order to find a set of metabolic fluxes that is consistent with the known metabolic constructions as well as with measured data.¹⁵⁰

Conclusions and future directions

Integrative analysis has the potential to provide insights into disease biology by illuminating the complex interplay between the host, microbiome, and the various biomolecules through which they interact. The integration of multiple -omics layers can elucidate interactions that may not be apparent when considering data layers individually, which improves the discovery of disease mechanisms.¹²⁰

An emerging application of multi-omics analysis is in precision medicine, where measurements from multiple -omics layers are used to inform treatment decisions, such that care is targeted toward the particular physiology of the patient. Due to the multifactorial role of the microbiota, such as in the production of bioactive compounds and the metabolism of pharmaceutical drugs, the microbiome provides a promising target for precision medicine. For example, it may be beneficial to adjust medications or dosages based on a patient’s microbiome composition or on other molecular phenotypes.¹⁵¹ Therapies may also use multi-omics information to target the microbiome itself, such as by selectively modulating bacterial abundances based on host physiology.¹⁵²

While a variety of approaches have been developed to support multi-omics integration, there is a current lack of standardization among them, which makes it difficult to interpret whether the results of multi-omic studies capture true signal or pipeline-related artifacts. A set of best practices should be established for integrating data layers together, which would help clarify the approaches that are best suited for each experimental design and improve the ability to compare results between studies. Standardization is especially relevant for microbiome analyses, which are sensitive to the pipeline parameters and bioinformatics software used.⁶⁵ Another important step in standardization is the development of readily accessible software packages that support multi-omics microbiome integration, as few currently exist.¹⁵³ Several software packages support the integration of multiple -omics layers in general.¹⁵⁴ However, these methods do not specifically address the challenges relevant to microbiome analysis, and it is necessary to benchmark how well they perform on microbiome datasets.

Another major difficulty in conducting a multi-omics study is the high cost of sample collection and data generation, which limits the ability of individual studies to detect associations with small effect sizes. Genetics research has benefited by the creation of large repositories, such as ENCODE and TCGA, that aggregate multi-omic genotypic and phenotypic data across different studies, and which researchers can use to perform large-scale computational analyses.^155,156 It would be beneficial for such repositories to be developed and expanded in the context of microbiome multi-omics, with an emphasis on taking into account differences in pipeline parameters and batch effects between studies.

Acknowledgments

We thank Rich Abdill for insightful comments on this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The work was supported by the National Institutes of Health [R35-GM128716 to RB].