Volume 5, Issue 1

As a powerful tool for analyzing full conditional (in-)dependencies between random variables, graphical models have become increasingly popular to infer genetic networks based on gene expression data. However, full (unconstrained) conditional relationships between random variables can be only estimated accurately if the number of observations is relatively large in comparison to the number of variables, which is usually not fulfilled for high-throughput genomic data.Recently, simplified graphical modeling approaches have been proposed to determine dependencies between gene expression profiles. For sparse graphical models such as genetic networks, it is assumed that the zero- and first-order conditional independencies still reflect reasonably well the full conditional independence structure between variables. Moreover, low-order conditional independencies have the advantage that they can be accurately estimated even when having only a small number of observations. Therefore, using only zero- and first-order conditional dependencies to infer the complete graphical model can be very useful. Here, we analyze the statistical and probabilistic properties of these low-order conditional independence graphs (called 0-1 graphs). We find that for faithful graphical models, the 0-1 graph contains at least all edges of the full conditional independence graph (concentration graph). For simple structures such as Markov trees, the 0-1 graph even coincides with the concentration graph. Furthermore, we present some asymptotic results and we demonstrate in a simulation study that despite their simplicity, 0-1 graphs are generally good estimators of sparse graphical models. Finally, the biological relevance of some applications is summarized.

We think of cluster analysis as class discovery. That is, we assume that there is an unknown mapping called clustering structure that assigns a class label to each observation, and the goal of cluster analysis is to estimate this clustering structure, that is, to estimate the number of clusters and cluster assignments. In traditional cluster analysis, it is assumed that such unknown mapping is unique. However, since the observations may cluster in more than one way depending on the variables used, it is natural to permit the existence of more than one clustering structure. This generalized clustering problem of estimating multiple clustering structures is the focus of this paper. We propose an algorithm for finding multiple clustering structures of observations which involves clustering both variables and observations. The number of clustering structures is determined by the number of variable clusters. The dissimilarity measure for clustering variables is based on nearest-neighbor graphs. The observations are clustered using weighted distances with weights determined by the clusters of the variables. The motivating application is to gene expression data.

This paper develops a Bayes regression model having change points for the analysis of array-CGH data by utilizing not only the underlying spatial structure of the genomic alterations but also the observation that the noise associated with the ratio of the fluorescence intensities is bigger when the intensities get smaller. We show that this Bayes regression approach is particularly suitable for the analysis of cDNA microarray-CGH data, which are generally noisier than those using genomic clones. A simulation study and a real data analysis are included to illustrate this approach.

Maintenance genes can be used for normalization in the comparison of gene expressions. Even though the absolute expression levels of maintenance genes may vary considerably among different tissues or cells, a set of maintenance genes may provide suitable normalization if their expression levels are relatively constant in the specific tissues or cells of interest. A statistical procedure is proposed to select maintenance genes for normalization of gene expression data from tissues or cells of interest. This procedure is based on simultaneous confidence intervals for practical equivalence of relative gene expressions in these tissues or cells. As an illustration, the procedure is applied to the maintenance gene expression data from Vandesompele et al. (2002).

Experiments to determine the complete 3-dimensional structures of protein complexes are difficult to perform and only a limited range of such structures are available.In contrast, large-scale screening experiments have identified thousands of pairwise interactions between proteins, but such experiments do not produce explicit structural information.In addition, the data produced by these high through-put experiments contain large numbers of false positive results, and can be biased against detection of certain types of interaction.Several methods exist that analyse such pairwise interaction data in terms of the constituent domains within proteins, scoring pairs of domain superfamilies according to their propensity to interact.These scores can be used to predict the strongest domain-domain contact (the contact with the largest surface area) between interacting proteins for which the domain-level structures of the individual proteins are known.We test this predictive approach on a set of pairwise protein interactions taken from the Protein Quaternary Structure (PQS) database for which the true domain-domain contacts are known.While the overall prediction success rate across the whole test data set is poor, we shown how interactions in the test data set for which the training data are not informative can be automatically excluded from the prediction process, giving improved prediction success rates at the expense of restricted coverage of the test data.

An important application of gene expression microarray data is classification of biological samples or prediction of clinical and other outcomes. One necessary part of multivariate statistical analysis in such applications is dimension reduction. This paper provides a comparison study of three dimension reduction techniques, namely partial least squares (PLS), sliced inverse regression (SIR) and principal component analysis (PCA), and evaluates the relative performance of classification procedures incorporating those methods. A five-step assessment procedure is designed for the purpose. Predictive accuracy and computational efficiency of the methods are examined. Two gene expression data sets for tumor classification are used in the study.

Modern methods of microarray data analysis are biased towards selecting those genes that display the most pronounced differential expression. The magnitude of differential expression does not necessarily indicate biological significance and other criteria are needed to supplement the information on differential expression. Three large sets of microarray data on childhood leukemia were analyzed by an original method introduced in this paper. A new type of stochastic dependence between expression levels in gene pairs was deciphered by our analysis. This modulation-like unidirectional dependence between expression signals arises when the expression of a ``gene-modulator'' is stochastically proportional to that of a ``gene-driver''. A total of more than 35% of all pairs formed from 12550 genes were conservatively estimated to belong to this type. There are genes that tend to form Type A relationships with the overwhelming majority of genes. However, this picture is not static: the composition of Type A gene pairs may undergo dramatic changes when comparing two phenotypes. The ability to identify genes that act as ``modulators'' provides a potential strategy of prioritizing candidate genes.

Recently Peres and Shields discovered a new method for estimating the order of a stationary fixed order Markov chain. They showed that the estimator is consistent by proving a threshold result. While this threshold is valid asymptotically in the limit, it is not very useful for DNA sequence analysis where data sizes are moderate. In this paper we give a novel interpretation of the Peres-Shields estimator as a sharp transition phenomenon. This yields a precise and powerful estimator that quickly identifies the core dependencies in data. We show that it compares favorably to other estimators, especially in the presence of variable dependencies. Motivated by this last point, we extend the Peres-Shields estimator to Variable Length Markov Chains. We compare it to a well-established estimator and show that it is superior in terms of the predictive likelihood. We give an application to the problem of detecting DNA sequence similarity in plasmids.

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.

A new data filtering method for SELDI-TOF MS proteomic spectra data is described. We examined technical repeats (2 per subject) of intensity versus m/z (mass/charge) of bone marrow cell lysate for two groups of childhood leukemia patients: acute myeloid leukemia (AML) and acute lymphoblastic leukemia (ALL). As others have noted, the type of data processing as well as experimental variability can have a disproportionate impact on the list of ``interesting'' proteins (see Baggerly et al. (2004)). We propose a list of processing and multiple testing techniques to correct for 1) background drift; 2) filtering using smooth regression and cross-validated bandwidth selection; 3) peak finding; and 4) methods to correct for multiple testing (van der Laan et al. (2005)). The result is a list of proteins (indexed by m/z) where average expression is significantly different among disease (or treatment, etc.) groups. The procedures are intended to provide a sensible and statistically driven algorithm, which we argue provides a list of proteins that have a significant difference in expression. Given no sources of unmeasured bias (such as confounding of experimental conditions with disease status), proteins found to be statistically significant using this technique have a low probability of being false positives.

In this article, we show how to apply our previously proposed Deletion/Substitution/Addition algorithm in the context of right-censoring for the prediction of survival. Furthermore, we introduce how to incorporate bagging into the algorithm to obtain a cross-validated bagged estimator. The method is used for predicting the survival time of patients with diffuse large B-cell lymphoma based on gene expression variables.

Using affected sibling pairs, the mean allele sharing statistic tests for linkage by testing if the mean proportion of alleles that are identical-by-descent (IBD) is equal to a half. The behavior of some versions of the mean allele sharing test statistic depends on whether or not families that are uninformative for their IBD status are included; the SIBPAL version provides less significant values when all families (informative and uninformative) are used than when only informative families are used. Here, we investigate this behavior both analytically and by simulation. Our investigation shows that the main issue is the choice of the variance estimator in the denominator of the statistic. The choice of the denominator is very important and is still not totally resolved. Our mathematical explanation supported by our simulation study might aid in the search for an optimum solution.

Simultaneously testing a collection of null hypotheses about a data generating distribution based on a sample of independent and identically distributed observations is a fundamental and important statistical problem involving many applications. Methods based on marginal null distributions (i.e., marginal p-values) are attractive since the marginal p-values can be based on a user supplied choice of marginal null distributions and they are computationally trivial, but they, by necessity, are known to either be conservative or to rely on assumptions about the dependence structure between the test-statistics. Re-sampling based multiple testing (Westfall and Young, 1993) involves sampling from a joint null distribution of the test-statistics, and controlling (possibly in a, for example, step-down fashion) the user supplied type-I error rate under this joint null distribution for the test-statistics. A generally asymptotically valid null distribution avoiding the need for the subset pivotality condition for the vector of test-statistics was proposed in Pollard, van der Laan (2003) for null hypotheses about general real valued parameters. This null distribution was generalized in Dudoit, vanderLaan, Pollard (2004) to general null hypotheses and test-statistics. In ongoing recent work van der Laan, Hubbard (2005), we propose a new generally asymptotically valid null distribution for the test-statistics and a corresponding bootstrap estimate, whose marginal distributions are user supplied, and can thus be set equal to the (most powerful) marginal null distributions one would use in univariate testing to obtain a p-value. Previous proposed null distributions either relied on a restrictive subset pivotality condition (Westfall and Young) or did not guarantee this latter property (Dudoit, vanderLaan, Pollard, 2004). It is argued and illustrated that the resulting new re-sampling based multiple testing methods provide more accurate control of the wished Type-I error in finite samples and are more powerful. We establish formal results and investigate the practical performance of this methodology in a simulation and data analysis.

As technology for microarray analysis becomes widespread, it is becoming increasingly important to be able to compare and combine the results of experiments that explore the same scientific question. In this article, we present a rank-aggregation approach for combining results from several microarray studies. The motivation for this approach is twofold; first, the final results of microarray studies are typically expressed as lists of genes, rank-ordered by a measure of the strength of evidence that they are functionally involved in the disease process, and second, using the information on this rank-ordered metric means that we do not have to concern ourselves with data on the actual expression levels, which may not be comparable across experiments. Our approach draws on methods for combining top-k lists from the computer science literature on meta-search. The meta-search problem shares several important features with that of combining microarray experiments, including the fact that there are typically few lists with many elements and the elements may not be common to all lists. We implement two meta-search algorithms, which use a Markov chain framework to convert pairwise preferences between list elements into a stationary distribution that represents an aggregate ranking (Dwork et al, 2001). We explore the behavior of the algorithms in hypothetical examples and a simulated dataset and compare their performance with that of an algorithm based on the order-statistics model of Thurstone (Thurstone, 1927). We apply all three algorithms to aggregate the results of five microarray studies of prostate cancer.

The evolution of drug resistance in HIV is characterized by the accumulation of resistance-associated mutations in the HIV genome. Mutagenetic trees, a family of restricted Bayesian tree models, have been applied to infer the order and rate of occurrence of these mutations. Understanding and predicting this evolutionary process is an important prerequisite for the rational design of antiretroviral therapies. In practice, mixtures models of K mutagenetic trees provide more flexibility and are often more appropriate for modelling observed mutational patterns.Here, we investigate the model selection problem for K-mutagenetic trees mixture models. We evaluate several classical model selection criteria including cross-validation, the Bayesian Information Criterion (BIC), and the Akaike Information Criterion. We also use the empirical Bayes method by constructing a prior probability distribution for the parameters of a mutagenetic trees mixture model and deriving the posterior probability of the model. In addition to the model dimension, we consider the redundancy of a mixture model, which is measured by comparing the topologies of trees within a mixture model. Based on the redundancy, we propose a new model selection criterion, which is a modification of the BIC.Experimental results on simulated and on real HIV data show that the classical criteria tend to select models with far too many tree components. Only cross-validation and the modified BIC recover the correct number of trees and the tree topologies most of the time. At the same optimal performance, the runtime of the new BIC modification is about one order of magnitude lower. Thus, this model selection criterion can also be used for large data sets for which cross-validation becomes computationally infeasible.

In the field of molecular evolution genome substitution models with neighbour dependent substitution rates have recently received much attention. It is well-known that substitution of nucleotides does not occur independently of neighbouring nucleotides, but there has been less focus on the phenomenon that this substitution process is also not time-reversible. In this paper I construct a pseudo-likelihood type method for inference in non-reversible substitution models with neighbour dependent substitution rates. I also construct an EM-algorithm for maximising the pseudo-likelihood. For human-mouse aligned sequence data a number of different models are investigated, where I show that strand-symmetric models are appropriate, and that overlapping di-nucleotide models do not fit the data well.

Consider the standard multiple testing problem where many hypotheses are to be tested, each hypothesis is associated with a test statistic, and large test statistics provide evidence against the null hypotheses. One proposal to provide probabilistic control of Type-I errors is the use of procedures ensuring that the expected number of false positives does not exceed a user-supplied threshold. Among such multiple testing procedures, we derive the most powerful method, meaning the test statistic cutoffs that maximize the expected number of true positives. Unfortunately, these optimal cutoffs depend on the true unknown data generating distribution, so could never be used in a practical setting. We instead consider splitting the sample so that the optimal cutoffs are estimated from a portion of the data, and then testing on the remaining data using these estimated cutoffs. When the null distributions for all test statistics are the same, the obvious way to control the expected number of false positives would be to use a common cutoff for all tests. In this work, we consider the common cutoff method as a benchmark multiple testing procedure. We show that in certain circumstances the use of estimated optimal cutoffs via sample splitting can dramatically outperform this benchmark method, resulting in increased true discoveries, while retaining Type-I error control. This paper is an updated version of the work presented in Rubin et al. (2005), later expanded upon by Wasserman and Roeder (2006).

Pixel saturation occurs when the pixel intensity exceeds the scanner upper threshold of detection and the recorded pixel intensity is then truncated at the threshold. Truncation of the pixel intensity causes the estimates of gene expression (i.e., intensity) to be biased. Microarray experiments are commonly affected by saturated pixels; as a result all higher level analyses are made on these biased gene expression estimates. In this paper, we propose a method for improving the quality of the signal for cDNA microarrays by making use of several scans at varying scanner sensitivities. For each spot, pixel level intensity readings are given as input to a Bayesian hierarchical model. The model uses the pixel intensities of the spot to provide a posterior distribution of the true expression level of the corresponding genes. The parameters of the hierarchical model are estimated jointly with these expression levels, thus performing an integrated analysis of the measurement data. The method improves in all ranges the accuracy with which intensities can be estimated and extends the dynamic range of measured gene expression at the high end. The method is generic and can be applied to data from any organism and for imaging with any scanner. Results from a real data set illustrate an improved precision in the estimation of the expression of genes compared to what can be achieved by applying standard methods and using only a single scan.

Statistical inference of graphical models has become an important tool in the reconstruction of biological networks of the type which model, for example, gene regulatory interactions. In particular, the construction of a score-based Bayesian posterior density over the space of models provides an intuitive and computationally feasible method of assessing model uncertainty and of assigning statistical confidence to structural features. One problem which frequently occurs with this approach is the tendency to overestimate the degree of model complexity. Spurious graphical features obtained in this way may affect the inference in unpredictable ways, even when using scoring techniques, such as the Bayesian Information Criterion (BIC), that are specifically designed to compensate for overfitting.In this article we propose a simple adjustment to a BIC-based scoring procedure. The method proceeds in two steps. In the first step we derive an independent estimate of the parametric complexity of the model. In the second we modify the BIC score so that the mean parametric complexity of the posterior density is equal to the estimated value. The method is applied to a set of test networks, and to a collection of genes from the yeast genome known to possess regulatory relationships. A Bayesian network model with binary responses is employed. In the examples considered, we find that the number of spurious graph edges inferred is reduced, while the effect on the identification of true edges is minimal.

Genetic epidemiology aims at identifying biological mechanisms responsible for human diseases. Genome-wide association studies, made possible by recent improvements in genotyping technologies, are now promisingly investigated. In these studies, common first-stage strategies focus on marginal effects but lead to multiple-testing and are unable to capture the possibly complex interplay between genetic factors.We have adapted the use of the local score statistic, already successfully applied to analyse long molecular sequences. Via sum statistics, this method captures local and possible distant dependences between markers. Dedicated to genome-wide association studies, it is fast to compute, able to handle large datasets, circumvents the multiple-testing problem and outlines a set of genomic regions (segments) for further analyses. Applied to simulated and real data, our approach outperforms classical Bonferroni and FDR corrections for multiple-testing. It is implemented in a software termed LHiSA for Local High-scoring Segments for Association and available at: http://stat.genopole.cnrs.fr/software/lhisa.

The spectral envelope, a frequency based technique for analyzing categorical time series, is applied to amino acid sequences to examine their periodicity. The periodic signatures of such sequences is related to the secondary structure of the folding patterns in the gene. For a pair of sequences, we define a spectral envelope covariance which emphasizes the common periodicities in the two sequences. This is used to give a similarity measure for the two sequences which can then be used in a neighbour joining algorithm to construct a phylogeny. We apply the spectral methods to myoglobin sequences from primates and cetaceans. The spectral envelope reflects the structure of this protein and the tree constructed using spectral methods shows strong agreement with published trees.The spectral envelope can be used to explore similarities between and within different protein families. Since we do not require aligned sequences, the spectral methods can be used to create phylogenies across different protein families. We apply the method to 11 protein families from PANDIT obtaining a tree where the families are separated and the relationship among the families is given.

There are many methods of correcting microarray data for non-biological sources of error. Authors routinely supply software or code so that interested analysts can implement their methods. Even with a thorough reading of associated references, it is not always clear how requisite parts of the method are calculated in the software packages. However, it is important to have an understanding of such details, as this understanding is necessary for proper use of the output, or for implementing extensions to the model.In this paper, the calculation of parameter estimates used in Robust Multichip Average (RMA), a popular preprocessing algorithm for Affymetrix GeneChip brand microarrays, is elucidated. The background correction method for RMA assumes that the perfect match (PM) intensities observed result from a convolution of the true signal, assumed to be exponentially distributed, and a background noise component, assumed to have a normal distribution. A conditional expectation is calculated to estimate signal. Estimates of the mean and variance of the normal distribution and the rate parameter of the exponential distribution are needed to calculate this expectation. Simulation studies show that the current estimates are flawed; therefore, new ones are suggested. We examine the performance of preprocessing under the exponential-normal convolution model using several different methods to estimate the parameters.

We outline a method of sample size calculation in microarray experiments on the basis of pilot data and illustrate its practical application with both simulated and real data. The method was shown to be consistent (as the number of `probed genes' tends to infinity) under general conditions in an earlier, more `theoretical' companion paper. Its implementation requires the values of test statistics, the sample size with which the statistics are computed, and the knowledge of their distribution under the null hypothesis.

We propose here a review of the methods available to compute pattern statistics on text generated by a Markov source. Theoretical, but also numerical aspects are detailed for a wide range of techniques (exact, Gaussian, large deviations, binomial and compound Poisson). The SPatt package (Statistics for Pattern, free software available at http://stat.genopole.cnrs.fr/spatt) implementing all these methods is then used to compare all these approaches in terms of computational time and reliability in the most complete pattern statistics benchmark available at the present time.

We propose a heuristic approach to the detection of evidence for recombination and gene conversion in multiple DNA sequence alignments. The proposed method consists of two stages. In the first stage, a sliding window is moved along the DNA sequence alignment, and phylogenetic trees are sampled from the conditional posterior distribution with MCMC. To reduce the noise intrinsic to inference from the limited amount of data available in the typically short sliding window, a clustering algorithm based on the Robinson-Foulds distance is applied to the trees thus sampled, and the posterior distribution over tree clusters is obtained for each window position. While changes in this posterior distribution are indicative of recombination or gene conversion events, it is difficult to decide when such a change is statistically significant. This problem is addressed in the second stage of the proposed algorithm, where the distributions obtained in the first stage are post-processed with a Bayesian hidden Markov model (HMM). The emission states of the HMM are associated with posterior distributions over phylogenetic tree topology clusters. The hidden states of the HMM indicate putative recombinant segments. Inference is done in a Bayesian sense, sampling parameters from the posterior distribution with MCMC. Of particular interest is the determination of the number of hidden states as an indication of the number of putative recombinant regions. To this end, we apply reversible jump MCMC, and sample the number of hidden states from the respective posterior distribution.

For situations where the number of tested hypotheses is increasingly large, the power to detect statistically significant multiple treatment effects decreases. As is the case with microarray technology, often researchers are interested in identifying differentially expressed genes for more than two types of cells or treatments. A two-step procedure is proposed for the purpose of increasing power to detect significant effects (i.e., to identify differentially expressed genes). Specifically, in the first step, the null hypothesis of equality across the mean expression levels for all treatments is tested for each gene. In the second step, only pairwise comparisons corresponding to the genes for which the treatment means are statistically different in the first step are tested. We propose an approach to estimate the overall FDR for both fixed rejection regions and fixed FDR significance levels. Also proposed is a procedure to find the FDR significance levels used in the first step and the second step such that the overall FDR can be controlled below a pre-specified FDR significance level. When compared via simulation the two-step approach has increased power over a one-step procedure, and controls the FDR at a desire significance level.

In a recent article in PLoS Genetics, Bock et al., (2006) undertake an extensive computational epigenetics analysis of the ability of DNA sequence-derived features, capturing attributes such as tetramer frequencies, repeats and predicted structure, to predict the methylation status of CpG islands. Their suite of analyses appears highly rigorous with regard to accompanying validation procedures, employing stringent Bonferroni corrections, stratified cross-validation, and follow-up experimental verification. Here, however, we showcase concerns with the validation steps, in part ascribable to the genome scale of the investigation, that serve as a cautionary note and indicate the heightened need for careful selection of analytic and companion validation methods. A series of new analyses of the same CpG island methylation data helps illustrate these issues, not just for this particular study, but also analogous investigations involving high-dimensional predictors with complex between-feature dependencies.

Previous nonparametric statistical methods on constructing the test and null statistics require having at least 4 arrays under each condition. In this paper, we provide an improved method of constructing the test and null statistics which only requires 2 arrays under one condition if the number of arrays under the other condition is at least 3. The conventional testing method defines the rejection region by controlling the probability of Type I error. In this paper, we propose to determine the critical values (or the cut-off points) of the rejection region by directly controlling the false discovery rate. Simulations were carried out to compare the performance of our proposed method with several existing methods. Finally, our proposed method is applied to the rat data of Pan et al. (2003). It is seen from both simulations and the rat data that our method has lower false discovery rates than those from the significance analysis of microarray (SAM) method of Tusher et al. (2001) and the mixture model method (MMM) of Pan et al. (2003).

One of the prevailing ideas in the literature on microarray data analysis is to pool the expression measures across genes and treat them as a sample drawn from some distribution. Several universal laws were proposed to analytically describe this distribution. This idea raises a number of concerns. The expression levels of genes are not identically distributed random variables so that treating them as a sample amounts to sampling from a mixture of equally weighted distributions, each being associated with a different gene. The expression levels of different genes are heavily dependent random variables so that the law of large numbers and statistical goodness-of-fit tests are normally inapplicable to this kind of data. This dependence represents a very serious pitfall in microarray data analysis.

This note is a comment on the article "Dimension Reduction for Classification with Gene Expression Microarray Data" that appeared in Statistical Applications in Genetics and Molecular Biology (Dai et al., 2006).