Volume 9, Issue 1

The term "epistasis" is sometimes used to describe some form of statistical interaction between genetic factors and is alternatively sometimes used to describe instances in which the effect of a particular genetic variant is masked by a variant at another locus. In general statistical tests for interaction are of limited use in detecting "epistasis" in the sense of masking. It is, however, shown that there are relations between empirical data patterns and epistasis that have not been previously noted. These relations can sometimes be exploited to empirically test for "epistatic interactions" in the sense of the masking of the effect of a particular genetic variant by a variant at another locus.

Studies have shown that many common diseases are influenced by multiple genes and their interactions. There is currently a strong interest in testing for association between combinations of these genes and disease, in particular because genes that affect the risk of disease only in the presence of another genetic variant may not be detected in marginal analysis. In this paper we propose the use of additive main effect and multiplicative interaction (AMMI) models to detect and to quantify gene-gene interaction effects for a quantitative trait. The objective of the present research is to demonstrate the practical advantages of these models to describe complex interaction between two unlinked loci. Although gene-gene interactions have often been defined as a deviance from additive genetic effects, the residual term has generally not been appropriately treated. The AMMI models allow for the analysis of a two way factorial data structure and combine the analysis of variance of the two main genotype effects with a principal component analysis of the residual multiplicative interaction. The AMMI models for gene-gene interaction presented here allow for the testing of non additivity between the two loci, and also describe how their interaction structure fits the existing non-additivity. Moreover, these models can be used to identify the specific two genotypes combinations that contribute to the significant gene-gene interaction. We describe the use of the biplot to display the structure of the interaction and evaluate the performance of the AMMI and the special cases of the AMMI previously described by Tukey and Mandel with simulated data sets. Our simulated study showed that the AMMI model is as powerful as general linear models when the interaction is not modeled in the presence of marginal effects. However, in the presence of pure epitasis, i.e. in the absence of marginal effects, the AMMI method was not found to be superior to other tested regression methods.

We present a Bayesian hierarchical model for quantitative real-time polymerase chain reaction (PCR) data, aiming at relative quantification of DNA copy number in different biological samples. The model is specified in terms of a hidden Markov model for fluorescence intensities measured at successive cycles of the polymerase chain reaction. The efficiency of the reaction is assumed to depend on the abundance of the target DNA through fluorescence intensities, and the relationship is specified based on the kinetics of the reaction. The model incorporates the intrinsic random nature of the process as well as measurement error. Taking a Bayesian inferential approach, marginal posterior distributions of the quantities of interest are estimated using Markov chain Monte Carlo. The method is applied to simulated data and an experimental data set.

The strength and weakness of microarray technology can be attributed to the enormous amount of information it is generating. To fully enhance the benefit of microarray technology for testing differentially expressed genes and classification, there is a need to minimize the amount of irrelevant genes present in microarray data. A major interest is to use probe-level data to call genes informative or noninformative based on the trade-off between the array-to-array variability and the measurement error. Existing works in this direction include filtering likely uninformative sets of hybridization (FLUSH; Calza et al., 2007) and I/NI calls for the exclusion of noninformative genes using FARMS (I/NI calls; Talloen et al., 2007; Hochreiter et al., 2006). In this paper, we propose a linear mixed model as a more flexible method that performs equally good as I/NI calls and outperforms FLUSH. We also introduce other criteria for gene filtering, such as, R2 and intra-cluster correlation. Additionally, we include some objective criteria based on likelihood ratio testing, the Akaike information criteria (AIC; Akaike, 1973) and the Bayesian information criterion (BIC; Schwarz, 1978 ).Based on the HGU-133A Spiked-in data set, it is shown that the linear mixed model approach outperforms FLUSH, a method that filters genes based on a quantile regression. The linear model is equivalent to a factor analysis model when either the factor loadings are set to a constant with the variance of the latent factor equal to one, or if the factor loadings are set to one together with unconstrained variance of the latent factor. Filtering based on conditional variance calls a probe set informative when the intensity of one or more probes is consistent across the arrays, while filtering using R2 or intra-cluster correlation calls a probe set informative only when average intensity of a probe set is consistent across the arrays. Filtering based on likelihood ratio test AIC and BIC are less stringent compared to the other criteria.

Tests for Hardy-Weinberg equilibrium (HWE) have been used to detect genotyping error, but those tests have low power unless the sample size is very large. We assessed the performance of measures of departure from HWE as an alternative way of screening for genotyping error. Three measures of the degree of disequilibrium (?, ,D, and F) were tested for their ability to detect genotyping error of 5% or more using simulations and a real dataset of 184 children with leukemia genotyped at 28 single nucleotide polymorphisms. The simulations indicate that all three disequilibrium coefficients can usefully detect genotyping error as judged by the area under the Receiver Operator Characteristic (ROC) curve. Their discriminative ability increases as the error rate increases, and is greater if the genotyping error is in the direction of the minor allele. Optimal thresholds for detecting genotyping error vary for different allele frequencies and patterns of genotyping error but allele frequency-specific thresholds can be nominated. Applying these thresholds would have picked up about 90% of genotyping errors in our actual dataset. Measures of departure from HWE may be useful for detecting genotyping error, but this needs to be confirmed in other real datasets.

Hidden Markov models (HMMs) play a major role in applications to unravel biomolecular functionality. Though HMMs are technically mature and widely applied in computational biology, there is a potential of methodical optimisation concerning its modelling of biological data sources with varying sequence lengths.Single building blocks of these models, the states, are associated with a certain holding time, being the link to the length distribution of represented sequence motifs. An adaptation of regular HMM topologies to bell-shaped sequence lengths is achieved by a serial chain-linking of hidden states, while residing in the class of conventional hidden Markov models. The factor of the repetition of states (r) and the parameter for state-specific duration of stay (p) are determined by fitting the distribution of sequence lengths with the method of moments (MM) and maximum likelihood (ML). Performance evaluations of differently adjusted HMM topologies underline the impact of an optimisation for HMMs based on sequence lengths. Secondary structure prediction on internal transcribed spacer 2 sequences demonstrates exemplarily the general impact of topological optimisations. In summary, we propose a general methodology to improve the modelling behaviour of HMMs by topological optimisation with ML and a fast and easily implementable moment estimator.

We wish to suggest the categorical analysis of variance as a means of quantifying the proportion of total genetic variation attributed to different sources of variation. This method potentially challenges researchers to rethink conclusions derived from a well-known method known as the analysis of molecular variance (AMOVA). The CATANOVA framework allows explicit definition, and estimation, of two measures of genetic differentiation. These parameters form the subject of interest in many research programmes, but are often confused with the correlation measures defined in AMOVA, which cannot be interpreted as relative contributions of particular sources of variation. Through a simulation approach, we show that under certain conditions, researchers who use AMOVA to estimate these measures of genetic differentiation may attribute more than justified amounts of total variation to population labels. Moreover, the two measures can also lead to incongruent conclusions regarding the genetic structure of the populations of interest. Fortunately, one of the two measures seems robust to variations in relative sample sizes used. Its merits are illustrated in this paper using mitochondrial haplotype and amplified fragment length polymorphism (AFLP) data.

In Genomewide association (GWA) studies investigating thousands of SNPs, large sample sizes are needed to obtain a reasonable power after correction for multiple testing. To obtain the necessary sample sizes, data from different populations/cohorts are combined. The problem of pooling evidence across cohorts bears some resemblance with meta-analysis of clinical trials, and in fact classical meta-analytic methodologies from that field are typically used in GWAs. However, in genetics, it can be expected that the cohorts show some amount of heterogeneity in the association measures that are used for significance testing. In this paper, we demonstrate how it is possible to exploit this heterogeneity to improve our ability to detect influential genetic variants. We also discuss how pathway analysis based on summary data can help resolve heterogeneity. The current standard method for testing SNPs across cohorts in GWAs will miss heterogeneous but important genetic variants affecting complex diseases. Our new testing strategy has the potential to detect them while maintaining sensitivity to variants with homogeneous effects.

Gene regulatory networks refer to the interactions that occur among genes and other cellular products. The topology of these networks can be inferred from measurements of changes in gene expression over time. However, because the measurement device (i.e., microarrays) typically yields information on thousands of genes over few biological replicates, these systems are quite difficult to elucidate. An approach with proven effectiveness for inferring networks is the Dynamic Bayesian Network. We have developed an iterative empirical Bayesian procedure with a Kalman filter that estimates the posterior distributions of network parameters. We compare our method to similar existing methods on simulated data and real microarray time series data. We find that the proposed method performs comparably on both model-based and data-based simulations in considerably less computational time. The R and C code used to implement the proposed method are publicly available in the R package ebdbNet.

In this work we deal with parameter estimation in a latent variable model, namely the multiple-hidden i.i.d. model, which is derived from multiple alignment algorithms. We first provide a rigorous formalism for the homology structure of k sequences related by a star-shaped phylogenetic tree in the context of multiple alignment based on indel evolution models. We discuss possible definitions of likelihoods and compare them to the criterion used in multiple alignment algorithms. Existence of two different Information divergence rates is established and a divergence property is shown under additional assumptions. This would yield consistency for the parameter in parametrization schemes for which the divergence property holds. We finally extend the definition of the multiple-hidden i.i.d. model and the results obtained to the case in which the sequences are related by an arbitrary phylogenetic tree. Simulations illustrate different cases which are not covered by our results.

Population structure is a recurrent problem for the detection of associations between a marker and a trait, because it can lead to an excess of false positives of the association tests. One popular way of circumventing this problem is the use of family based tests, which consider the transmission of the genotype from the parents to the offspring. Here we focus on quantitative traits and study the Abecasis “orthogonal" quantitative transmission disequilibrium test, which is commonly used in family based association studies. We derive the probability distribution of this test under a general model of structured population. Our derivations show that this test leads to a small excess of false positives due to population structure. They also illustrate and quantify how the heterogeneity in genotypes and phenotypes between populations affect the power of the test. We finally show that the excess of false positives observed for the Abecasis “orthogonal" test may also be found for the Allison “linear" test, though at a lower extent.

Recent applications of principal components analysis (PCA) and multidimensional scaling (MDS) in human population genetics have found that "statistical maps" based on the genotypes in population-genetic samples often resemble geographic maps of the underlying sampling locations. To provide formal tests of these qualitative observations, we describe a Procrustes analysis approach for quantitatively assessing the similarity of population-genetic and geographic maps. We confirm in two scenarios, one using single-nucleotide polymorphism (SNP) data from Europe and one using SNP data worldwide, that a measurably high level of concordance exists between statistical maps of population-genetic variation and geographic maps of sampling locations. Two other examples illustrate the versatility of the Procrustes approach in population-genetic applications, verifying the concordance of SNP analyses using PCA and MDS, and showing that statistical maps of worldwide copy-number variants (CNVs) accord with statistical maps of SNP variation, especially when CNV analysis is limited to samples with the highest-quality data. As statistical maps with PCA and MDS have become increasingly common for use in summarizing population relationships, our examples highlight the potential of Procrustes-based quantitative comparisons for interpreting the results in these maps.

Nuisance factors in a protein-array study add obfuscating variation to spot intensity measurements, diminishing the accuracy and precision of protein concentration predictions. The effects of nuisance factors may be reduced by design of experiments, and by estimating and then subtracting nuisance effects. Estimated nuisance effects also inform about the quality of the study and suggest refinements for future studies.We demonstrate a method to reduce nuisance effects by incorporating a non-interfering internal calibration in the study design and its complemental analysis of variance. We illustrate this method by applying a chip-level internal calibration in a biomarker discovery study.The variability of sample intensity estimates was reduced 16% to 92% with a median of 58%; confidence interval widths were reduced 8% to 70% with a median of 35%. Calibration diagnostics revealed processing nuisance trends potentially related to spot print order and chip location on a slide.The accuracy and precision of a protein-array study may be increased by incorporating a non-interfering internal calibration. Internal calibration modeling diagnostics improve confidence in study results and suggest process steps that may need refinement. Though developed for our protein-array studies, this internal calibration method is applicable to other targeted array-based studies.

We present a weighted-LASSO method to infer the parameters of a first-order vector auto-regressive model that describes time course expression data generated by directed gene-to-gene regulation networks. These networks are assumed to own prior internal structures of connectivity which drive the inference method. This prior structure can be either derived from prior biological knowledge or inferred by the method itself. We illustrate the performance of this structure-based penalization both on synthetic data and on two canonical regulatory networks (the yeast cell cycle regulation network and the E. coli S.O.S. DNA repair network).

Multiallelic markers, such as microsatellites, provide a powerful tool for studying the genetic structure and organization of an outcrossing population. However, statistical methods of analyzing multiallelic markers in current literature are limited in scope due to the complexity of the multiple alleles. We present a closed-form EM algorithm framework to estimate trigenic linkage disequilibria coefficients of three multiallelic markers and present joint and separate statistical hypothesis tests of different linkage disequilibria. Linkage disequilibria analysis with three multiallelic markers is shown to be considerably more powerful than a two marker analysis or a three marker analysis that treats the multiallelic markers as biallelic markers. A three multiallelic marker model was used to analyze marker data from Lycoris longituba, a tulip-like ornamental plant in China, where each marker consisted of two to four distinct alleles. This algorithm will be useful for studying the pattern of genetic variation for outcrossing populations.

Partial least squares (PLS) is a well known dimension reduction method which has been recently adapted for high dimensional classification problems in genome biology. We develop sparse versions of the recently proposed two PLS-based classification methods using sparse partial least squares (SPLS). These sparse versions aim to achieve variable selection and dimension reduction simultaneously. We consider both binary and multicategory classification. We provide analytical and simulation-based insights about the variable selection properties of these approaches and benchmark them on well known publicly available datasets that involve tumor classification with high dimensional gene expression data. We show that incorporation of SPLS into a generalized linear model (GLM) framework provides higher sensitivity in variable selection for multicategory classification with unbalanced sample sizes between classes. As the sample size increases, the two-stage approach provides comparable sensitivity with better specificity in variable selection. In binary classification and multicategory classification with balanced sample sizes, the two-stage approach provides comparable variable selection and prediction accuracy as the GLM version and is computationally more efficient.

There are a number of common human diseases for which the genetic component may include an epistatic interaction of multiple genes. Detecting these interactions with standard statistical tools is difficult because there may be an interaction effect, but minimal or no main effect. Reconstructability analysis (RA) uses Shannon's information theory to detect relationships between variables in categorical datasets. We applied RA to simulated data for five different models of gene-gene interaction, and find that even with heritability levels as low as 0.008, and with the inclusion of 50 non-associated genes in the dataset, we can identify the interacting gene pairs with an accuracy of ?80%. We applied RA to a real dataset of type 2 non-insulin-dependent diabetes (NIDDM) cases and controls, and closely approximated the results of more conventional single SNP disease association studies. In addition, we replicated prior evidence for epistatic interactions between SNPs on chromosomes 2 and 15.

Inferring the structure of networks usually involves the attempt of retrieving their modular organization and knowing its possible interpretation, while quantifying the involved computational complexity through the methods and algorithms to be used. In protein interactomics, it is assumed that even the most recently created interactomes are known only up to a certain degree of coverage and accuracy, due to both experimental and computational limitations. Therefore, we need to infer from the measured interactomes about real interactomes as much as we infer from samples relative to a reference population. In order to exploit additional information sources, it is common to integrate multiple omic data and pursue method fusion. Particularly after the advent of high-throughput technologies, the increased complexity of data-intensive applications has determined an important role for network inference. Consequently, advances in spectral clustering, community detection algorithms and modularity optimization methods have been proposed, according to both deterministic and probabilistic solutions. We have considered the two kinds of approaches, and applied some of the available methods to two human interactomes obtained from high-throughput small-scale experiments and mass spectrometry measurements. The main motivation of this study is refining the resolution spectrum at which protein modularity maps can be studied. First, we started by a coarse-grained interactome decomposition through core and community structures, and by applying sub-sampling to the interactome adjacency matrix. Then, we switched to stochastic methods to uncover fine-grained interactome components, and applied both variational and mixture statistical models. Lastly, we integrated our analysis with the biological validation of the retrieved modules. Overall, the proposed approach shows potential for calibrating modularity detection in protein interactomes at different resolutions.

Integration of data from multiple omics platforms has become a major challenge in studying complex systems and traits. For integrating data from multiple platforms, the underlying spaces from which the top ranked elements come from are likely to be different. Thus, taking the underlying spaces into consideration explicitly is important, as failure to do so would lead to inefficient use of data and might render biases and/or sub-optimal results. We propose two space oriented classes of heuristic algorithms for integrating ranked lists from omic scale data. These algorithms are either Borda inspired or Markov chain based that take the underlying spaces of the individual ranked lists into account explicitly. We applied this set of algorithms to a number of problems, including one that aims at aggregating results from three cDNA and two Affymetrix gene expression studies in which the underlying spaces between Affymetrix and cDNA platforms are clearly different.

The significance of risk effects in genetic association studies is assessed using the odds ratio for various genetic models (dominant, recessive and co-dominant) by merging genotypes. These models are not independent and there is no a priori biological justification for their choice. Consequently, the interpretation of their results can be problematic, especially when multiallelic variants and disease progression are investigated. The introduction of the generalized odds ratio (ORG) may be a remedy. The ORG utilizes the complete genotype distribution and it provides an estimate of the magnitude of the association, given that the mutational load and/or the phenotype are treated as a graded exposure and/or outcome. The performance of the ORG was tested in 13 meta-analyses with binary outcomes (12 with biallelic and one 3-allelic variants) and in one meta-analysis that investigated disease progression. Six biallelic meta-analyses produced a significant ORG, indicating higher risk of disease given that the diseased subjects have a higher mutational load compared to the non-diseased ones. Four of the six meta-analyses showed significance for all genetic models. The multiallelic meta-analysis produced a significant ORG, indicating that the mutational load is implicated in disease susceptibility; on the contrary, the multiple genetic models produced diverse results. In the disease progression meta-analysis, the risk of progression was related to mutational load of the variant whereas the conventional analysis did not reveal this association. Application of the ORG could overcome the shortcomings of multiple model testing or erroneous model specification and provides an alternative and robust way for genetic association testing.

Cellular functions of living organisms are carried out through complex systems of interacting components. Including such interactions in the analysis, and considering sub-systems defined by biological pathways instead of individual components (e.g. genes), can lead to new findings about complex biological mechanisms. Networks are often used to capture such interactions and can be incorporated in models to improve the efficiency in estimation and inference. In this paper, we propose a model for incorporating external information about interactions among genes (proteins/metabolites) in differential analysis of gene sets. We exploit the framework of mixed linear models and propose a flexible inference procedure for analysis of changes in biological pathways. The proposed method facilitates the analysis of complex experiments, including multiple experimental conditions and temporal correlations among observations. We propose an efficient iterative algorithm for estimation of the model parameters and show that the proposed framework is asymptotically robust to the presence of noise in the network information. The performance of the proposed model is illustrated through the analysis of gene expression data for environmental stress response (ESR) in yeast, as well as simulated data sets.

Research on analyzing microarray data has focused on the problem of identifying differentially expressed genes to the neglect of the problem of how to integrate evidence that a gene is differentially expressed with information on the extent of its differential expression. Consequently, researchers currently prioritize genes for further study either on the basis of volcano plots or, more commonly, according to simple estimates of the fold change after filtering the genes with an arbitrary statistical significance threshold. While the subjective and informal nature of the former practice precludes quantification of its reliability, the latter practice is equivalent to using a hard-threshold estimator of the expression ratio that is not known to perform well in terms of mean-squared error, the sum of estimator variance and squared estimator bias. On the basis of two distinct simulation studies and data from different microarray studies, we systematically compared the performance of several estimators representing both current practice and shrinkage. We find that the threshold-based estimators usually perform worse than the maximum-likelihood estimator (MLE) and they often perform far worse as quantified by estimated mean-squared risk. By contrast, the shrinkage estimators tend to perform as well as or better than the MLE and never much worse than the MLE, as expected from what is known about shrinkage. However, a Bayesian measure of performance based on the prior information that few genes are differentially expressed indicates that hard-threshold estimators perform about as well as the local false discovery rate (FDR), the best of the shrinkage estimators studied. Based on the ability of the latter to leverage information across genes, we conclude that the use of the local-FDR estimator of the fold change instead of informal or threshold-based combinations of statistical tests and non-shrinkage estimators can be expected to substantially improve the reliability of gene prioritization at very little risk of doing so less reliably. Since the proposed replacement of post-selection estimates with shrunken estimates applies as well to other types of high-dimensional data, it could also improve the analysis of SNP data from genome-wide association studies.

There has been increasing interest in predicting patients' survival after therapy by investigating gene expression microarray data. In the regression and classification models with high-dimensional genomic data, boosting has been successfully applied to build accurate predictive models and conduct variable selection simultaneously. We propose the Buckley-James boosting for the semiparametric accelerated failure time models with right censored survival data, which can be used to predict survival of future patients using the high-dimensional genomic data. In the spirit of adaptive LASSO, twin boosting is also incorporated to fit more sparse models. The proposed methods have a unified approach to fit linear models, non-linear effects models with possible interactions. The methods can perform variable selection and parameter estimation simultaneously. The proposed methods are evaluated by simulations and applied to a recent microarray gene expression data set for patients with diffuse large B-cell lymphoma under the current gold standard therapy.

Regional co-expression refers to the phenomenon of contiguous genes exhibiting similar expression patterns. Among others, DNA copy number aberrations may be causally involved in regional co-expression. We propose a random coefficients model to explain regional co-expression from DNA copy number information, while modeling residual co-expression due to other causes by a correlated error structure. We show how the model parameters may be estimated (computationally efficient and consistently) from high-dimensional data, and suggest several robustifications of the estimation procedure. From the model we are able to assess whether there is a shared effect on expression levels due to the DNA copy number aberrations, but also whether this effect is homogeneous across genes. In two examples we use the proposed methodology to investigate the association between DNA copy number aberrations and regional co-expression.

We consider the problem of locating multiple interacting quantitative trait loci (QTL) influencing traits measured in counts. In many applications the distribution of the count variable has a spike at zero. Zero-inflated generalized Poisson regression (ZIGPR) allows for an additional probability mass at zero and hence an improvement in the detection of significant loci. Classical model selection criteria often overestimate the QTL number. Therefore, modified versions of the Bayesian Information Criterion (mBIC and EBIC) were successfully used for QTL mapping. We apply these criteria based on ZIGPR as well as simpler models. An extensive simulation study shows their good power detecting QTL while controlling the false discovery rate. We illustrate how the inability of the Poisson distribution to account for over-dispersion leads to an overestimation of the QTL number and hence strongly discourages its application for identifying factors influencing count data. The proposed method is used to analyze the mice gallstone data of Lyons et al. (2003). Our results suggest the existence of a novel QTL on chromosome 4 interacting with another QTL previously identified on chromosome 5. We provide the corresponding code in R.

We present a new methodology for discriminating genomic symbolic sequences, which combines wavelet analysis and a self-organizing map algorithm. Wavelets are used to extract variation across various scales in the oligonucleotide patterns of a sequence. The variation is quantified by the estimated wavelet variance, which yields a feature vector. Feature vectors obtained from many genomic sequences, possibly of different lengths, are then classified with a nonparametric self-organizing map scheme. When applied to nearly 200 entire mitochondrial DNA sequences, or their fragments, the method predicts species taxonomic group membership very well, and allows the results to be visualized. When only thousands of nucleotides are available, wavelet-based feature vectors of short oligonucleotide patterns are more efficient in discrimination than frequency-based feature vectors of long patterns. This new data analysis strategy could be extended to numeric genomic data. The routines needed to perform the computations are readily available in two packages of software R.

We present a method for estimating and providing a confidence interval for the number of DNA replication origins in the genome of the yeast Kluyveromyces lactis. The method requires an initial set of verified sites from which a position specific frequency matrix (PSFM) can be constructed. We further assume that we have access to a sparingly used experimental procedure which can verify the functionality of a few, but not all, computationally predicted sites. While our motivation comes from estimating the number of autonomously replicating sequences (ARSs), our method can also be applied to estimating the genome-wide number of “functional” transcription factor binding sites, where functionality is determined by experimental verification of the transcription factor binding event using, for example, ChIP data. The reliability of our method is demonstrated by correctly predicting the known number of Saccharomyces cerevisiae ARSs as well as the number of S. cerevisiae probes that bind to the transcription factor ABF1.

High density tiling arrays are an effective strategy for genome-wide identification of transcription factor binding regions. Sliding window methods that calculate moving averages of log ratios or t-statistics have been useful for the analysis of tiling array data. Here, we present a method that generalizes the moving average approach to evaluate sliding windows of p-values by using combined p-value statistics. In particular, the combined p-value framework can be useful in situations when taking averages of the corresponding test-statistic for the hypothesis may not be appropriate or when it is difficult to assess the significance of these averages. We exhibit the strengths of the combined p-values methods on Drosophila tiling array data and assess their ability to predict genomic regions enriched for transcription factor binding. The predictions are evaluated based on their proximity to target genes and their enrichment of known transcription factor binding sites. We also present an application for the generalization of the moving average based on integrating two different tiling array experiments.

Statistical methods generating sparse models are of great value in the gene expression field, where the number of covariates (genes) under study moves about the thousands while the sample sizes seldom reach a hundred of individuals. For phenotype classification, we propose different lasso logistic regression approaches with specific penalizations for each gene. These methods are based on a generalized soft-threshold (GSoft) estimator. We also show that a recent algorithm for convex optimization, namely, the cyclic coordinate descent (CCD) algorithm, provides with a way to solve the optimization problem significantly faster than with other competing methods. Viewing GSoft as an iterative thresholding procedure allows us to get the asymptotic properties of the resulting estimates in a straightforward manner. Results are obtained for simulated and real data. The leukemia and colon datasets are commonly used to evaluate new statistical approaches, so they come in useful to establish comparisons with similar methods. Furthermore, biological meaning is extracted from the leukemia results, and compared with previous studies. In summary, the approaches presented here give rise to sparse, interpretable models that are competitive with similar methods developed in the field.

Granger causality (GC) tests are ideally suited to investigate time series data generated by bivariate vector autoregressive (VAR) processes. Recent studies have applied GC analysis and its extensions for modeling functional relationships and network structure from temporal gene expression profiles. The present study investigates GC analysis of human cell-cycle gene expression profiles that can be modeled as a first-order bivariate VAR. Analytical results presented establish the contribution of the VAR process parameters, including auto-regulatory feedback and noise variance to the mean-squared forecast error, as a critical component in identifying statistically significant GC relationships. These results in turn discourage blind inference of functional relationship between a given pair of genes solely based on the result of the statistical tests for GC. The presence of significant auto-regulatory feedback and discrepancy in noise variance is demonstrated across the cell-cycle gene expression profiles by VAR parameter estimation. It is emphasized that discrepancies in noise variance can be due to artifacts and can lead to spurious existence of functional relationship between a given pair of genes. VAR parameter estimation is encouraged for better of GC interpretation of the results. Published case studies on GC analysis of the same publicly available cell-cycle gene expression data are reinvestigated for transparency.

The genetic control of a complex trait can be studied by testing and mapping the genotypes of the underlying quantitative trait loci (QTLs) through their associations with observable marker genotypes. All existing statistical methods for QTL mapping assume an equilibrium population, allowing marker-QTL associations to be simply described at the gametic level. However, many mapping populations in practice may deviate from equilibrium; thus, gametic associations cannot reflect marker-QTL associations at the genotype level. We develop a robust model for mapping QTLs in a non-equilibrium natural population in which individuals are not necessarily randomly mating due to various evolutionary forces. Without use of Hardy-Weinberg equilibrium, the new model founds marker-QTL associations directly on the genotypes, specified by a group of disequilibrium parameters. Simulation studies were performed to test the statistical properties of the model, which suggests that the new model covers current mapping models and can be safely used for any data set.

The recent emergence of massively parallel sequencing technologies has enabled an increasing number of human genome re-sequencing studies, notable among them being the 1000 Genomes Project. The main aim of these studies is to identify the yet unknown genetic variants in a genomic region, mostly low frequency variants (frequency less than 5%). We propose here a set of statistical tools that address how to optimally design such studies in order to increase the number of genetic variants we expect to discover. Within this framework, the tradeoff between lower coverage for more individuals and higher coverage for fewer individuals can be naturally solved.The methods here are also useful for estimating the number of genetic variants missed in a discovery study performed at low coverage.We show applications to simulated data based on coalescent models and to sequence data from the ENCODE project. In particular, we show the extent to which combining data from multiple populations in a discovery study may increase the number of genetic variants identified relative to studies on single populations.

How best to summarize large and complex datasets is a problem that arises in many areas of science. We approach it from the point of view of seeking data summaries that minimize the average squared error of the posterior distribution for a parameter of interest under approximate Bayesian computation (ABC). In ABC, simulation under the model replaces computation of the likelihood, which is convenient for many complex models. Simulated and observed datasets are usually compared using summary statistics, typically in practice chosen on the basis of the investigator's intuition and established practice in the field. We propose two algorithms for automated choice of efficient data summaries. Firstly, we motivate minimisation of the estimated entropy of the posterior approximation as a heuristic for the selection of summary statistics. Secondly, we propose a two-stage procedure: the minimum-entropy algorithm is used to identify simulated datasets close to that observed, and these are each successively regarded as observed datasets for which the mean root integrated squared error of the ABC posterior approximation is minimized over sets of summary statistics. In a simulation study, we both singly and jointly inferred the scaled mutation and recombination parameters from a population sample of DNA sequences. The computationally-fast minimum entropy algorithm showed a modest improvement over existing methods while our two-stage procedure showed substantial and highly-significant further improvement for both univariate and bivariate inferences. We found that the optimal set of summary statistics was highly dataset specific, suggesting that more generally there may be no globally-optimal choice, which argues for a new selection for each dataset even if the model and target of inference are unchanged.

Complex diseases will have multiple functional sites, and it will be invaluable to understand the cross-locus interaction in terms of linkage disequilibrium (LD) between those sites (epistasis) in addition to the haplotype-LD effects. We investigated the statistical properties of a class of matrix-based statistics to assess this epistasis. These statistical methods include two LD contrast tests (Zaykin et al., 2006) and partial least squares regression (Wang et al., 2008). To estimate Type 1 error rates and power, we simulated multiple two-variant disease models using the SIMLA software package. SIMLA allows for the joint action of up to two disease genes in the simulated data with all possible multiplicative interaction effects between them. Our goal was to detect an interaction between multiple disease-causing variants by means of their linkage disequilibrium (LD) patterns with other markers. We measured the effects of marginal disease effect size, haplotype LD, disease prevalence and minor allele frequency have on cross-locus interaction (epistasis).In the setting of strong allele effects and strong interaction, the correlation between the two disease genes was weak (r = 0.2). In a complex system with multiple correlations (both marginal and interaction), it was difficult to determine the source of a significant result. Despite these complications, the partial least squares and modified LD contrast methods maintained adequate power to detect the epistatic effects; however, for many of the analyses we often could not separate interaction from a strong marginal effect. While we did not exhaust the entire parameter space of possible models, we do provide guidance on the effects that population parameters have on cross-locus interaction.

As a consequence of the “large p small n” characteristic for microarray data, hypothesis tests based on individual genes often result in low average power. There are several proposed tests that attempt to improve power. Among these, the FS test that was developed using the concept of James-Stein shrinkage to estimate the variances showed a striking average power improvement. In this paper, we establish a framework in which we model the key parameters with a distribution to find an optimal Bayes test which we call the MAP test (where MAP stands for Maximum Average Power). Under this framework, the FS test can be derived as an empirical Bayes test approximating the MAP test corresponding to modeling the variances. By modeling both the means and the variances with a distribution, a MAP statistic is derived which is optimal in terms of average power but is computationally intensive. An empirical Bayes test called the FSS test is derived as an approximation to the MAP tests and can be computed instantaneously. The FSS statistic shrinks both the means and the variances and has numerically identical average power to the MAP tests. Much numerical evidence is presented in this paper that shows that the proposed test performs uniformly better in average power than the other tests in the literature, including the classical F test, the FS test, the test of Wright and Simon, the moderated t-test, SAM, Efron's t test, the B-statistic and Storey's optimal discovery procedure. A theory is established which indicates that the proposed test is optimal in power when controlling the false discovery rate (FDR).

High correlations were obtained between probes in seemingly unrelated probe sets, following an examination of the data from thousands of Affymetrix GeneChips. Investigation revealed that these unexpected correlations were between probes that were adjacent to high-valued probes. Using carefully selected probes, together with simple linear models, the extent of blur has been measured for each CEL file. The cause is shown to be attributable to poorly performing scanners. Blur can result in the doubling of the values of thousands of probes. This in turn can lead to the doubling of the expression level for hundreds of probe sets.

Quantitative trait loci (QTL) mapping often results in data on a number of traits that have well-established causal relationships. Many multi-trait QTL mapping methods that account for the correlation among multiple traits have been developed to improve the statistical power and the precision of QTL parameter estimation. However, none of these methods are capable of incorporating the causal structure among the traits. Consequently, genetic functions of the QTL may not be fully understood. Structural equation modeling (SEM) allows researchers to explicitly characterize the causal structure among the variables and to decompose effects into direct, indirect, and total effects. In this paper, we developed a multi-trait SEM method of QTL mapping that takes into account the causal relationships among traits related to grain yield. Performance of the proposed method is evaluated by simulation study and applied to data from a wheat experiment. Compared with single trait analysis and the multi-trait least-squares analysis, our multi-trait SEM improves statistical power of QTL detection and provides important insight into how QTLs regulate traits by investigating the direct, indirect, and total QTL effects. The approach also helps build biological models that more realistically reflect the complex relationships among QTL and traits and is more precise and efficient in QTL mapping than single trait analysis.

We propose a new approach for clustering DNA features using array CGH data from multiple tumor samples. We distinguish data-collapsing (joining contiguous DNA clones or probes with extremely similar data into regions) from clustering (joining contiguous, correlated regions based on a maximum likelihood principle). The model-based clustering algorithm accounts for the apparent spatial patterns in the data. We evaluate the randomness of the clustering result by a cluster stability score in combination with cross-validation. Moreover, we argue that the clustering really captures spatial genomic dependency by showing that coincidental clustering of independent regions is very unlikely.Using the region and cluster information, we combine testing of these for association with a clinical variable in a hierarchical multiple testing approach. This allows for interpreting the significance of both regions and clusters while controlling the Family-Wise Error Rate simultaneously. We prove that in the context of permutation tests and permutation-invariant clusters it is allowed to perform clustering and testing on the same data set. Our procedures are illustrated on two cancer data sets.

Characterizing dynamic gene expression pattern and predicting patient outcome is now significant and will be of more interest in the future with large scale clinical investigation of microarrays. However, there is currently no method that has been developed for prediction of patient outcome using longitudinal gene expression, where gene expression of patients is being monitored across time. Here, we propose a novel prediction approach for patient survival time that makes use of time course structure of gene expression. This method is applied to a burn study. The genes involved in the final predictors are enriched in the inflammatory response and immune system related pathways. Moreover, our method is consistently better than prediction methods using individual time point gene expression or simply pooling gene expression from each time point.

Probabilistic mixture models provide a popular approach to cluster noisy gene expression data for exploring gene function. Since gene expression data obtained from microarray experiments are often associated with significant sources of technical and biological noise, replicated experiments are typically used to deal with data variability, and internal replication (e.g. from multiple probes per gene in an experiment) provides valuable information about technical sources of noise. However, current implementations of mixture models either do not consider the correlation between the replicated measurements for the same experimental condition, or ignore the probe-level measurement error, and thus overlook the rich information about technical noise. Moreover, most current methods use non-robust Gaussian components to describe the data, and these methods are therefore sensitive to non-Gaussian clusters and outliers. In many cases, this will lead to over-estimation of the number of model components as multiple Gaussian components are used to fit a non-Gaussian cluster. We propose a robust Student's t-mixture model, which explicitly handles replicated gene expression data, includes the consideration of probe-level measurement error when available and automatically selects the appropriate number of model components using a minimum message length criterion. We apply the model to gene expression data using probe-level measurements from an Affymetrix probe-level model, multi-mgMOS, which provides uncertainty estimates. The proposed Student's t-mixture model shows robust performance on synthetic data sets with realistic noise characteristics in comparison to a standard Gaussian mixture model and two other previously published methods. We also compare performance with these methods on two yeast time-course data sets and show that the new method obtains more biologically meaningful clusters in terms of enrichment statistics for GO categories and interactions between transcription factors and genes. Automatically selecting the number of components is more computationally efficient than using a model selection approach and allows the methods to be applied to larger data sets.

Klebanov et al. (2006) proposed a new type of stochastic dependence, Type A dependence, between gene expression levels. They estimated the abundance of Type A pairs by testing the correlation coefficients of gene pairs. We propose a new model, hidden regulator dependence, as an alternative to Type A dependence. We show that the correlation based procedure proposed by Klebanov et al. (2006) fails to differentiate hidden regulator dependence from Type A dependence, although their probabilistic structures are quite different.

Permutation tests are amongst the most commonly used statistical tools in modern genomic research, a process by which p-values are attached to a test statistic by randomly permuting the sample or gene labels. Yet permutation p-values published in the genomic literature are often computed incorrectly, understated by about 1/m, where m is the number of permutations. The same is often true in the more general situation when Monte Carlo simulation is used to assign p-values. Although the p-value understatement is usually small in absolute terms, the implications can be serious in a multiple testing context. The understatement arises from the intuitive but mistaken idea of using permutation to estimate the tail probability of the test statistic. We argue instead that permutation should be viewed as generating an exact discrete null distribution. The relevant literature, some of which is likely to have been relatively inaccessible to the genomic community, is reviewed and summarized. A computation strategy is developed for exact p-values when permutations are randomly drawn. The strategy is valid for any number of permutations and samples. Some simple recommendations are made for the implementation of permutation tests in practice.