Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis

Abstract

Background: Experimental designs that take advantage of high-throughput sequencing to generate datasets include RNA sequencing (RNA-seq), chromatin immunoprecipitation sequencing (ChIP-seq), sequencing of 16S rRNA gene fragments, metagenomic analysis and selective growth experiments. In each case the underlying data are similar and are composed of counts of sequencing reads mapped to a large number of features in each sample. Despite this underlying similarity, the data analysis methods used for these experimental designs are all different, and do not translate across experiments. Alternative methods have been developed in the physical and geological sciences that treat similar data as compositions. Compositional data analysis methods transform the data to relative abundances with the result that the analyses are more robust and reproducible.

Results: Data from an in vitro selective growth experiment, an RNA-seq experiment and the Human Microbiome Project 16S rRNA gene abundance dataset were examined by ALDEx2, a compositional data analysis tool that uses Bayesian methods to infer technical and statistical error. The ALDEx2 approach is shown to be suitable for all three types of data: it correctly identifies both the direction and differential abundance of features in the differential growth experiment, it identifies a substantially similar set of differentially expressed genes in the RNA-seq dataset as the leading tools and it identifies as differential the taxa that distinguish the tongue dorsum and buccal mucosa in the Human Microbiome Project dataset. The design of ALDEx2 reduces the number of false positive identifications that result from datasets composed of many features in few samples.

Conclusion: Statistical analysis of high-throughput sequencing datasets composed of per feature counts showed that the ALDEx2 R package is a simple and robust tool, which can be applied to RNA-seq, 16S rRNA gene sequencing and differential growth datasets, and by extension to other techniques that use a similar approach.

Keywords: 16S rRNA gene sequencing; Dirichlet distribution; Monte Carlo sampling; RNA-seq; centered log-ratio transformation; compositional data; differential abundance; high-throughput sequencing; microbiome.

Figures

Figure 1

Outline of the approach for one feature in three control and three experimental samples. The count values for feature i, sample j are converted to probabilities by Monte Carlo sampling from the Dirichlet distribution with the addition of a uniform prior. Each count value is now represented by a vector of probabilities 1:n, where n is the number of Monte Carlo instances sampled: three instances are shown in the example, but 128 are used by default. Each probability in the vector is consistent with the number of counts in feature i given the total number of reads observed for sample j. Each Monte Carlo Dirichlet instance is center log-ratio transformed giving a vector of transformed values. These values are the base 2 logarithm of the abundance of the feature in each Dirichlet instance in each sample divided by the geometric mean abundance of the Dirichlet instance of the sample. Significance tests for control samples (C1 : C3) vs experimental samples (E1 : E3) are performed on each element in the vector of clr values. Each resulting P value is corrected using the Benjamini–Hochberg procedure. The expected values are reported for both the distribution of P values and for the distribution of Benjamini–Hochberg corrected values. clr, centered log-ratio; FDR, false discovery rate.

Figure 2

Effect of DMC sampling on the selex dataset. The first column shows the results when the data is clr transformed without DMC sampling, the next three show the effect of 1, 16 and 128 DMC samples followed by the clr transformation. Features that pass a threshold P<0.05 are shown in cyan and those where the fdr statistic is <0.05, are shown in red. Features where the median clr value is below the geometric mean are highlighted in black if they are not significant. Those where the median clr value is greater than the geometric mean are shown in gray. clr, centered log-ratio; DMC, Dirichlet Monte Carlo.

Figure 3

MA plot for DESeq. The base 2 logarithm of average expression across all samples for a feature is plotted vs the base 2 logarithm of fold-change. Points that are significantly different with a fdr less than 0.05 are in red, all others are in gray.

Figure 4

Differential features in common between ALDEx2, DESeq and baySeq. Genes are colored in light yellow if ALDEx and at least one of the other two tools identified them as significantly different with an fdr <0.05, black if they were identified by both baySeq and DESeq, magenta if only by baySeq, cyan if only by DESeq, and orange if only by ALDEx2. Small gray dots are non-differential genes. The Venn diagram illustrates the number of differentially abundant genes identified by each method. MA, mean difference between conditions vs average expression; MW, mean difference between conditions vs maximum within-condition variance.

Figure 5

OTUs with different relative abundances between tongue dorsum and buccal mucosa. Each OTU is colored by membership in the taxonomic level indicated. OTU abundance values are median relative abundance values derived from ALDEx2. OTU, operational taxonomic unit.