Using substitution matrices to estimate probability distributions for biological sequences.

Accurately estimating probabilities from observations is important for probabilistic-based approaches to problems in computational biology. In this paper we present a biologically-motivated method for estimating probability distributions over discrete alphabets from observations using a mixture model of common ancestors. The method is an extension of substitution matrix-based probability estimation methods. In contrast to previous such methods, our method has a simple Bayesian interpretation and has the advantage over Dirichlet mixtures that it is both effective and simple to compute for large alphabets. The method is applied to estimate amino acid probabilities based on observed counts in an alignment and is shown to perform comparably to previous methods. The method is also applied to estimate probability distributions over protein families and improves protein classification accuracy.     

FastMG: a simple,fast, and accurate maximum likelihood procedure to estimate amino acid replacement rate matrices from large data sets

Background

Amino acid replacement rate matrices are a crucial component of many protein analysis systems such as sequence similarity search, sequence alignment, and phylogenetic inference. Ideally, the rate matrix reflects the mutational behavior of the actual data under study; however, estimating amino acid replacement rate matrices requires large protein alignments and is computationally expensive and complex. As a compromise, sub-optimal pre-calculated generic matrices are typically used for protein-based phylogeny. Sequence availability has now grown to a point where problem-specific rate matrices can often be calculated if the computational cost can be controlled.

Results

The most time consuming step in estimating rate matrices by maximum likelihood is building maximum likelihood phylogenetic trees from protein alignments. We propose a new procedure, called FastMG, to overcome this obstacle. The key innovation is the alignment-splitting algorithm that splits alignments with many sequences into non-overlapping sub-alignments prior to estimating amino acid replacement rates. Experiments with different large data sets showed that the FastMG procedure was an order of magnitude faster than without splitting. Importantly, there was no apparent loss in matrix quality if an appropriate splitting procedure is used.

Conclusions

FastMG is a simple, fast and accurate procedure to estimate amino acid replacement rate matrices from large data sets. It enables researchers to study the evolutionary relationships for specific groups of proteins or taxa with optimized, data-specific amino acid replacement rate matrices. The programs, data sets, and the new mammalian mitochondrial protein rate matrix are available at http://fastmg.codeplex.com.     

A reliable and rapid procedure to estimate drug partitioning in biomembranes

A direct method using derivative spectrophotometry was developed for determining membrane-water molar partition coefficients (Kp) of the anticancer drugs tamoxifen (TAM) and 4-hydroxytamoxifen (OHTAM). This method explores a shift in the absorption spectra of the drugs when removed from the aqueous phase to a hydrophobic environment. Partition of TAM and OHTAM depends on membrane composition and on drug concentration, temperature and presence of cholesterol. Unlike OHTAM, partition of TAM in DMPC bilayers, liposomes of sarcoplasmic reticulum (SR) lipids and native membranes of SR and mitochondria decreases linearly with drug concentration. Additionally, the partition of these drugs is higher in SR native membranes than in liposomes of SR lipids. The partition also depends on membrane type, being higher in mitochondria than in SR membranes. Maximal partitionings in DMPC are observed at temperatures in the range of the main phase transition. Cholesterol strongly affects the incorporation of drugs and maximal inhibition was observed in DMPC bilayers.     

A stingless bee (Melipona seminigra) uses optic flow to estimate flight distances

Foragers of a stingless bee, Melipona seminigra, are able to use the optic flow experienced en route to estimate flight distance. After training the bees to collect food inside a flight tunnel with black-and-white stripes covering the side walls and the floor, their search behavior was observed in tunnels lacking a reward. Like honeybees, the bees accurately estimated the distance to the previously offered food source as seen from the sections of the tunnel where they turned around in search of the food. Changing the visual flow by decreasing the width of the flight tunnel resulted in the underestimation of the distance flown. The removal of image motion cues either in the ventral or lateral field of view reduced the bees' ability to gauge distances. When the feeder inside the tunnel was displaced together with the bees feeding on it while preventing the bee from seeing any image motion during the displacement the bees experienced different distances on their way to the food source and during their return to the nest. In the subsequent test the bees searched for the food predominantly at the distance associated with their return flight.     

A nonsingular optimization approach to the feed rate profile optimization of fedbatch cultures

A new method to calculate the optimal feed rate profile for fedbatch culture is proposed. Instead of the usual singular control approach of taking the feed rate as the control variable, the substrate concentration profile is used as the transformed control variable to avoid the computational difficulty associated with the singular control. Thus, the problem is converted into a nonsingular optimization problem of determining the optimal substrate concentration profile subject to a constraint. The equivalent feed rate profile to match the optimal substrate concentration profile is then generated. With this method the computational difficulty associated with singular controls for high-order systems is circumvented. The proposed method is illustrated by a number of examples.     

A model to estimate insulin sensitivity in dairy cows

Impairment of the insulin regulation of energy metabolism is considered to be an etiologic key component for metabolic disturbances. Methods for studies of insulin sensitivity thus are highly topical. There are clear indications that reduced insulin sensitivity contributes to the metabolic disturbances that occurs especially among obese lactating cows. Direct measurements of insulin sensitivity are laborious and not suitable for epidemiological studies. We have therefore adopted an indirect method originally developed for humans to estimate insulin sensitivity in dairy cows. The method, "Revised Quantitative Insulin Sensitivity Check Index" (RQUICKI) is based on plasma concentrations of glucose, insulin and free fatty acids (FFA) and it generates good and linear correlations with different estimates of insulin sensitivity in human populations. We hypothesized that the RQUICKI method could be used as an index of insulin function in lactating dairy cows. We calculated RQUICKI in 237 apparently healthy dairy cows from 20 commercial herds. All cows included were in their first 15 weeks of lactation. RQUICKI was not affected by the homeorhetic adaptations in energy metabolism that occurred during the first 15 weeks of lactation. In a cohort of 24 experimental cows fed in order to obtain different body condition at parturition RQUICKI was lower in early lactation in cows with a high body condition score suggesting disturbed insulin function in obese cows. The results indicate that RQUICKI might be used to identify lactating cows with disturbed insulin function.     

Simple stereological procedure to estimate the number and dimensions of root hairs

We demonstrate a simple procedure for systematic, uniformly random sampling of a root system of known length to obtain practically unbiased estimates of the total number and dimensions of root hairs. Irrespective of the length of the root system, only 100 root hairs need to be counted to estimate the total number with sufficient precision. Numbers and dimensions of root hairs were estimated for five crested wheatgrass (Agropyron cristatum L.) root systems that had been grown for one month in a gel. Less than one hour was required to obtain estimates of root hair parameters for a single plant. There was low variability of spatial density of root hairs within a given branching order (CV < 15%); however, because of large variation in the total length of laterals, the total number of root hairs varied greatly (CV ∼70%). On average, root hairs provided half of the total surface area of a root system and a total length 20 times that of the roots. This revised version was published online in June 2006 with corrections to the Cover Date.     

A proprioception based regulation model to estimate the trunk muscle forces

Evaluation of loads acting on the spine requires the knowledge of the muscular forces acting on it, but muscles redundancy necessitates developing a muscle forces attribution strategy. Optimisation, EMG, or hybrid models allow evaluating muscle force patterns, yielding a unique muscular arrangement or/and requiring EMG data collection. This paper presents a regulation model of the trunk muscles based on a proprioception hypothesis, which searches to avoid the spinal joint overloading. The model is also compared to other existing models for evaluation. Compared to an optimisation model, the proposed alternative muscle pattern yielded a significant spine postero-anterior shear decrease. Compared to a model based on combination of optimisation criteria, present model better fits muscle activation observed using EMG (38% improvement). Such results suggest that the proposed model, based on regulation of all spinal components, may be more relevant from a physiologic point of view.     

Conjugate momentum estimate using non-linear dynamic model of the sit-to-stand correlates well with accelerometric surface data

The purpose of this study was the development of a non-linear double inverted constrained pendulum model for the analysis of the movement of sit-to-stand (STS) transition. Ten able-bodied subjects perform five trials in their natural speed. Kinematics, kinetics as well as body worn accelerometer data were collected during the STS task using optoelectronic motion capture, force plate and inertial measurement unit, respectively. The conjugate momentum for the whole body which includes linear and angular motion correlates well with the accelerometric surface spanned by the accelerometer data. The partitioning of the conjugate momentum indicates a clear coordination between upper and lower limb after seat-off period. Moreover, the normalization procedure indicates a clear minimal and somehow invariant threshold value of the conjugate momentum to approximately 0.3 (body mass×body length) to perform the sit-to-stand for able-bodied subject. This threshold correlates well with the data obtained from accelerometeric index. The proposed accelerometric index is relevant to assess STS performance and to detect failed STS in clinics and outside a laboratory for patients with reduced mobility.     

A procedure to estimate okadaic acid in whole dinoflagellate cells using immunological techniques

A single procedure to detect and estimate okadaic acid in isolated whole cells was developed based on immunofluorescence and microscope photometry. This procedure allows the study of variations in okadaic acid concentration per cell although it is no substitute for HPLC procedures. Cells from mid-log exponential and stationary phase from two different clonal cultures of the okadaic-acid-producing dinoflagellate Prorocentrum lima (PI 5V and PI 7V) were analyzed. The results showed that: (1) cells from saturated phase cultures contain more okadaic acid than those from exponentially-growing mid-log phase; (2) genetic differences exist in okadaic acid production between the clones used; (3) okadaic acid is synthesized continuously during the whole cell cycle.     

A simplified procedure to determine the optimal rate of freezing biological systems

The effect of several cell-level parameters on the predicted optimal cooling rate B(opt) of an arbitrary biological system has been studied using a well-defined water transport model. An extensive investigation of the water transport model revealed three key cell level parameters: reference permeability of the membrane to water L(pg), apparent activation energy E(Lp), and the ratio of the available surface area for water transport to the initial volume of intracellular water (SA/WV). We defined B(opt) as the "highest" cooling rate at which a predefined percent of the initial water volume is trapped inside the cell (values ranging from 5% to 80%) at a predefined end temperature (values ranging from -5 degrees C to -40 degrees C). Irrespective of the choice of the percent of initial water volume trapped and the end temperature, an exact and linear relationship exists between L(pg), SA/WV, and B(opt0. However, a nonlinear and inverse relationship is found between E(Lp) and B(opt). Remarkably, for a variety of biological systems a comparison of the published experimentally determined values of B(opt) agreed quite closely with numerically predicted B(opt) values when the model assumed 5% of initial water is trapped inside the cell at a temperature of -15 degrees C. This close agreement between the experimental and model predicted optimal cooling rates is used to develop a generic optimal cooling rate chart and a generic optimal cooling rate equation that greatly simplifies the prediction of the optimal rate of freezing of biological systems.     

Scaffold filling under the breakpoint and related distances

Motivated by the trend of genome sequencing without completing the sequence of the whole genomes, a problem on filling an incomplete multichromosomal genome (or scaffold) I with respect to a complete target genome G was studied. The objective is to minimize the resulting genomic distance between I' and G, where I' is the corresponding filled scaffold. We call this problem the onesided scaffold filling problem. In this paper, we conduct a systematic study for the scaffold filling problem under the breakpoint distance and its variants, for both unichromosomal and multichromosomal genomes (with and without gene repetitions). When the input genome contains no gene repetition (i.e., is a fragment of a permutation), we show that the two-sided scaffold filling problem (i.e., G is also incomplete) is polynomially solvable for unichromosomal genomes under the breakpoint distance and for multichromosomal genomes under the genomic (or DCJ--Double-Cut-and-Join) distance. However, when the input genome contains some repeated genes, even the one-sided scaffold filling problem becomes NP-complete when the similarity measure is the maximum number of adjacencies between two sequences. For this problem, we also present efficient constant-factor approximation algorithms: factor-2 for the general case and factor 1.33 for the one-sided case.     

An improved multi-joint EMG-assisted optimization approach to estimate joint and muscle forces in a musculoskeletal model of the lumbar spine

Muscle force partitioning methods and musculoskeletal system simplifications are key modeling issues that can alter outcomes, and thus change conclusions and recommendations addressed to health and safety professionals. A critical modeling concern is the use of single-joint equilibrium to estimate muscle forces and joint loads in a multi-joint system, an unjustified simplification made by most lumbar spine biomechanical models. In the context of common occupational tasks, an EMG-assisted optimization method (EMGAO) is modified in this study to simultaneously account for the equilibrium at all lumbar joints (M-EMGAO). The results of this improved approach were compared to those of its conventional single-joint equivalent (S-EMGAO) counterpart, the latter method being applied to the same lumbar joints but one at a time. Despite identical geometrical configurations and passive contributions used in both models, computed outcomes clearly differed between single- and multi-joint methods, especially at larger trunk flexed postures and during asymmetric lifting. Moreover, muscle forces predicted by L5-S1 single-joint analyses do not maintain mechanical equilibrium at other spine joints crossed by the same muscles. Assuming that the central nervous system does not attempt to balance the external moments one joint at a time and that a given muscle cannot exert different forces at different joints, the proposed multi-joint method represents a substantial improvement over its single-joint counterpart. This improved approach, hence, resolves trunk muscle forces with biological integrity but without compromising mechanical equilibrium at the lumbar joints.     

Limits to the rate of adaptive substitution in sexual populations

In large populations, many beneficial mutations may be simultaneously available and may compete with one another, slowing adaptation. By finding the probability of fixation of a favorable allele in a simple model of a haploid sexual population, we find limits to the rate of adaptive substitution, [Formula: see text], that depend on simple parameter combinations. When variance in fitness is low and linkage is loose, the baseline rate of substitution is [Formula: see text], where [Formula: see text] is the population size, [Formula: see text] is the rate of beneficial mutations per genome, and [Formula: see text] is their mean selective advantage. Heritable variance [Formula: see text] in log fitness due to unlinked loci reduces [Formula: see text] by [Formula: see text] under polygamy and [Formula: see text] under monogamy. With a linear genetic map of length [Formula: see text] Morgans, interference is yet stronger. We use a scaling argument to show that the density of adaptive substitutions depends on [Formula: see text], [Formula: see text], [Formula: see text], and [Formula: see text] only through the baseline density: [Formula: see text]. Under the approximation that the interference due to different sweeps adds up, we show that [Formula: see text], implying that interference prevents the rate of adaptive substitution from exceeding one per centimorgan per 200 generations. Simulations and numerical calculations confirm the scaling argument and confirm the additive approximation for [Formula: see text]; for higher [Formula: see text], the rate of adaptation grows above [Formula: see text], but only very slowly. We also consider the effect of sweeps on neutral diversity and show that, while even occasional sweeps can greatly reduce neutral diversity, this effect saturates as sweeps become more common-diversity can be maintained even in populations experiencing very strong interference. Our results indicate that for some organisms the rate of adaptive substitution may be primarily recombination-limited, depending only weakly on the mutation supply and the strength of selection.     

A cooperative quay crane-based stochastic model to estimate vessel handling time

Having a good estimate of a vessel’s handling time is essential for planning and scheduling container terminal resources, such as berth positions, quay cranes (QCs) and transport vehicles. However, estimating the expected vessel handling time is not straightforward, because it depends on vessel characteristics, resource allocation decisions, and uncertainties in terminal processes. To estimate the expected vessel handling time, we propose a two-level stochastic model. The higher level model consists of a continuous-time Markov chain (CTMC) that captures the effect of QC assignment and scheduling on vessel handling time. The lower level model is a multi-class closed queuing network that models the dynamic interactions among the terminal resources and provides an estimate of the transition rate input parameters to the higher level CTMC model. We estimate the expected vessel handling times for several container load and unload profiles and discuss the effect of terminal layout parameters and crane service time variabilities on vessel handling times. From numerical experiments, we find that by having QCs cooperate, the vessel handling times are reduced by up to 15 %. The vessel handling time is strongly dependent on the variation in the QC service time and on the vehicle travel path topology.     

A simple procedure for estimating the false discovery rate

MOTIVATION: The most used criterion in microarray data analysis is nowadays the false discovery rate (FDR). In the framework of estimating procedures based on the marginal distribution of the P-values without any assumption on gene expression changes, estimators of the FDR are necessarily conservatively biased. Indeed, only an upper bound estimate can be obtained for the key quantity pi0, which is the probability for a gene to be unmodified. In this paper, we propose a novel family of estimators for pi0 that allows the calculation of FDR. RESULTS: The very simple method for estimating pi0 called LBE (Location Based Estimator) is presented together with results on its variability. Simulation results indicate that the proposed estimator performs well in finite sample and has the best mean square error in most of the cases as compared with the procedures QVALUE, BUM and SPLOSH. The different procedures are then applied to real datasets. AVAILABILITY: The R function LBE is available at http://ifr69.vjf.inserm.fr/lbe CONTACT: broet@vjf.inserm.fr.     

A gradient of silent substitution rate in the human pseudoautosomal region

It has been demonstrated that recombination in the human p-arm pseudoautosomal region (p-PAR) is at least twenty times more frequent than the genomic average of approximately 1 cM/Mb, which may affect substitution patterns and rates in this region. Here I report the analysis of substitution patterns and rates in 10 human, chimpanzee, gorilla, and orangutan genes across the p-PAR. Between species silent divergence in the p-PAR forms a gradient, increasing toward the telomere. The correlation of silent divergence with distance from the p-PAR boundary is highly significant (rho = 0.911, P < 0.001). After exclusion of the CpG dinucleotides this correlation is still significant (rho = 0.89, P < 0.01), thus the substitution rate gradient cannot be explained solely by the differences in the extent of methylation across the p-PAR. Frequent recombination in the PAR may result in a relatively strong effect of biased gene conversion (BGC), which, because of the increased probability of fixation of the G or C nucleotides at (A or T)/(G or C) segregating sites, may affect substitution rates. BGC, however, does not seem to be the factor creating the substitution rate gradient in the p-PAR, because the only gradient is still detactable if only A<-->T and G<-->C substitutions are taken into account (rho = 0.82, P < 0.01). I hypothesize that the substitution rate gradient in the p-PAR is due to the mutagenic effect of recombination, which is very frequent in the distal human p-PAR and might be lower near the p-PAR boundary.     

Using the nucleotide substitution rate matrix to detect horizontal gene transfer

Background  

Horizontal gene transfer (HGT) has allowed bacteria to evolve many new capabilities. Because transferred genes perform many medically important functions, such as conferring antibiotic resistance, improved detection of horizontally transferred genes from sequence data would be an important advance. Existing sequence-based methods for detecting HGT focus on changes in nucleotide composition or on differences between gene and genome phylogenies; these methods have high error rates.