Extracting replicable associations across multiple studies: Empirical Bayes algorithms for controlling the false discovery rate.

In almost every field in genomics, large-scale biomedical datasets are used to report associations. Extracting associations that recur across multiple studies while controlling the false discovery rate is a fundamental challenge. Here, we propose a new method to allow joint analysis of multiple studies. Given a set of p-values obtained from each study, the goal is to identify associations that recur in at least k > 1 studies while controlling the false discovery rate. We propose several new algorithms that differ in how the study dependencies are modeled, and compare them and extant methods under various simulated scenarios. The top algorithm, SCREEN (Scalable Cluster-based REplicability ENhancement), is our new algorithm that works in three stages: (1) clustering an estimated correlation network of the studies, (2) learning replicability (e.g., of genes) within clusters, and (3) merging the results across the clusters. When we applied SCREEN to two real datasets it greatly outperformed the results obtained via standard meta-analysis. First, on a collection of 29 case-control gene expression cancer studies, we detected a large set of consistently up-regulated genes related to proliferation and cell cycle regulation. These genes are both consistently up-regulated across many cancer studies, and are well connected in known gene networks. Second, on a recent pan-cancer study that examined the expression profiles of patients with and without mutations in the HLA complex, we detected a large active module of up-regulated genes that are both related to immune responses and are well connected in known gene networks. This module covers thrice more genes as compared to the original study at a similar false discovery rate, demonstrating the high power of SCREEN. An implementation of SCREEN is available in the supplement.

A hierarchical Bayesian model for flexible module discovery in three-way time-series data.

MOTIVATION: Detecting modules of co-ordinated activity is fundamental in the analysis of large biological studies. For two-dimensional data (e.g. genes × patients), this is often done via clustering or biclustering. More recently, studies monitoring patients over time have added another dimension. Analysis is much more challenging in this case, especially when time measurements are not synchronized. New methods that can analyze three-way data are thus needed. RESULTS: We present a new algorithm for finding coherent and flexible modules in three-way data. Our method can identify both core modules that appear in multiple patients and patient-specific augmentations of these core modules that contain additional genes. Our algorithm is based on a hierarchical Bayesian data model and Gibbs sampling. The algorithm outperforms extant methods on simulated and on real data. The method successfully dissected key components of septic shock response from time series measurements of gene expression. Detected patient-specific module augmentations were informative for disease outcome. In analyzing brain functional magnetic resonance imaging time series of subjects at rest, it detected the pertinent brain regions involved. AVAILABILITY AND IMPLEMENTATION: R code and data are available at http://acgt.cs.tau.ac.il/twigs/.

Is this the right normalization? A diagnostic tool for ChIP-seq normalization.

BACKGROUND: Chip-seq experiments are becoming a standard approach for genome-wide profiling protein-DNA interactions, such as detecting transcription factor binding sites, histone modification marks and RNA Polymerase II occupancy. However, when comparing a ChIP sample versus a control sample, such as Input DNA, normalization procedures have to be applied in order to remove experimental source of biases. Despite the substantial impact that the choice of the normalization method can have on the results of a ChIP-seq data analysis, their assessment is not fully explored in the literature. In particular, there are no diagnostic tools that show whether the applied normalization is indeed appropriate for the data being analyzed. RESULTS: In this work we propose a novel diagnostic tool to examine the appropriateness of the estimated normalization procedure. By plotting the empirical densities of log relative risks in bins of equal read count, along with the estimated normalization constant, after logarithmic transformation, the researcher is able to assess the appropriateness of the estimated normalization constant. We use the diagnostic plot to evaluate the appropriateness of the estimates obtained by CisGenome, NCIS and CCAT on several real data examples. Moreover, we show the impact that the choice of the normalization constant can have on standard tools for peak calling such as MACS or SICER. Finally, we propose a novel procedure for controlling the FDR using sample swapping. This procedure makes use of the estimated normalization constant in order to gain power over the naive choice of constant (used in MACS and SICER), which is the ratio of the total number of reads in the ChIP and Input samples. CONCLUSIONS: Linear normalization approaches aim to estimate a scale factor, r, to adjust for different sequencing depths when comparing ChIP versus Input samples. The estimated scaling factor can easily be incorporated in many peak caller algorithms to improve the accuracy of the peak identification. The diagnostic plot proposed in this paper can be used to assess how adequate ChIP/Input normalization constants are, and thus it allows the user to choose the most adequate estimate for the analysis.

repfdr: a tool for replicability analysis for genome-wide association studies.

MOTIVATION: Identification of single nucleotide polymorphisms that are associated with a phenotype in more than one study is of great scientific interest in the genome-wide association studies (GWAS) research. The empirical Bayes approach for discovering whether results have been replicated across studies was shown to be a reliable method, and close to optimal in terms of power. RESULTS: The R package repfdr provides a flexible implementation of the empirical Bayes approach for replicability analysis and meta-analysis, to be used when several studies examine the same set of null hypotheses. The usefulness of the package for the GWAS community is discussed. AVAILABILITY AND IMPLEMENTATION: The R package repfdr can be downloaded from CRAN. CONTACT: ruheller@gmail.com SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Opposing effects of thyroid hormones on cancer risk: a population-based study.

OBJECTIVE: The association between dysregulated thyroid hormone function and cancer risk is inconclusive, especially among different age groups and uncommon malignancies. We sought to determine the relation of TSH and free T4 levels with overall cancer risk as well as risk of specific cancer types. DESIGN AND METHODS: Data on thyroid hormone profile was collected from 375 635 Israeli patients with no prior history of cancer. Cancer cases were identified via the Israel National Cancer Registry. Cox proportional hazards model was used to assess hazard ratios for overall cancer as well as 20 cancer subgroups. RESULTS: In this study, 23 808 cases of cancer were detected over median follow up of 10.9 years. Among patients younger than 50 at inclusion, TSH in the hyperthyroid range, elevated free T4 and subclinical hyperthyroidism were associated with increased cancer risk (HR: 1.3, 1.28 and 1.31, respectively). In contrast, patients 50 or older with clinical hyperthyroidism were at lower cancer risk (HR: 0.64). Elevated TSH was associated with decreased risk of prostate cancer (HR: 0.67). Log-TSH elevation was associated with decreased risk of thyroid cancer (HR: 0.82) and increased risk of melanoma (HR: 1.11) and uterine cancer (HR: 1.27). Elevated free T4 was associated with increased lung cancer risk (HR: 1.54), while free T4 levels above the normal range and clinical hyperthyroidism were related to lower colorectal cancer risk (HR: 0.59 and 0.08, respectively). CONCLUSIONS: Thyroid hormones display opposing effects on cancer risk, based on patient age and cancer type.

Pre-diagnosis thyroid hormone dysfunction is associated with cancer mortality.

Research on the association between thyroid hormone levels and cancer mortality remains limited and inconclusive. We determined the relation of thyroid stimulating hormone (TSH), free T4 (FT4), and free T3 (FT3) levels with mortality in overall cancer and specific tumor types. Thyroid hormone levels 1-5 years prior to cancer diagnosis, as well as multiple clinical and demographic parameters, were retrospectively collected for 10,325 Israeli cancer patients, diagnosed between 2000 and 2016. Patients treated with thyroid altering medications were excluded. Cancer diagnosis was determined via the Israel National Cancer Registry. Multivariate-adjusted Cox proportional hazards model was used to assess the hazard ratios (HRs) based on thyroid hormone function for cancer mortality. A total of 5265 patients died during the follow-up period (median of 4.4 years). TSH, FT4, and FT3 levels in the hypothyroid range were associated with increase in overall mortality (adjusted HR 1.20, 1.74, 1.87, respectively). We further analyzed the association between TSH and mortality in 14 cancer subgroups. Specifically, TSH in both the hyperthyroid and hypothyroid range was associated with melanoma mortality (adjusted HR 2.20, 4.47, respectively). In conclusion, pre-diagnosis of thyroid dysfunction is associated with increased cancer mortality, a relation likely driven by specific cancer types. These findings suggest that thyroid hormones may potentially serve as prognostic markers in cancer.

Computational characterization of B-cell epitopes.

Characterizing B-cell epitopes is a fundamental step for understanding the immunological basis of bio-recognition. To date, epitope analyses have either been based on limited structural data, or sequence data alone. In this study, our null hypothesis was that the surface of the antigen is homogeneously antigenic. To test this hypothesis, a large dataset of antibody-antigen complex structures, together with crystal structures of the native antigens, has been compiled. Computational methods were developed and applied to detect and extract physico-chemical, structural, and geometrical properties that may distinguish an epitope from the remaining antigen surface. Rigorous statistical inference was able to clearly reject the null hypothesis showing that epitopes are distinguished from the remaining antigen surface in properties such as amino acid preference, secondary structure composition, geometrical shape, and evolutionary conservation. Specifically, epitopes were found to be significantly enriched with tyrosine and tryptophan, and to show a general preference for charged and polar amino acids. Additionally, epitopes were found to show clear preference for residing on planar parts of the antigen that protrude from the surface, yet with a rugged surface shape at the atom level. The effects of complex formation on the structural properties of the antigen were also computationally characterized and it is shown that epitopes undergo compression upon antibody binding. This correlates with the finding that epitopes are enriched with unorganized secondary structure elements that render them flexible. Thus, this study extends the understanding of the underlying processes required for antibody binding, and reveals new aspects of the antibody-antigen interaction.

Getting more out of meta-analyses: a new approach to meta-analysis in light of unexplained heterogeneity.

BACKGROUND AND OBJECTIVES: Meta-analyses sometimes summarize results in the presence of substantial unexplained between-study heterogeneity. As GRADE criteria highlight, unexplained heterogeneity reduces certainty in the evidence, resulting in limited confidence in average effect estimates. The aim of this paper is to provide a new clinically useful approach to estimating an intervention effect in light of unexplained heterogeneity. METHODS: We used a random-effects model to estimate the distribution of an intervention-effect across various groups of patients given data derived from meta-analysis. The model provides a distribution of the probabilities of various possible effects in a new group of patients. We examined how our method influenced the conclusions of two meta-analyses. RESULTS: In one example, our method illustrated that evidence from a meta-analysis did not support authors' highly publicized conclusion that hypericum is as effective as other antidepressants. In the second example, our method provided insight into a subgroup analysis of the effect of ribavirin in hepatitis C, demonstrating clear important benefit in one subgroup but not in others. CONCLUSION: Analysing the distribution of an intervention-effect in random-effects models may enable clinicians to improve their understanding of the probability of particular-intervention effects in a new population.

Associating quantitative behavioral traits with gene expression in the brain: searching for diamonds in the hay.

UNLABELLED: Gene expression and phenotypic functionality can best be associated when they are measured quantitatively within the same experiment. The analysis of such a complex experiment is presented, searching for associations between measures of exploratory behavior in mice and gene expression in brain regions. The analysis of such experiments raises several methodological problems. First and foremost, the size of the pool of potential discoveries being screened is enormous yet only few biologically relevant findings are expected, making the problem of multiple testing especially severe. We present solutions based on screening by testing related hypotheses, then testing the hypotheses of interest. In one variant the subset is selected directly, in the other one a tree of hypotheses is tested hierarchical; both variants control the False Discovery Rate (FDR). Other problems in such experiments are in the fact that the level of data aggregation may be different for the quantitative traits (one per animal) and gene expression measurements (pooled across animals); in that the association may not be linear; and in the resolution of interest only few replications exist. We offer solutions to these problems as well. The hierarchical FDR testing strategies presented here can serve beyond the structure of our motivating example study to any complex microarray study. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Genomic analysis of COP9 signalosome function in Drosophila melanogaster reveals a role in temporal regulation of gene expression.

The COP9 signalosome (CSN), an eight-subunit protein complex, is conserved in all higher eukaryotes. CSN intersects the ubiquitin-proteasome pathway, modulating signaling pathways controlling various aspects of development. We are using Drosophila as a model system to elucidate the function of this important complex. Transcriptome data were generated for four csn mutants, sampled at three developmental time points. Our results are highly reproducible, being confirmed using two different experimental setups that entail different microarrays and different controls. Our results indicate that the CSN acts as a transcriptional repressor during development of Drosophila, resulting in achronic gene expression in the csn mutants. 'Time shift' analysis with the publicly available Drosophila transcriptome data indicates that genes repressed by the CSN are normally induced primarily during late embryogenesis or during metamorphosis. These temporal shifts are likely due to the roles of the CSN in regulating transcription factors. A null mutation in CSN subunit 4 and hypomorphic mutations in csn5 lead to more severe defects than seen in the csn5-null mutants strain, suggesting that CSN5 carries only some of the CSN function.

Data mining in a behavioral test detects early symptoms in a model of amyotrophic lateral sclerosis.

"What's wrong with my genetically engineered animal?" is a common yet often difficult to answer question in behavioral phenotyping. We present here a method termed Pattern Array for mining movement patterns and isolating those that best capture an effect of a genetic manipulation. We demonstrate the method by searching for early motor symptoms in the open-field behavior of SOD1 mutant rats, an animal model of amyotrophic lateral sclerosis. Pattern Array was able to identify a unique motor pattern that differentiated the SOD1 mutants from the wild-type controls 2 months before disease onset. This pattern included heavy braking while moving near the arena wall but turning away from it. SOD1 mutants performed this pattern significantly less than wild-type controls in 2 independent data sets. At such early age the SOD1 mutants could not be differentiated from the controls by standard behavioral measures or by subjective observation. The early discovered symptom may enable investigators to test therapies aimed for intervention rather than remediation. Our results demonstrate the feasibility and potential of detecting subtle behavioral effects using data mining strategies.

SKY analysis of childhood neural tumors and cell lines demonstrates a susceptibility of aberrant chromosomes to further rearrangements.

Malignant solid tumors are commonly characterized by a large number of complex structural and numerical chromosomal alterations, which often reflect the level of genomic instability and can be associated with disease progression. The aim of this study was to evaluate whether chromosomes that harbor primary aberrations have a higher susceptibility to accumulate further alterations. We used spectral karyotyping (SKY), to compare the individual chromosomal instability of two chromosome types: chromosomes that have a primary aberration and chromosomes without an aberration, in 13 primary childhood neural tumors and seven cell lines. We found that chromosomes that contain a primary aberration are significantly (p-value<0.001) more likely to gain further structural rearrangements or to undergo numerical changes (22.6%, 36 of 159 chromosomes) than chromosomes with no initial aberration (4.9%, 54 of 1099 chromosomes). These results are highly suggestive that aberrant chromosomes in solid tumors have a higher susceptibility to accumulate further rearrangements than "normal" chromosomes.

Prior exercise training improves the outcome of acute myocardial infarction in the rat. Heart structure, function, and gene expression.

OBJECTIVES: The aim of this research was to investigate the structural, functional, and molecular features of the remodeling heart in prior swim-trained infarcted rats. BACKGROUND: Physical exercise training is a known protective factor against cardiovascular morbidity and mortality. The structural and molecular aspects underlying this protection in the remodeling heart have not been investigated. METHODS: After seven weeks of swimming exercise training, rats underwent surgical ligation of the left coronary artery followed by a four-week sedentary period. Untrained control rats underwent the same surgical protocol. Left ventricular function was assessed by echocardiography four weeks after infarction, and hearts were sampled for histological and molecular analysis. Ribonucleic acid from the surviving left ventricle was analyzed by complementary deoxyribonucleic acid arrays followed by Northern blotting or quantitative reverse transcription polymerase chain reaction of selected messenger ribonucleic acids (mRNAs). RESULTS: Scar area was 1.6-fold smaller (p = 0.0002), arteriolar density was 1.7-fold higher (p = 0.0002), and left ventricular shortening fraction was 1.9-fold higher (p = 0.003) in the exercise-trained compared with sedentary hearts. Eleven genes whose expression level varied by at least +/-1.5-fold distinguished the prior exercised rats from their sedentary counterparts. Compared with sedentary, the exercised hearts displayed 9- and 2.4-times lower levels of atrial natriuretic peptide and aldolase mRNA (p = 0.03 and 0.04, respectively), and a 2.7- and 1.9-fold higher abundance of cytochrome c-oxidase and fatty acid binding protein, respectively (p < 0.03, each). CONCLUSIONS: Swimming exercise training before acute myocardial infarction reduces scar size, increases arteriole density, and manifests adaptation of stress- and energy-metabolism-related genes that may contribute to the improved heart function observed during remodeling.

Quantitative trait Loci analysis using the false discovery rate.

False discovery rate control has become an essential tool in any study that has a very large multiplicity problem. False discovery rate-controlling procedures have also been found to be very effective in QTL analysis, ensuring reproducible results with few falsely discovered linkages and offering increased power to discover QTL, although their acceptance has been slower than in microarray analysis, for example. The reason is partly because the methodological aspects of applying the false discovery rate to QTL mapping are not well developed. Our aim in this work is to lay a solid foundation for the use of the false discovery rate in QTL mapping. We review the false discovery rate criterion, the appropriate interpretation of the FDR, and alternative formulations of the FDR that appeared in the statistical and genetics literature. We discuss important features of the FDR approach, some stemming from new developments in FDR theory and methodology, which deem it especially useful in linkage analysis. We review false discovery rate-controlling procedures--the BH, the resampling procedure, and the adaptive two-stage procedure-and discuss the validity of these procedures in single- and multiple-trait QTL mapping. Finally we argue that the control of the false discovery rate has an important role in suggesting, indicating the significance of, and confirming QTL and present guidelines for its use.

The accuracy of Gram stain of respiratory specimens in excluding Staphylococcus aureus in ventilator-associated pneumonia.

OBJECTIVE: To evaluate the Gram stain of deep tracheal aspirate as a tool to direct empiric antibiotic therapy, and more specifically as a tool to exclude the need for empiric antibiotic coverage against Staphylococcus aureus in ventilator-associated pneumonia (VAP). DESIGN: A prospective, single-center, observational, cohort study. SETTING: All wards at a community hospital. PATIENTS: Adult patients requiring mechanical ventilation, identified as having VAP in a 54-month prospective surveillance database. INTERVENTIONS: Sampling of lower airway secretions by deep endotracheal aspiration was taken from each patient who developed VAP. Samples were sent immediately for Gram stain and qualitative bacterial cultures. Demographic and relevant clinical data were collected; Gram stain, culture, and antibiotic susceptibility results were documented; and outcome was followed prospectively. MEASUREMENTS AND MAIN RESULTS: The analysis included 114 consecutive patients with 115 episodes of VAP from June 2007 to January 2012. Sensitivity of Gram stain compared with culture was 90.47% for gram-positive cocci, 69.6% for gram-negative rods, and 50% for sterile cultures. Specificity was 82.5%, 77.8%, and 79%, respectively. Negative predictive value was high for gram-positive cocci (97%) and sterile cultures (96%) but low for gram-negative rods (20%). Acinetobacter baumanii (45%) and Pseudomonas aeruginosa (38 %) were the prevailing isolates. S aureus was found in 18.3% of the patients. Most isolates were multiresistant. CONCLUSIONS: Absence of gram-positive bacteria on Gram stain had a high negative predictive value. These data can be used to narrow the initial empiric antibiotic regimen and to avoid unnecessary exposure of patients to vancomycin and other antistaphyloccocal agents.

Embryo quality and implantation rates are not influenced by total motile count values in an ICSI programme: a novel point of view.

Total motile count (TMC) is a useful tool for sperm evaluation, comprising both quantitative and motility parameters. Although frequently used, TMC has not yet been evaluated as a contributory variable for intracytoplasmic sperm injection (ICSI) cycles. In this study we evaluate the possible role of TMC as a prognostic parameter in cycles designated for ICSI. We also test the existence of a possible TMC-threshold value that might be predictive for ICSI cycle outcome in the everyday practice. This is a retrospective cohort study in which the research question is addressed by a locally weighted regression (LOESS) analysis. Primary outcome measures are fertilization rate, good quality embryos rate and implantation rate. A total of 666 patients were included, contributing 1456 cycles. The effect of TMC over the fertilization rate was significant, depicting an inverted U-shaped curve: with up to approximately 10 million motile sperm, fertilization rates increased as TMC increased, but from this point on decreased. A slight increment in the rate of good embryo formation with increasing value of TMC was noted, but this did not reach a statistical significance. TMC values demonstrated no effect in the case of implantation rates. ICSI may offer an advantage related to fertilization rates for the sub-fertile male population, with a motile sperm count up to 10 million.

Selective inference in complex research.

We explain the problem of selective inference in complex research using a recently published study: a replicability study of the associations in order to reveal and establish risk loci for type 2 diabetes. The false discovery rate approach to such problems will be reviewed, and we further address two problems: (i) setting confidence intervals on the size of the risk at the selected locations and (ii) selecting the replicable results.

A data mining approach to in vivo classification of psychopharmacological drugs.

Data mining is a powerful bioinformatics strategy that has been successfully applied in vitro to screen for gene-expression profiles predicting toxicological or carcinogenic response ('class predictors'). In this report we used a data mining algorithm named Pattern Array (PA) in vivo to analyze mouse open-field behavior and characterize the psychopharmacological effects of three drug classes--psychomotor stimulant, opioid, and psychotomimetic. PA represents rodent movement with approximately 100,000 complex patterns, defined as multiple combinations of several ethologically relevant variables, and mines them for those that maximize any effect of interest, such as the difference between drug classes. We show that PA can discover behavioral predictors of all three drug classes, thus developing a reliable drug-classification scheme in small group sizes. The discovered predictors showed orderly dose dependency despite being explicitly mined only for class differences, with the high doses scoring 4-10 standard deviations from the vehicle group. Furthermore, these predictors correctly classified in a dose-dependent manner four 'unknown' drugs (ie that were not used in the training process), and scored a mixture of a psychomotor stimulant and an opioid as being intermediate between these two classes. The isolated behaviors were highly heritable (h(2)>50%) and replicable as determined in 10 inbred strains across three laboratories. PA can in principle be applied for mining behaviors predicting additional properties, such as within-class differences between drugs and within-drug dose-response, all of which can be measured automatically in a single session per animal in an open-field arena, suggesting a high potential as a tool in psychotherapeutic drug discovery.

Cluster-based analysis of FMRI data.

We propose a method for the statistical analysis of fMRI data that tests cluster units rather than voxel units for activation. The advantages of this analysis over previous ones are both conceptual and statistical. Recognizing that the fundamental units of interest are the spatially contiguous clusters of voxels that are activated together, we set out to approximate these cluster units from the data by a clustering algorithm especially tailored for fMRI data. Testing the cluster units has a two-fold statistical advantage over testing each voxel separately: the signal to noise ratio within the unit tested is higher, and the number of hypotheses tests compared is smaller. We suggest controlling FDR on clusters, i.e., the proportion of clusters rejected erroneously out of all clusters rejected and explain the meaning of controlling this error rate. We introduce the powerful adaptive procedure to control the FDR on clusters. We apply our cluster-based analysis (CBA) to both an event-related and a block design fMRI vision experiment and demonstrate its increased power over voxel-by-voxel analysis in these examples as well as in simulations.

Identifying differentially expressed genes using false discovery rate controlling procedures.

MOTIVATION: DNA microarrays have recently been used for the purpose of monitoring expression levels of thousands of genes simultaneously and identifying those genes that are differentially expressed. The probability that a false identification (type I error) is committed can increase sharply when the number of tested genes gets large. Correlation between the test statistics attributed to gene co-regulation and dependency in the measurement errors of the gene expression levels further complicates the problem. In this paper we address this very large multiplicity problem by adopting the false discovery rate (FDR) controlling approach. In order to address the dependency problem, we present three resampling-based FDR controlling procedures, that account for the test statistics distribution, and compare their performance to that of the naïve application of the linear step-up procedure in Benjamini and Hochberg (1995). The procedures are studied using simulated microarray data, and their performance is examined relative to their ease of implementation. RESULTS: Comparative simulation analysis shows that all four FDR controlling procedures control the FDR at the desired level, and retain substantially more power then the family-wise error rate controlling procedures. In terms of power, using resampling of the marginal distribution of each test statistics substantially improves the performance over the naïve one. The highest power is achieved, at the expense of a more sophisticated algorithm, by the resampling-based procedures that resample the joint distribution of the test statistics and estimate the level of FDR control. AVAILABILITY: An R program that adjusts p-values using FDR controlling procedures is freely available over the Internet at www.math.tau.ac.il/~ybenja.