Statistical genetics in and out of quasi-linkage equilibrium.

This review is about statistical genetics, an interdisciplinary topic between statistical physics and population biology. The focus is on the phase ofquasi-linkage equilibrium(QLE). Our goals here are to clarify under which conditions the QLE phase can be expected to hold in population biology and how the stability of the QLE phase is lost. The QLE state, which has many similarities to a thermal equilibrium state in statistical mechanics, was discovered by M Kimura for a two-locus two-allele model, and was extended and generalized to the global genome scale byNeher&Shraiman (2011). What we will refer to as the Kimura-Neher-Shraiman theory describes a population evolving due to the mutations, recombination, natural selection and possibly genetic drift. A QLE phase exists at sufficiently high recombination rate (r) and/or mutation ratesµwith respect to selection strength. We show how in QLE it is possible to infer the epistatic parameters of the fitness function from the knowledge of the (dynamical) distribution of genotypes in a population. We further consider the breakdown of the QLE regime for high enough selection strength. We review recent results for the selection-mutation and selection-recombination dynamics. Finally, we identify and characterize a new phase which we call the non-random coexistence where variability persists in the population without either fixating or disappearing.

Global analysis of more than 50,000 SARS-CoV-2 genomes reveals epistasis between eight viral genes.

Genome-wide epistasis analysis is a powerful tool to infer gene interactions, which can guide drug and vaccine development and lead to deeper understanding of microbial pathogenesis. We have considered all complete severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) genomes deposited in the Global Initiative on Sharing All Influenza Data (GISAID) repository until four different cutoff dates, and used direct coupling analysis together with an assumption of quasi-linkage equilibrium to infer epistatic contributions to fitness from polymorphic loci. We find eight interactions, of which three are between pairs where one locus lies in gene ORF3a, both loci holding nonsynonymous mutations. We also find interactions between two loci in gene nsp13, both holding nonsynonymous mutations, and four interactions involving one locus holding a synonymous mutation. Altogether, we infer interactions between loci in viral genes ORF3a and nsp2, nsp12, and nsp6, between ORF8 and nsp4, and between loci in genes nsp2, nsp13, and nsp14. The paper opens the prospect to use prominent epistatically linked pairs as a starting point to search for combinatorial weaknesses of recombinant viral pathogens.

Dynamic insights on transcription initiation and RNA processing during bacterial adaptation.

Transcription initiation and RNA processing govern gene expression and enable bacterial adaptation by reshaping the RNA landscape. The aim of this study was to simultaneously observe these two fundamental processes in a transcriptome responding to an environmental signal. A controlled σE system in E. coli was coupled to our previously described tagRNA-seq method to yield process kinetics information. Changes in transcription initiation frequencies (TIF) and RNA processing frequencies (PF) were followed using 5' RNA tags. Changes in TIF showed a binary increased/decreased pattern that alternated between transcriptionally activated and repressed promoters, providing the bacterial population with transcriptional oscillation. PF variation fell into three categories of cleavage activity: (i) constant and independent of RNA levels, (ii) increased once RNA has accumulated, and (iii) positively correlated to changes in TIF. This work provides a comprehensive and dynamic view of major events leading to transcriptomic reshaping during bacterial adaptation. It unveils an interplay between transcription initiation and the activity of specific RNA cleavage sites. This study utilized a well-known genetic system to analyze fundamental processes and can serve as a blueprint for comprehensive studies that exploit the RNA metabolism to decipher and understand bacterial gene expression control.

Interacting networks of resistance, virulence and core machinery genes identified by genome-wide epistasis analysis.

Recent advances in the scale and diversity of population genomic datasets for bacteria now provide the potential for genome-wide patterns of co-evolution to be studied at the resolution of individual bases. Here we describe a new statistical method, genomeDCA, which uses recent advances in computational structural biology to identify the polymorphic loci under the strongest co-evolutionary pressures. We apply genomeDCA to two large population data sets representing the major human pathogens Streptococcus pneumoniae (pneumococcus) and Streptococcus pyogenes (group A Streptococcus). For pneumococcus we identified 5,199 putative epistatic interactions between 1,936 sites. Over three-quarters of the links were between sites within the pbp2x, pbp1a and pbp2b genes, the sequences of which are critical in determining non-susceptibility to beta-lactam antibiotics. A network-based analysis found these genes were also coupled to that encoding dihydrofolate reductase, changes to which underlie trimethoprim resistance. Distinct from these antibiotic resistance genes, a large network component of 384 protein coding sequences encompassed many genes critical in basic cellular functions, while another distinct component included genes associated with virulence. The group A Streptococcus (GAS) data set population represents a clonal population with relatively little genetic variation and a high level of linkage disequilibrium across the genome. Despite this, we were able to pinpoint two RNA pseudouridine synthases, which were each strongly linked to a separate set of loci across the chromosome, representing biologically plausible targets of co-selection. The population genomic analysis method applied here identifies statistically significantly co-evolving locus pairs, potentially arising from fitness selection interdependence reflecting underlying protein-protein interactions, or genes whose product activities contribute to the same phenotype. This discovery approach greatly enhances the future potential of epistasis analysis for systems biology, and can complement genome-wide association studies as a means of formulating hypotheses for targeted experimental work.

DCA for genome-wide epistasis analysis: the statistical genetics perspective.

Direct coupling analysis (DCA) is a now widely used method to leverage statistical information from many similar biological systems to draw meaningful conclusions on each system separately. DCA has been applied with great success to sequences of homologous proteins, and also more recently to whole-genome population-wide sequencing data. We here argue that the use of DCA on the genome scale is contingent on fundamental issues of population genetics. DCA can be expected to yield meaningful results when a population is in the quasi-linkage equilibrium (QLE) phase studied by Kimura and others, but not, for instance, in a phase of clonal competition. We discuss how the exponential (Potts model) distributions emerge in QLE, and compare couplings to correlations obtained in a study of about 3000 genomes of the human pathogen Streptococcus pneumoniae.

Theory of Nonequilibrium Local Search on Random Satisfaction Problems.

We study local search algorithms to solve instances of the random k-satisfiability problem, equivalent to finding (if they exist) zero-energy ground states of statistical models with disorder on random hypergraphs. It is well known that the best such algorithms are akin to nonequilibrium processes in a high-dimensional space. In particular, algorithms known as focused, and which do not obey detailed balance, outperform simulated annealing and related methods in the task of finding the solution to a complex satisfiability problem, that is to find (exactly or approximately) the minimum in a complex energy landscape. A physical question of interest is if the dynamics of these processes can be well predicted by the well-developed theory of equilibrium Gibbs states. While it has been known empirically for some time that this is not the case, an alternative systematic theory that does so has been lacking. In this Letter we introduce such a theory based on the recently developed technique of cavity master equations and test it on the paradigmatic random 3-satisfiability problem. Our theory predicts the qualitative form of the phase boundary between the satisfiable (SAT) and unsatisfiable (UNSAT) region of the phase diagram where the numerics of a focused Metropolis search and cavity master equation cannot be distinguished.

Landscape of RNA polyadenylation in E. coli.

Polyadenylation is thought to be involved in the degradation and quality control of bacterial RNAs but relatively few examples have been investigated. We used a combination of 5΄-tagRACE and RNA-seq to analyze the total RNA content from a wild-type strain and from a poly(A)polymerase deleted mutant. A total of 178 transcripts were either up- or down-regulated in the mutant when compared to the wild-type strain. Poly(A)polymerase up-regulates the expression of all genes related to the FliA regulon and several previously unknown transcripts, including numerous transporters. Notable down-regulation of genes in the expression of antigen 43 and components of the type 1 fimbriae was detected. The major consequence of the absence of poly(A)polymerase was the accumulation of numerous sRNAs, antisense transcripts, REP sequences and RNA fragments resulting from the processing of entire transcripts. A new algorithm to analyze the position and composition of post-transcriptional modifications based on the sequence of unencoded 3΄-ends, was developed to identify polyadenylated molecules. Overall our results shed new light on the broad spectrum of action of polyadenylation on gene expression and demonstrate the importance of poly(A) dependent degradation to remove structured RNA fragments.

Dynamics inside the cancer cell attractor reveal cell heterogeneity, limits of stability, and escape.

The observed intercellular heterogeneity within a clonal cell population can be mapped as dynamical states clustered around an attractor point in gene expression space, owing to a balance between homeostatic forces and stochastic fluctuations. These dynamics have led to the cancer cell attractor conceptual model, with implications for both carcinogenesis and new therapeutic concepts. Immortalized and malignant EBV-carrying B-cell lines were used to explore this model and characterize the detailed structure of cell attractors. Any subpopulation selected from a population of cells repopulated the whole original basin of attraction within days to weeks. Cells at the basin edges were unstable and prone to apoptosis. Cells continuously changed states within their own attractor, thus driving the repopulation, as shown by fluorescent dye tracing. Perturbations of key regulatory genes induced a jump to a nearby attractor. Using the Fokker-Planck equation, this cell population behavior could be described as two virtual, opposing influences on the cells: one attracting toward the center and the other promoting diffusion in state space (noise). Transcriptome analysis suggests that these forces result from high-dimensional dynamics of the gene regulatory network. We propose that they can be generalized to all cancer cell populations and represent intrinsic behaviors of tumors, offering a previously unidentified characteristic for studying cancer.

Whole-genome mapping of 5' RNA ends in bacteria by tagged sequencing: a comprehensive view in Enterococcus faecalis.

Enterococcus faecalis is the third cause of nosocomial infections. To obtain the first snapshot of transcriptional organizations in this bacterium, we used a modified RNA-seq approach enabling to discriminate primary from processed 5' RNA ends. We also validated our approach by confirming known features in Escherichia coli. We mapped 559 transcription start sites (TSSs) and 352 processing sites (PSSs) in E. faecalis. A blind motif search retrieved canonical features of SigA- and SigN-dependent promoters preceding transcription start sites mapped. We discovered 85 novel putative regulatory RNAs, small- and antisense RNAs, and 72 transcriptional antisense organizations. Presented data constitute a significant insight into bacterial RNA landscapes and a step toward the inference of regulatory processes at transcriptional and post-transcriptional levels in a comprehensive manner.

The bulk and the tail of minimal absent words in genome sequences.

Minimal absent words (MAW) of a genomic sequence are subsequences that are absent themselves but the subwords of which are all present in the sequence. The characteristic distribution of genomic MAWs as a function of their length has been observed to be qualitatively similar for all living organisms, the bulk being rather short, and only relatively few being long. It has been an open issue whether the reason behind this phenomenon is statistical or reflects a biological mechanism, and what biological information is contained in absent words. In this work we demonstrate that the bulk can be described by a probabilistic model of sampling words from random sequences, while the tail of long MAWs is of biological origin. We introduce the concept of a core of a MAW, which are sequences present in the genome and closest to a given MAW. We show that in E. faecalis, E. coli and yeast the cores of the longest MAWs, which exist in two or more copies, are located in highly conserved regions the most prominent example being ribosomal RNAs. We also show that while the distribution of the cores of long MAWs is roughly uniform over these genomes on a coarse-grained level, on a more detailed level it is strongly enhanced in 3' untranslated regions (UTRs) and, to a lesser extent, also in 5' UTRs. This indicates that MAWs and associated MAW cores correspond to fine-tuned evolutionary relationships, and suggest that they can be more widely used as markers for genomic complexity.

SEK: sparsity exploiting k-mer-based estimation of bacterial community composition.

MOTIVATION: Estimation of bacterial community composition from a high-throughput sequenced sample is an important task in metagenomics applications. As the sample sequence data typically harbors reads of variable lengths and different levels of biological and technical noise, accurate statistical analysis of such data is challenging. Currently popular estimation methods are typically time-consuming in a desktop computing environment. RESULTS: Using sparsity enforcing methods from the general sparse signal processing field (such as compressed sensing), we derive a solution to the community composition estimation problem by a simultaneous assignment of all sample reads to a pre-processed reference database. A general statistical model based on kernel density estimation techniques is introduced for the assignment task, and the model solution is obtained using convex optimization tools. Further, we design a greedy algorithm solution for a fast solution. Our approach offers a reasonably fast community composition estimation method, which is shown to be more robust to input data variation than a recently introduced related method. AVAILABILITY AND IMPLEMENTATION: A platform-independent Matlab implementation of the method is freely available at http://www.ee.kth.se/ctsoftware; source code that does not require access to Matlab is currently being tested and will be made available later through the above Web site.

Improving contact prediction along three dimensions.

Correlation patterns in multiple sequence alignments of homologous proteins can be exploited to infer information on the three-dimensional structure of their members. The typical pipeline to address this task, which we in this paper refer to as the three dimensions of contact prediction, is to (i) filter and align the raw sequence data representing the evolutionarily related proteins; (ii) choose a predictive model to describe a sequence alignment; (iii) infer the model parameters and interpret them in terms of structural properties, such as an accurate contact map. We show here that all three dimensions are important for overall prediction success. In particular, we show that it is possible to improve significantly along the second dimension by going beyond the pair-wise Potts models from statistical physics, which have hitherto been the focus of the field. These (simple) extensions are motivated by multiple sequence alignments often containing long stretches of gaps which, as a data feature, would be rather untypical for independent samples drawn from a Potts model. Using a large test set of proteins we show that the combined improvements along the three dimensions are as large as any reported to date.

Detection and quantitative estimation of spurious double stranded DNA formation during reverse transcription in bacteria using tagRNA-seq.

Standard RNA-seq has a well know tendency to generate "ghost" antisense reads due to formation of spurious second strand cDNA in the sequencing process. We recently reported on a novel variant of RNA-seq coined "tagRNA-seq" introduced for the purpose of distinguishing primary from processed transcripts in bacteria. Incidentally, the additional information provided by the tags is also very suitable for detection of true anti-sense RNA transcripts and quantification of spurious antisense signals in a sample. We briefly explain how to perform such a detection and illustrate on previously published datasets.

Maximum likelihood reconstruction for Ising models with asynchronous updates.

We describe how the couplings in an asynchronous kinetic Ising model can be inferred. We consider two cases: one in which we know both the spin history and the update times and one in which we know only the spin history. For the first case, we show that one can average over all possible choices of update times to obtain a learning rule that depends only on spin correlations and can also be derived from the equations of motion for the correlations. For the second case, the same rule can be derived within a further decoupling approximation. We study all methods numerically for fully asymmetric Sherrington-Kirkpatrick models, varying the data length, system size, temperature, and external field. Good convergence is observed in accordance with the theoretical expectations.

A simple and efficient method to search for selected primary transcripts: non-coding and antisense RNAs in the human pathogen Enterococcus faecalis.

Enterococcus faecalis is a commensal bacterium and a major opportunistic human pathogen. In this study, we combined in silico predictions with a novel 5'RACE-derivative method coined '5'tagRACE', to perform the first search for non-coding RNAs (ncRNAs) encoded on the E. faecalis chromosome. We used the 5'tagRACE to simultaneously probe and characterize primary transcripts, and demonstrate here the simplicity, the reliability and the sensitivity of the method. The 5'tagRACE is complementary to tiling arrays or RNA-sequencing methods, and is also directly applicable to deep RNA sequencing and should significantly improve functional studies of bacterial RNA landscapes. From 45 selected loci of the E. faecalis chromosome, we discovered and mapped 29 novel ncRNAs, 10 putative novel mRNAs and 16 antisense transcriptional organizations. We describe in more detail the oxygen-dependent expression of one ncRNA located in an E. faecalis pathogenicity island, the existence of an ncRNA that is antisense to the ncRNA modulator of the RNA polymerase, SsrS and provide evidences for the functional interplay between two distinct toxin-antitoxin modules.

Inverse Ising inference using all the data.

We show that a method based on logistic regression, using all the data, solves the inverse Ising problem far better than mean-field calculations relying only on sample pairwise correlation functions, while still computationally feasible for hundreds of nodes. The largest improvement in reconstruction occurs for strong interactions. Using two examples, a diluted Sherrington-Kirkpatrick model and a two-dimensional lattice, we also show that interaction topologies can be recovered from few samples with good accuracy and that the use of l(1) regularization is beneficial in this process, pushing inference abilities further into low-temperature regimes.

Anomalous thermodynamics at the microscale.

Particle motion at the microscale is an incessant tug-of-war between thermal fluctuations and applied forces on one side and the strong resistance exerted by fluid viscosity on the other. Friction is so strong that completely neglecting inertia--the overdamped approximation--gives an excellent effective description of the actual particle mechanics. In sharp contrast to this result, here we show that the overdamped approximation dramatically fails when thermodynamic quantities such as the entropy production in the environment are considered, in the presence of temperature gradients. In the limit of vanishingly small, yet finite, inertia, we find that the entropy production is dominated by a contribution that is anomalous, i.e., has no counterpart in the overdamped approximation. This phenomenon, which we call an entropic anomaly, is due to a symmetry breaking that occurs when moving to the small, finite inertia limit. Anomalous entropy production is traced back to futile phase-space cyclic trajectories displaying a fast downgradient sweep followed by a slow upgradient return to the original position.

Temporal epistasis inference from more than 3 500 000 SARS-CoV-2 genomic sequences.

We use direct coupling analysis (DCA) to determine epistatic interactions between loci of variability of the SARS-CoV-2 virus, segmenting genomes by month of sampling. We use full-length, high-quality genomes from the GISAID repository up to October 2021 for a total of over 3 500 000 genomes. We find that DCA terms are more stable over time than correlations but nevertheless change over time as mutations disappear from the global population or reach fixation. Correlations are enriched for phylogenetic effects, and in particularly statistical dependencies at short genomic distances, while DCA brings out links at longer genomic distance. We discuss the validity of a DCA analysis under these conditions in terms of a transient auasilinkage equilibrium state. We identify putative epistatic interaction mutations involving loci in spike.

Optimal protocols and optimal transport in stochastic thermodynamics.

Thermodynamics of small systems has become an important field of statistical physics. Such systems are driven out of equilibrium by a control, and the question is naturally posed how such a control can be optimized. We show that optimization problems in small system thermodynamics are solved by (deterministic) optimal transport, for which very efficient numerical methods have been developed, and of which there are applications in cosmology, fluid mechanics, logistics, and many other fields. We show, in particular, that minimizing expected heat released or work done during a nonequilibrium transition in finite time is solved by the Burgers equation and mass transport by the Burgers velocity field. Our contribution hence considerably extends the range of solvable optimization problems in small system thermodynamics.