Journal of Computational Biology

The study of evolution between species (phylogenetics) and the study of evolution within a species (population genetics) are highly related, as the same biological mechanisms are fundamental to both fields. Although both have been studied for a long time, their joint study in a unified setting has been prevented by the different time scales they consider and the different data types they employ. A similar discrepancy holds for their whole-genome specializations, comparative genomics and pangenomics. Two active areas in these fields are genome rearrangement studies and graphical pangenomics, respectively. Since the emergence of graphical pangenomics, these have existed as separate fields, despite observations that central data structures representing genomic variants in both fields are highly similar. While there exists a wealth of theoretical results for various rearrangement models in comparative genomics, the application to pangenomic data is hampered by the limitations of rearrangement problem formulations. On the practical side, pangenomes typically contain too many individual genomes for classical problems, such as the often NP-hard parsimony problems, to be solved, or for all-vs-all comparisons using rearrangement distances to be performed. On the theoretical side, some assumptions in the formulation of rearrangement problems, such as the assumption of an underlying tree, are inadequate for many pangenomes. In this work, we propose the Complete Ancestral Reconstruction for Pangenomes (CARP) problem, which overcomes these limitations while retaining intuitive relationships to both classical rearrangement problems and pangenome graphs.

Distance-guided tree construction with unknown tree topology and branch lengths has been a long studied problem. In contrast, distance-guided branch lengths assignment with fixed tree topology has not yet been systematically investigated, despite having significant applications. In this paper, we provide a formal mathematical formulation of this problem and propose two representative methods for solving this problem, each with its own strength. We evaluate the performance of these two methods under various settings using simulated data, providing guidance for the choice of methods in respective cases. We demonstrate a practical application of this operation through an extension we termed FunUniFrac, which quantifies the differences in functional units between metagenomic samples over a functional tree with assigned branch lengths, allowing clustering of metagenomic samples by functional similarity instead of taxonomic similarity in traditional methods, thus expanding the realm of comparative studies in metagenomics.

MotivationIntra-sample heterogeneity describes the phenomenon where a genomic sample contains a diverse set of genomic sequences. In practice, the true string sets in a sample are often unknown due to limitations in sequencing technology. In order to compare heterogeneous samples, genome graphs can be used to represent such sets of strings. However, a genome graph is generally able to represent a string set universe that contains multiple sets of strings in addition to the true string set. This difference between genome graphs and string sets is not well characterized. As a result, a distance metric between genome graphs may not match the distance between true string sets. ResultsWe extend a genome graph distance metric, Graph Traversal Edit Distance (GTED) proposed by Ebrahimpour Boroojeny et al., to FGTED to model the distance between heterogeneous string sets and show that GTED and FGTED always underestimate the Earth Movers Edit Distance (EMED) between string sets. We introduce the notion of string set universe diameter of a genome graph. Using the diameter, we are able to upper-bound the deviation of FGTED from EMED and to improve FGTED so that it reduces the average error in empirically estimating the similarity between true string sets. On simulated TCR sequences and Hepatitis B virus genomes, we show that the diameter-corrected FGTED reduces the average deviation of the estimated distance from the true string set distances by more than 250%. AvailabilityData and source code for reproducing the experiments are available at: https://github.com/Kingsford-Group/gtedemedtest/ Contactcarlk@cs.cmu.edu

The transition towards graph pangenomes is posing several new challenging questions, most notably how to extend the classical notion of read alignment from a sequence-to-sequence to a sequence-to-graph setting. Especially on variation graphs, where paths corresponding to individual genomes are labeled, notions of alignments that are strongly inspired by the classical ones are usually able to capture only variations that can be expressed by mismatches or gaps, such as SNPs or short insertions and deletions. On the other hand the recent investigation of pangenomes at bacterial scale (Colquhoun et al, 2021) shows that most tools are tailored for human pangenomes and are not suited to bacteria which exhibit, among other characteristics, a larger variability. Such variability leads to the need for incorporating a greater flexibility when computing an alignment. In this paper, we extend the usual notion of sequence-to-graph alignment by including recombinations among the variations that explicitly represented and evaluated in an alignment. From a computational modeling point of view, a recombination corresponds to identifying a new path of the variation graph which is a mosaic of two different paths, possibly joined by a new arc. We provide a dynamic programming algorithm for computing an optimal alignment that allows recombinations with an affine penalty. We have implemented our approach with the tool RecGraph and we have analyzed its accuracy over some over some bacterial pangenome graphs.

In this paper, we present diploid sequence graphs, graphs whose paths encode pairs of haplotypes. We describe an efficient algorithm for creating a diploid graph from a directed acyclic (haploid) sequence graph, such that the diploid graph represents all the possible pairings of haplotypes present in the sequence graph and their similarity relationships. Starting with the sequence graph, our method uses a graph decomposition approach based on an extension of the SPQR-tree to systematically identify structural patterns that reduce redundancy while preserving genetic variation. We develop a polynomial-time algorithm that parsimoniously enumerates all disjoint paths with shared endpoints in two-terminal directed acyclic graphs. In the future, we envisage that diploid graphs may enable more accurate modeling of recombination, phasing, and variation-aware alignment in diploid genomes.

Co-linear chaining is a widely used technique in sequence alignment tools that follow seed-filter-extend methodology. It is a mathematically rigorous approach to combine short exact matches. For colinear chaining between two sequences, efficient subquadratic-time chaining algorithms are well-known for linear, concave and convex gap cost functions [Eppstein et al. JACM92]. However, developing extensions of chaining algorithms for directed acyclic graphs (DAGs) has been challenging. Recently, a new sparse dynamic programming framework was introduced that exploits small path cover of pangenome reference DAGs, and enables efficient chaining [Makinen et al. TALG19, RECOMB18]. However, the underlying problem formulation did not consider gap cost which makes chaining less effective in practice. To address this, we develop novel problem formulations and optimal chaining algorithms that support a variety of gap cost functions. We demonstrate empirically the ability of our provably-good chaining implementation to align long reads more precisely in comparison to existing aligners. For mapping simulated long reads from human genome to a pangenome DAG of 95 human haplotypes, we achieve 98.7% precision while leaving < 2% reads unmapped. Implementationhttps://github.com/at-cg/minichain

Fixed-length k-mers have long been the standard in sequence analysis. However, they impose a uniform resolution across heterogeneous genomes, often resulting in significant redundancy and a loss of contextual sensitivity. To address these limitations, we introduce Minimum Unique Substrings (MUSs), which are variable-length sequence units that adapt to the local complexity of the genome. MUSs function as context-aware markers that naturally define repeat boundaries by extending only until uniqueness is achieved. We build upon the theoretical relationship between MUSs and maximal repeats, extending this framework to sequencing reads by establishing a read-consistent definition of uniqueness. We present a linear-time ([O] (n)) algorithm based on a generalized suffix tree and introduce the concept of "outposts." These outposts act as anchors for uniqueness, enabling precise localization of MUS boundaries within the sequencing data. Empirical studies of E. coli K-12 and human HiFi reads reveal distinct distributions in MUS lengths that reflect their respective genomic architectures. The compact bacterial genome produces a highly dense set of MUSs with a narrow length distribution (averaging 30.44 bp). In contrast, the repeat-rich human genome requires longer substrings to resolve uniqueness, resulting in an increased mean length (36.08 bp) and a broader distribution that delineates complex repetitive elements. The MUS framework achieves 100% unique coverage with an average length of 36.08 bp, surpassing the 69% coverage of k = 61. By reducing the total number of tokens by over 99%, it provides higher resolution and superior data compression compared to fixed-length k-mer sampling. These results demonstrate that MUSs provide a biologically meaningful, context-sensitive alternative to k-mers, with direct applications in genome assembly, repeat characterization, and comparative genomics.

A well-studied approximate version of the graph matching problem is directly relevant for the study of protein-protein interaction networks. Called by the computational biology community Global Network Alignment, the two networks to be matched are derived from the protein-protein interaction (PPI) networks from organisms of two different species. If these two species evolved recently from a common ancestor, we can view the two PPI networks as a single network that evolved over time. It is the two versions of this network that we want to align using approximate graph matching. The first spectral method for the PPI global alignment problem proposed by the biological community was the IsoRank method of Singh et al. This method for global biological network alignment is still used today. However, with the advent of many more experiments, the size of the networks available to match has grown considerably, making running IsoRank unfeasible on these networks without access to state of the art computational resources. In this paper, we develop a new IsoRank approximation, which exploits the mathematical properties of IsoRanks linear system to solve the problem in quadratic time with respect to the maximum size of the two PPI networks. We further propose a computationally cheaper refinement to this initial approximation so that the updated result is even closer to the original IsoRank formulation. In experiments on synthetic and real PPI networks, we find that the results of our approximate IsoRank are not only nearly as accurate as the original IsoRank results but are also much faster, which makes the global alignment of large-scale biological networks feasible and scalable.

Graph based non-linear reference structures such as variation graphs and colored de Bruijn graphs enable incorporation of full genomic diversity within a population. However, transitioning from a simple string-based reference to graphs requires addressing many computational challenges, one of which concerns accurately mapping sequencing read sets to graphs. Paired-end Illumina sequencing is a commonly used sequencing platform in genomics, where the paired-end distance constraints allow disambiguation of repeats. Many recent works have explored provably good index-based and alignment-based strategies for mapping individual reads to graphs. However, validating distance constraints efficiently over graphs is not trivial, and existing sequence to graph mappers rely on heuristics. We introduce a mathematical formulation of the problem, and provide a new algorithm to solve it exactly. We take advantage of the high sparsity of reference graphs, and use sparse matrix-matrix multiplications (SpGEMM) to build an index which can be queried efficiently by a mapping algorithm for validating the distance constraints. Effectiveness of the algorithm is demonstrated using real reference graphs, including a human MHC variation graph, and a pan-genome de-Bruijn graph built using genomes of 20 B. anthracis strains. While the one-time indexing time can vary from a few minutes to a few hours using our algorithm, answering a million distance queries takes less than a second.\n\n2012 ACM Subject ClassificationMathematics of computing [-&gt;] Paths and connectivity problems; Applied computing [-&gt;] Computational genomics

Metric multidimensional scaling is one of the classical methods for embedding data into low-dimensional Euclidean space. It creates the low-dimensional embedding by approximately preserving the pairwise distances between the input points. However, current state-of-the-art approaches only scale to a few thousand data points. For larger data sets such as those occurring in single-cell RNA sequencing experiments, the running time becomes prohibitively large and thus alternative methods such as PCA are widely used instead. Here, we propose a neural network based approach for solving the metric multidimensional scaling problem that is orders of magnitude faster than previous state-of-the-art approaches, and hence scales to data sets with up to a few million cells. At the same time, it provides a non-linear mapping between high- and low-dimensional space that can place previously unseen cells in the same embedding.

Multidomain proteins are mosaics of domains, protein modules that are associated with a specific structure or function and are found in diverse combinations. This modular organization facilitates the evolution of novel protein functions, but the principles that govern the relationship between the domain content of a protein and its function is poorly understood. In particular, do domains always contribute the same function, or does the functional contribution of a domain depend on the neighboring domains in the protein? To answer this question, we used vector embeddings, which account for local contextual signals, to model the protein domain content of multidomain proteins. We observe that multidomain architectures that are semantically similar share more functional attributes than multidomain architectures selected based on domain content similarity, alone, suggesting that context is important for understanding the relationship between domain content and protein function. Surprisingly, vector semantics also identified multidomain architecture pairs with significantly high functional similarity, despite having no domains in common at all, suggesting that vector semantics may be discovering domain "synonyms". Taken together, our results underscore the importance of contextual models for understanding the interplay between domain architecture evolution and functional innovation in multidomain proteins.

In many bioinformatics applications the task is to identify biologically significant locations in an individual genome. In our work, we are interested in finding high-density clusters of such biologically meaningful locations in a graph representation of a pangenome, which is a collection of related genomes. Different formulations of finding such clusters were previously studied for sequences. In this work, we study an extension of this problem for graphs, which we formalize as finding a set of vertex-disjoint paths with a maximum score in a weighted directed graph. We provide a linear-time algorithm for a special class of graphs corresponding to elastic-degenerate strings, one of pangenome representations. We also provide a fixed-parameter tractable algorithm for directed acyclic graphs with a special path decomposition of a limited width.

In computational biology, mapping a sequence s onto a sequence graph G poses a significant challenge. One possible approach to tackling this problem is to find a walk p in G that spells a sequence most similar to s. This challenge is formally known as the Graph Sequence Mapping Problem (GSMP). In this paper, we delve into an alternative problem formulation known as the De Bruijn Graph Sequence Mapping Problem (BSMP). Both problems have three variants: changes only in the sequence, changes in the graph, and changes in both the sequence and the graph. We concentrate on addressing the variant involving changes in the graph. In the literature, when this problem does not allow the De Bruijn graph to induce new arcs after changes, it becomes NP-complete, as proven by Gibney et. al [4]. However, we reformulate the problem by considering the characteristics of the arcs induced in the De Bruijn graph. This reformulation alters the problem definition, thereby enabling the application of a polynomial-time algorithm for its resolution. Approaching the problem with this arc-inducing characteristic is new, and the algorithm proposed in this work is new in the literature.

We present a method to provide a biologically meaningful representation of the space of protein sequences. While billions of protein sequences are available, organizing this vast amount of information into functional categories is daunting, time-consuming and incomplete. We present our unsupervised approach that combines Transformer protein language models, UMAP graphs, and spectral clustering to create meaningful clusters in the protein spaces. To demonstrate the meaningfulness of the clusters, we show that they preserve most of the signal present in a dataset of manually curated enzyme protein families.

Phylogenetic inference has become a crucial tool for interpreting cancer genomic data, but continuing advances in our understanding of somatic mutability in cancer, genomic technologies for profiling it, and the scale of data available have created a persistent need for new algorithms able to deal with these challenges. One particular need has been for new forms of consensus tree algorithms, which present special challenges in the cancer space for dealing with heterogeneous data, short evolutionary time scales, and rapid mutation by a wide variety of somatic mutability mechanisms. We develop a new consensus tree method for clonal phylogenetics, ConTreeDP, based on a formulation of the Maximum Directed Partition Support Consensus Tree (MDPSCT) problem. We demonstrate theoretically and empirically that our approach can efficiently and accurately compute clonal consensus trees from cancer genomic data. Availability: https://github.com/CMUSchwartzLab/ConTreeDP

A pangenome is a collection of taxonomically related genomes, often from the same species, serving as a representation of their genomic diversity. The study of pangenomes, or pangenomics, aims to quantify and compare this diversity, which has significant relevance in fields such as medicine and biology. Originally conceptualized as sets of genes, pangenomes are now commonly represented as pangenome graphs. These graphs consist of nodes representing genomic sequences and edges connecting consecutive sequences within a genome. Among possible pangenome graphs, a common option is the compacted de Bruijn graph. In our work, we focus on the colored compacted de Bruijn graph, where each node is associated with a set of colors that indicate the genomes traversing it. In response to the evolution of pangenome representation, we introduce a novel method for comparing pangenomes by their node counts, addressing two main challenges: the variability in node counts arising from graphs constructed with different numbers of genomes, and the large influence of rare genomic sequences. We propose an approach for interpolating and extrapolating node counts in colored compacted de Bruijn graphs, adjusting for the number of genomes. To tackle the influence of rare genomic sequences, we apply Hill numbers, a well-established diversity index previously utilized in ecology and metagenomics for similar purposes, to proportionally weight both rare and common nodes according to the frequency of genomes traversing them.

Single-cell sequencing (SCS) enables investigating tumor evolution at a single cell resolution. A common type of analysis to investigate evolutionary structure from an SCS experiment is to determine a phylogenetic tree structure from the data. This problem has been well-studied under the assumption that mutations only accumulate in the evolution of cancer and there is a simple characterization of when the data is compatible with a perfect phylogeny based on the absence of a special "conflict" submatrix. SCS data can be represented as a binary matrix, where the ij-th entry indicates whether cell i has mutation j. In practice, SCS data is noisy, so a natural question is what is the minimum number of entries to flip in the data matrix, in order that the matrix becomes "conflict-free" and thus compatible with a perfect phylogeny. Furthermore, the false positive rate is orders of magnitude smaller than the false negative rate, so that at most a few false positives occur with high probability. We consider a variation of the minimum-flip problem in which the number of false positives in the solution is a parameter. Often, there can be multiple optimal solutions, so a natural question is what relations are true among all optimal solutions for a small range of possible false positives values; we call such relations essential. In this work, we propose an efficient algorithm based on integer linear programming to determine all essential relations in the data. We test our software tool, EssentCell, on several data sets and discuss the results found.

MotivationUnderstanding the phylogenetic relationship among organisms is the key in contemporary evolutionary study and sequence analysis is the workhorse towards this goal. Conventional approaches to sequence analysis are based on sequence alignment, which is neither scalable to large-scale datasets due to computational inefficiency nor adaptive to next-generation sequencing (NGS) data. Alignment-free approaches are typically used as computationally effective alternatives yet still suffering the high demand of memory consumption. One desirable sequence comparison method at large-scale requires succinctly-organized sequence data management, as well as prompt sequence retrieval given a never-before-seen sequence as query. ResultsIn this paper, we proposed a novel approach, referred to as SAINT, for efficient and accurate alignment-free sequence comparison. Compared to existing alignment-free sequence comparison methods, SAINT offers advantages in two aspects: (1) SAINT is a weakly-supervised learning method where the embedding function is learned automatically from the easily-acquired data; (2) SAINT utilizes the non-linear deep learning-based model which potentially better captures the complicated relationship among genome sequences. We have applied SAINT to real-world datasets to demonstrate its empirical utility, both qualitatively and quantitatively. Considering the extensive applicability of alignment-free sequence comparison methods, we expect SAINT to motivate a more extensive set of applications in sequence comparison at large scale. AvailabilityThe open source, Apache licensed, python-implemented code will be available upon acceptance. Supplementary informationSupplementary data are available at Bioinformatics online.

MotivationSeveral popular methods exist to predict function from multiple protein-protein association networks. For example, both the Mashup algorithm, introduced by Cho, Peng and Berger, and deepNF, introduced by Gligorijevi[c], Barotand, and Bonneau, analyze the diffusion in each network first, to characterize the topological context of each node. In Mashup, the high-dimensional topological patterns in individual networks are canonically represented using low-dimensional vectors, one per gene or protein, to yield the multi-network embedding. In deepNF, a multimodal autoencoder is trained to extract common network features across networks that yield a low-dimensional embedding. Neither embedding takes into account known functional labels; rather, these are then used by the machine learning methods applied after embedding. ResultsWe introduce MELISSA (MultiNetwork Embedding with Label Integrated Semi-Supervised Augmentation) which incorporates functional labels in the embedding stage. The function labels induce sets of "must link" and "cannot link" constraints which guide a further semi-supervised dimension reduction to yield an embedding that captures both the network topology and the information contained in the annotations. We find that the MELISSA embedding improves on both the Mashup and deepNF embeddings in creating more functionally enriched neighborhoods for predicting GO labels for multiplex association networks in both yeast and humans. AvailabilityMELISSA is available at https://github.com/XiaozheHu/melissa

AO_SCPLOWBSTRACTC_SCPLOWThere is a considerable body of work in the field of computer science on the topic of sparse graph recovery, particularly with regards to the innovative deep learning approaches that have been recently introduced. Despite this abundance of research, however, these methods are often not applied to the recovery of Gene Regulatory Networks (GRNs). This work aims to initiate this trend by highlighting the potential benefits of using these computational techniques in the recovery of GRNs from single cell RNA sequencing or bulk sequencing based gene expression data. GRNs are directed graphs that capture the direct dependence between transcription factors (TFs) and their target genes. Understanding these interactions is vital for studying the mechanisms in cell differentiation, growth and development. We categorize graph recovery methods into four main types based on the underlying formulations: Regression-based, Graphical Lasso, Markov Networks and Directed Acyclic Graphs. We selected representative methods from each category and made modifications to incorporate transcription factor information as a prior to ensure successful reconstruction of GRNs.