Unpaired Single-Cell Dataset Alignment with Wavelet Optimal Transport

We thank the reviewer for their time, effort, and constructive feedback. We address the reviewer’s concerns and questions below:

---

Comparing at multiple scales is the primary motivation for this approach. I think it would be helpful to show (e.g., in experiments) which scales are most relevant and how much they contribute to performance.

Please view our response to this question in the “Response to Common Questions & Concerns” comment.

---

There are limited discussions of when to use which aggregation method (sum, max, potentially other weighting schemes) and how this choice affects results. The overall impression is that there are many tunable components (scale selection, filtering approach, aggregation method, kernel choice) without sufficient guidance or understanding of their interactions.

Please view our response to this question in the “Response to Common Questions & Concerns” comment.

---

L-WOT may require some discussion on when or whether it converges at all.

Since L-WOT is a bilevel optimization problem where both the inner and outer loops are nonconvex, it is difficult to make any rigorous statements about the convergence of L-WOT. However, in practice, we find that limiting the number of outer loops to 100 is sufficient to obtain good performance in experiments.

---

Is entropy always a good heuristic to emphasize informative scales?

We added a new section (Appendix 5.2.2), which explores the different wavelet scales and the informative ones. Specifically, in subsection “Filter Scales” of Appendix 5.2.2, we find that entropy closely matches the ideal filter for both single-cell datasets, demonstrating that is a good heuristic at least for scGEM and SNARE-seq.

However, heuristics by definition can not always be good in all scenarios; the entropy heuristic is no exception. For instance, if one scale is uniformly distributed at random, the entropy heuristic would emphasize this scale even though it is complete noise. In practice, as shown in the experiments, we still obtain good results with this heuristic even if there are potential cases in which it may fail.

---

Using KDE with tunable bandwidth seems to introduce further complexity (in the Appendix, the authors report fixing Gaussian KDE bw at 0.4; have the authors considered adaptive bandwidth like Scott's or Silverman's?)

We did not try adaptive bandwidths for the Gaussian KDE since each dataset is unit normalized.

---

It is not clear whether the experiments exclusively use the Chebyshev polynomial approximation. The paper would be strengthened if there is some analysis/comparison between the Chebyshev approximation vs the exact computation to provide some understanding of the trade-off between speed and performance.

Our implementation exclusively uses the Chebyshev polynomial approximation, as stated in Section 3.1. This approximation is well-established in spectral graph wavelets literature [1] providing extensive theoretical guarantees and empirical validation. Since our work builds upon this foundation, we believe that re-validating the Chebyshev approximation would not substantially strengthen our contribution or provide new insights beyond what is already established in the literature. Our experimental results across multiple settings demonstrate that using the Chebyshev approximation effectively serves our method objective in handling noise and non-isometry in single-cell dataset alignment.

---

I think that runtime/memory benchmarking for all methods should be discussed more thoroughly (beside what is provided in Appendix E). What are the runtime for each method in each experiment? Moreover, constructing large, fully-connected graphs is expensive. It helps to understand the computational/memory cost when evaluating different methods for practical applications.

We believe the runtime/memory benchmarking is comprehensively covered in Appendix E where we provide

1. A detailed comparative table showing how each method scales with respect to feature dimensionality, number of samples, number of scales, and choice of wavelet kernels
2. Explicit timing benchmarks comparing our methods against GW-OT across different dataset sizes (from n=100 to n=10,000)

Regarding fully connected graphs, this is a preprocessing step for computing intra-dataset distances and is not part of the core optimal transport algorithm - both our method and baselines require this distance computation.

However, if there are specific aspects of the computational analysis you would like us to elaborate on, we would be happy to clarify those sections.

---

Why is this work not positioned as a general framework for handling noisy, non-isometric spaces (many of which the authors leave as future directions)? Are there inherent limitations to it when applied to other domains?

Please view our response to this question in the “Response to Common Questions & Concerns” comment.

---

We thank the reviewer for their time, effort, and constructive feedback. We address the reviewer’s concerns and questions below:

---

Few misspellings in the main text (e.g., in line 012, 'throughput' should be 'throughout').

This is not a misspelling. High throughput technologies are instruments that can analyze and create data for a large batch of samples.

---

Figure 1 lacks clarity in illustrating both the task and framework.

Can you be more specific about what part of the figure lacks clarity?

The task is visualized through two graphs (X and Y) with nodes that need to be matched across datasets, specifically highlighting how nodes C and D in X correspond to nodes 3 and 7 in Y. The framework is explicitly shown through:

• Multiple layers representing different scales (S0, S1, S2...)
• The spectral graph wavelet coefficients (ψ) that define relationships between nodes
• A visual representation of how these multi-scale views contribute to the total cost C
• The filter F that removes uninformative scales and wavelets The figure directly communicates our core concept: WOT finds matches between points by considering their relationships at multiple scales, illustrated through the layered representation and the mathematical formulation above.

If you have more specific points for improvement, we would be happy to include them.

---

Lack of comparisons with more recent single-cell paired alignment methods such as Harmony and scDML.

In unpaired dataset alignment, the methods referenced in Table 2 are current state-of-the-art methods.

Note that while many newer works may fall into the category of data alignment, they often have very strong assumptions between the data spaces and incorporate those assumptions into their method. For example, Harmony [1] assumes knowledge of cell type labeling or batch information. Other examples include assuming a prior knowledge graph between features [3], utilizing weakly linked features [2], and more [4,5]. In short, these methods are not truly "unpaired alignment" methods and thus would not be a fair and meaningful baseline. Additionally, these methods cannot operate on modalities such as brightfield imaging/scRNA-seq, etc., or other modalities outside of single-cell biology where these assumptions may not hold.

Furthermore, scDML is a data clustering technique, not a cross dataset alignment technique. Specifically, it corrects for batch effect for datasets in the same space (i.e. given two RNA-seq datasets, cluster the two). This is a different task altogether than ours because we’re trying to align datasets not in the same space (i.e. given an RNA-seq dataset and an ATAC-seq dataset, align the two).

---

The description of the weakness of related works in Section 2.1 remains for more detailed analysis about these methods.

Our paper provides an analysis of related works' limitations through both general discussion and experiments:

• Section 2.1 introduces existing methods and their general limitations
• Section 4.1 explicitly demonstrates how baseline methods like GW fail in high-noise scenarios
• Section 4.2 shows baseline limitations in handling non-isometric relationships
• Real single-cell experiments in Section 4.3 validate these limitations in practice

Rather than just stating weaknesses, we systematically demonstrate them through carefully designed experiments that isolate failure modes. The progression from controlled experiments to real data provides clear, empirical evidence of where and why existing methods fall short (noise and non-isometry). If the reviewer has specific aspects of the baseline methods they feel warrant deeper analysis, we welcome more detailed feedback.

---

The study lacks a broader range of experimental scenarios involving real single-cell multi-omic datasets. Results would be more convincing with diverse single-cell multi-omic datasets from different tissues and sequencing technologies, to demonstrate the method's effectiveness across various empirical conditions.

We agree that aligning samples across different tissues would be useful. However, our current choice of the scGEM and SNARE-seq datasets was deliberate and comprehensive for a couple of reasons:

1. These datasets represent distinct biological scenarios and technical challenges:
   • scGEM captures a dynamic process (cell reprogramming) with gradual state transitions
   • SNARE-seq represents discrete cell types with clear cluster boundaries
   • They use different measurement technologies (gene expression/DNA methylation vs RNA/chromatin accessibility)
2. These datasets are well-established benchmarks in the field, enabling direct comparison with multiple baseline methods using standardized evaluation metrics.

While additional datasets could provide further validation, our current experiments already demonstrate WOT's performance across significantly different biological contexts and technical conditions.

---

We thank the reviewer for their time, effort, and constructive feedback. We address the reviewer’s concerns and questions below:

---

I read it through multiple times and it was still quite confusing and not sure who the audience is. This is clearly not going to be the computational biology practitioners as the paper is too technical without providing connections and insights sufficiently to the single-cell applications.

We respectfully disagree that the paper lacks clear motivation and connection to single-cell biology. Our method was specifically designed to address two critical challenges unique to single-cell data alignment: high technical noise and non-isometric relationships between modalities. These are not arbitrary technical innovations looking for applications - they directly respond to fundamental challenges in the single-cell field as detailed in our introduction.

Unlike traditional GW methods that were made for isometric matching, WOT was intentionally made to handle the complex realities of single-cell data. The progression of our experiments deliberately and consistently demonstrates this connection back to single-cell: we first (1) isolate and validate WOT's ability to handle noise (bifurcation experiment) and non-isometry (shape correspondence experiment) before (2) showing its effectiveness in real single-cell datasets.

We also respectfully disagree that the paper's technical depth makes it inaccessible to computational biology practitioners. Computational biology, particularly single-cell analysis, regularly uses technically challenging mathematical methods from other fields - from manifold learning in trajectory inference to probabilistic modeling. Our target audience is precisely these practitioners.

Additionally, much like how GW was initially developed for object matching but has now been applied to even the field of single-cell [1], our method is initially developed for single-cell but can also be applied in future works to other domains; this general applicability is a benefit rather than a drawback.

---

It was not clear what alignment of "unpaired" single-cell dataset would actually mean. What would be the use case for this, how would biologists benefit from such alignment for scientific research?

As stated throughout the introduction and mathematically explicit in lines 158-161, unpaired alignment between single cell datasets means that we want to define the mappings between cells of one modality (like scRNA-seq) to cells of another modality (like ATAC-seq) without any apriori knowledge on the joint distribution of these modalities.

Since each modality provides a different measurement (i.e. view) of the cell state, biologists often want to have data from multiple modalities from the same cell. However, this case is often impossible since many measurements are destructive, meaning you cannot conduct multiple measurements on the same cell. Hence, having unpaired single-cell alignment methods like WOT is important because it effectively provides data from multiple modalities on the “same” cell even when it is biologically impossible to measure these modalities on the same cell.

---

The authors need to keep referring back & connecting to the single-cell example/application throughout the presentation of the algorithm. Section 3 is almost entirely separated from the single-cell application, and it is unclear how these two are related. For instance, the "scale" is discussed a lot in the algorithm, but what does scale actually mean in single-cells?

We agree that some terms like “scale” in the methods section may be ambiguous in how they relate to single-cell. However, since Section 3 is our Methods section, we deliberately structured the paper to first establish the mathematical foundation of WOT before demonstrating its practical relevance to single-cell biology in the Experiment section. This is standard practice for methods papers, where a clean technical presentation enables readers to fully understand the approach before seeing its application.

To clarify the meaning of scale, it intuitively represents the high or low frequency patterns of the single cells dataset. While there is not a direct a one-to-one correspondence with a specific biological variation, one could imagine the low-frequency scales (higher valued scales) corresponding to global, slow-moving patterns like cell type while high-frequency scales (lower valued scales) represent local, fast-moving patterns like noise. Each scale represents a different level of pattern resolution from the more global scale information to local scale information. The goal of WOT is to separate the various scales of each single-cell dataset and align the most common scales, which directly reduces noise and non-isometry between different datasets.

To further address these concerns, we added a new section (Appendix D.2.2) to visualize and attempt to explain the scales of wavelets for single-cell datasets.

---

We thank the reviewer for their time, effort, and constructive feedback. We address the reviewer’s concerns and questions below:

---

Weakness 1

While we acknowledge that the range of flexibility with our framework can be initially daunting, we want to make clear that most hyperparameters can be fixed to default values as specified in Table 3 of Appendix C. There are only three hyperparameters that are variable between each experiment: (1) the entropic $\epsilon$ regularization, which is a hyperparameter common to any entropic OT method, (2) the aggregation operation, and (3) the weight normalization of the graph edges.

We added Algorithm 2 to Appendix C, which summarizes the unsupervised hyperparameter procedure originally described at the beginning of Appendix C in our initial submission. We do not use result metrics from experiments to guide hyperparameter selection.

This unsupervised procedure explains the different parameter values observed across datasets. Importantly, we default to sum aggregation and RBF normalization in most cases, only deviating when necessary according to Algorithm 2. The $\epsilon$ parameter shows the most variation, which aligns with previous work showing that fixed defaults often yield uninformative transport plans [1]

Further, to address your concerns about how we initially selected the defaults for the hyperparameters, we added Table 4, describing where each value came from.

---

Question 1

We have written a new section (Appendix C.2) with details on hyperparameter selection for baselines.

---

​​Question 2

Multiple refers to evaluating multiple kernels in the experiment. For instance, in Experiment 3, we evaluated both a simple tight kernel and the heat kernel. We have added this clarification in Table 5. The specific scales are determined by the wavelet kernel, which we default to the values in [3]

---

Question 3

We directly use the reduced-dimension single-cell datasets provided by [1] and [5]. For completeness, we summarize their procedure:

• scGEM: both datasets are reduced to their corresponding dimensions through normalized PCA
• SNARE-seq: RNA-seq was reduced through PCA and ATAC-seq was reduced through the topic modeling framework cisTopic [6]

---

​​Question 4

Thank you for pointing this out. We use the quadratic loss $L(a,b) := \frac{1}{2}(a - b)^2$. We have revised the beginning of the experiment section to include this important detail.

---

​​Question 5

For the bifurcation experiment, we use the Euclidean distance as the intra-domain distances for both GW and our method. We focused on comparing WOT with only GW-OT in the bifurcation experiment since GW-OT forms the theoretical foundation of WOT (as shown in Remark 1). This controlled comparison allowed us to specifically demonstrate the benefits of incorporating wavelets into the optimal transport framework, particularly in high-noise scenarios.

---

Question 6

Our work specifically focuses on the unpaired setting since it reflects the reality of single-cell data collection - paired samples are often impossible or prohibitively expensive to obtain. While we haven't evaluated performance with validation data, this was a deliberate choice aligned with our goal of developing methods that work effectively "out of the box" in real-world scenarios where paired samples are unavailable. Evaluating performance with validation data would be an interesting direction for future work, but it would not address the core unpaired challenge.

---

[1] SCOT: single-cell multi-omics alignment with optimal transport. Journal of computational biology, 2022.

[2] Large sample analysis of the median heuristic. 2017.

[3] PyGSP: Graph Signal Processing in Python. 2014.

[4] Unsupervised topological alignment for single-cell multi-omics integration. Bioinformatics, 2020.

[5] MATCHER: manifold alignment reveals correspondence between single cell transcriptome and epigenome dynamics. Genome Biology, 2017.

[6] cisTopic: cis-regulatory topic modeling on single-cell ATAC-seq data. Nature Methods, 2019.

Edit: I previously posted this under the wrong comment, deleted it and moved it here.

Dear authors,

Thank you very much for responding to my questions and updating the document. I do believe single-cell multi-omic integration in a fully unpaired setting, as described in the paper, is an important task for scientific applications, so I think the details of self-tuning procedure are important. My initial worry was that the hyperparameters may have been chosen with some validation data (since the descriptions on this were missing), leading to unfair benchmarking, but it appears this is not the case, so I will incline towards raising my score. However, I have some additional questions and concerns regarding the self-tuning procedure in Algorithm 2 in Supplementary Section C. My understanding is that you are doing the following:

1. Start out by a very small $\epsilon$ value, $10^{-4}$, and use sum for aggregation function and RBF for norm. If this already gives a valid coupling, just use these hyperparameters.
2. But small $\epsilon$ values could give degenerate couplings, e.g. due to NaNs. In that case you gradually increase $\epsilon$ until you get a valid coupling (no NaNs) thats sufficiently different than a uniform coupling, "sufficiently different" described by $\eta$.
3. If you've raised $\epsilon$ up to $\epsilon$ $\geq 1$ and still don't have a valid coupling, you cycle through this procedure trying out different aggregation functions first, and then different norm values.

Based on this here are my questions:
1. Tending to small $\epsilon$s (i.e. stopping hyperparameter tuning at the first point of sufficiently high enough $\epsilon$ to get a valid coupling) makes sense for aligning datasets that are already paired (e.g. SNARE-seq, scGEM experiments from the paper). These datasets have underlying 1-to-1 mapping of cells since they were jointly measured so the ground-truth coupling is super sparse (I am aware you aren't using this info when solving the coupling, just using it for benchmarking). However, this likely won't be the case for real-world separately sequenced datasets since the underlying biological manifold won't have 1-to-1 matches between cells. In 1-to-many settings, higher $\epsilon$s would likely be more favorable, so this self-tuning procedure may not be as effective. I don't know if this is practically possible in the next 5 days of the rebuttal phase but: Could you test the method out on a couple of real-world datasets, where we may have some confidence in cell-type annotations so they could be used for benchmarking (e.g. ideally through FACS sorting but that could be hard to find, expert annotated data from a bio publication could be fine too) ? If possible, It would be additionally nice to see how large the discrepancy is between the quality of alignment with the self-tuned procedure vs a small grid of varying hyperparameters to see where in the range the self-tuned hyperparameters are performance-wise.

2. More of a concern than a question: It appears that you don't consider $\rho_1$ and $\rho_2$ in the self-tuning procedure and fix them to 1.0. This works well in your practice because all the dataset you showcase your algorithm on are essentially balanced dataset; they have the same number of cells and same proportion of cell types because they came from experiments where measurements were jointly taken. This will not be the case for real-world applications and when datasets are unbalanced, $\rho$ values really matter based on my practical experience with unbalanced OT. I think testing your algorithm on some real-world separately sequenced datasets will likely show this.

3. Shouldn't the threshold value eta be dependent on size of T matrix (number of samples)?

4. Are the runtimes reported in the table on page 22 for a single combination of hyperparameter? I ask because $\epsilon$ updates in Algorithm 2 are going to be super small initially. For a starting value of $\epsilon=10^{-4}$, the next 10 $\epsilon$ values would be: 0.00010000005, 0.0001000001, 0.00010000015, 0.0001000002, 0.00010000025, 0.0001000003, 0.00010000035, 0.0001000004, 0.00010000045, 0.0001000005, 0.00010000055. These are tiny increments and for each update, the whole spectral wavelet GWOT is run again. This seems like quite a computationally heavy procedure and could prohibit the adoption of the method in real-world applications. If you run your algorithm on a couple of real-world datasets, could you also report how long the self-tuning procedure takes for these?