Transcription factors (TFs) are essential players in the regulation of gene expression and thus have been the subject of interest in the context of diseases and the manipulation of specific cell functions or pathways. The complex interplay of TFs and genes can be modeled as a network, in which each edge represents a regulatory influence. A challenge in such networks is the identification of a set of TFs with the maximum regulatory influence, which is of particular importance in the context of perturbation studies. We model the task of finding a set of TFs to maximize the influence on a set of genes as a probing problem in a bipartite graph. We test different adaptive and non-adaptive algorithms on simulated graphs and discuss their properties and adaptivity gap. Then, we apply the algorithms on real-life data to find TFs regulating genes involved in T-cell mediated immunity and lymphoid leukemia. The approach has minimal data requirements and can be readily applied to all other types of bipartite networks.
Cancer phylogenies are key to understanding tumor evolution. However, due to the uncertainty in phylogenetic estimation, one typically infers many, equally-plausible phylogenies from bulk DNA sequencing data of tumors, hindering downstream analysis that relies on correct phylogenies. To resolve this challenge, we introduce Sapling, a method to solve two variants of the Backbone Tree Inference from Reads problem, which seeks a small set of backbone trees on a subset of mutations that collectively summarize the space of plausible cancer phylogenies. We prove that the problems are NP-hard. On simulated and real data, we demonstrate that Sapling is capable of inferring high-quality backbone trees that adequately summarize the space of plausible cancer phylogenies. In addition, we demonstrate that Sapling is able to infer full-size trees with higher likelihoods than state-of-the-art methods.
BackgroundAtom-to-atom maps play an important role in many applications. However, they are often difficult to obtain. The KEGG reaction database does not provide atom-to-atom maps for its reactions and instead offers a description of local changes for pairs of reactant and product molecules in terms of so-called RCLASSes. Developed for classification purposes, RCLASS data are difficult to use for purposes such as the construction of atom-to-atom maps or reaction rules. DPO graph transformation rules, on the other hand, work as a convenient and efficient representation, particularly for these applications. The RCLASS data can be understood as collections of local graph patterns in the reactants and products of a reaction, together with partial correspondences of atoms. The problem of converting RCLASS data into DPO rules, therefore, is a special case of the graph reconstruction problem, which consists of inferring a graph from a collection of subgraphs.ResultsWe developed laveau, a tool that computes explicit DPO rules from KEGG reactions and RCLASS data. The algorithm proceeds stepwise, starting with a translation of individual RDM codes, specifically developed by the KEGG database, into equivalent RDM pattern graphs. Multiple RDM pattern graphs for the same RCLASS are then combined based on their embeddings into the reactant and product molecules, observing certain consistency conditions. In the final step, these combined pairwise patterns are merged into a pair of subgraphs of reactants and products, respectively. If RCLASSes connecting all pairs of reactant and product molecules are available, the complete reaction center(s) is/are contained in the union of these subgraphs. The atom-to-atom map inherited from the RDM codes then defines a DPO transformation rule. Application of these rules to the reactants then yields complete atom-to-atom maps (AAMs). Starting from 3195 RCLASSes, laveau generates a total of 1232 DPO rules and 1594 AAMs.ConclusionsThe laveau software makes it possible to extract local atom-to-atom maps from the RCLASSes of the KEGG database, covering a large set of enzyme-catalyzed reactions. The results are made available in the form of DPO rules for use in atom-level models of metabolic networks, filling a crucial gap in the available data.
Background: In the Nearest-Neighbor Thermodynamic Model, a standard approach for RNA secondary structure prediction, the energy of the multiloops is modeled using a linear entropic penalty governed by three branching parameters. Although these parameters are typically fixed, recent work has shown that reparametrizing the multiloop score and considering alternative branching conformations can lead to significantly better structure predictions. However, prior approaches for exploring the alternative branching structures were computationally inefficient for long sequences. Results: We present a novel algorithm that partitions the parameter space, identifying all distinct branching structures (optimal under different branching parameters) for a given RNA sequence using the fewest possible minimum free energy computations. Our method efficiently computes the full parameter-space partition and the associated optimal structures, enabling a comprehensive evaluation of the structural landscape across parameter choices. We apply this algorithm to the Archive II benchmarking dataset, assessing the maximum attainable prediction accuracy for each sequence under the reparameterized multiloop model. We find that the potential for improvement over default predictions is substantial in many cases, and that the optimal prediction accuracy is highly sensitive to auxiliary modeling decisions, such as the treatment of lonely base pairs and dangling ends. Conclusion: Our results support the hypothesis that the conventional choice of multiloop parameters may limit prediction accuracy and that exploring alternative parameterizations is both tractable and worthwhile. The efficient partitioning algorithm we introduce makes this exploration feasible for longer sequences and larger datasets. Furthermore, we identify several open challenges in identifying the optimal structure.
Chemical reaction databases typically report the molecular structures of reactant and product compounds, as well as their stoichiometry, but lack information, in particular, on the correspondence of reactant and product atoms. These atom-to-atom maps (AAM), however, are crucial for applications including chemical synthesis planning in organic chemistry and the analysis of isotope labeling experiments in modern metabolomics. AAMs therefore need to be reconstructed computationally. This situation is aggravated, furthermore, by the fact that chemically correct AAMs are, fundamentally, determined by quantum-mechanical phenomena and thus cannot be reliably computed by solving graph-theoretical optimization problems defined by the reactant and product structures. A viable solution for this problem is to shift the focus into first identifying a partial AAM containing the reaction center, i.e., covering the atoms incident with all bonds that change during a reaction. This then leads to the problem of extending the partial map to the full reaction. The AAM of a reaction is faithfully represented by the Imaginary Transition State (ITS) graph, providing a convenient graph-theoretic framework to address the questions of when and how a partial AAM can be extended. We show that an unique extension exists whenever, and only if, these partial AAMs cover the reaction center. Moreover, uniqueness results are generalized to partial AAMs in situations where hydrogen atoms are not represented explicitly. In this case their extension can be computed by solving a constrained graph-isomorphism search between specific subgraphs of ITS graphs. We close by benchmarking different tools for this task.
Abstract Background Several recent cancer phylogeny inference methods have used the k -Dollo evolutionary model for single-nucleotide variants, which requires a phylogeny T on binary sequencing data matrix B such that each variant is gained once and lost at most k times. The 1-Dollo variant has been studied extensively but its hardness remains open. Results We prove that the 1-Dollo Linear Phylogeny (1DLP) problem, where we additionally require the resulting 1-Dollo phylogeny T to be linear, is equivalent to verifying whether matrix B has the Consecutive Ones Property, which can be determined in polynomial time. We also show that some practical extensions of 1DLP, such as the minimization of false negatives, are NP-hard. We then show how to recursively decompose any 1-Dollo phylogeny T , not necessarily linear, into several 1-Dollo linear phylogenies and extend this characterization to all matrices B that admit 1-Dollo phylogenies. We use this characterization to develop Dolphyin, a new exponential-time algorithm for inferring 1-Dollo phylogenies. Dolphyin is runtime-competitive with integer linear programming-based algorithm SPhyR (El-Kebir 2018) on simulated datasets and infers 1-Dollo phylogenies with false negative sequencing error rates at or below simulated ground truth rates. We apply Dolphyin to acute myeloid leukemia datasets and find that the majority of the cancers can be explained by 1-Dollo phylogenies with error rates in line with the used sequencing technology. Conclusion Our work develops a novel, combinatorial algorithm for practical inference of 1-Dollo phylogenies.
Mutational processes shape cancer genomes, leaving characteristic marks that are termed signatures. The level of activity of each such process, or its signature exposure, provides important information on the disease, improving patient stratification and the prediction of drug response. Thus, there is growing interest in developing refitting methods that accurately decipher those exposures. Previous work in this domain was unsupervised in nature, employing algebraic decomposition and probabilistic inference methods. We present SuRe, a supervised approach to signature refitting that demonstrates superiority over current methods. SuRe leverages a neural network model to capture correlations between signature exposures in real data. We show that SuRe outperforms previous methods on sparse mutation data from both tumor-type-specific and pan-cancer data sets, with an increasing performance advantage as the data become sparser. We further demonstrate the model’s utility in clinical settings by predicting homologous recombination deficiency in breast cancer from sparse data. Furthermore, SuRe outperforms standard methods in the unsupervised stratification of over 13,000 patients from large-scale panel sequencing cohorts, highlighting its potential for analyzing targeted sequencing data.
Error correction of long reads is an important initial step in genome assembly workflows. For organisms with ploidy greater than one, it is important to preserve haplotype-specific variation during read correction. This challenge has driven the development of several haplotype-aware correction methods. However, existing methods are based on either ad-hoc heuristics or deep learning approaches. In this paper, we introduce a rigorous formulation for this problem. Our approach builds on the minimum error correction framework used in reference-based haplotype phasing. We prove that the proposed formulation for error correction of reads in de novo context, i.e., without using a reference genome, is NP-hard. To make our exact algorithm scale to large datasets, we introduce practical heuristics. Experiments using PacBio HiFi sequencing datasets from human and plant genomes show that our approach achieves accuracy comparable to state-of-the-art methods. Implementation:https://github.com/at-cg/HALE.
Motivation: Spliced alignment refers to the alignment of messenger RNA (mRNA) or protein sequences to eukaryotic genomes. It plays a critical role in gene annotation and the study of gene functions. Accurate spliced alignment demands sophisticated modeling of splice sites, but current aligners use simple models, which may affect their accuracy given dissimilar sequences. Results: We implemented minisplice to learn splice signals with a one-dimensional convolutional neural network (1D-CNN) and trained a model with 7,026 parameters for vertebrate and insect genomes. It captures conserved splice signals across phyla and reveals GC-rich introns specific to mammals and birds. We used this model to estimate the empirical splicing probability for every GT and AG in genomes, and modified minimap2 and miniprot to leverage pre-computed splicing probability during alignment. Evaluation on human long-read RNA-seq data and cross-species protein datasets showed our method greatly improves the junction accuracy especially for noisy long RNA-seq reads and proteins of distant homology. Availability and implementation: https://github.com/lh3/minisplice
K-mer-based analysis of genomic data is ubiquitous, but the presence of repetitive k-mers continues to pose problems for the accuracy of many methods. For example, the Mash tool (Ondov et al 2016) can accurately estimate the substitution rate between two low-repetitive sequences from their k-mer sketches; however, it is inaccurate on repetitive sequences such as the centromere of a human chromosome. Follow-up work by Blanca et al. (2021) has attempted to model how mutations affect k-mer sets based on strong assumptions that the sequence is non-repetitive and that mutations do not create spurious k-mer matches. However, the theoretical foundations for extending an estimator like Mash to work in the presence of repeat sequences have been lacking. In this work, we relax the non-repetitive assumption and propose a novel estimator for the mutation rate. We derive theoretical bounds on our estimator's bias. Our experiments show that it remains accurate for repetitive genomic sequences, such as the alpha satellite higher order repeats in centromeres. We demonstrate our estimator's robustness across diverse datasets and various ranges of the substitution rate and k-mer size. Finally, we show how sketching can be used to avoid dealing with large k-mer sets while retaining accuracy. Our software is available at https://github.com/medvedevgroup/Repeat-Aware_Substitution_Rate_Estimator.
The 𝒞 -Orientation problem asks whether it is possible to orient an undirected graph to a directed phylogenetic network of a desired network class 𝒞 . This problem arises, for example, when visualising evolutionary data, as popular methods such as Neighbor-Net are distance-based and inevitably produce undirected graphs. The complexity of 𝒞 -Orientation remains open for many classes 𝒞 , including binary tree-child networks, and practical methods are still lacking. In this paper, we propose (1) an exact FPT algorithm for 𝒞 -Orientation, applicable to any class 𝒞 admitting a tractable membership test, and parameterised by the reticulation number and the maximum size of minimal basic cycles, and (2) a very fast heuristic for Tree-Child Orientation. While the state-of-the-art for 𝒞 -Orientation is a simple exponential time algorithm whose computational bottleneck lies in searching for appropriate reticulation vertex placements, our methods significantly reduce this search space. Experiments show that, although our FPT algorithm is still exponential, it significantly outperforms the existing method. The heuristic runs even faster but with increasing false negatives as the reticulation number grows. Given this trade-off, we also discuss theoretical directions for improvement and biological applicability of the heuristic approach.
Adding rank support to strings over a fixed-sized alphabet has numerous applications. Prominent among those is the (bidirectional) FM-Index which is commonly utilized to index and analyze genomic data. At its core lies the rank operation on the Burrows-Wheeler-Transform (BWT) which, given a position in the BWT and a character, answers how often the specified character appears from the start to that position. Implementing those rank queries is usually based on bit vectors with rank support. In this work, we discuss three implementation improvements. First, a novel approach named paired-blocks that reduces the space overhead of the support structure by half to a total of only 1.6% . Second, a method for masking bits for the population count (also known as popcount) which greatly improves the runtime of 512-bit wide blocks in conjunction with AVX512 SIMD extensions. Third, a revised method for EPR-dictionaries (Pockrandt et al. in International conference on research in computational molecular biology. Springer, New York, 2017) called flattened bit vectors (fBV) with less space consumption and faster rank operations on strings, which is competitive in size and depending on the parameters between 2× and 9× faster than Wavelet Trees (Gog et al. in 13th International Symposium on Experimental Algorithms. Springer, New York, 2014).
A biological study can produce a limited number of marker genes, not large enough to be used in gene set enrichment analysis. Here we suggest VOL-Gene, a graph-based algorithm that partitions all genes into non-overlapping classes of functionally related genes, thus assigning a single function to each gene. To this end, many functional signatures are combined into a single weighted graph, which is partitioned into cliques. For a poorly annotated marker gene, our approach fetches a number of genes that belong to the same class, some of which can be well annotated and are likely to take part in the same biological process.
Finding shortest unique substrings (SUS) is a fundamental problem in string processing with applications in bioinformatics. In this paper, we present an algorithm for solving a variant of the SUS problem, the left-bounded shortest unique substrings (LSUS). This variant is particularly important in applications such as PCR primer design. Our algorithm runs in O(n) time using 2n memory words plus n bytes for an input string of length n. Experimental results with real and artificial datasets show that our algorithm is the fastest alternative in practice, being two times faster (on the average) than related works, while using a similar peak memory footprint.
The increasing number and volume of genomic and metagenomic data necessitates scalable and robust computational models for precise analysis. Sketching techniques utilizing k -mers from a biological sample have proven to be useful for large-scale analyses. In recent years, FracMinHash has emerged as a popular sketching technique and has been used in several useful applications. Recent studies on FracMinHash proved unbiased estimators for the containment and Jaccard indices. However, theoretical investigations for other metrics are still lacking. In this paper, we present a theoretical framework for estimating similarity/distance metrics by using FracMinHash sketches, when the metric is expressible in a certain form. We establish conditions under which such an estimation is sound and recommend a minimum scale factor s for accurate results. Experimental evidence supports our theoretical findings. We also present frac-kmc, a fast and efficient FracMinHash sketch generator program. frac-kmc is the fastest known FracMinHash sketch generator, delivering accurate and precise results for cosine similarity estimation on real data. frac-kmc is also the first parallel tool for this task, allowing for speeding up sketch generation using multiple CPU cores – an option lacking in existing serialized tools. We show that by computing FracMinHash sketches using frac-kmc, we can estimate pairwise similarity speedily and accurately on real data. frac-kmc is freely available here: https://github.com/KoslickiLab/frac-kmc/
Analyzing and comparing sequences of symbols is among the most fundamental problems in computer science, possibly even more so in bioinformatics. Maximal Common Subsequences (MCSs), i.e., inclusion-maximal sequences of non-contiguous symbols common to two or more strings, have only recently received attention in this area, despite being a basic notion and a natural generalization of more common tools like Longest Common Substrings/Subsequences. In this paper we simplify and engineer recent advancements in MCSs into a practical tool called M C D A G , the first publicly available tool that can index MCSs of real genomic data, and show that its definition can be generalized to multiple strings. We demonstrate that our tool can index pairs of sequences exceeding 10,000 base pairs within minutes, utilizing only 4-7% more than the minimum required nodes. For three or more sequences, we observe experimentally that the minimum index may exhibit a significant increase in the number of nodes.
Orthology inference lies at the foundation of comparative genomics research. The correct identification of loci which descended from a common ancestral sequence is not only complicated by sequence divergence but also duplication and other genome rearrangements. The conservation of gene order, i.e. synteny, is used in conjunction with sequence similarity as an additional factor for orthology determination. Current approaches, however, rely on genome annotations and are therefore limited. Here we present an annotation-free approach and compare it to synteny analysis with annotations. We find that our approach works better in closely related genomes whereas there is a better performance with annotations for more distantly related genomes. Overall, the presented algorithm offers a useful alternative to annotation-based methods and can outperform them in many cases.
Orthologous genes, which arise through speciation, play a key role in comparative genomics and functional inference. In particular, graph-based methods allow for the inference of orthology estimates without prior knowledge of the underlying gene or species trees. This results in orthology graphs, where each vertex represents a gene, and an edge exists between two vertices if the corresponding genes are estimated to be orthologs. Orthology graphs inferred under a tree-like evolutionary model must be cographs. However, real-world data often deviate from this property, either due to noise in the data, errors in inference methods or, simply, because evolution follows a network-like rather than a tree-like process. The latter, in particular, raises the question of whether and how orthology graphs can be derived from or, equivalently, are explained by phylogenetic networks. In this work, we study the constraints imposed on orthology graphs when the underlying evolutionary history follows a phylogenetic network instead of a tree. We show that any orthology graph can be represented by a sufficiently complex level-k network. However, such networks lack biologically meaningful constraints. In contrast, level-1 networks provide a simpler explanation, and we establish characterizations for level-1 explainable orthology graphs, i.e., those derived from level-1 evolutionary histories. To this end, we employ modular decomposition, a classical technique for studying graph structures. Specifically, an arbitrary graph is level-1 explainable if and only if each primitive subgraph is a near-cograph (a graph in which the removal of a single vertex results in a cograph). Additionally, we present a linear-time algorithm to recognize level-1 explainable orthology graphs and to construct a level-1 network that explains them, if such a network exists. Finally, we demonstrate the close relationship of level-1 explainable orthology graphs to the substitution operation, weakly chordal and perfect graphs, as well as graphs with twin-width at most 2.
Due to the increasing availability of high-quality genome sequences, pan-genomes are gradually replacing single consensus reference genomes in many bioinformatics pipelines to better capture genetic diversity. Traditional bioinformatics tools using the FM-index face memory limitations with such large genome collections. Recent advancements in run-length compressed indices like Gagie et al.’s r-index and Nishimoto and Tabei’s move structure, alleviate memory constraints but focus primarily on backward search for MEM-finding. Arakawa et al.’s br-index initiates complete approximate pattern matching using bidirectional search in run-length compressed space, but with significant computational overhead due to complex memory access patterns. We introduce b-move, a novel bidirectional extension of the move structure, enabling fast, cache-efficient, lossless approximate pattern matching in run-length compressed space. It achieves bidirectional character extensions up to 7 times faster than the br-index, closing the performance gap with FM-index-based alternatives. For locating occurrences, b-move performs ϕ and ϕ ^-1 operations up to 7 times faster than the br-index. At the same time, it maintains the favorable memory characteristics of the br-index, for example, all available complete E. coli genomes on NCBI’s RefSeq collection can be compiled into a b-move index that fits into the RAM of a typical laptop. b-move proves practical and scalable for pan-genome indexing and querying. We provide a C++ implementation of b-move, supporting efficient lossless approximate pattern matching including locate functionality, available at https://github.com/biointec/b-move under the AGPL-3.0 license.
Bayesian phylogenetics typically estimates a posterior distribution, or aspects thereof, using Markov chain Monte Carlo methods. These methods integrate over tree space by applying local rearrangements to move a tree through its space as a random walk. Previous work explored the possibility of replacing this random walk with a systematic search, but was quickly overwhelmed by the large number of probable trees in the posterior distribution. In this paper we develop methods to sidestep this problem using a recently introduced structure called the subsplit directed acyclic graph (sDAG). This structure can represent many trees at once, and local rearrangements of trees translate to methods of enlarging the sDAG. Here we propose two methods of introducing, ranking, and selecting local rearrangements on sDAGs to produce a collection of trees with high posterior density. One of these methods successfully recovers the set of high posterior density trees across a range of data sets. However, we find that a simpler strategy of aggregating trees into an sDAG in fact is computationally faster and returns a higher fraction of probable trees.
