A review of common methods used in the analysis of human microbiome sequencing data

Search strategy

The following electronic databases were searched using the relevant keywords and search syntax appropriate to each database. The bibliographies of the articles were hand searched to identify additional studies.

Only primary articles in the English language were considered. The search used a combination of the following search terms:

A total of 6100 studies were extracted from electronic databases. Of the 6100 abstracts included, 381 were randomly selected to be assessed by the author to determine the study population, type of metagenomic sequencing used and analysis methods. To further assess the types of methods used in metagenomic studies an R script was written to extract keywords from all abstracts returned by the search strategy in order to determine sequencing methods and possible analysis methods for the studies that were not selected to be reviewed. All text extraction and analyses were performed in R (version 4.1.2) (R Core Team, 2020). The full text of an additional five randomly selected abstracts was reviewed for each section of the review (normalisation methods, diversity metrics, and, inferential statistics for microbiome data). For full text review, 18 papers that discussed methods for compositional data in their abstracts were included. In total the full text of 38 randomly selected articles was reviewed.

Normalisation of metagenomic data

Comparison between samples within microbiome sequencing experiments is challenging because it is often very difficult to collect the same number of reads for each sample even if the same sample is sequenced multiple times. Differences in sample read count can arise from differences in the sequencing machine, difficulties in loading the same molar amounts of the sequencing libraries, or because of random variation (Gloor, Macklaim et al., 2017). The total read depth is a strong confounder for any distance or diversity metric used and is thus also a strong confounder for any downstream analyses that make use of these distance or diversity metrics (McMurdie and Holmes, 2014).

Initially “rarefying” was a common method of normalising read counts between samples, however McMurdie and Holmes (2014) question the validity of rarefying because “although rarefying does equalise variances, it does so only by inflating the variances in all samples to the largest (worst) value among them at the cost of discriminating power (increased uncertainty)”. The prevalence of rarefying in microbiome sequencing studies is likely to due to software pipelines for analysing microbiome data that implement rarefying as the default method for dealing with large differences in read counts between samples. A variety of alternative normalisation strategies for metagenomic data have been proposed and are found in the literature. Commonly used alternatives to rarefying include the trimmed mean of M-values (TMM) (McCarthy et al., 2012; Robinson and Oshlack, 2010), and the normalisation scaling factor calculated by DESeq2 (Love et al., 2014). The utility of count normalisation strategies in general has also been questioned, as the number of counts returned by sequencing instruments do not contain information on the absolute number of molecules in the original environment due to the already described constraint to a fixed total number of reads. Researchers using a compositional data analysis approach normalise the data prior to analysis when sequence count data are converted to log-ratios (Gloor, Macklaim et al., 2017).

The choice of normalisation method can dramatically alter the ability of methods to distinguish signal from noise depending on the methods used (Lin et al., 2016; Weiss et al., 2017). For example, the choice of normalisation methods can alter conclusions drawn from community comparisons using beta diversity metrics (McKnight et al., 2019). Some normalisation methods, such as TMM normalisation, can give different results depending on how low abundance reads get removed from the dataset.

Microbiome diversity metrics

A key scientific objective of many metagenomic sequencing experiments is to investigate if a particular microbial composition is associated with a given outcome. For example, investigating if the diversity in a sample from an individual with a particular disease differs from that of a comparable healthy control. Microbiome sequencing studies have been used to investigate questions in a wide range of research areas including epidemiology (Willner et al., 2009), clinical diagnostics, and as potential treatments for certain diseases (Rossen et al., 2015).

One common way this is done in the literature is to calculate a summary measure of the diversity present in a given sample, typically an alpha diversity measure, and then compare the mean species diversity between the two groups representing those samples. These diversity metrics are then often used in further downstream analyses, such as ordination or regression. Beta diversity, a measure of the dissimilarity of microbial communities between samples, is also commonly calculated in microbiome studies. Beta diversity dissimilarity matrices are used for ordination, such as principal coordinates analysis (PCoA) or principal component analysis (PCA), or classification depending on the research question. These methods allow for the visualization of complex high-dimensional dissimilarity matrices in lower dimensional spaces. Ordination plots can then be supplemented by the inclusion of sample metadata to visualize any possible sample clustering based on the metadata, for example PCA plots can be generated and coloured by sample disease status. Although the total number of reads in a sample is a strong confounder for dissimilarity-based methods (McMurdie and Holmes, 2014), dissimilarity matrices generated after the data has been normalised may still provide useful information about sample dissimilarity. Three such dissimilarity matrices are commonly used in the literature; Bray-Curtis dissimilarity, Jensen-Shannon distance, and UniFrac distance (Lozupone et al., 2011).

The use of these diversity metrics in metagenomic studies has been criticised as inappropriate for compositional data, as the dissimilarity metrics do not take into account the compositional nature of microbiome data. Beta diversity metrics are sensitive to the total sequencing read depth of a sample, and many distance or dissimilarity methods discriminate between samples largely based on the most relatively abundant features in the samples, not on the features that are necessarily the most variable between samples (Gorvitovskaia et al., 2016; Wong et al., 2016). This can lead to a substanial shift in the location of included samples in an ordination depending on included and excluded dataset features (Gloor, Macklaim et al., 2017; Wong et al., 2016). Alternatives for compositional data such as Aitchison distance and adjustments to UniFrac distances have been proposed (Wong et al., 2016) that do account for the compositional nature of the data.

Inferential statistics for microbiome data

Alpha and beta diversity metrics provide useful descriptive summaries of the microbial diversity present in a given sample. However, in many microbiome studies the aim is not merely to determine if samples are more or less diverse than one another but rather to determine which microbes are differentially abundant between the comparison groups of interest. In compositional data structures, when the proportion of one microorganism increases, the proportions of the other microorganisms in the sample must by definition decrease to allow for the standardised proportion of all OTUs to sum of 1 (Aitchison, 1986; Gloor, Wu et al., 2016). Furthermore, microbiome data is sparse and high-dimensional, typically containing many more features than there are samples and many OTUs are either not detected in all samples or are removed during data quality control and pre-processing pipelines (Paulson et al., 2013). Count distributions arising from microbiomes sequencing data are also known to be overdispersed (McMurdie and Holmes, et al., 2014; Richards, 2008). The combination of these features of microbiome sequencing data makes determining true differences in microbial composition challenging and calls for the use of specialized analytic methods.

Identifying differentially abundant features in microbiome data is a common goal of microbiome sequencing studies. As such there are a number of methods that have been applied to detecting differentially abundant features. Although not appropriate for compositional data, classical statistical methods such as parametric statistical tests (for example, two sample t-tests, Wilcoxon rank sum tests, and ANOVA) and measures of correlation, including Spearman’s rank correlation, and Pearson correlation are frequently used to investigate differentially abundant organisms and microbial networks. More sophisticated methods have been developed to assess differentially abundant features between samples by taking into account the nature of microbiome data as well as the number of hypothesis tests performed during differential abundance analysis. For example, one of the most common methods for performing differential abundance analyses on an entire microbiome sequencing dataset at once is implemented in the R package DESeq2 (Love et al., 2014). Given the challenges associated with metagenomic data and multiple testing, methods for detection of differentially abundant organisms in microbiome sequencing studies are an active area of research. The two most common software pipelines for bioinformatic processing of metagenomic sequencing data are QIIME2 (Bolyen et al., 2019) and mothur (Schloss et al., 2009), both of these software pipelines include tools to prepare, quality control, and analyse data from metagenomic experiments. Failing to account for compositionality of metagenomic data means that findings for microbiome analysis carried out with these methods are particularly sensitive to negative correlation bias. Furthermore, many current methods for differential abundance analysis in microbiome studies have been shown to have unacceptably high false discovery rates. In some cases, depending on the method and data pre-processing pipeline the false discovery rates exceeded 80% (Hawinkel et al., 2019; Weiss et al., 2017).

Compositional data analysis of microbiome

An important step in microbiome sequencing projects is visualization of the high-dimensional abundance tables. Typically, this is done by reducing the dimensionality of the data via PCoA, and then plotting the resulting principal coordinates. In the compositional data analysis framework, the equivalent of the PCoA plot is the variance-based compositional principal component (PCA) biplot (Aitchison and Greenacre, 2002). Compositional alternatives for investigating the relative abundance of OTUs between groups have also been developed to replace traditional differential abundance analysis methods. For example, the ALDEx2 R package provides compositional methods for estimation of differential expression by performing statistical tests on the centered log-ratio (CLR) values from a modelled probability distribution of the dataset and reports the expected values of parametric and non-parametric statistical tests along with estimate of effect-size (Fernandes, Macklaim et al., 2013; Fernandes, Reid et al., 2014). The ANCOM R package offers an alternative to ALDEx2 for differential abundance testing while taking into account compositionality. ANCOM performs statistical tests on point estimates of data transformed by an additive log-ratio (ALR) transformation, where invariant taxa are chosen as the reference for the ALR transformation (Mandal et al., 2015).

Although not yet widely used in the literature, there exist statistical methods for the analysis of compositional microbial data from cross-sectional microbiome sequencing experiments (Gloor, Macklaim et al., 2017). For example, sequencing read counts are explicitly normalized across samples when analysed as compositional, because compositional analyses make use of log-ratio transformations of the abundance tables (Aitchison, 1986). Log-ratios have the important mathematical property that their sample space is the whole real number line, and thus many standard statistical analyses that have been developed for real random variables can be used in conjunction with log-ratios while appropriately taking into account the compositionality of the data. Similarly in some applications the Aitchison distance, defined as the Euclidean distance between two log-ratio points, can be used as the compositional alternative to common ecological distance measures, such as Bray-Curtis, UniFrac, and Shannon distances (Gloor, Macklaim et al., 2017). There are also proposed compositional alternatives to determine correlation and other measures of association for compositional data, such as the proportionality (Φ) (Lovell et al., 2015) and proportionality coefficient (p) (Erb and Notredame, 2016) metrics.

The 12 papers that analysed microbiome data using methods for compositional data can be broadly classified into two groups; studies that analysed data from longitudinal microbiome experiments (three papers) and studies that analysed data from cross-sectional microbiome sequencing experiments (nine papers). Of the papers that analysed cross-sectional microbiome data three papers did not fully utilise a compositional data analysis pipeline, instead only making use of ANCOM to evaluate differentially abundant OTUs (Chung et al., 2019; Iszatt et al., 2019; Lee et al., 2014). All other papers that were included in this review transformed the data via the CLR transform prior to performing statistical tests on the OTU tables. Although the utility of ecological diversity measures for metagenomic data is debated, three papers that utilised compositional methods reported conventional ecological diversity metrics, such as Shannon and UniFrac distances (Genco et al., 2019; Peters et al., 2017; San-Juan-Vergara et al., 2018), instead of the compositional alternatives such as Aitchison distance.

Only three studies used compositional methods to investigate the change in human microbial composition over time, a specific focus of this review, however all three of the studies evaluated time points independently from one another. For example, Dahl et al. (2018) investigated the gut microbiome in preterm infants at ten days, four months, and one year. The investigators performed a standard analysis of metagenomic data without accounting for compositionality, calculated traditional diversity metrics and evaluated differences in beta diversity using PERMANOVA. The authors then used ANCOM to test for differentially abundant OTUs at each time point separately. The other two papers that investigated the microbiome at two timepoints transformed all OTU counts via the centered log-ratio transformation. In the study by Dijkhuizen et al. (2019) zero values were dealt with by imputing the missing values using the robCompositions R package. Log ratio-transformed OTU relative abundances were then analysed using principal components analysis (PCA), and group difference in log-ratio transformed relative abundances were assessed with logistic regression. Notably, baseline, inactive disease, and persistent activity samples were analysed independently. Finally, Stanaway et al. (2016) measured exposure to seasonal pesticides during the spring-to-summer fruit-thinning season and the offseason. Prior to transformation into log-ratios alpha diversity measures were calculated, then OTU tables were then transformed using the CLR transform. The PCA on CLR transformed data was used as an exploratory examination of the beta diversity in spring/summer and winter samples separately. Group differences in PC1 scores from the PCA were compared using t-tests. The PCA data was then used to determine clusters in the spring/summer and winter season data respectively, in order to determine if a cluster with exposure to azinphos-methyl could be distinguished from non-exposure based on OTU table data.

An area of open research in compositional methods for microbiome data is the application of compositional analysis methods to data from longitudinal microbiome sequencing studies where we must assume there exist correlation within individuals. Theoretical methods for the analysis of some forms of longitudinal dependence in compositional data have been developed (Pawlowsky-Glahn et al., 2015), however these methods have not yet been applied or tested for use in microbiome sequencing studies.