Hyperbolic Genome Embeddings

We thank the reviewer for this positive assessment of our work. We appreciate your recognition of the paper's many strengths, including its significant impact on the genomics research community, the rigorous and high-quality execution of experiments, and the comprehensive benchmarking efforts. We are also grateful for the thoughtful questions you raised which prompted discerning experiments that have strengthened the paper and unified the motivating arguments.

I believe the analysis could be improved by investigating exactly how the HCNN represents sequences related through evolution. A common intuition I have seen for hyperbolic embeddings in computer vision is that datum closer to the origin tend to represent high level features and vice versa. It would be convincing if some analogy could be made for the genomic sequences through the experiments. For example, while Figure 2 demonstrates that hierarchical structure exists in the HCNN representation, it would be interesting to see whether the sequences placed at the "top" of the learned hierarchy reconcile with biological intuition.

Thanks for bringing up an interesting consideration. Overall, we report two findings:

1. Sequences embeddings near the top of the hierarchical structure represent ambiguous sequences (visualized in Figure 15 of the updated manuscript).
2. Ablating conserved regions of sequence (corresponding to regions with strong evolutionary signal) results in the loss of definitive features for sequences, shifting them closer to the top (more ambiguous part) of the hierarchy (visualized in Figure 16 of the updated manuscript). Thus, sequences with definitive features with respect to the class label lie closer to the Poincaré disk boundary.

In exploring the sequence representations of HCNNs, we started with the intuition built in [1], where hyperbolic image embeddings of MNIST (Figure 1 of [1]) near the center of the Poincaré disk represent the most ambiguous looking digits, while clear images lie near the boundary. Similarly, in the new Figure 15, we observe that in the processed pseudogene dataset in TEB, the sequence embeddings that lie close to the center of the Poincaré disk (the top of the hierarchy) correspond to low confidence embeddings for HCNNs, while the embeddings near the disk boundaries show the highest classification confidence. This is consistent with the idea that well defined sequences are at the bottom of the hierarchy where there is more space to separate out nuanced differences between sequences based on distinctive features.

Next, since we cannot readily identify features from DNA sequence as one can with images, we conduct an experiment to dissect the sequence features informing the hyperbolic genome embedding. Given our dataset of processed pseudogenes, we examine the changes made to the HCNN representation by perturbing a fixed pseudogene sequence. For a fixed sequence, we:

1. Compute the Genomic Evolutionary Rate Profiling (GERP) [2] score for each nucleotide along the sequence. GERP scores quantify evolutionary constraints at specific genomic positions, identifying which positions are functionally important based on selective pressure.
2. Mutate a fraction of the nucleotides under the highest selective pressure.
3. Use our HCNN to generate an embedding for this perturbed instance of our original sequence.

Figure 16 visualizes this experiment using a processed pseudogene sequence and a background sequence. As the evolutionary signal under strong selection is eroded by the introduced mutations, it is likely that the features that make the pseudogene more “gene-like” are degraded. This degradation ultimately makes the sequence more ambiguous to the HCNN, and we observe that the perturbed representations move closer to the top of the hierarchy (near the center of the Poincaré disk), where the low confidence sequences lie. Removal of these evolutionary features actively hinders the model in identifying pseudogenes. However, perturbing conserved regions from the noisy background sequences does not appear to have this effect, as the model focuses on learning features common to the pseudogene class. Please see section A.10 of the updated paper for visualizations and writeup.

[1] Khrulkov, V., Mirvakhabova, L., Ustinova, E., Oseledets, I., & Lempitsky, V. (2020). Hyperbolic image embeddings. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (pp. 6418-6428).

[2] Cooper, G. M., Stone, E. A., Asimenos, G., Green, E. D., Batzoglou, S., & Sidow, A. (2005). Distribution and intensity of constraint in mammalian genomic sequence. Genome research, 15(7), 901-913.

Thank you for conducting a thorough rebuttal, it was genuinely exciting to see the results of the additional analysis. I found the perturbation experiments particularly creative and compelling. I believe the set of new experiments successfully identifies the hypothetical use-cases and representational advantages of the HCNN in genomics, and consequently will be of high practical value to the community. I will raise my score to reflect the improved manuscript.

Thank you for taking the time to write your thoughtful review. We are glad to see that you recognize the novelty of our approach and the comprehensiveness of our evaluation. We are very appreciative of the discerning comments you have provided, and believe that the questions you have raised have led to illuminating experiments that have strengthened the paper and unified the motivating arguments. We provide specifics in our inline responses.

Authors propose a new dataset which can evaluate ability of models to classify transposable elements. This class of sequences are widely abundant however it is unclear as to their biological significance compared to sequences encoding mature RNAs.

Thanks for raising this point - we realize that we did not provide adequate background on the biological significance of transposable elements (TEs). We have clarified this in the paper by adding the following explanation (with citations):

TEs can influence gene expression and regulation by acting as alternative promoters, providing transcription factor binding sites, introducing alternative splicing, mediating epigenetic modifications, and being co-opted as regulatory elements throughout evolution. As such, TEs have been also implicated in disease pathogenesis.

W1: This work can benefit from additional justification for why the hierarchical tree structure is a good inductive bias for genomic data. Although the original data generating process can take on a tree structure is it necessarily the case that we would want to add this as an inductive bias to the model?

Ultimately, most of the inductive biases in deep learning genomics models have been inherited from the language and image deep learning domains. While some of these inductive biases are applicable to genomics, many lack strong theoretical or empirical justification and persist primarily due to the absence of alternative frameworks. Our work functions as a challenge to the status quo reliance on Euclidean models in the field of genomics.

Given that evolution is the guiding principle behind all sequence organization [1, 2], incorporating an inductive bias that mirrors evolutionary processes is essential. We demonstrate this necessity through two key approaches:

1. Empirical Evidence: We show strong, statistically significant performance improvements across a wide variety of genomics problems pertinent to the field.

2. Theoretical Framework: We reformulate a well-established measure of hyperbolicity to illustrate the alignment between our approach and the unique, inherent characteristics of genomic data.

[1] Shapiro, J. A. (2002). Genome organization and reorganization in evolution: Formatting for computation and function. Annals of the New York Academy of Sciences, 981(1), 111–134. https://doi.org/10.1111/j.1749-6632.2002.tb04917.x

[2] Lynch, M., & Conery, J. S. (2006). The origins of eukaryotic gene structure. Molecular Biology and Evolution, 23(2), 450–468. https://doi.org/10.1093/molbev/msj050

Thank you for taking the time to write this detailed review. We are glad to see that you enjoyed the writing and presentation of our work, as well as the background motivation. We also appreciate your recognition of the Transposable Elements Benchmark as a valuable contribution to the research community.

Our initial paper was insufficiently clear in that our main goal is to highlight the compatibility of the HCNN framework as specifically well suited to the genomics domain. While we ourselves do not present a GFM in this paper, we provide multiple lines of evidence that future GFMs/GMs may benefit from the inductive biases inherent in hyperbolic architectures. We believe that most of your concerns stem from comparisons made to existing GFMs, however our intention was to open up a new class of modeling approaches for genomics with our work. As such, if we have better recontextualized the core contribution of the paper and dispelled your concerns regarding model efficiency, we hope that you might be willing to raise your score. We hope that by starting a dialogue early, we will be able to address any responses you may have.

We’ll now respond inline to specific comments:

Comparing the number of parameters between CNN-based models and Transformer-based models is not convincing. Due to the underlying difference between these two model architectures, fewer parameters do not guarantee efficiency. A more suitable comparison would be the time and memory usage of training the models on the same dataset.

We completely agree that fewer parameters do not always guarantee efficiency. You raise an important point in terms of examining model efficiency across multiple criteria, which carries practical importance for the wider research community. In light of your comment, we compared the time, memory, and computational burden of training the four classes of models that achieved SOTA on GUE: NT-2500M multi, DNABERT-2, HCNN-S, and HCNN-M. Our results are summarized in the following table:

Memory: We report the minimum required VRAM necessary to fit each model to the same GUE dataset. For HCNN-S and HCNN-M, we report the VRAM required to train the models from scratch, and for DNABERT-2 and NT-2500M-multi we report the VRAM necessary to finetune the model on the dataset. We use an asterisk to note that for NT-2500M-multi, the authors reported fine-tuning on eight A100 GPUs (640 GB), which is likely closer to the recommended VRAM [1].

Time: We compute the wall-clock time by fitting our models on the same GUE dataset. We compare times across two machines - a relatively smaller machine 1 (two RTX 3090s), and a larger machine 2 (four L40s) that can support NT-2500M-multi. We use the recommended 3 epochs for fine-tuning DNABERT-2/NT-2500M-multi per the DNABERT-2 repo, and use the default number of training epochs (100) for HCNN-S/HCNN-M. We additionally report the relative increase in time compared to the fastest model (HCNN-S) in parentheses next to the absolute time.

FLOPs: We measure each model’s computational cost by considering the relative Floating Point Operations (FLOPs)—which is the total number of multiplication and addition operations during a forward and backward pass—compared to HCNN-S, which performs the lowest number of FLOPs.

From the following experiments, we can conclude HCNNs maintain supremacy per all relevant criteria: parameters, runtime, computation cost, and memory. Across most criteria, the difference is often several orders of magnitude in improvement, cementing HCNNs as a viable resource for the genomics research community.

There is a lack of baselines. As a CNN-based model, the model should also be compared with HyenaDNA and Caduceus.

Thanks for the suggestion - we agree that it may be helpful to include comparisons to diverse, non-Transformer/CNN architectures such as HyenaDNA and Caduceus. We have added HyenaDNA and Caduceus to Table 4 (which is now Table 5 in the revised paper). HCNNs outperform HyenaDNA on 27/28 datasets, and Caduceus on 19/28 datasets.

Results are reported here:

We are pleased to report that the results of these experiments provide a clearer understanding of how phylogenetic tree structure is learned. Our models primarily operate by uncovering the underlying phylogenetic tree structure amidst noise. While this learning mechanism may also offer advantages in disentangling distinct evolutionary patterns (Scenario B), the results consistently indicate that excelling in Scenario C is somewhat orthogonal to excelling in Scenario A. This is likely because sequence differentiation is performed at the tree level rather than at the clade level. Notably, this finding aligns with results observed on real-world datasets: the most significant performance gains are observed in scenarios where an evolutionary signal (e.g., transcription factor binding sites, epigenetic marks, or transposable elements) is effectively distinguished from background noise (e.g., background sequences).

Overall, we find these results to support our original hypothesis discussed in the introduction of our paper - the inductive biases of our models extend past approximating explicit edit distance mutations to instead focus on latent hierarchical structure that likely occupies a lower dimensional manifold.

For experiment 2, we observe that in the processed pseudogene dataset in TEB, the sequence embeddings that lie close to the center of the Poincaré disk (the top of the hierarchy) correspond to low confidence embeddings for HCNNs (Figure 15), while the embeddings near the disk boundaries show the highest classification confidence. This is consistent with the idea that well defined sequences are at the bottom of the hierarchy where there is more space to separate out nuanced differences between sequences based on distinctive features. Since we cannot readily identify features from DNA sequence, we conduct an experiment to dissect the sequence features informing the hyperbolic genome embedding. Given our dataset of processed pseudogenes, we examine the changes made to the HCNN representation by perturbing a fixed pseudogene sequence. For a fixed sequence, we:

1. Compute the Genomic Evolutionary Rate Profiling (GERP) [1] score for each nucleotide along the sequence. GERP scores quantify evolutionary constraints at specific genomic positions, identifying which positions are functionally important based on selective pressure.
2. Mutate a fraction of the nucleotides under the highest selective pressure
3. Use our HCNN to generate an embedding for this perturbed instance of our original sequence.

Figure 16 visualizes this experiment using a processed pseudogene sequence and a background sequence. As the evolutionary signal under strong selection is eroded by the introduced mutations, it is likely that the features that make the pseudogene more “gene-like” are degraded. This degradation ultimately makes the sequence more ambiguous to the HCNN, and we observe that the perturbed representations move closer to the top of the hierarchy (near the center of the Poincaré disk), where the low confidence sequences lie. Removal of these evolutionary features actively hinders the model in identifying pseudogenes. However, perturbing conserved regions from the noisy background sequences does not appear to have this effect, as the model focuses on learning features common to the pseudogene class.

[1] Cooper, G. M., Stone, E. A., Asimenos, G., Green, E. D., Batzoglou, S., & Sidow, A. (2005). Distribution and intensity of constraint in mammalian genomic sequence. Genome research, 15(7), 901-913.

Thank you for taking the time to review our work. We appreciate your recognition of hyperbolic architectures as a promising avenue for genome sequence modeling.

How does the CNN/HCNN performance scale with number of parameters? Figure 5 seems to indicate that performance gets worse with increasing hidden dim?

Thanks for raising this point - we would like to clarify that Figure 5 shows the average net increase in performance is higher at lower dimensions when using HCNNs compared to CNNs. Given the importance of learning compact representations in all deep learning domains, we intended this figure to demonstrate that HCNNs offer an added benefit in this area. We have added a second plot to Figure 5 that shows that HCNN/CNN performance improves with increasing hidden dimension size.

Additionally figure 5 should include error bars. Also the cumulative improvement y-axis is unclear. It should just be the average MCC value.

We agree that this change would add clarity to the figure, so we have updated it to show the average MCC value with error bars.

Figure 1 shows that sequences which obey a phylogenetic tree structure are embedded into a hyperbolic space which learns this structure. However I do not see any results which indicate that the phylogenetic structure is being learned.

Thank you for raising an intriguing and relevant point. We have conducted a number of additional experiments to elucidate the ways in which our models learn phylogenetic structure (experiment 1), provide interpretability for the learned phylogenetic embeddings (experiment 2), and demonstrate a concrete example of how our learning process confers advantages for phylogenetic data over existing models (experiment 3). We summarize the contributions by experiment here:

Exp 1). Creating synthetic datasets using an evolutionary process simulator (section 4.1): we extend our synthetic dataset experiment to investigate three potential scenarios underlying the data-generating processes of sequences. We determine that our models primarily operate by uncovering underlying phylogenetic tree structure amidst noise.

Exp 2). In silico mutagenesis (section A.10): We investigate how the phylogenetic structure of sequences are represented by HCNNs through a perturbation experiment. We find that sequences near the top of the learned hierarchical structure represent ambiguous sequences, and ablating conserved regions of sequences results in the loss of definitive features, shifting them closer to the top of the hierarchy.

Exp 3). Homology splitting (section A.7): We show that the inductive biases of HCNNs enable improved generalizability to unseen homology branches, which is an important learning priority for all biological sequence models.

We further elaborate on each experiment to provide necessary details. For experiment 1, we reconsidered the various plausible data-generating processes for biological sequences. We defined three potential cases of biological signal transmission that might be learned by the models, and constructed 3 classes of synthetic datasets (visualized in the new Figure 2):

(A). Intra-tree differentiation: this was our original synthetic experiment. In this scenario, we generated sequences from a single phylogenetic tree, with labels derived from clade membership. This experiment tested whether a hyperbolic model would better differentiate parts of a single phylogenetic tree from each other.

(B) Inter-tree differentiation: in this scenario sequences are generated from different phylogenetic trees, with labels derived from phylogeny membership. This experiment tests whether the model is differentiating different evolutionary processes from one another.

(C) Tree identification: In this scenario, sequences for one label are generated from a phylogenetic tree, whereas the other label identifies sequences that are either completely randomly generated or sampled from background sequence (and therefore without any guiding evolutionary relationship). This experiment tests whether the model is differentiating evolutionary processes from noise.

We additionally tested two variations of synthetic datasets - one in which sequences were completely randomly generated, and one in which sequences are sampled from existing genomes to see if the fully artificial sequences we originally used were not generalizable to our real world dataset results.