Volume 11, Issue 4

The goal of association mapping is to identify genetic variants that predict disease, and as the field of human genetics matures, the number of successful association studies is increasing. Many such studies have shown that for many diseases, risk is explained by a reasonably large number of variants that each explains a very small amount of disease risk. This is prompting the use of genetic risk scores in building predictive models, where information across several variants is combined for predictive modeling. In the current study, we compare the performance of four previously proposed genetic risk score methods and present a new method for constructing genetic risk score that incorporates explained variance information. The methods compared include: a simple count Genetic Risk Score, an odds ratio weighted Genetic Risk Score, a direct logistic regression Genetic Risk Score, a polygenic Genetic Risk Score, and the new explained variance weighted Genetic Risk Score. We compare the methods using a wide range of simulations in two steps, with a range of the number of deleterious single nucleotide polymorphisms (SNPs) explaining disease risk, genetic modes, baseline penetrances, sample sizes, relative risks (RR) and minor allele frequencies (MAF). Several measures of model performance were compared including overall power, C-statistic and Akaike’s Information Criterion. Our results show the relative performance of methods differs significantly, with the new explained variance weighted GRS (EV-GRS) generally performing favorably to the other methods.

We analytically derive the first and second derivatives of the likelihood in maximum likelihood methods for phylogeny. These results enable the Newton-Raphson method to be used for maximising likelihood, which is important because there is a need for faster methods for optimisation of parameters in maximum likelihood methods. Furthermore, the calculation of the Hessian matrix also opens up possibilities for standard likelihood theory to be applied, for inference in phylogeny and for model selection problems. Another application of the Hessian matrix is local influence analysis, which can be used for detecting a number of biologically interesting phenomena. The pruning algorithm has been used to speed up computation of likelihoods for a tree. We explain how it can be used to speed up the computation for the first and second derivatives of the likelihood with respect to branch lengths and other parameters. The results in this paper apply not only to bifurcating trees, but also to general multifurcating trees. We demonstrate the use of our Hessian calculation for the three applications listed above, and compare with existing methods for those applications.

The task of analyzing high-dimensional single nucleotide polymorphism (SNP) data in a case-control design using multivariable techniques has only recently been tackled. While many available approaches investigate only main effects in a high-dimensional setting, we propose a more flexible technique, cluster-localized regression (CLR), based on localized logistic regression models, that allows different SNPs to have an effect for different groups of individuals. Separate multivariable regression models are fitted for the different groups of individuals by incorporating weights into componentwise boosting, which provides simultaneous variable selection, hence sparse fits. For model fitting, these groups of individuals are identified using a clustering approach, where each group may be defined via different SNPs. This allows for representing complex interaction patterns, such as compositional epistasis, that might not be detected by a single main effects model. In a simulation study, the CLR approach results in improved prediction performance, compared to the main effects approach, and identification of important SNPs in several scenarios. Improved prediction performance is also obtained for an application example considering urinary bladder cancer. Some of the identified SNPs are predictive for all individuals, while others are only relevant for a specific group. Together with the sets of SNPs that define the groups, potential interaction patterns are uncovered.

We study exact tests for (2 x 2) and (2 x 3) contingency tables, in particular exact chi-squared tests and exact tests of Fisher type. In practice, these tests are typically carried out without randomization, leading to reproducible results but not exhausting the significance level. We discuss that this can lead to methodological and practical issues in a multiple testing framework when many tables are simultaneously under consideration as in genetic association studies.Realized randomized p-values are proposed as a solution which is especially useful for data-adaptive (plug-in) procedures. These p-values allow to estimate the proportion of true null hypotheses much more accurately than their non-randomized counterparts. Moreover, we address the problem of positively correlated p-values for association by considering techniques to reduce multiplicity by estimating the "effective number of tests" from the correlation structure.An algorithm is provided that bundles all these aspects, efficient computer implementations are made available, a small-scale simulation study is presented and two real data examples are shown.

For the problem of multiple testing, the Benjamini-Hochberg (B-H) procedure has become a very popular method in applications. We show how the B-H procedure can be interpreted as a test based on the spacings corresponding to the p-value distributions. This interpretation leads to the incorporation of the empirical null hypothesis, a term coined by Efron (2004). We develop a mixture modelling approach for the empirical null hypothesis for the B-H procedure and demonstrate some theoretical results regarding both finite-sample as well as asymptotic control of the false discovery rate. The methodology is illustrated with application to two high-throughput datasets as well as to simulated data.

Mutations that confer a selective advantage to an organism are the raw material upon which natural selection acts. The number of such mutations that are available is a central quantity of interest for understanding the tempo and trajectory of adaptive evolution. While this quantity is typically unknown, it can be estimated with varying levels of accuracy based on data obtained experimentally. We propose a method for estimating the number of beneficial mutations that accounts for the evolutionary forces that generate the data. Our model-based parametric approach is compared to an adjusted nonparametric abundance-based coverage estimator. We show that, in general, our estimator performs better. When the number of mutations is small, however, the performances of the two estimators are similar.

Cancer patients often develop multiple malignancies that may be either metastatic spread of a previous cancer (clonal tumors) or new primary cancers (independent tumors). If diagnosis cannot be easily made on the basis of the pathology review, the patterns of somatic mutations in the tumors can be compared. Previously we have developed statistical methods for testing clonality of two tumors using their loss of heterozygosity (LOH) profiles at several candidate markers. These methods can be applied to all possible pairs of tumors when multiple tumors are analyzed, but this strategy can lead to inconsistent results and loss of statistical power. In this work we will extend clonality tests to three and more malignancies from the same patient. A non-parametric test can be performed using any possible subset of tumors, with the subsequent adjustment for multiple testing. A parametric likelihood model is developed for 3 or 4 tumors, and it can be used to estimate the phylogenetic tree of tumors. The proposed tests are more powerful than combination of all possible pairwise tests.

Viruses such as HIV and Hepatitis C (HCV) replicate rapidly and with high transcription error rates, which may facilitate their escape from immune detection through the encoding of mutations at key positions within human leukocyte antigen (HLA)-specific peptides, thus impeding T-cell recognition. Large-scale population-based host-viral association studies are conducted as hypothesis-generating analyses which aim to determine the positions within the viral sequence at which host HLA immune pressure may have led to these viral escape mutations. When transmission of the virus to the host is HLA-associated, however, standard tests of association can be confounded by the viral relatedness of contemporarily circulating viral sequences, as viral sequences descended from a common ancestor may share inherited patterns of polymorphisms, termed ‘founder effects’. Recognizing the correspondence between this problem and the confounding of case-control genome-wide association studies by population stratification, we adapt methods taken from that field to the analysis of host-viral associations. In particular, we consider methods based on principal components analysis within a logistic regression framework motivated by alternative formulations in the Frisch-Waugh-Lovell Theorem. We demonstrate via simulation their utility in detecting true host-viral associations whilst minimizing confounding by associations generated by founder effects. The proposed methods incorporate relatively robust, standard statistical procedures which can be easily implemented using widely available software, and provide alternatives to the more complex computer intensive methods often implemented in this area.

An important and challenging problem in systems biology is the inference of gene regulatory networks from short non-stationary time series of transcriptional profiles. A popular approach that has been widely applied to this end is based on dynamic Bayesian networks (DBNs), although traditional homogeneous DBNs fail to model the non-stationarity and time-varying nature of the gene regulatory processes. Various authors have therefore recently proposed combining DBNs with multiple changepoint processes to obtain time varying dynamic Bayesian networks (TV-DBNs). However, TV-DBNs are not without problems. Gene expression time series are typically short, which leaves the model over-flexible, leading to over-fitting or inflated inference uncertainty. In the present paper, we introduce a Bayesian regularization scheme that addresses this difficulty. Our approach is based on the rationale that changes in gene regulatory processes appear gradually during an organism's life cycle or in response to a changing environment, and we have integrated this notion in the prior distribution of the TV-DBN parameters. We have extensively tested our regularized TV-DBN model on synthetic data, in which we have simulated short non-homogeneous time series produced from a system subject to gradual change. We have then applied our method to real-world gene expression time series, measured during the life cycle of Drosophila melanogaster, under artificially generated constant light condition in Arabidopsis thaliana, and from a synthetically designed strain of Saccharomyces cerevisiae exposed to a changing environment.

Recently, expression quantitative loci (eQTL) mapping studies, where expression levels of thousands of genes are viewed as quantitative traits, have been used to provide greater insight into the biology of gene regulation. Originally, eQTLs were detected by applying standard QTL detection tools (using a “one gene at-a-time” approach), but this method ignores many possible interactions between genes. Several other methods have proposed to overcome these limitations, but each of them has some specific disadvantages. In this paper, we present an integrated hierarchical Bayesian model that jointly models all genes and SNPs to detect eQTLs. We propose a model (named iBMQ) that is specifically designed to handle a large number G of gene expressions, a large number S of regressors (genetic markers) and a small number n of individuals in what we call a ``large G, large S, small n'' paradigm. This method incorporates genotypic and gene expression data into a single model while 1) specifically coping with the high dimensionality of eQTL data (large number of genes), 2) borrowing strength from all gene expression data for the mapping procedures, and 3) controlling the number of false positives to a desirable level. To validate our model, we have performed simulation studies and showed that it outperforms other popular methods for eQTL detection, including QTLBIM, R-QTL, remMap and M-SPLS. Finally, we used our model to analyze a real expression dataset obtained in a panel of mice BXD Recombinant Inbred (RI) strains. Analysis of these data with iBMQ revealed the presence of multiple hotspots showing significant enrichment in genes belonging to one or more annotation categories.

Background: Among the most commonly applied microarray normalization methods are intensity-dependent normalization methods such as lowess or loess algorithms. Their computational complexity makes them slow and thus less suitable for normalization of large datasets. Current implementations try to circumvent this problem by using a random subset of the data for normalization, but the impact of this modification has not been previously assessed. We developed a novel intensity-dependent normalization method for microarrays that is fast, simple and can include weighing of observations. Results: Our normalization method is based on the P-spline scatterplot smoother using all data points for normalization. We show that using a random subset of the data for normalization should be avoided as unstable results can be produced. However, in certain cases normalization based on an invariant subset is desirable, for example, when groups of samples before and after intervention are compared. We show in the context of DNA methylation arrays that a constant weighted P-spline normalization yields a more reliable normalization curve than the one obtained by normalization on the invariant subset only. Conclusions: Our novel intensity-dependent normalization method is simpler and faster than current loess algorithms, and can be applied to one- and two-colour array data, similar to normalization based on loess. Availability: An implementation of the method is currently available as an R package called TurboNorm from www.bioconductor.org .

We evaluate power performance to detect the correct mode of inheritance in gene-disease associations of two different approaches: the MAX test and the degree of dominance index or h-index. The MAX test is a special case of the conditional independence tests that simultaneously test for association and select the most likely genetic model based on a three-dimensional normal distribution. The h-index is based on the philosophy of using orthogonal contrasts to infer the mode of inheritance quantitatively. A population genetic model is developed where the real mode of inheritance is known a priori and power performance can be accurately determined. The simulations showed that none of the two approaches generally outperforms the other, nor each of them provides a panacea to estimate efficiently the mode of inheritance in all parameter space. However, the simultaneous application of both approaches can provide insights in determining the underlying mode of inheritance.

When modeling time course microarray data special interest may reside in identifying time frames in which gene expression levels follow a monotonic (increasing or decreasing) trend. A trajectory may change its regime because of the reaction to treatment or of a natural developmental phase, as in our motivating example about identification of genes involved in embryo development of mice with the 22q11 deletion. To this aim we propose a new flexible Bayesian autoregressive hidden Markov model based on three latent states, corresponding to stationarity, to an increasing and to a decreasing trend. In order to select a list of genes, we propose decision criteria based on the posterior distribution of the parameters of interest, taking into account the uncertainty in parameter estimates. We also compare the proposed model with two simpler models based on constrained formulations of the probability transition matrix.

The problem of locating quantitative trait loci (QTL) for experimental populations can be approached by multiple regression analysis. In this context variable selection using a modification of the Bayesian Information Criterion (mBIC) has been well established in the past. In this article a memetic algorithm (MA) is introduced to find the model which minimizes the selection criterion. Apart from mBIC also a second modification (mBIC2) is considered, which has the property of controlling the false discovery rate. Given the Bayesian nature of our selection criteria, we are not only interested in finding the best model, but also in computing marker posterior probabilities using all models visited by MA. In a simulation study MA (with mBIC and mBIC2) is compared with a parallel genetic algorithm (PGA) which has been previously suggested for QTL mapping. It turns out that MA in combination with mBIC2 performs best, where determining QTL positions based on marker posterior probabilities yields even better results than using the best model selected by MA. Finally we consider a real data set from the literature and show that MA can also be extended to multiple interval mapping, which potentially increases the precision with which the exact location of QTLs can be estimated.

Score-based pairwise alignments are widely used in bioinformatics in particular with molecular database search tools, such as the BLAST family. Due to sophisticated heuristics, such algorithms are usually fast but the underlying scoring model unfortunately lacks a statistical description of the reliability of the reported alignments. In particular, close to gaps, in low-score or low-complexity regions, a huge number of alternative alignments arise which results in a decrease of the certainty of the alignment. ppALIGN is a software package that uses hidden Markov Model techniques to compute position-wise reliability of score-based pairwise alignments of DNA or protein sequences. The design of the model allows for a direct connection between the scoring function and the parameters of the probabilistic model. For this reason it is suitable to analyze the outcomes of popular score based aligners and search tools without having to choose a complicated set of parameters. By contrast, our program only requires the classical score parameters (the scoring function and gap costs). The package comes along with a library written in C++, a standalone program for user defined alignments (ppALIGN) and another program (ppBLAST) which can process a complete result set of BLAST. The main algorithms essentially exhibit a linear time complexity (in the alignment lengths), and they are hence suitable for on-line computations. We have also included alternative decoding algorithms to provide alternative alignments. ppALIGN is a fast program/library that helps detect and quantify questionable regions in pairwise alignments. Due to its structure, the input/output interface it can to be connected to other post-processing tools. Empirically, we illustrate its usefulness in terms of correctly predicted reliable regions for sequences generated using the ROSE model for sequence evolution, and identify sensor-specific regions in the denitrifying betaproteobacterium Aromatoleum aromaticum.