The activities of most enzymes and drugs depend on interactions between proteins and small molecules. Accurate prediction of these interactions could greatly accelerate pharmaceutical and biotechnological research. Current machine learning models designed for this task have a limited ability to generalize beyond the proteins used for training. This limitation is likely due to a lack of information exchange between the protein and the small molecule during the generation of the required numerical representations. Here, we introduce ProSmith, a machine learning framework that employs a multimodal Transformer Network to simultaneously process protein amino acid sequences and small molecule strings in the same input. This approach facilitates the exchange of all relevant information between the two molecule types during the computation of their numerical representations, allowing the model to account for their structural and functional interactions. Our final model combines gradient boosting predictions based on the resulting multimodal Transformer Network with independent predictions based on separate deep learning representations of the proteins and small molecules. The resulting predictions outperform recently published state-of-the-art models for predicting protein-small molecule interactions across three diverse tasks: predicting kinase inhibitions; inferring potential substrates for enzymes; and predicting Michaelis constants KM. The Python code provided can be used to easily implement and improve machine learning predictions involving arbitrary protein-small molecule interactions.

The physiology of biological cells evolved under physical and chemical constraints, such as mass conservation across the network of biochemical reactions, nonlinear reaction kinetics, and limits on cell density. For unicellular organisms, the fitness that governs this evolution is mainly determined by the balanced cellular growth rate. We previously introduced growth balance analysis (GBA) as a general framework to model and analyze such nonlinear systems, revealing important analytical properties of optimal balanced growth states. It has been shown that at optimality, only a minimal subset of reactions can have nonzero flux. However, no general principles have been established to determine if a specific reaction is active at optimality. Here, we extend the GBA framework to study the optimality of each biochemical reaction, and we identify the mathematical conditions determining whether a reaction is active or not at optimal growth in a given environment. We reformulate the mathematical problem in terms of a minimal number of dimensionless variables and use the Karush-Kuhn-Tucker (KKT) conditions to identify fundamental principles of optimal resource allocation in GBA models of any size and complexity. Our approach helps to identify from first principles the economic values of biochemical reactions, expressed as marginal changes in cellular growth rate; these economic values can be related to the costs and benefits of proteome allocation into the reactions' catalysts. Our formulation also generalizes the concepts of Metabolic Control Analysis to models of growing cells. We show how the extended GBA framework unifies and extends previous approaches of cellular modeling and analysis, putting forward a program to analyze cellular growth through the stationarity conditions of a Lagrangian function. GBA thereby provides a general theoretical toolbox for the study of fundamental mathematical properties of balanced cellular growth.

A substantial fraction of the bacterial cytosol is occupied by catalysts and their substrates. While a higher volume density of catalysts and substrates might boost biochemical fluxes, the resulting molecular crowding can slow down diffusion, perturb the reactions' Gibbs free energies, and reduce the catalytic efficiency of proteins. Due to these tradeoffs, dry mass density likely possesses an optimum that facilitates maximal cellular growth and that is interdependent on the cytosolic molecule size distribution. Here, we analyze the balanced growth of a model cell, accounting systematically for crowding effects on reaction kinetics. Its optimal cytosolic volume occupancy depends on the nutrient-dependent resource allocation into large ribosomal vs. small metabolic macromolecules, reflecting a tradeoff between the saturation of metabolic enzymes, favoring larger occupancies with higher encounter rates, and the inhibition of the ribosomes, favoring lower occupancies with unhindered diffusion of tRNAs. Our predictions across growth rates are quantitatively consistent with the experimentally observed reduction in volume occupancy on rich media compared to minimal media in E. coli. Strong deviations from optimal cytosolic occupancy only lead to minute reductions in growth rate, which are nevertheless evolutionarily relevant due to large bacterial population sizes. In sum, cytosolic density variation in bacterial cells appears to be consistent with an optimality principle of cellular efficiency.

Summary: The fixation index FST can be used to identify non-neutrally evolving loci from genome-scale SNP data across two or more populations. Recent years have seen the development of sophisticated approaches to estimate FST based on Markov-Chain Monte-Carlo simulations. Here, we present a vectorized R implementation of an extension of the widely used BayeScan software for codominant markers, adding the option to group individual SNPs into pre-defined blocks. A typical application of this new approach is the identification of genomic regions, genes, or gene sets containing SNPs that evolved under directional selection. Availability and implementation: The R implementation of our method, which builds on the powerful population genetics and genomics software PopGenome, is available freely from CRAN. Supplementary information: Supplementary data are available at Bioinformatics online.

Energy metabolism is central to cellular biology. Thus, genome-scale models of heterotrophic unicellular species must account appropriately for the utilization of external nutrients to synthesize energy metabolites such as ATP. However, metabolic models designed for flux-balance analysis (FBA) may contain thermodynamically impossible energy-generating cycles: without nutrient consumption, these models are still capable of charging energy metabolites (such as ADP→ATP or NADP+→NADPH). Here, we show that energy-generating cycles occur in over 85% of metabolic models without extensive manual curation, such as those contained in the ModelSEED and MetaNetX databases; in contrast, such cycles are rare in the manually curated models of the BiGG database. Energy generating cycles may represent model errors, e.g., erroneous assumptions on reaction reversibilities. Alternatively, part of the cycle may be thermodynamically feasible in one environment, while the remainder is thermodynamically feasible in another environment; as standard FBA does not account for thermodynamics, combining these into an FBA model allows erroneous energy generation. The presence of energy-generating cycles typically inflates maximal biomass production rates by 25%, and may lead to biases in evolutionary simulations. We present efficient computational methods (i) to identify energy generating cycles, using FBA, and (ii) to identify minimal sets of model changes that eliminate them, using a variant of the GlobalFit algorithm.

Constraint-based metabolic modeling methods such as Flux Balance Analysis (FBA) are routinely used to predict the effects of genetic changes and to design strains with desired metabolic properties. The major bottleneck in modeling genome-scale metabolic systems is the establishment and manual curation of reliable stoichiometric models. Initial reconstructions are typically refined through comparisons to experimental growth data from gene knockouts or nutrient environments. Existing methods iteratively correct one erroneous model prediction at a time, resulting in accumulating network changes that are often not globally optimal. We present GlobalFit, a bi-level optimization method that finds a globally optimal network, by identifying the minimal set of network changes needed to correctly predict all experimentally observed growth and non-growth cases simultaneously. When applied to the genome-scale metabolic model of Mycoplasma genitalium, GlobalFit decreases unexplained gene knockout phenotypes by 79%, increasing accuracy from 87.3% (according to the current state-of-the-art) to 97.3%. While currently available computers do not allow a global optimization of the much larger metabolic network of E. coli, the main strengths of GlobalFit are already played out when considering only one growth and one non-growth case simultaneously. Application of a corresponding strategy halves the number of unexplained cases for the already highly curated E. coli model, increasing accuracy from 90.8% to 95.4%.

MOTIVATION: Constraint-based metabolic modeling methods such as Flux Balance Analysis (FBA) are routinely used to predict metabolic phenotypes, e.g. growth rates, ATP yield or the fitness of gene knockouts. One frequent difficulty of constraint-based solutions is the inclusion of thermodynamically infeasible loops (or internal cycles), which add nonbiological fluxes to the predictions. RESULTS: We propose a simple postprocessing of constraint-based solutions, which removes internal cycles from any given flux distribution [Formula: see text] without disturbing other fluxes not involved in the loops. This new algorithm, termed CycleFreeFlux, works by minimizing the sum of absolute fluxes [Formula: see text] while (i) conserving the exchange fluxes and (ii) using the fluxes of the original solution to bound the new flux distribution. This strategy reduces internal fluxes until at least one reaction of every possible internal cycle is inactive, a necessary and sufficient condition for the thermodynamic feasibility of a flux distribution. If alternative representations of the input flux distribution in terms of elementary flux modes exist that differ in their inclusion of internal cycles, then CycleFreeFlux is biased towards solutions that maintain the direction given by [Formula: see text] and towards solutions with lower total flux [Formula: see text]. Our method requires only one additional linear optimization, making it computationally very efficient compared to alternative strategies. AVAILABILITY AND IMPLEMENTATION: We provide freely available R implementations for the enumeration of thermodynamically infeasible cycles as well as for cycle-free FBA solutions, flux variability calculations and random sampling of solution spaces. CONTACT: [email protected].

SUMMARY: The statistical programming language R has become a de facto standard for the analysis of many types of biological data, and is well suited for the rapid development of new algorithms. However, variant call data from population-scale resequencing projects are typically too large to be read and processed efficiently with R's built-in I/O capabilities. WhopGenome can efficiently read whole-genome variation data stored in the widely used variant call format (VCF) file format into several R data types. VCF files can be accessed either on local hard drives or on remote servers. WhopGenome can associate variants with annotations such as those available from the UCSC genome browser, and can accelerate the reading process by filtering loci according to user-defined criteria. WhopGenome can also read other Tabix-indexed files and create indices to allow fast selective access to FASTA-formatted sequence files. AVAILABILITY AND IMPLEMENTATION: The WhopGenome R package is available on CRAN at http://cran.r-project.org/web/packages/WhopGenome/. A Bioconductor package has been submitted. CONTACT: [email protected].

BACKGROUND: Figures of phylogenetic trees are widely used to illustrate the result of evolutionary analyses. However, one cannot easily extract a machine-readable representation from such images. Therefore, new software emerges that helps to preserve phylogenies digitally for future research. RESULTS: TreeSnatcher Plus is a GUI-driven JAVA application that semi-automatically generates a Newick format for multifurcating, arbitrarily shaped, phylogenetic trees contained in pixel images. It offers a range of image pre-processing methods and detects the topology of a depicted tree with adequate user assistance. The user supervises the recognition process, makes corrections to the image and to the topology and repeats steps if necessary. At the end TreeSnatcher Plus produces a Newick tree code optionally including branch lengths for rectangular and freeform trees. CONCLUSIONS: Although illustrations of phylogenies exist in a vast number of styles, TreeSnatcher Plus imposes no limitations on the images it can process with adequate user assistance. Given that a fully automated digitization of all figures of phylogenetic trees is desirable but currently unrealistic, TreeSnatcher Plus is the only program that reliably facilitates at least a semi-automatic conversion from such figures into a machine-readable format.

Peter D. Keightley, Martin J. Lercher, Adam Eyre-Walker Recently, two groups have examined the level of sequence constraint in noncoding DNA flanking mammalian genes, and appear to have found conflicting results. By comparing 500-bp blocks in mice and rats, we found that mean nucleotide divergence within 2 kb of the start and stop codons of protein-coding genes is substantially lower than that of introns, and decreases when approaching the coding sequence [1]. If nucleotide changes within introns are largely free from selection, this implies that noncoding blocks close to genes evolved under selective constraints, presumably because they contain gene expression control regions. In contrast, we find that upstream sequences in hominids do not evolve slower than introns, while downstream regions are under about half of the constraint seen in murids [1]. By analysing a similar set of noncoding DNA sequences, Bush and Lahn also found that the mean level of selective constraints in upstream regions between humans and chimpanzees is very low. However, their slightly more complex main analysis was to search for 16-bp sequences within upstream regions that are strongly conserved between humans, mice, and either dogs or chickens. They then examined the divergence between humans and chimpanzees at the flanking nucleotides, finding substantially reduced divergence compared with the genomic mean. This demonstrated selective constraints at certain upstream sequences in hominids. An analogous analysis of mouse–rat sequences showed that the selective constraints are about twice as strong in murids as in hominids [2]. These two findings—on one hand, a near absence of selective constraints in blocks upstream of hominid genes [1], and on the other, evidence for strong selective constraints in these regions [2]—appear to contradict each other. How can we square the two sets of results? The answer is rather simple—windows with high conservation scores are relatively rare, and they contribute little to the mean calculated over 500-bp windows (unfortunately, Bush and Lahn do not tell us the fraction of 59 alignments within high conservation scores). Bush and Lahn also suggest that the apparent discrepancy ‘‘likely results from the fact that in large 500-bp blocks, functional elements that are under constraint are mixed with large sections of nonfunctional DNA, which are not under constraint’’ [2]. We believe that this interpretation, while formally correct, obscures important and interesting information that can be gained from combining the two studies. Some sequences outside the conserved 16-mers identified by Bush and Lahn are also likely to be functional, since the same 500-bp regions (largely ‘‘nonfunctional’’ according to Bush and Lahn) show strong evidence of evolutionary constraints between mice and rats [1]. Bush and Lahn also note that constraint, on either side of conserved windows, is greater in murids than in hominids. However, they observe a much smaller difference than that seen in our analysis. This is deceptive because by concentrating attention on regions that are conserved between humans, mice, and dogs, they ignore the fact that there might be many more highly conserved regions in murids than there are in hominids. In summary, there is no conflict between our results and those of Bush and Lahn; they concentrate their attention on a preselected subset of the sites we considered and so have a different perspective on the problem. What is clear from both studies is that there is a qualitative difference in the level of conservation in the 59 flanking sequences between murids and hominids. We have argued that this is likely to be due to the fixation of slightly deleterious mutations in hominids that are otherwise selectively eliminated in rodents. Differences in constraints between hominids and murids demonstrate that the overwhelming majority of changes at upstream regulatory sites have only small effects on fitness. This has counterintuitive consequences: to obtain a comprehensive list of human regulatory sites, it might be better to examine conservation in murid rather than hominid genomes. “

Modern multiring positron emission tomographs allow the acquisition of 3-D data sets to increase their sensitivity. A substantial part of this data is due to scattered radiation. The authors describe the experimental dependence of point source scatter distributions on energy window setting, source location, and scatter volume in geometries relevant for brain studies. The point source scatter distribution was parametrized accurately by a broad, 2-D Gaussian, which included a shift parameter to account for asymmetry of the scatter medium relative to the source. This parametrization was used to formulate two fast scatter correction algorithms suitable for brain scans. In both algorithms, a 2-D subset of the measured projections was transformed into a scatter projection. An image of the 3-D scatter distribution was reconstructed using 2-D algorithms. It was then subtracted from the total (true+scattered) 3-D image. Both algorithms were implemented in different combinations with the additional attenuation correction and were tested on point source and phantom measurements. It was shown that, for the situation typical for brain scans, reconstructed scatter fractions could be reduced to 5% or less.