Statistical evaluation of methods for identification of differentially abundant genes in comparative metagenomics.

BACKGROUND: Metagenomics is the study of microbial communities by sequencing of genetic material directly from environmental or clinical samples. The genes present in the metagenomes are quantified by annotating and counting the generated DNA fragments. Identification of differentially abundant genes between metagenomes can provide important information about differences in community structure, diversity and biological function. Metagenomic data is however high-dimensional, contain high levels of biological and technical noise and have typically few biological replicates. The statistical analysis is therefore challenging and many approaches have been suggested to date. RESULTS: In this article we perform a comprehensive evaluation of 14 methods for identification of differentially abundant genes between metagenomes. The methods are compared based on the power to detect differentially abundant genes and their ability to correctly estimate the type I error rate and the false discovery rate. We show that sample size, effect size, and gene abundance greatly affect the performance of all methods. Several of the methods also show non-optimal model assumptions and biased false discovery rate estimates, which can result in too large numbers of false positives. We also demonstrate that the performance of several of the methods differs substantially between metagenomic data sequenced by different technologies. CONCLUSIONS: Two methods, primarily designed for the analysis of RNA sequencing data (edgeR and DESeq2) together with a generalized linear model based on an overdispersed Poisson distribution were found to have best overall performance. The results presented in this study may serve as a guide for selecting suitable statistical methods for identification of differentially abundant genes in metagenomes.

Upper Respiratory Infections and MRI Activity in Relapsing-Remitting Multiple Sclerosis.

BACKGROUND: Although clinical reports have suggested a relationship between systemic infections and multiple sclerosis (MS) relapses, MRI evidence supporting an association is conflicting. Here we evaluated the temporal relationship between upper respiratory infections (URIs) and MRI activity in relapsing-remitting (RR) MS. METHODS: We combined individual data on URI with data on active lesions in pre-scheduled MRI examinations performed every 4 weeks for 28 weeks in 69 patients. A 4-week at-risk (AR) period started, by definition, 1 week before the onset of a URI. We recorded the relationship between the number of active lesions in each MRI with (1) the number of days of AR time in the immediately preceding 4-week period and (2) the number of days passed since the onset of a preceding URI. RESULTS: Average MRI lesions/day showed no difference between AR (0.0764) and not-AR (0.0774) periods. The number of lesions in 483 pre-scheduled MRI examinations did not correlate with the AR proportion in the prior 4-week period (rho = -0.03), and time from URI onset did not correlate with lesion number on the next MRI examination (rho = 0.003). CONCLUSION: The occurrence of a URI did not increase the risk of MRI activity evaluated in an adjacent 4-week window in RRMS.

Changes in HbA1c and frequency of measuring HbA1c and adjusting glucose-lowering medications in the 10 years following diagnosis of type 2 diabetes: a population-based study in the UK.

AIMS/HYPOTHESIS: The aim of this work was to study levels of HbA1c and patterns of adjusting glucose-lowering drugs in patients with impaired glycaemic control over 10 years after diagnosis of type 2 diabetes. METHODS: We studied 4,529 individuals in The Health Improvement Network Database newly diagnosed with type 2 diabetes in the year 2000. RESULTS: From 6 months to 10 years after diagnosis, the HbA1c increased from 7.04% (53.4 mmol/mol) to 7.49% (58.3 mmol/mol) (average annual change: 0.047% [0.51 mmol/mol]). The greatest annual change occurred between 6 months and 2 years (0.21% [2.30 mmol/mol] increase per year, p < 0.001), followed by the 2-5 year time period (0.033% [0.36 mmol/mol] increase per year, p < 0.001). No significant increase in HbA1c occurred between 5 and 10 years (p = 0.20). In multivariable analyses, patients who were younger (p < 0.001), with higher BMI (p = 0.033) and who were current insulin users (p = 0.024) at diagnosis had greater increases in HbA1c between 6 months and 2 years. For individuals with HbA1c above 7.0% (53 mmol/mol) the mean time to next measurement of HbA1c was 0.53 years and increase in doses or changes to other glucose-lowering medications were performed in 26% of cases. CONCLUSIONS/INTERPRETATION: HbA1c increases by approximately 0.5% (5 mmol/mol) over 10 years after diagnosis of type 2 diabetes, with the main increase appearing in the first years after diagnosis. More frequent monitoring of HbA1c and adjustments of glucose-lowering drugs may be essential to prevent the decline.

A novel method for cross-species gene expression analysis.

BACKGROUND: Analysis of gene expression from different species is a powerful way to identify evolutionarily conserved transcriptional responses. However, due to evolutionary events such as gene duplication, there is no one-to-one correspondence between genes from different species which makes comparison of their expression profiles complex. RESULTS: In this paper we describe a new method for cross-species meta-analysis of gene expression. The method takes the homology structure between compared species into account and can therefore compare expression data from genes with any number of orthologs and paralogs. A simulation study shows that the proposed method results in a substantial increase in statistical power compared to previously suggested procedures. As a proof of concept, we analyzed microarray data from heat stress experiments performed in eight species and identified several well-known evolutionarily conserved transcriptional responses. The method was also applied to gene expression profiles from five studies of estrogen exposed fish and both known and potentially novel responses were identified. CONCLUSIONS: The method described in this paper will further increase the potential and reliability of meta-analysis of gene expression profiles from evolutionarily distant species. The method has been implemented in R and is freely available at http://bioinformatics.math.chalmers.se/Xspecies/.

Transcriptional and metabolic data integration and modeling for identification of active pathways.

With the growing availability of omics data generated to describe different cells and tissues, the modeling and interpretation of such data has become increasingly important. Pathways are sets of reactions involving genes, metabolites, and proteins highlighting functional modules in the cell. Therefore, to discover activated or perturbed pathways when comparing two conditions, for example two different tissues, it is beneficial to use several types of omics data. We present a model that integrates transcriptomic and metabolomic data in order to make an informed pathway-level decision. Since metabolites can be seen as end-points of perturbations happening at the gene level, the gene expression data constitute the explanatory variables in a sparse regression model for the metabolite data. Sophisticated model selection procedures are developed to determine an appropriate model. We demonstrate that the transcript profiles can be used to informatively explain the metabolite data from cancer cell lines. Simulation studies further show that the proposed model offers a better performance in identifying active pathways than, for example, enrichment methods performed separately on the transcript and metabolite data.

Expression profiling of GIST: CD133 is associated with KIT exon 11 mutations, gastric location and poor prognosis.

In gastrointestinal stromal tumors (GISTs), KIT exon 11 deletions are associated with poor prognosis. The aim of this study was to determine the gene expression profiles of GISTs carrying KIT exon 11 deletions and to identify genes associated with poor prognosis. Expression profiling was performed on nine tumors with KIT exon 11 deletions and 7 without KIT exon 11 mutations using oligonucleotide microarrays. In addition, gene expression profiles for 35 GISTs were analyzed by meta-analysis. Expression of CD133 (prominin-1) protein was examined by tissue microarray (TMA) analysis of 204 GISTs from a population-based study in western Sweden. Survival analysis was performed on patients subjected to R0 resection (n=180) using the Cox proportional hazards model. Gene expression profiling, meta-analysis, and qPCR showed up regulation of CD133 in GISTs carrying KIT exon 11 deletions. Immunohistochemical analysis on TMA confirmed CD133 expression in 28% of all tumors. CD133 positivity was more frequent in gastric GISTs (48%) than in small intestinal GISTs (4%). CD133 positivity was also more frequent in GISTs with KIT exon 11 mutations (41%) than in tumors with mutations in KIT exon 9, platelet-derived growth factor receptor α (PDGFRA), or wild-type tumors (0-17%). Univariate survival analysis showed a significant correlation between the presence of CD133 protein and shorter overall survival (hazard ratio=2.23, p=0.027). Multivariate analysis showed that CD133 provided additional information on patient survival compared to age, sex, National Institutes of Health (NIH) risk group and mutational status. CD133 is expressed in a subset of predominantly gastric GISTs with KIT exon 11 mutations and poor prognosis.

mRNA stability changes precede changes in steady-state mRNA amounts during hyperosmotic stress.

Under stress, cells need to optimize the activity of a wide range of gene products during the response phases: shock, adaptation, and recovery. This requires coordination of several levels of regulation, including turnover and translation efficiencies of mRNAs. Mitogen-activated protein (MAP) kinase pathways are implicated in many aspects of the environmental stress response, including initiation of transcription, translation efficiency, and mRNA turnover. In this study, we analyze mRNA turnover rates and mRNA steady-state levels at different time points following mild hyperosmotic shock in Saccharomyces cerevisiae cells. The regulation of mRNA stability is transient and affects most genes for which there is a change in transcript level. These changes precede and prepare for the changes in steady-state levels, both regarding the initial increase and the later decline of stress-induced mRNAs. The inverse is true for stress-repressed genes, which become stabilized during hyperosmotic stress in preparation of an increase as the cells recover. The MAP kinase Hog1 affects both steady-state levels and stability of stress-responsive transcripts, whereas the Hog1-activated kinase Rck2 influences steady-state levels without a major effect on stability. Regulation of mRNA stability is a wide-spread, but not universal, effect on stress-responsive transcripts during transient hyperosmotic stress. By destabilizing stress-induced mRNAs when their steady-state levels have reached a maximum, the cell prepares for the subsequent recovery phase when these transcripts are to return to normal levels. Conversely, stabilization of stress-repressed mRNAs permits their rapid accumulation in the recovery phase. Our results show that mRNA turnover is coordinated with transcriptional induction.

Historical HbA1c Values May Explain the Type 2 Diabetes Legacy Effect: UKPDS 88.

OBJECTIVE: Type 2 diabetes all-cause mortality (ACM) and myocardial infarction (MI) glycemic legacy effects have not been explained. We examined their relationships with prior individual HbA1c values and explored the potential impact of instituting earlier, compared with delayed, glucose-lowering therapy. RESEARCH DESIGN AND METHODS: Twenty-year ACM and MI hazard functions were estimated from diagnosis of type 2 diabetes in 3,802 UK Prospective Diabetes Study participants. Impact of HbA1c values over time was analyzed by weighting them according to their influence on downstream ACM and MI risks. RESULTS: Hazard ratios for a one percentage unit higher HbA1c for ACM were 1.08 (95% CI 1.07-1.09), 1.18 (1.15-1.21), and 1.36 (1.30-1.42) at 5, 10, and 20 years, respectively, and for MI was 1.13 (1.11-1.15) at 5 years, increasing to 1.31 (1.25-1.36) at 20 years. Imposing a one percentage unit lower HbA1c from diagnosis generated an 18.8% (95% CI 21.1-16.0) ACM risk reduction 10-15 years later, whereas delaying this reduction until 10 years after diagnosis showed a sevenfold lower 2.7% (3.1-2.3) risk reduction. Corresponding MI risk reductions were 19.7% (22.4-16.5) when lowering HbA1c at diagnosis, and threefold lower 6.5% (7.4-5.3%) when imposed 10 years later. CONCLUSIONS: The glycemic legacy effects seen in type 2 diabetes are explained largely by historical HbA1c values having a greater impact than recent values on clinical outcomes. Early detection of diabetes and intensive glucose control from the time of diagnosis is essential to maximize reduction of the long-term risk of glycemic complications.

Evolutionary forces act on promoter length: identification of enriched cis-regulatory elements.

Transcription factors govern gene expression by binding to short DNA sequences called cis-regulatory elements. These sequences are typically located in promoters, which are regions of variable length upstream of the open reading frames of genes. Here, we report that promoter length and gene function are related in yeast, fungi, and plants. In particular, the promoters for stress-responsive genes are in general longer than those of other genes. Essential genes have, on the other hand, relatively short promoters. We utilize these findings in a novel method for identifying relevant cis-regulatory elements in a set of coexpressed genes. The method is shown to generate more accurate results and fewer false positives compared with other common procedures. Our results suggest that genes with complex transcriptional regulation tend to have longer promoters than genes responding to few signals. This phenomenon is present in all investigated species, indicating that evolution adjust promoter length according to gene function. Identification of cis-regulatory elements in Saccharomyces cerevisiae can be done with the web service located at http://enricher.zool.gu.se.

Testing of chromosomal clumping of gene properties.

Clumping of gene properties like expression or mutant phenotypes along chromosomes is commonly detected using completely random null-models where their location is equally likely across the chromosomes. Interpretation of statistical tests based on these assumptions may be misleading if dependencies exist that are unequal between chromosomes or in different chromosomal parts. One such regional dependency is the telomeric effect, observed in several studies of Saccharomyces cerevisiae, under which e.g. essential genes are less likely to reside near the chromosomal ends. In this study we demonstrate that standard randomisation test procedures are of limited applicability in the presence of telomeric effects. Several extensions of such standard tests are here suggested for handling clumping simultaneously with regional differences in essentiality frequencies in sub-telomeric and central gene positions. Furthermore, a general non-homogeneous discrete Markov approach for combining parametrically modelled position dependent probabilities of a dichotomous property with a simple single parameter clumping is suggested. This Markov model is adapted to the observed telomeric effects and then simulations are used to demonstrate properties of the suggested modified randomisation tests. The model is also applied as a direct alternative tool for statistical analysis of the S. cerevisiae genome for clumping of phenotypes.

Genetic basis of arsenite and cadmium tolerance in Saccharomyces cerevisiae.

BACKGROUND: Arsenic and cadmium are widely distributed in nature and pose serious threats to the environment and human health. Exposure to these nonessential toxic metals may result in a variety of human diseases including cancer. However, arsenic and cadmium toxicity targets and the cellular systems contributing to tolerance acquisition are not fully known. RESULTS: To gain insight into metal action and cellular tolerance mechanisms, we carried out genome-wide screening of the Saccharomyces cerevisiae haploid and homozygous diploid deletion mutant collections and scored for reduced growth in the presence of arsenite or cadmium. Processes found to be required for tolerance to both metals included sulphur and glutathione biosynthesis, environmental sensing, mRNA synthesis and transcription, and vacuolar/endosomal transport and sorting. We also identified metal-specific defence processes. Arsenite-specific defence functions were related to cell cycle regulation, lipid and fatty acid metabolism, mitochondrial biogenesis, and the cytoskeleton whereas cadmium-specific defence functions were mainly related to sugar/carbohydrate metabolism, and metal-ion homeostasis and transport. Molecular evidence indicated that the cytoskeleton is targeted by arsenite and that phosphorylation of the Snf1p kinase is required for cadmium tolerance. CONCLUSION: This study has pin-pointed core functions that protect cells from arsenite and cadmium toxicity. It also emphasizes the existence of both common and specific defence systems. Since many of the yeast genes that confer tolerance to these agents have homologues in humans, similar biological processes may act in yeast and humans to prevent metal toxicity and carcinogenesis.

PROPHECY--a yeast phenome database, update 2006.

Connecting genotype to phenotype is fundamental in biomedical research and in our understanding of disease. Phenomics--the large-scale quantitative phenotypic analysis of genotypes on a genome-wide scale--connects automated data generation with the development of novel tools for phenotype data integration, mining and visualization. Our yeast phenomics database PROPHECY is available at http://prophecy.lundberg.gu.se. Via phenotyping of 984 heterozygous diploids for all essential genes the genotypes analysed and presented in PROPHECY have been extended and now include all genes in the yeast genome. Further, phenotypic data from gene overexpression of 574 membrane spanning proteins has recently been included. To facilitate the interpretation of quantitative phenotypic data we have developed a new phenotype display option, the Comparative Growth Curve Display, where growth curve differences for a large number of mutants compared with the wild type are easily revealed. In addition, PROPHECY now offers a more informative and intuitive first-sight display of its phenotypic data via its new summary page. We have also extended the arsenal of data analysis tools to include dynamic visualization of phenotypes along individual chromosomes. PROPHECY is an initiative to enhance the growing field of phenome bioinformatics.

Modelling of zero-inflation improves inference of metagenomic gene count data.

Metagenomics enables the study of gene abundances in complex mixtures of microorganisms and has become a standard methodology for the analysis of the human microbiome. However, gene abundance data is inherently noisy and contains high levels of biological and technical variability as well as an excess of zeros due to non-detected genes. This makes the statistical analysis challenging. In this study, we present a new hierarchical Bayesian model for inference of metagenomic gene abundance data. The model uses a zero-inflated overdispersed Poisson distribution which is able to simultaneously capture the high gene-specific variability as well as zero observations in the data. By analysis of three comprehensive datasets, we show that zero-inflation is common in metagenomic data from the human gut and, if not correctly modelled, it can lead to substantial reductions in statistical power. We also show, by using resampled metagenomic data, that our model has, compared to other methods, a higher and more stable performance for detecting differentially abundant genes. We conclude that proper modelling of the gene-specific variability, including the excess of zeros, is necessary to accurately describe gene abundances in metagenomic data. The proposed model will thus pave the way for new biological insights into the structure of microbial communities.

Searching for neurodegeneration in multiple sclerosis at clinical onset: Diagnostic value of biomarkers.

BACKGROUND: Neurodegeneration occurs during the early stages of multiple sclerosis. It is an essential, devastating part of the pathophysiology. Tools for measuring the degree of neurodegeneration could improve diagnostics and patient characterization. OBJECTIVE: This study aimed to determine the diagnostic value of biomarkers of degeneration in patients with recent clinical onset of suspected multiple sclerosis, and to evaluate these biomarkers for characterizing disease course. METHODS: This cross-sectional study included 271 patients with clinical features of suspected multiple sclerosis onset and was the baseline of a prospective study. After diagnostic investigations, the patients were classified into the following disease groups: patients with clinically isolated syndrome (n = 4) or early relapsing remitting multiple sclerosis (early RRMS; n = 93); patients with relapsing remitting multiple sclerosis with disease durations ≥2 years (established RRMS; n = 39); patients without multiple sclerosis, but showing symptoms (symptomatic controls; n = 89); and patients diagnosed with other diseases (n = 46). In addition, we included healthy controls (n = 51) and patients with progressive multiple sclerosis (n = 23). We analyzed six biomarkers of neurodegeneration: cerebrospinal fluid neurofilament light chain levels; cerebral spinal fluid glial fibrillary acidic protein; cerebral spinal fluid tau; retinal nerve fiber layer thickness; macula volume; and the brain parenchymal fraction. RESULTS: Except for increased cerebral spinal fluid neurofilament light chain levels, median 670 ng/L (IQR 400-2110), we could not find signs of early degeneration in the early disease group with recent clinical onset. However, the intrathecal immunoglobin G production and cerebral spinal fluid neurofilament light chain levels showed diagnostic value. Moreover, elevated levels of cerebral spinal fluid glial fibrillary acidic protein, thin retinal nerve fiber layers, and low brain parenchymal fractions were associated with progressive disease, but not with the other phenotypes. Thin retinal nerve fiber layers and low brain parenchymal fractions, which indicated neurodegeneration, were associated with longer disease duration. CONCLUSIONS: In clinically suspected multiple sclerosis, intrathecal immunoglobin G production and neurofilament light chain levels had diagnostic value. Therefore, these biomarkers could be included in diagnostic work-ups for multiple sclerosis. We found that the thickness of the retinal nerve fiber layer and the brain parenchymal fraction were not different between individuals that were healthy, symptomatic, or newly diagnosed with multiple sclerosis. This finding suggested that neurodegeneration had not reached a significant magnitude in patients with a recent clinical onset of multiple sclerosis.

Weighted analysis of general microarray experiments.

BACKGROUND: In DNA microarray experiments, measurements from different biological samples are often assumed to be independent and to have identical variance. For many datasets these assumptions have been shown to be invalid and typically lead to too optimistic p-values. A method called WAME has been proposed where a variance is estimated for each sample and a covariance is estimated for each pair of samples. The current version of WAME is, however, limited to experiments with paired design, e.g. two-channel microarrays. RESULTS: The WAME procedure is extended to general microarray experiments, making it capable of handling both one- and two-channel datasets. Two public one-channel datasets are analysed and WAME detects both unequal variances and correlations. WAME is compared to other common methods: fold-change ranking, ordinary linear model with t-tests, LIMMA and weighted LIMMA. The p-value distributions are shown to differ greatly between the examined methods. In a resampling-based simulation study, the p-values generated by WAME are found to be substantially more correct than the alternatives when a relatively small proportion of the genes is regulated. WAME is also shown to have higher power than the other methods. WAME is available as an R-package. CONCLUSION: The WAME procedure is generalized and the limitation to paired-design microarray datasets is removed. The examined other methods produce invalid p-values in many cases, while WAME is shown to produce essentially valid p-values when a relatively small proportion of genes is regulated. WAME is also shown to have higher power than the examined alternative methods.

Quantitative transcriptome, proteome, and sulfur metabolite profiling of the Saccharomyces cerevisiae response to arsenite.

Arsenic is ubiquitously present in nature, and various mechanisms have evolved enabling cells to evade toxicity and acquire tolerance. Herein, we explored how Saccharomyces cerevisiae (budding yeast) respond to trivalent arsenic (arsenite) by quantitative transcriptome, proteome, and sulfur metabolite profiling. Arsenite exposure affected transcription of genes encoding functions related to protein biosynthesis, arsenic detoxification, oxidative stress defense, redox maintenance, and proteolytic activity. Importantly, we observed that nearly all components of the sulfate assimilation and glutathione biosynthesis pathways were induced at both gene and protein levels. Kinetic metabolic profiling evidenced a significant increase in the pools of sulfur metabolites as well as elevated cellular glutathione levels. Moreover, the flux in the sulfur assimilation pathway as well as the glutathione synthesis rate strongly increased with a concomitant reduction of sulfur incorporation into proteins. By combining comparative genomics and molecular analyses, we pinpointed transcription factors that mediate the core of the transcriptional response to arsenite. Taken together, our data reveal that arsenite-exposed cells channel a large part of assimilated sulfur into glutathione biosynthesis, and we provide evidence that the transcriptional regulators Yap1p and Met4p control this response in concert.

Sensitive and robust gene expression changes in fish exposed to estrogen--a microarray approach.

BACKGROUND: Vitellogenin is a well established biomarker for estrogenic exposure in fish. However, effects on gonadal differentiation at concentrations of estrogen not sufficient to give rise to a measurable vitellogenin response suggest that more sensitive biomarkers would be useful. Induction of zona pellucida genes may be more sensitive but their specificities are not as clear. The objective of this study was to find additional sensitive and robust candidate biomarkers of estrogenic exposure. RESULTS: Hepatic mRNA expression profiles were characterized in juvenile rainbow trout exposed to a measured concentration of 0.87 and 10 ng ethinylestradiol/L using a salmonid cDNA microarray. The higher concentration was used to guide the subsequent identification of generally more subtle responses at the low concentration not sufficient to induce vitellogenin. A meta-analysis was performed with data from the present study and three similar microarray studies using different fish species and platforms. Within the generated list of presumably robust responses, several well-known estrogen-regulated genes were identified. Two genes, confirmed by quantitative RT-PCR (qPCR), fulfilled both the criteria of high sensitivity and robustness; the induction of the genes encoding zona pellucida protein 3 and a nucleoside diphosphate kinase (nm23). CONCLUSION: The cross-species, cross-platform meta-analysis correctly identified several robust responses. This adds confidence to our approach used for identifying candidate biomarkers. Specifically, we propose that analyses of an nm23 gene together with zona pellucida genes may increase the possibilities to detect an exposure to low levels of estrogenic compounds in fish.

Quality optimised analysis of general paired microarray experiments.

In microarray experiments, several steps may cause sub-optimal quality and the need for quality control is strong. Often the experiments are complex, with several conditions studied simultaneously. A linear model for paired microarray experiments is proposed as a generalisation of the paired two-sample method by Kristiansson et al. (2005). Quality variation is modelled by different variance scales for different (pairs of) arrays, and shared sources of variation are modelled by covariances between arrays. The gene-wise variance estimates are moderated in an empirical Bayes approach. Due to correlations all data is typically used in the inference of any linear combination of parameters. Both real and simulated data are analysed. Unequal variances and strong correlations are found in real data, leading to further examination of the fit of the model and of the nature of the datasets in general. The empirical distributions of the test-statistics are found to have a considerably improved match to the null distribution compared to previous methods, which implies more correct p-values provided that most genes are non-differentially expressed. In fact, assuming independent observations with identical variances typically leads to optimistic p-values. The method is shown to perform better than the alternatives in the simulation study.

PROPHECY--a database for high-resolution phenomics.

The rapid recent evolution of the field phenomics--the genome-wide study of gene dispensability by quantitative analysis of phenotypes--has resulted in an increasing demand for new data analysis and visualization tools. Following the introduction of a novel approach for precise, genome-wide quantification of gene dispensability in Saccharomyces cerevisiae we here announce a public resource for mining, filtering and visualizing phenotypic data--the PROPHECY database. PROPHECY is designed to allow easy and flexible access to physiologically relevant quantitative data for the growth behaviour of mutant strains in the yeast deletion collection during conditions of environmental challenges. PROPHECY is publicly accessible at http://prophecy.lundberg.gu.se.

Variability in Metagenomic Count Data and Its Influence on the Identification of Differentially Abundant Genes.

Metagenomics is the study of microorganisms in environmental and clinical samples using high-throughput sequencing of random fragments of their DNA. Since metagenomics does not require any prior culturing of isolates, entire microbial communities can be studied directly in their natural state. In metagenomics, the abundance of genes is quantified by sorting and counting the DNA fragments. The resulting count data are high-dimensional and affected by high levels of technical and biological noise that make the statistical analysis challenging. In this article, we introduce an hierarchical overdispersed Poisson model to explore the variability in metagenomic data. By analyzing three comprehensive data sets, we show that the gene-specific variability varies substantially between genes and is dependent on biological function. We also assess the power of identifying differentially abundant genes and show that incorrect assumptions about the gene-specific variability can lead to unacceptable high rates of false positives. Finally, we evaluate shrinkage approaches to improve the variance estimation and show that the prior choice significantly affects the statistical power. The results presented in this study further elucidate the complex variance structure of metagenomic data and provide suggestions for accurate and reliable identification of differentially abundant genes.