Volume 7, Issue 1

Based on previous studies related to the yeast cell cycle, it is well known that the underlying cellular network in yeast consists of many interactions between genes that have periodic expression patterns during the cell division cycle. In this study, it is proposed that cell cycle-specific gene expression can be understood as a phenomenon of collective synchronization or, in other words, an ensemble of non-identical oscillating response signals from different systems. Therefore, we aimed to apply the theory of statistical multivariate phase synchronization to understand the cell's cyclic transcriptome as a phenomenon of collective synchronization. To this end, a novel algorithm called Self-Organizing Maps with statistical Phase Synchronization (SOMPS) is proposed and evaluated using yeast cell cycle-specific gene expression data. From the evaluation experiments, we draw the following conclusions: 1) It is possible to find groups of genes that have biological interactions with each other and significantly share gene ontology slim terms of biological processes using the theory of multivariate phase synchronization with cell cycle-specific gene expression signals; 2) Among all output clusters of SOMPS, a relatively large cluster with high periodicity with respect to its trained mean field can be considered a prominent cluster; 3) For each gene, it is possible to identify the degree of the strength of its biological interactions with other genes using the coupling strength of synchronization with its trained mean field; and 4) It is feasible to understand cell cycle-specific expression patterns as a phenomenon of collective synchronization.

The relationship between speciation times and the corresponding times of gene divergence is of interest in phylogenetic inference as a means of understanding the past evolutionary dynamics of populations and of estimating the timing of speciation events. It has long been recognized that gene divergence times might substantially pre-date speciation events. Although the distribution of the difference between these has previously been studied for the case of two populations, this distribution has not been explicitly computed for larger species phylogenies. Here we derive a simple method for computing this distribution for trees of arbitrary size. A two-stage procedure is proposed which (i) considers the probability distribution of the time from the speciation event at the root of the species tree to the gene coalescent time conditionally on the number of gene lineages available at the root; and (ii) calculates the probability mass function for the number of gene lineages at the root. This two-stage approach dramatically simplifies numerical analysis, because in the first step the conditional distribution does not depend on an underlying species tree, while in the second step the pattern of gene coalescence prior to the species tree root is irrelevant. In addition, the algorithm provides intuition concerning the properties of the distribution with respect to the various features of the underlying species tree. The methodology is complemented by developing probabilistic formulae and software, written in R. The method and software are tested on five-taxon species trees with varying levels of symmetry. The examples demonstrate that more symmetric species trees tend to have larger mean coalescent times and are more likely to have a unimodal gamma-like distribution with a long right tail, while asymmetric trees tend to have smaller mean coalescent times with an exponential-like distribution. In addition, species trees with longer branches generally have shorter mean coalescent times, with branches closest to the root of the tree being most influential.

Multiple changes at the DNA level are at the basis of complex diseases. Identifying the genetic networks that are influenced by these changes might help in understanding the development of these diseases. Canonical correlation analysis is used to associate gene expressions with DNA-markers and thus reveals sets of co-expressed and co-regulated genes and their associating DNA-markers. However, when the number of variables gets high, e.g. in the case of microarray studies, interpretation of these results can be difficult. By adapting the elastic net to canonical correlation analysis the number of variables reduces, and interpretation becomes easier, moreover, due to the grouping effect of the elastic net co-regulated and co-expressed genes cluster. Additionally, our adaptation works well in situations where the number of variables exceeds by far the number of subjects.

The development of an organism represents a complex dynamic process, which is controlled by a network of genes and multiple environmental factors. Programmed cell death (PCD), a physiological cell suicide process, occurs during the development of most organisms and is, typically, a complex dynamic trait. Understanding how genes control this complex developmental process has been a long-standing topic in PCD studies. In this article, we propose a nonparametric model, based on orthogonal Legendre polynomials, to map genes or quantitative trait loci (QTLs) that govern the dynamic features of the PCD process. The model is built under the maximum likelihood-based functional mapping framework and is implemented with the EM algorithm. A general information criterion is proposed for selecting the optimal Legendre order that best fits the dynamic pattern of the PCD process. The consistency of the order selection criterion is established. A nonstationary structured antedependence model (SAD) is applied to model the covariance structure among the phenotypes measured at different time points. The developed model generates a number of hypothesis tests regarding the genetic control mechanism of the PCD process. Extensive simulation studies are conducted to investigate the statistical behavior of the model. Finally, we apply the model to a rice tiller number data set in which several QTLs are identified. The developed model provides a quantitative and testable framework for assessing the interplay between genes and the developmental PCD process, and will have great implications for elucidating the genetic architecture of the PCD process.

Multiple branching trees have been used to model the acquisition of HIV drug resistance mutations, and several different algorithms have been developed to construct the tree set that best describes the data. These algorithms have mainly focused on the structure of the tree set. The focal point of this paper is estimation of functions of the tree set parameters that incorporate uncertainty in the tree set. The functions of interest are the state probabilities, the co-occurrence of mutations and the order of acquisition. Such functions are of interest because they help characterize the genetic pathways that lead to multi-drug resistance. We propose a bootstrap technique to account for the additional variability in estimates due to uncertainty in the tree set. The methods are applied to genetic sequences of patients from a database compiled by the Forum for Collaborative HIV Research in an effort to characterize genetic pathways to resistance to drugs from the nucleoside reverse transcriptase inhibitor (NRTI) class. The main results were that patients with a 211K mutation in the RT region of the viral genome were more likely to have a 215Y mutation and less likely to have a 70R mutation compared to patients without a 211K mutation.

In this article, we introduce the drifting Markov models (DMMs) which are inhomogeneous Markov models designed for modeling the heterogeneities of sequences (in our case DNA or protein sequences) in a more flexible way than homogeneous Markov chains or even hidden Markov models (HMMs). We focus here on the polynomial drift: the transition matrix varies in a polynomial way. To show the reliability of our models on DNA, we exhibit high similarities between the probability distributions of nucleotides obtained by our models and the frequencies of these nucleotides computed by using a sliding window. In a further step, these DMMs can be used as the states of an HMM: on each of its segments, the observed process can be modeled by a drifting Markov model. Search of rare words in DNA sequences remains possible with DMMs and according to the fits provided, DMMs turn out to be a powerful tool for this purpose. The software is available on request from the author. It will soon be integrated on seq++ library (http://stat.genopole.cnrs.fr/seqpp/).

One application of gene expression arrays is to derive molecular profiles, i.e., sets of genes, which discriminate well between two classes of samples, for example between tumour types. Users are confronted with a multitude of classification methods of varying complexity that can be applied to this task. To help decide which method to use in a given situation, we compare important characteristics of a range of classification methods, including simple univariate filtering, penalised likelihood methods and the random forest.Classification accuracy is an important characteristic, but the biological interpretability of molecular profiles is also important. This implies both parsimony and stability, in the sense that profiles should not vary much when there are slight changes in the training data. We perform a random resampling study to compare these characteristics between the methods and across a range of profile sizes. We measure stability by adopting the Jaccard index to assess the similarity of resampled molecular profiles.We carry out a case study on five well-established cancer microarray data sets, for two of which we have the benefit of being able to validate the results in an independent data set. The study shows that those methods which produce parsimonious profiles generally result in better prediction accuracy than methods which don't include variable selection. For very small profile sizes, the sparse penalised likelihood methods tend to result in more stable profiles than univariate filtering while maintaining similar predictive performance.

Cross-validation based point estimates of prediction accuracy are frequently reported in microarray class prediction problems. However these point estimates can be highly variable, particularly for small sample numbers, and it would be useful to provide confidence intervals of prediction accuracy.We performed an extensive study of existing confidence interval methods and compared their performance in terms of empirical coverage and width. We developed a bootstrap case cross-validation (BCCV) resampling scheme and defined several confidence interval methods using BCCV with and without bias-correction.The widely used approach of basing confidence intervals on an independent binomial assumption of the leave-one-out cross-validation errors results in serious under-coverage of the true prediction error. Two split-sample based methods previously proposed in the literature tend to give overly conservative confidence intervals. Using BCCV resampling, the percentile confidence interval method was also found to be overly conservative without bias-correction, while the bias corrected accelerated (BCa) interval method of Efron returns substantially anti-conservative confidence intervals. We propose a simple bias reduction on the BCCV percentile interval. The method provides mildly conservative inference under all circumstances studied and outperforms the other methods in microarray applications with small to moderate sample sizes.

Nested Effects Models (NEMs) are a class of graphical models introduced to analyze the results of gene perturbation screens. NEMs explore noisy subset relations between the high-dimensional outputs of phenotyping studies, e.g., the effects showing in gene expression profiles or as morphological features of the perturbed cell.In this paper we expand the statistical basis of NEMs in four directions. First, we derive a new formula for the likelihood function of a NEM, which generalizes previous results for binary data. Second, we prove model identifiability under mild assumptions. Third, we show that the new formulation of the likelihood allows efficiency in traversing model space. Fourth, we incorporate prior knowledge and an automated variable selection criterion to decrease the influence of noise in the data.

The level of differential gene expression may be defined as a fold change, a frequency of upregulation, or some other measure of the degree or extent of a difference in expression across groups of interest. On the basis of expression data for hundreds or thousands of genes, inferring which genes are differentially expressed or ranking genes in order of priority introduces a bias in estimates of their differential expression levels. A previous correction of this feature selection bias suffers from a lack of generality in the method of ranking genes, from requiring many biological replicates, and from unnecessarily overcompensating for the bias.For any method of ranking genes on the basis of gene expression measured for as few as three biological replicates, a simple leave-one-out algorithm corrects, with less overcompensation, the bias in estimates of the level of differential gene expression. In a microarray data set, the bias correction reduces estimates of the probability of upregulation or downregulation from 100% to as low as 60%, even for genes with estimated local false discovery rates close to 0. A simulation study quantifies both the advantage of smoothing estimates of bias before correction and the degree of overcompensation.

The bootstrap is a tool that allows for efficient evaluation of prediction performance of statistical techniques without having to set aside data for validation. This is especially important for high-dimensional data, e.g., arising from microarrays, because there the number of observations is often limited. For avoiding overoptimism the statistical technique to be evaluated has to be applied to every bootstrap sample in the same manner it would be used on new data. This includes a selection of complexity, e.g., the number of boosting steps for gradient boosting algorithms. Using the latter, we demonstrate in a simulation study that complexity selection in conventional bootstrap samples, drawn with replacement, is severely biased in many scenarios. This translates into a considerable bias of prediction error estimates, often underestimating the amount of information that can be extracted from high-dimensional data. Potential remedies for this complexity selection bias, such as alternatively using a fixed level of complexity or of using sampling without replacement are investigated and it is shown that the latter works well in many settings. We focus on high-dimensional binary response data, with bootstrap .632+ estimates of the Brier score for performance evaluation, and censored time-to-event data with .632+ prediction error curve estimates. The latter, with the modified bootstrap procedure, is then applied to an example with microarray data from patients with diffuse large B-cell lymphoma.

It is a common practice to use resampling methods such as the bootstrap for calculating the p-value for each test when performing large scale multiple testing. The precision of the bootstrap p-values and that of the false discovery rate (FDR) relies on the number of bootstraps used for testing each hypothesis. Clearly, the larger the number of bootstraps the better the precision. However, the required number of bootstraps can be computationally burdensome, and it multiplies the number of tests to be performed. Further adding to the computational challenge is that in some applications the calculation of the test statistic itself may require considerable computation time. As technology improves one can expect the dimension of the problem to increase as well. For instance, during the early days of microarray technology, the number of probes on a cDNA chip was less than 10,000. Now the Affymetrix chips come with over 50,000 probes per chip. Motivated by this important need, we developed a simple adaptive bootstrap methodology for large scale multiple testing, which reduces the total number of bootstrap calculations while ensuring the control of the FDR. The proposed algorithm results in a substantial reduction in the number of bootstrap samples. Based on a simulation study we found that, relative to the number of bootstraps required for the Benjamini-Hochberg (BH) procedure, the standard FDR methodology which was the proposed methodology achieved a very substantial reduction in the number of bootstraps. In some cases the new algorithm required as little as 1/6th the number of bootstraps as the conventional BH procedure. Thus, if the conventional BH procedure used 1,000 bootstraps, then the proposed method required only 160 bootstraps. This methodology has been implemented for time-course/dose-response data in our software, ORIOGEN, which is available from the authors upon request.

Nucleosomes, the fundamental repeating subunits of all eukaryotic chromatin, are responsible for packaging DNA into chromosomes inside the cell nucleus and controlling gene expression. While it has been well established that nucleosomes exhibit higher affinity for select DNA sequences, until recently it was unclear whether such preferences exerted a significant, genome-wide effect on nucleosome positioning in vivo. This question was seemingly and recently resolved in the affirmative: a wide-ranging series of experimental and computational analyses provided extensive evidence that the instructions for wrapping DNA around nucleosomes are contained in the DNA itself. This subsequently labeled second genetic code was based on data-driven, structural, and biophysical considerations. It was subjected to an extensive suite of validation procedures, with one conclusion being that intrinsic, genome-encoded, nucleosome organization explains approximately 50% of in vivo nucleosome positioning. Here, we revisit both the nature of the underlying sequence preferences, and the performance of the proposed code. A series of new analyses, employing spectral envelope (Fourier transform) methods for assessing key sequence periodicities, classification techniques for evaluating predictive performance, and discriminatory motif finding methods for devising alternate models, are applied. The findings from the respective analyses indicate that signature dinucleotide periodicities are absent from the bulk of the high affinity nucleosome-bound sequences, and that the predictive performance of the code is modest. We conclude that further exploration of the role of sequence-based preferences in genome-wide nucleosome positioning is warranted. This work offers a methodologic counterpart to a recent, high resolution determination of nucleosome positioning that also questions the accuracy of the proposed code and, further, provides illustrations of techniques useful in assessing sequence periodicity and predictive performance.

We develop an approach for microarray differential expression analysis, i.e. identifying genes whose expression levels differ between two or more groups. Current approaches to inference rely either on full parametric assumptions or on permutation-based techniques for sampling under the null distribution. In some situations, however, a full parametric model cannot be justified, or the sample size per group is too small for permutation methods to be valid.We propose a semi-parametric framework based on partial mixture estimation which only requires a parametric assumption for the null (equally expressed) distribution and can handle small sample sizes where permutation methods break down. We develop two novel improvements of Scott's minimum integrated square error criterion for partial mixture estimation [Scott, 2004a,b]. As a side benefit, we obtain interpretable and closed-form estimates for the proportion of EE genes. Pseudo-Bayesian and frequentist procedures for controlling the false discovery rate are given. Results from simulations and real datasets indicate that our approach can provide substantial advantages for small sample sizes over the SAM method of Tusher et al. [2001], the empirical Bayes procedure of Efron and Tibshirani [2002], the mixture of normals of Pan et al. [2003] and a t-test with p-value adjustment [Dudoit et al., 2003] to control the FDR [Benjamini and Hochberg, 1995].

The availability of very dense genetic maps is changing in a fundamental way the methods used to identify the genetic basis of both rare and common inherited traits. The ability to directly compare the genomes of two related individuals and quickly identify those regions that are inherited identical-by-descent (IBD) from a recent common ancestor would be of utility in a wide range of genetic mapping methods. Here, we describe a simple method for using dense SNP maps to identify regions of the genome likely to be inherited IBD by family members. This method is based on identifying obligate recombination events and examining the pattern of distribution of such events along the genetic map. Specifically, we use the length of a consecutive set of biallelic markers that have a high probability of having avoided such obligate recombination events. This ``SNP streak" is derived from subsets of samples within a pedigree and allows us to make statistical inferences about the ancestry of the region(s) containing stretches of markers with these properties. We show that the use of subsets of more than two samples has the advantage of identifying shorter shared subsegments as significant. This mitigates the effects of errors in SNP calls. We provide specific examples of microarray-based SNP data, using a family with a complex pedigree and with a rare form of inherited kidney disease, to illustrate this approach.

As more genetic loci are genotyped simultaneously and as the interest in effects of combinations of loci increases, the need for more powerful analysis methods is increased. In the present paper we present a method aimed at increasing the power of likelihood ratio tests for case-control studies investigating possible two-locus effects. The method is based on the notion that the expected effect pattern of one locus, as well as the expected pattern of a penetrance matrix representing the effect of two loci, is a monotone one. By using an algorithm for making the estimated penetrance matrix monotone, the alternative hypothesis is restricted to monotone penetrance matrices only. The evaluation of the likelihood ratio tests for several underlying monotone models shows that the power is substantially increased by using a monotone alternative as compared to when an unrestricted alternative is used.

Little is known about how human cancers grow because direct observations are impractical. Cancers are clonal populations and the billions of cancer cells present in a visible tumor are progeny of a single transformed cell. Therefore, human cancers can be represented by somatic cell ancestral trees that start from a single transformed cell and end with billions of present day cancer cells. We use a genealogical approach to infer tumor growth from somatic trees, employing haplotype DNA methylation pattern variation, or differences between specific CpG sites or "tags," in the cancer genome. DNA methylation is an epigenetic mark that is copied, with error, during genome replication. At our tags, neutral copy errors in DNA methylation appear to occur at random, and much more frequently than sequence copy errors. To reconstruct a cancer tree, we sample and compare human colorectal genomes within small geographic regions (a cancer fragment), between fragments on the same side of the tumor, and between fragments from opposite tumor halves. The combined information on both physical distance and epigenetic distance informs our model for tumor ancestry. We use approximate Bayesian computation, a simulation-based method, to model tumor growth under a variety of evolutionary scenarios, estimating parameters that fit observed DNA methylation patterns. We conclude that methylation patterns sampled from human cancers are consistent with replication errors and certain simple cancer growth models. The inferred cancer trees are consistent with Gompertzian growth, a well-known cancer growth pattern.

As the field of genomics matures, more complex genotypes and phenotypes are being studied. Fanconi anemia (FA), for example, is an inherited chromosome instability syndrome with a complex array of variable disease phenotypes including congenital malformations, hematological manifestations, and cancer. To better understand specific aspects of the genetic etiology of FA and other rare diseases with complex phenotypes, it is often necessary to reduce the dimensions of the disease phenotype information. Towards this end, we extend a novel non-parametric approach to include information about a hierarchical structure among disease phenotypes. The proposed extension increases information content of the phenotype scores obtained and, thereby, the power of genotype-phenotype relationships studies.

The allelic association or linkage disequilibrium between two loci is a parameter of fundamental interest in modern population genetics for evolutionary inference and association mapping studies. Among the many measures available, the optimal measure of allelic association ? presents a strong evolutionary theory basis and is modeled on the physical distance along the chromosome with the Malécot equation for isolation by distance. Moreover, ? is equal to the absolute value of D', the standardized measure of gametic disequilibrium. We studied here the statistical properties of the ? sample estimator. We derived its asymptotic probability distribution and showed that it is neither asymptotically normal nor unbiased when ?=0 or when allelic frequencies are equal at both loci, in contrast to previous claims. This asymptotic study leads to propose a new test for absence of linkage disequilibrium. We compared it to Pearson's ?2 test for independence in a contingency table and showed by simulations that the range in power of these two tests depends on the sign of D'. The new test outperformed slightly the ?2 test, when D', polarized with respect to major alleles, is negative. Finally, we derived the asymptotic bias and information of the ? estimator that are due to the experimental sampling and showed by simulation that its bias is large in small samples. The consequences of these findings on applications using the ? measure are then discussed in particular for constructing LD unit maps, and call for a revised statistical treatment.

Making sound proteomic inferences using ELISA microarray assay requires both an accurate prediction of protein concentration and a credible estimate of its error. We present a method using monotonic spline statistical models (MS), penalized constrained least squares fitting (PCLS) and Monte Carlo simulation (MC) to predict ELISA microarray protein concentrations and estimate their prediction errors. We contrast the MSMC (monotone spline Monte Carlo) method with a LNLS (logistic nonlinear least squares) method using simulated and real ELISA microarray data sets.MSMC rendered good fits in almost all tests, including those with left and/or right clipped standard curves. MS predictions were nominally more accurate; especially at the extremes of the prediction curve. MC provided credible asymmetric prediction intervals for both MS and LN fits that were superior to LNLS propagation-of-error intervals in achieving the target statistical confidence. MSMC was more reliable when automated prediction across simultaneous assays was applied routinely with minimal user guidance.

Normalization of expression levels applied to microarray data can help in reducing measurement error. Different methods, including cyclic loess, quantile normalization and median or mean normalization, have been utilized to normalize microarray data. Although there is considerable literature regarding normalization techniques for mRNA microarray data, there are no publications comparing normalization techniques for microRNA (miRNA) microarray data, which are subject to similar sources of measurement error. In this paper, we compare the performance of cyclic loess, quantile normalization, median normalization and no normalization for a single-color microRNA microarray dataset. We show that the quantile normalization method works best in reducing differences in miRNA expression values for replicate tissue samples. By showing that the total mean squared error are lowest across almost all 36 investigated tissue samples, we are assured that the bias correction provided by quantile normalization is not outweighed by additional error variance that can arise from a more complex normalization method. Furthermore, we show that quantile normalization does not achieve these results by compression of scale.

Genome-wide association studies are now widely used tools to identify genes and/or regions which may contribute to the development of various diseases. With case-control data a 2x3 contingency table can be constructed for each SNP to perform genotype-based tests of association. An increasingly common technique to increase the power to detect an association is to collapse each 2x3 table into a table assuming either a dominant or recessive mode of inheritance (2x2 table). We consider three different methods of determining which genetic model to choose and show that each of these methods of collapsing genotypes increases the type I error rate (i.e., the rate of false positives). However, one of these methods does lead to an increase in power compared with the usual genotype- and allele-based tests for most genetic models.

Many clustering methods require that the number of clusters believed present in a given data set be specified a priori, and a number of methods for estimating the number of clusters have been developed. However, the selection of the number of clusters is well recognized as a difficult and open problem and there is a need for methods which can shed light on specific aspects of the data. This paper adopts a model for clustering based on a specific structure for a similarity matrix. Publicly available gene expression data sets are analyzed to illustrate the method and the performance of our method is assessed by simulation.

Information obtained by microarray studies frequently cause discrepancies between projects, laboratories and/or even repeated experiments. Such inconsistencies may be due to the lack of generality in the intellectual frameworks that form the basis for understanding the data. This article proposes a parametric framework that can handle a wide range of experimental data obtained by GeneChip expression arrays. The framework is based on a parsimonious model, which has been developed according to thermodynamic estimations of the process of hybridization. Using the model, probe data were normalized and summarized into gene expression levels. Verification of the appropriateness of the model is demonstrated statistically by the use of real data obtained from several project series. Furthermore, improved stabilities in changes in expression are observed in comparison with other currently used methods. Estimations of transcriptome differences between organs coincided between various projects.

The analysis of high-dimensional data sets is often forced to rely upon well-chosen summary statistics. A systematic approach to choosing such statistics, which is based upon a sound theoretical framework, is currently lacking. In this paper we develop a sequential scheme for scoring statistics according to whether their inclusion in the analysis will substantially improve the quality of inference. Our method can be applied to high-dimensional data sets for which exact likelihood equations are not possible. We illustrate the potential of our approach with a series of examples drawn from genetics. In summary, in a context in which well-chosen summary statistics are of high importance, we attempt to put the `well' into `chosen.'

A composite-conditional-likelihood (CCL) approach is proposed to map the position of a trait-influencing mutation (TIM) using the ancestral recombination graph (ARG) and importance sampling to reconstruct the genealogy of DNA sequences with respect to windows of marker loci and predict the linkage disequilibrium pattern observed in a sample of cases and controls. The method is designed to fine-map the location of a disease mutation, not as an association study. The CCL function proposed for the position of the TIM is a weighted product of conditional likelihood functions for windows of a given number of marker loci that encompass the TIM locus, given the sample configuration at the marker loci in those windows. A rare recessive allele is assumed for the TIM and single nucleotide polymorphisms (SNPs) are considered as markers. The method is applied to a range of simulated data sets. Not only do the CCL profiles converge more rapidly with smaller window sizes as the number of simulated histories of the sampled sequences increases, but the maximum-likelihood estimates for the position of the TIM remain as satisfactory, while requiring significantly less computing time. The simulations also suggest that non-random samples, more precisely, a non-proportional number of controls versus the number of cases, has little effect on the estimation procedure as well as sample size and marker density beyond some threshold values. Moreover, when compared with some other recent methods under the same assumptions, the CCL approach proves to be competitive.

We propose a suite of statistical methods for inferring a cis-regulatory module, which is a combination of several transcription factors binding in the promoter regions to regulate gene expression. The approach is an integrative analysis that combines information from multiple types of biological data, including genomic DNA sequences, genome-wide location analysis (ChIP-chip experiments), and gene expression microarray. More specifically, we use a hidden Markov model to first predict a cluster of transcription factor binding sites in DNA sequences. The predictions are refined by regression analysis on gene expression microarray data and/or ChIP-chip binding experiments. In regression analysis, we particularly apply factor analysis, whose statistical model characterizes the modular structure of cis-regulation. When groups of coexpressed genes are available, we further apply canonical correlation analysis to infer relationships between a group of genes and their common set of transcription factors. Our approach is validated on the well-studied yeast cell cycle gene regulation. It is then used to study condition-specific regulators for a set of Ste12 target genes. The multiple data sources provide information of transcriptional regulation from different aspects. Therefore, the integrative analysis offers a fine prediction on transcriptional regulatory code and infers potential regulatory networks.

We address the problem of detection and correction of spatial flaws in oligonucleotide microarrays. We present two similar procedures, of which one is intended solely for use with replicates and the other has wider applicability. By constructing a set of replicates, with one realistically flawed, we are able to examine the extent to which our procedures are capable of repairing the flaw. We find that, for this purpose, our procedures are superior to the existing `Harshlight' procedure.

In this paper we propose the minimum entropy clustering (MEC) method for clustering genes based on their phylogenetic signals. This entropy based method will cluster two genes together when their concatenation can decrease the entropy. An integral feature of MEC is that it chooses the number of clusters automatically, which is a major advantage over the other methods. Our simulation results show that this method is quite successful in clustering genes with a common phylogeny.

Methods for multivariate meta-analysis of genetic association studies are reviewed, summarized and presented in a unified framework. Modifications of standard models are described in detail in order to be applied in genetic association studies. The model based on summary data is uniformly defined for both discrete and continuous outcomes and analytical expressions for the covariance of the two jointly modeled outcomes are derived for both cases. The models based on the binary nature of the data are fitted using both prospective and retrospective likelihood. Furthermore, formal tests for assessing the genetic model of inheritance are developed based on standard normal theory. The general model is compared to the recently proposed genetic model-free bivariate approach (either using summary or binary data), and it is clearly shown that the estimates provided by this approach are nearly identical to the estimates derived by the general bivariate model using the aforementioned tests for the genetic model. The methods developed here as well as the tests, are easily implemented in all major statistical packages, escaping the need of self written software. The methods are applied in several already published meta-analyses of genetic association studies (with both discrete and continuous outcomes) and the results are compared against the widely used univariate approach as well as against the genetic model free approaches. Illustrative examples of code in Stata are given in the appendix. It is anticipated that the methods developed in this work will be widely applied in the meta-analysis of genetic association studies.

Importance sampling or Markov Chain Monte Carlo sampling is required for state-of-the-art statistical analysis of population genetics data. The applicability of these sampling-based inference techniques depends crucially on the proposal distribution. In this paper, we discuss importance sampling for the infinite sites model. The infinite sites assumption is attractive because it constraints the number of possible genealogies, thereby allowing for the analysis of larger data sets. We recall the Griffiths-Tavaré and Stephens-Donnelly proposals and emphasize the relation between the latter proposal and exact sampling from the infinite alleles model. We also introduce a new proposal that takes knowledge of the ancestral state into account. The new proposal is derived from a new result on exact sampling from a single site. The methods are illustrated on simulated data sets and the data considered in Griffiths and Tavaré (1994).

We introduce a novel statistical concept, called a supervised distance matrix, which quantifies pairwise similarity between variables in terms of their association with an outcome. Supervised distance matrices are derived in two stages. First, the observed data is transformed based on particular working models for association. Examples of transformations include residuals or influence curves from regression models. In the second stage, a choice of distance measure is used to compute all pairwise distances between variables in the transformed data. We present consistent estimators of the resulting distance matrix, including an inverse probability of censoring weighted estimator for use with right-censored outcomes. Supervised distance matrices can be used with standard (unsupervised) clustering algorithms to identify groups of similarly predictive variables and to discover subpopulations of related samples. This approach is illustrated using simulations and an analysis of gene expression data with a censored survival outcome. The proposed methods are widely applicable in genomics and other fields where high-dimensional data is collected on each subject.

We address a potential shortcoming of three probabilistic models for detecting interspecific recombination in DNA sequence alignments: the multiple change-point model (MCP) of Suchard et al. (2003), the dual multiple change-point model (DMCP) of Minin et al. (2005), and the phylogenetic factorial hidden Markov model (PFHMM) of Husmeier (2005). These models are based on the Bayesian paradigm, which requires the solution of an integral over the space of branch lengths. To render this integration analytically tractable, all three models make the same assumption that the vectors of branch lengths of the phylogenetic tree are independent among sites. While this approximation reduces the computational complexity considerably, we show that it leads to the systematic prediction of spurious topology changes in the Felsenstein zone, that is, the area in the branch lengths configuration space where maximum parsimony consistently infers the wrong topology due to long-branch attraction. We apply two Bayesian hypothesis tests, based on an inter- and an intra-model approach to estimating the marginal likelihood. We then propose a revised model that addresses these shortcomings, and compare it with the aforementioned models on a set of synthetic DNA sequence alignments systematically generated around the Felsenstein zone.

Recent biotechnology advances allow for multiple types of omics data, such as transcriptomic, proteomic or metabolomic data sets to be integrated. The problem of feature selection has been addressed several times in the context of classification, but needs to be handled in a specific manner when integrating data. In this study, we focus on the integration of two-block data that are measured on the same samples. Our goal is to combine integration and simultaneous variable selection of the two data sets in a one-step procedure using a Partial Least Squares regression (PLS) variant to facilitate the biologists' interpretation. A novel computational methodology called ``sparse PLS" is introduced for a predictive analysis to deal with these newly arisen problems. The sparsity of our approach is achieved with a Lasso penalization of the PLS loading vectors when computing the Singular Value Decomposition.Sparse PLS is shown to be effective and biologically meaningful. Comparisons with classical PLS are performed on a simulated data set and on real data sets. On one data set, a thorough biological interpretation of the obtained results is provided. We show that sparse PLS provides a valuable variable selection tool for highly dimensional data sets.

To rigorously determine whether a gene or a set of genes have alterations that are involved in carcinogenesis requires a comparison of the prevalence of identified changes to a control mutation frequency present in tumor DNA. To facilitate this task, we develop a testing approach and the associated R library, called TRAB, that evaluates whether the frequency of somatic mutation in a given gene is higher than that observed in a control group of genes. Specifically, we test the null hypothesis that the frequency belongs to a control population of frequencies, against the alternative hypothesis that the frequency is higher. Mutation frequencies in the control group are themselves allowed to be variable. TRAB computes the a posteriori probability and the Bayes factor for the hypothesis using a hierarchical Bayesian approach.