Continual integration of single-cell multimodal data with MIRACLE

Single-cell continual integration with CLThe scenario of single-cell continual integration

Single-cell continual integration aims to continuously incorporate new data into the atlas. Let \({{\mathcal{X}}}_{t}={\{{{\mathbf{x}}}_{t,n}\}}_{n = 1}^{{N}_{t}}\) be the batch of Nt cells arriving at time t, with xt,n the observation of cell n, and let \({{\mathcal{Y}}}_{t}={\{{{\mathbf{y}}}_{t,n}\}}_{n = 1}^{{N}_{t}}\) be its integration results, where yt,n may include the low-dimensional embedding or batch-corrected data of cell n. The goal is to update all results \({{\mathcal{Y}}}_{1:t}=\mathop{\bigcup }\nolimits_{\tau = 1}^{t}{{\mathcal{Y}}}_{\tau }\) whenever new data \({{\mathcal{X}}}_{t}\) arrives.

Offline strategy

To achieve continual integration, the common approach is an offline strategy6,7,8,9

Whenever new data \({{\mathcal{X}}}_{t}\) arrives, all data \({{\mathcal{X}}}_{1:t}=\mathop{\bigcup }\nolimits_{\tau = 1}^{t}{{\mathcal{X}}}_{\tau }\) are used to train parameters θt via the offline function foff, and the inference function g predicts the results \({{\mathcal{Y}}}_{1:t}\). Although effective, retraining on all data each time becomes increasingly inefficient.

Online strategy 1 using generalization

To avoid repeated retraining, online strategies can be adopted. One efficient approach is generalization11,13,14

The model is first pretrained offline on large external datasets \({{\mathcal{X}}}_{0}\) to obtain fixed parameters θ0, which are then generalized to new data by inferring \({{\mathcal{Y}}}_{t}\) from each \({{\mathcal{X}}}_{t}\) via g, so past results need no update. This avoids training during integration but relies on external datasets that often fail to capture unknown biological and technical variations in new data, giving poor adaptability.

Online strategy 2 using fine-tuning

To adapt to unknown changes in new data, online strategies based on fine-tuning have emerged10,12,15

At step t, fFT fine-tunes the previous parameters θt−1 on new data \({{\mathcal{X}}}_{t}\) to obtain θt, which g uses to predict \({{\mathcal{Y}}}_{1:t}\) (θ0 is pretrained or randomly initialized). However, as t grows, θt gradually forgets past information, degrading results, especially for earlier data.

Online strategy 3 (ours) using CL

To adapt to new data while preventing forgetting, we propose an online integration strategy based on CL17,18

Unlike fine-tuning, CL introduces a rehearsal memory \({{\mathcal{R}}}_{t}\subseteq {{\mathcal{X}}}_{1:t}\) that retains key data to prevent forgetting. Combining new data \({{\mathcal{X}}}_{t}\) with the previous memory \({{\mathcal{R}}}_{t-1}\), fCL updates θt−1 to θt; the memory is then updated by the sampling function h, and predictions \({{\mathcal{Y}}}_{1:t}\) are made by g. Initially, \({{\mathcal{R}}}_{0}={\uptheta}\) and θ0 is pretrained or randomly initialized.

Dynamic architecture adaptation

In continual integration, new data \({{\mathcal{X}}}_{t}\) often originate from various sequencing technologies and may introduce new modalities or features. To seamlessly incorporate these, we propose dynamically expanding the model architecture in the parameter update equation (equation (4a)). Let \({\tilde{{\mathbf{\uptheta}}}}_{t}\) be the model parameters associated with the new modalities or features. We randomly initialize \({\tilde{{\mathbf{\uptheta}}}}_{t}\) and combine it with the previous parameters θt−1 to get the initial value \({\hat{{\mathbf{\uptheta}}}}_{t}\):

On the basis of \({\hat{{\mathbf{\uptheta}}}}_{t}\), we can then obtain θt by training on \({{\mathcal{X}}}_{t}\) and \({{\mathcal{R}}}_{t-1}\).

Distribution-preserving reservoir sampling

In equation (4b), the rehearsal memory \({{\mathcal{R}}}_{t}\) grows with more iterations, harming efficiency and knowledge sharing. To bound it, we propose a batched reservoir sampling algorithm that limits the capacity of \({{\mathcal{R}}}_{t}\) when a new batch \({{\mathcal{X}}}_{t}\) arrives. To keep integration unbiased, we further combine stratified sampling to preserve the interbatch distribution with a ball tree sampling method to preserve the intrabatch distribution. We call this combined scheme DPRS.

Batched reservoir sampling for batch data

Classical reservoir sampling gives each cell an equal probability of being retained under a fixed capacity but handles only one new cell at a time. To handle new batches with many cells, we propose batched reservoir sampling.

Assume that the rehearsal memory has a capacity limit of M cells. After each integration, we add new data \({{\mathcal{X}}}_{t}\) of size Nt to the previous rehearsal memory \({{\mathcal{R}}}_{t-1}\) of size Mt−1, resulting in \({{\mathcal{R}}}_{t}\). However, if Nt + Mt−1 > M, we first sample \({{\mathcal{R}}}_{t-1}\) and \({{\mathcal{X}}}_{t}\) with sampling ratios α and β, respectively, where α and β satisfy

where Sτ = N1 + N2 + ⋯ + Nτ is the total number of cells in the first τ batches. Equation (6) ensures that \({{\mathcal{R}}}_{t}\) has size M, with new and old data proportions consistent with their original batches. From equation (6), we get α = St−1M/(StMt−1) and β = M/St. Note that, when Nt = 1, batched reservoir sampling reduces to classical reservoir sampling.

Stratified sampling for interbatch distribution

To maintain the cell distribution across batches, in addition to batched reservoir sampling, we further introduce batch-based stratified sampling on past data. Specifically, in memory \({{\mathcal{R}}}_{t-1}\), assuming that \({{\mathcal{R}}}_{t-1}^{(\tau )}\) is the data from the τth batch and is of size \({M}_{t-1}^{(\tau )}\), we sample each \({{\mathcal{R}}}_{t-1}^{(\tau )}\) with a ratio of α for τ = 1, 2, …, t − 1. This ensures that the cell proportions in \({{\mathcal{R}}}_{t}\) match those in the original batches.

Ball tree sampling for intrabatch distribution

Within a batch, the common SRS poorly preserves the data distribution at low sampling ratios54; we therefore propose a ball tree sampling algorithm.

Assume that the data before sampling are generated from the underlying distribution p(x). We divide the data space \({\mathcal{D}}\) into infinitesimal intervals \({\{{{\mathcal{D}}}_{i}\}}_{i = 1}^{\infty }\), each with a Lebesgue measure \(\lambda ({{\mathcal{D}}}_{i})\) close to 0. Then, for any x, assuming it lies in the j(x)th interval, p(x) can be approximated as

where N is the total number of cells to be sampled, Nj(x) is the number of cells in \({{\mathcal{D}}}_{j({\mathbf{x}})}\) and P(·) denotes probability. Similarly, after sampling with ratio γ, the underlying distribution \(\tilde{p}({\mathbf{x}})\) is approximated as

where γN is the total number of cells after sampling, and \({\tilde{N}}_{j({\mathbf{x}})}\) is the number of cells in \({{\mathcal{D}}}_{j({\mathbf{x}})}\) after sampling. To ensure the sampled distribution matches the original, \(\tilde{p}({\mathbf{x}})=p({\mathbf{x}})\) is required. Combining equations (7) and (8), we get \({\tilde{N}}_{j({\mathbf{x}})}\approx \gamma {N}_{j({\mathbf{x}})}\), which means that, for any interval \({{\mathcal{D}}}_{i}\), we should sample at a ratio of γ.

Partitioning \({\mathcal{D}}\) by equal intervals (equal \(\lambda ({{\mathcal{D}}}_{i})\)) handles nonuniform data poorly—lacking detail in dense regions and missing sparse ones—and the number of intervals grows exponentially with dimensionality (the ‘curse of dimensionality’). We instead use equal-frequency partitioning, with the same number of cells Ni per interval, which adapts interval size to local density and grows only linearly with the total sample size N.

For fine-grained preservation, we set Ni = 1/γ (for example, Ni = 3 when γ = 1/3) and draw one sample per interval. Because 1/γ is not always an integer, we partition at multiple resolutions and select samples across resolutions to achieve any ratio. For the partitioning, we use the ball tree55, which handles nonuniform, high-dimensional data better than the k-dimensional tree. We term this intrabatch algorithm BTS (see Supplementary Information section 3 for details).

Mosaic data handling

A key challenge in continual integration is mosaic data, where batches contain varying modality combinations. MIRACLE addresses this as a general CL framework wrapped around a base model that processes the mosaic structure at each incremental step. We selected MIDAS as the base model for its effectiveness in static mosaic integration56: MIDAS is a multimodal VAE that uses self-supervised learning to align modalities into a shared embedding space and information-theoretic principles to disentangle the biological cell state from technical noise, performing dimensionality reduction, batch correction and imputation in a single, offline training pass. MIRACLE extends MIDAS to online, incremental atlas construction through two mechanisms: (1) a dynamic architecture adaptation mechanism, to incorporate new modalities or features from incoming data, and (2) a data rehearsal mechanism, based on our DPRS approach, to mitigate catastrophic forgetting in a memory-efficient manner.

Suppose that, at time t, the set of modalities present in the observed data \({{\mathcal{X}}}_{1:t}\) is \({\mathcal{M}}\), for example, \({\mathcal{M}}=\{\,\text{ATAC},\text{RNA},\text{ADT}\,\}\). For the nth cell observation xτ,n in batch τ ∈ {1, …, t}, we denote it as

where \({{\mathcal{M}}}_{\tau }\subseteq {\mathcal{M}}\) is the set of modalities in batch τ, \({{\mathbf{o}}}_{\tau ,n}^{m}\in {{\mathbb{N}}}^{{D}_{\tau ,n}^{m}}\) is the observation vector of modality m with size \({D}_{\tau ,n}^{m}\) and sτ,n = τ is the batch ID.

Using the CL framework of MIRACLE with MIDAS as the base model, we obtain the integration result yτ,n for each cell observation xτ,n

where \({\hat{{\mathbf{o}}}}_{\tau ,n}^{m}\in {{\mathbb{N}}}^{{D}_{t}^{m}}\), with size \({D}_{t}^{m}\), is the imputed and batch-corrected count vector for modality m and \({{\mathbf{z}}}_{\tau ,n}\in {{\mathbb{R}}}^{{D}^{z}}\) is the batch-corrected low-dimensional cell embedding that represents the biological state of the cell.

MIRACLE training

MIRACLE uses MIDAS as the base model to handle single-cell multimodal data with diverse modality combinations. During the tth integration step, the model’s parameters θt are composed of those for the encoders (\({{\mathbf{\uptheta}}}_{t}^{\,\text{enc}\,}\)), decoders (\({{\mathbf{\uptheta}}}_{t}^{\,\text{dec}\,}\)) and classifiers (\({{\mathbf{\uptheta}}}_{t}^{\,\text{cls}\,}\)). The optimization of these parameters is achieved by iteratively minimizing two loss functions inherited from the MIDAS base model: the encoder–decoder loss \({l}^{f,g}({{\mathbf{\uptheta}}}_{t}^{\,\text{enc}\,},{{\mathbf{\uptheta}}}_{t}^{\,\text{dec}\,};{{\mathbf{x}}}_{\tau ,n},{{\mathbf{\uptheta}}}_{t}^{\,\text{cls}\,})\) and the classifier loss \({l}^{r}({{\mathbf{\uptheta}}}_{t}^{\,\text{cls}\,};{{\mathbf{x}}}_{\tau ,n},{{\mathbf{\uptheta}}}_{t}^{\,\text{enc}\,})\). In contrast to MIDAS, this process in MIRACLE is conducted over both new data \({{\mathcal{X}}}_{t}\) and rehearsal data \({{\mathcal{R}}}_{t-1}\), where \({{\mathbf{x}}}_{\tau ,n}\in {{\mathcal{X}}}_{t}\) for τ = t and \({{\mathbf{x}}}_{\tau ,n}\in {{\mathcal{R}}}_{t-1}^{(\tau )}\) for 1 ≤τ < t. The complete MIRACLE training process, corresponding to equation (4a), is presented in Supplementary Information section 4, with the constituent loss functions detailed in Supplementary Information section 5.

Integration of multiple batches at once

In practical applications, new data often arrive in multiple batches, denoted as \(\{{{\mathcal{X}}}_{t},{{\mathcal{X}}}_{t+1},\ldots \}\), which may also consist of mosaic data. In this case, integrating these data batches simultaneously, rather than sequentially, can enhance the efficiency of continuous integration. Leveraging MIDAS’s capability to support mosaic data, MIRACLE can be seamlessly extended to accommodate this scenario.

Datasets

All datasets used in this study are publicly available; their per-batch protocols, accessions and filtered cell numbers are detailed in Supplementary Table 12. We group them below by their role in our experiments and, for each, downloaded count matrices for gene unique molecular identifiers (UMIs), ATAC fragments and ADTs as applicable.

Large single-modality datasets

Two large RNA sequencing datasets were used to benchmark distribution-preserving sampling and scalability. DHCM23 is an snRNA-seq dataset of human dilated and hypertrophic cardiomyopathy (523,369 cells across 42 batches; Broad Single Cell Portal SCP1303). HLCA4 is a single-cell RNA sequencing (scRNA-seq) dataset of human lung (2.3 million cells from 48 studies, each treated as one batch), known for strong batch effects57.

Multimodal PBMC datasets

Ten human PBMC datasets spanning diverse modality combinations were used for bimodal, mosaic and atlas-level integration: DOGMA29 (DOGMA-seq; ATAC+RNA+ADT; four batches; GEO GSE166188), TEA30 (TEA-seq; ATAC+RNA+ADT; five batches; GEO GSE158013), TEA Multiome30 (10x Multiome; ATAC+RNA; two batches; GEO GSE158013), 10X Multiome58,59,60,61 (10x Multiome; ATAC+RNA; four batches62), WNN28 (CITE-seq; RNA+ADT; eight batches collected before HIV-vaccine administration63), DOGMA+CITE64 (one CITE-seq batch with RNA+ADT and three DOGMA-seq batches with ATAC+RNA+ADT; GEO GSE200417), ISSAAC30 (ISSAAC-seq; ATAC+RNA; one batch; ArrayExpress E-MTAB-11264), NEAT31 (NEAT-seq; ATAC+RNA; two batches; GEO GSE178707), ASAP29 (ASAP-seq; ATAC+ADT; two batches; GEO GSE156473) and ASAP-CITE29 (CITE-seq; RNA+ADT; two batches; GEO GSE156473). Subsets of these datasets were combined to form the DOTEA mosaic dataset and the PBMC reference atlas (Supplementary Tables 9 and 11).

Cross-tissue datasets

Three human immune tissues were used to assess cross-tissue integration: tonsil32 (RNA; one batch; GEO GSE165860), the BMMC29 (ASAP-seq; ATAC+ADT; one batch; GEO GSE156477) and spleen33 (RNA; one batch; GEO GSE159929).

Respiratory infection datasets

Three PBMC datasets from respiratory infections were used for cross-disease integration. COVID-1936 (CITE-seq; RNA+ADT; ArrayExpress E-MTAB-10026) spans severity levels (asymptomatic, mild, moderate, severe and critical); after quality control, we kept the two largest batches per level (ten batches total). Flu A37 (scRNA-seq whole blood; GEO GSE243629) provided seven flu A virus-infected batches (adults, pregnant women and children), from which we retained only PBMC cell types. TB38 (CITE-seq; RNA+ADT; GEO GSE158769) comprises memory T cells from active or latent infection; after quality control, we kept the 10 largest batches per type (20 batches total).

Data preprocessing

For multimodal data spanning multiple batches, we performed quality control and feature selection per batch and then unified features across batches, using the Seurat package (version 4.1.0)28 for the RNA and ADT count matrices.

For RNA, we filtered low-quality cells by the number of detected genes per cell, total UMI count and mitochondrial read percentage, normalized and log-transformed the counts with NormalizeData, removed low-frequency genes and selected the top 4,000 highly variable genes per batch with FindVariableFeatures; their union via SelectionIntegrationFeatures yielded 4,000 final features.

For ADT, we excluded low-quality cells by total protein tag count, applied centered log-ratio normalization with NormalizeData and retained all ADT features without further selection.

For ATAC fragments, we used the Signac package (version 1.6.0)65 with MACS2 peak calling66, performed quality control by transcription start site enrichment score and nucleosome signal, merged peaks with Signac’s reduce function and recounted fragments within the merged peaks. RNA UMI counts, ADT tag counts and binarized ATAC fragment counts served as model inputs.

Newly arrived data followed the same procedure; to prevent excessive expansion of ATAC features during continual integration, we aligned peaks with those from previous batches.

Third-party cell-type labels

To comprehensively evaluate PBMC datasets both qualitatively and quantitatively, we used the third-party tool Seurat for cell type annotation via label transfer. The reference dataset was the CITE-seq PBMC atlas from ref. 28. Label transfer was performed using Seurat’s FindTransferAnchors and TransferData functions with the ‘cca’ reduction method for reference mapping. When raw RNA expression data were unavailable, we generated a gene activity matrix from ATAC data using Signac’s GeneActivity function, which was then used for label transfer. However, owing to the predominance of CD4+ T cells in the NEAT dataset, finding a suitable reference for accurate annotation was challenging. Therefore, after the initial annotation using Seurat’s label transfer, we applied its FindClusters function to further cluster the data and refine the results. For the flu A dataset, originally derived from whole blood, we performed clustering in Seurat, integrated Seurat’s label transfer results with SingleR67 annotations for manual refinement, and retained only PBMC-associated cell types.

For the tonsil dataset, we used Azimuth (https://azimuth.hubmapconsortium.org), a reference-based single-cell analysis tool, to automate cell annotation on the basis of gene expression data. We then manually refined the annotations by integrating Seurat’s clustering results. Similarly, for the BMMC dataset, we employed Azimuth’s ATAC application for ATAC-based cell annotation. For the spleen dataset, we used SingleR for unbiased scRNA-seq annotation, followed by manual refinement using Seurat clustering.

Sampling methods for distribution preservation

To evaluate the distribution preservation capability of our BTS algorithm, we compared its performance with SRS, Sketch, and scSampler. Because DPRS samples in the embedding space, for each comparison we used MIRACLE to convert the dataset into low-dimensional embeddings, randomly selected 10,000 cells as the reference sample and repeated the experiment 50 times across different sampling ratios.

SRS

We used the subsample function in the Scanpy package68 to perform SRS. To ensure diverse sampling results, we set different random seeds for each run of the experiment.

Sketch

The Sketch method24 is designed to preserve rare populations during sampling. We utilized the SketchData function in the Seurat package (version 5.0.0) for sample selection69.

scSampler

The python package scSampler25 is also designed to preserve rare populations and can be downloaded from ref. 70. Since our input data were already low-dimensional embeddings, we skipped the default reduction process.

Evaluation of distribution preservation

We used the MMD metric to evaluate the consistency of data distributions before and after sampling. A smaller MMD value indicates better sampling performance. Given two sets of samples \({\mathcal{X}}=\{{{\mathbf{x}}}_{1},\ldots ,{{\mathbf{x}}}_{{N}_{1}}\}\) and \({{\mathcal{X}}}^{{\prime} }=\{{{\mathbf{x}}}_{1}^{{\prime} },\ldots ,{{\mathbf{x}}}_{{N}_{2}}^{{\prime} }\}\), MMD is defined as

where k( ⋅ , ⋅) is the Gaussian kernel function.

Implementation of comparative integration methods

We compared horizontal, rectangular and mosaic integration methods (Supplementary Table 13) in both offline and online scenarios. In the offline scenario, all batches were available from the start; in the online scenario, batches arrived incrementally. For methods that originally lacked online support, we added it through CL or fine-tuning-based transfer learning. Below, each method is described with its repository and key settings, with minor offline/online variants grouped under their base method.

Horizontal integration methods
Geneformer

Geneformer27,71 is a large-scale model pretrained on 30 million cells; we averaged its per-gene embeddings to obtain cell embeddings.

PCA

PCA was applied after selecting highly variable genes, retaining the top 32 principal components as embeddings.

Rectangular integration methods
SCALEX

SCALEX13,72 is an online scRNA-/scATAC-seq method; lacking native ADT support, we concatenated RNA and ADT. SCALEX-projection trained on the first batch and projected the remaining batches (embedding dimension 32).

trVAE

trVAE73,74 supports a single modality, so we concatenated RNA and ADT. trVAE+scArches added scArches transfer learning, training a reference on the first batch and adaptors on the remaining batches.

Concerto

Concerto12,75 integrates RNA and ADT. Concerto-transfer pretrained on the first batch and fine-tuned sequentially on the remaining batches.

Symphony

Symphony14,76 handles new-data queries; we concatenated RNA and ADT, building a reference on the first batch and querying the rest. Symphony-offline instead trained on all batches together.

online iNMF

Online iNMF10,77 iteratively updates RNA models; we concatenated RNA and ADT, integrated by batch and used the final model for all embeddings. Online iNMF-offline trained all batches together.

CCA+WNN and RPCA+WNN

Seurat’s CCA6 removes batch effects in feature space (FindIntegrationAnchors with reduction ‘cca’, then IntegrateData), followed by PCA and WNN omics fusion; RPCA+WNN is identical with reduction ‘rpca’. Their ‘-continual’ variants integrate each new batch against previously corrected counts, resembling CL with unbounded rehearsal memory and thus declining in efficiency as integrations accumulate.

sciPENN

sciPENN34,78 integrates RNA and ADT using the intersection of RNA features and the union of ADT features. sciPENN-transfer trained on the initial batch and fine-tuned on subsequent ones.

Mosaic integration methods
MIRACLE

MIRACLE employs MIDAS19,79 as its base model, inheriting its neural network architecture and training hyperparameters. Its rehearsal memory was set to 200K cells for all bimodal and trimodal tasks—substantially above the 20K sufficient for single-modality RNA integration on DHCM—to better preserve biological fidelity in more complex multimodal data.

MIRACLE-offline and MIRACLE-transfer

MIRACLE-offline, the upper benchmark, trains the same base model on all batches simultaneously, avoiding catastrophic forgetting and requiring no rehearsal memory. MIRACLE-transfer, the lower benchmark, fine-tunes the same model on only the new data at each step without rehearsal, making it prone to forgetting.

Multigrate

Multigrate80,81 integrates RNA, ADT and ATAC (KL and integ set to 0.1 and 3000). Multigrate+scArches11,82 added scArches transfer learning, using the first batch (or the first two for the DOTEA mosaic, to cover all modalities) as reference and adaptors for the rest.

scVAEIT

scVAEIT35,83 integrates RNA, ADT and ATAC. scVAEIT-transfer (our extension) trained on the first batch and fine-tuned on subsequent ones.

totalVI

totalVI7 (ref. 84, incorporated into scArches) integrates RNA and ADT. totalVI+scArches built the reference on the first batch and adaptors on the rest.

Evaluation metrics for integration

To comprehensively evaluate the integration performance of MIRACLE and various state-of-the-art methods, we used the scIB and scMIB metrics. For horizontal and rectangular integration tasks, we used the scIB metric to assess batch correction and biological conservation performance. For mosaic integration tasks, we used the scMIB metric to evaluate batch correction, modality alignment and biological conservation performance. Modality alignment is crucial for mosaic integration as it reflects the model’s capabilities in applications such as modality imputation and cross-modal label transfer. Note that, for other mosaic integration methods compared, batch correction can only be performed in the embedding space rather than the feature space, so we excluded feature space-based metrics for batch correction and biological conservation in scMIB.

Specifically, the batch correction metrics include the graph integration local inverse Simpson’s index (iLISI) yiLISI, graph connectivity ygc and kNN batch effect test (kBET) ykBET. The modality alignment metrics include the modality averaged silhouette width (ASW) yASW, fraction of samples closer than the true match (FOSCTTM) yFOSCTTM, label transfer F1 yltF1, ATAC area under the receiver operating characteristic (AUROC) yAUROC, RNA Pearson’s ryRNAr and ADT Pearson’s ryADTr. The biological conservation metrics include the normalized MI (NMI) yNMI, adjusted Rand index (ARI) yARI, isolated label F1 yilF1 and graph cell-type LISI (cLISI) ycLISI. The definitions of these metrics are detailed in Supplementary Information section 6.

We average each type of metrics to obtain its comprehensive score, including the batch correction score ybatch, modality alignment score ymod and biological conservation score ybio

Following ref. 26, the scIB overall score yscIB is the sum of ybatch weighted by 0.4 and ybio weighted by 0.6

As an extension of scIB, following ref. 19, the scMIB overall score yscMIB is the sum of ybatch weighted by 0.3, ymod weighted by 0.3, and ybio weighted by 0.4

CL metrics for catastrophic forgetting

To quantitatively assess MIRACLE’s ability to mitigate catastrophic forgetting, we adopted two classic metrics from the CL literature: average accuracy (ACC) and backward transfer (BWT)85. These metrics are established benchmarks in the field for evaluating a model’s retention of previously learned knowledge. While originally designed for classification tasks, we have adapted them for the context of single-cell data integration.

Let T be the total number of integration steps (that is, the total number of batches to be integrated sequentially). We define yt,τ as the overall integration performance score (for example, the scIB or scMIB overall score) on the cumulative data from the first τ batches, \({{\mathcal{X}}}_{1:\tau }=\mathop{\bigcup }\nolimits_{i = 1}^{\tau }{{\mathcal{X}}}_{i}\), evaluated using the model parameters θt obtained after integrating up to batch t (where 1 ≤ τ ≤ t ≤ T).

ACC

ACC measures the average performance across all cumulative tasks at the end of the entire CL process, providing a summary of how well the final model performs. It is defined as

In this formula, yT,τ denotes the integration performance on the first τ batches as evaluated by the final model, θT, which has been trained on all T batches. A higher ACC indicates better overall performance of the final model across all historical integration stages.

BWT

BWT directly quantifies the influence of learning new tasks on the performance of previously learned tasks. A negative BWT is a hallmark of catastrophic forgetting, indicating that the model’s performance on older tasks has degraded after learning new ones. BWT is calculated as

Here, yτ,τ represents the integration performance on the first τ batches evaluated immediately after the τth batch was integrated (that is, using model θτ). This score can be considered the ‘peak’ performance for that cumulative task. The difference, yT,τ − yτ,τ, therefore directly measures the performance drop (or gain) on an earlier task set after subsequent learning has occurred. A BWT value close to zero signifies successful knowledge retention, while a large negative value indicates substantial catastrophic forgetting.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.