Multi-species integration, alignment and annotation of single-cell RNA-seq data with CAMEX

Domain: chemistry
Published: 2026-02-21
Editorial outcome: major_revision

Abstract
Single-cell RNA-seq (scRNA-seq) data from multiple species present remarkable opportunities to explore cellular origins and evolution. However, integrating and annotating scRNA-seq data across different species remains challenging due to the variations in sequencing techniques, ambiguity of homologous relationships, and limited biological knowledge. To tackle the above challenges, we introduce CAMEX, a heterogeneous Graph Neural Network (GNN) tool that leverages many-to-many homologous relationships for multi-species integration, alignment, and annotation of scRNA-seq data from multiple species. Notably, CAMEX outperforms state-of-the-art methods integration on various cross-species benchmarking datasets (ranging from one to eleven species). Besides, CAMEX facilitates the alignment of diverse species across different developmental stages, significantly enhancing our understanding of organ and organism origins. Furthermore, CAMEX enables the detection of species-specific cell types and marker genes through cell and gene embedding. In short, CAMEX holds the potential to provide invaluable insights into how evolutionary forces operate across different species at single-cell resolution.

Introduction
Single-cell omics techniques enable the profiling of thousands of cells in one experiment, revealing a comprehensive landscape of heterogeneous cells1. These techniques facilitate the discovery of rare cell states, infer gene regulatory networks, and enhance our understanding of development and differentiation processes at a finer resolution2. With advancements in sequencing methods, scRNA-seq data from different species, including humans3, monkeys4, mice5, zebrafish6, and more, have been generated by various studies. These datasets provide clues on the origin of cellular diversity and evolutionary drivers of cellular morphology7. Diverse studies have revealed many cryptic cell types across multiple species8. Comparing gene expression of individual cells across species uncovers species-specific transcriptomic regulation9, provides novel insights into the phylogeny of cells and organs, and elucidates how cells function in different developmental stages10. Integrating scRNA-seq data across species advances our understanding of molecular and cellular biology and becomes an essential first step toward comparative biology11,12. These scRNA-seq data enable cross-species comparisons and analyses, shedding light on cellular functions and evolutionary dynamics. Only through these comparisons can we eventually decipher and reconstruct species evolution13.

Single-cell RNA sequencing experiments can only process a limited number of cells every time. Data from various batches exhibit significant variations due to differences in capture time, operating personnel, and sequencing platforms14. For instance, batch effects make it challenging to align the developmental stages of different species, hinder cross-species cell annotation and the identification of novel cell types and subtypes. Overall, these batch effects can obscure biological signals and interfere with our ability to analyze specific biological changes of interest15. Researchers have developed many methods, such as Harmony16, scVI17, and SingleCellNet (SCN)18, for integrating and annotating scRNA-seq data within a single species. Furthermore, inspired by the rapid advancement of large models in natural language processing, analogous foundation models have emerged in the single-cell field, such as scGPT19, which facilitates downstream tasks like batch effect removal. However, existing methods align features in the expression profile matrix across different species solely through one-to-one homologous genes, resulting in the loss of much biological information20.

Several methods have made initial explorations into the integration and annotation of cross-species scRNA-seq data, such as BPA (Biological Process Activity)21, SAMap22, SATURN23 and CAME24. However, the limited biological knowledge available for non-model organisms poses challenges for implementing BPA. SAMAP relies on sequence alignment, a computationally intensive process, and is primarily used for visualization purposes rather than generating embeddings for downstream tasks19. SATURN measures inter-species gene similarity using protein-language models, but faces challenges from proteome-transcriptome mismatches (notably from alternative splicing events) and requires cell annotation inputs that might limit its application. While our previous annotation method, CAME, utilizes many-to-many homologous gene relationships, it is restricted to pairwise species alignment (two species) within a shared embedding space and cannot perform integration without cell labels across multiple species datasets.

In a word, integrating scRNA-seq data obtained from diverse sources and non-model organisms faces challenges due to technical and biological factors. First, the data from different sequencing platforms vary greatly, including the difference in detected UMIs and genes. Second, aligning the features of expression profile matrices from other species through the relationships between orthologous genes poses a significant difficulty. Last, processing scRNA-seq data from non-model organisms with limited biological knowledge, e.g., specific regulatory networks, distinctive communication mechanisms, and individual diversity, is a complex challenge.

To this end, we introduce CAMEX, a heterogeneous graph neural network (GNN) tool that leverages many-to-many homologous relationships for integration, alignment, and annotation of scRNA-seq data from multiple species. Specifically, we construct a heterogeneous graph consisting of cells and genes based on expression matrices and many-to-many homology relationships. We design a heterogeneous graph neural network as an encoder to encode cells and genes from different species into the common embedding space. For the batch correction task, we design two decoders, including a dot-product decoder to reconstruct the edges in the heterogeneous graph and a fully connected decoder to reconstruct the cell expression. For the cell annotation task, we feed the cell embedding into the graph attention network (GAT) for classification. CAMEX achieved competitive performance compared to state-of-the-art (SOTA) methods in integrating scRNA-seq data in different scenarios, including the liver (involving four species with two technologies), ovary (encompassing three species with rare populations) and pancreas (comprising three species with multiple batches). Furthermore, CAMEX successfully integrated a testis dataset across 11 different species, preserving the sperm differentiation trajectory. Additionally, CAMEX correctly integrated RNA-seq data from seven species, spanning multiple organs and various developmental stages. Notably, CAMEX could align and predict developmental time beyond the Carnegie stages across different species. Finally, we demonstrated that CAMEX effectively annotates multi-species cell labels, even for distantly related species, and can discover new subgroups and marker genes.

Results
Overview of CAMEX
CAMEX takes diverse single-cell datasets from various species generated using the same or different technologies as input (Fig. 1). The expression matrix can be transformed into a bipartite graph A, with cells and genes as vertices. An edge A is established between a cell n and a gene m if the cell n expresses the gene m. Then, we construct a k-nearest neighbor (kNN) graph to connect cells according to their expression similarity. In addition, self-loops are added for each cell and gene to prevent excessive smoothing in the GNN’s message passing. Finally, we connect the genes of different species based on many-to-many homologous relationships. The initial feature of each cell node corresponds to its preprocessed expression, while the initial feature of each gene node is set to empty. CAMEX employs a heterogeneous GNN encoder to embed each node into a low-dimensional common space by nonlinearly propagating features from neighboring nodes to center nodes.

As to the integration task, CAMEX utilizes an encoder-decoder architecture with no cell type labels in any dataset. Specifically, it employs a heterogeneous GNN as an encoder, a dot product decoder to reconstruct the edges in the cell-gene heterogeneous graph, and a fully connected decoder to reconstruct the cell expression. The loss function comprises both cross-entropy loss between the predicted and raw labels of edges and mean square error (MSE) loss between the predicted and input expression of cells. CAMEX minimizes them to optimize the model parameters by the back-propagation algorithm in an unsupervised manner.

As to the annotation task, CAMEX adopts an encoder-classifier architecture with at least one given reference dataset containing cell labels. It employs a heterogeneous GNN as an encoder and a GAT to predict cell labels. The loss function computes the cross-entropy loss between the predicted and ground truth labels of cell annotations in the reference dataset. CAMEX minimizes it to optimize the model parameters by the back-propagation algorithm in a semi-supervised way.

We collected various cross-species scRNA-seq datasets to benchmark CAMEX and other baseline methods. In the integration task, we applied CAMEX to the datasets and compared its performance against other SOTA methods, i.e., Harmony16, PyLiger25,26, scVI17, Scanorama27, scGEN28, scGPT19, Seurat29, SAMap22, and SATURN30. For SAMap and SATURN, we conducted our experiments following their respective tutorials (Methods). For CAMEX, we leveraged many-to-many homologous gene relationships to construct heterogeneous networks. For the other methods, we employed one-to-one homologous genes to align expression matrices for the baseline methods. We conducted a quantitative evaluation of batch correction and biological variation conservation. Specifically, the batch correction score using metrics including kBET (k-nearest-neighbor Batch Effect Test)31, PCR batch (Principal Component Regression of the batch covariate)31, ASW batch (Average Silhouette Width across batch)31, and graph connectivity14. The biological variation conservation score using metrics including ARI (Adjusted Rand Index)32, NMI (Normalized Mutual Information)33, Cell-type ASW14, ISO label F1 (Isolated label F1)14, and ISO label silhouette (Isolated label silhouette)14. Higher values for each metric indicate better performance (Methods). Finally, we computed an overall score as a 1:1 weighted average of the batch correction score and the biological variation conservation score to indicate the overall model performance. Notably, some metrics for SAMap were not calculable. Although we computed its Batch Correction, Biological Variation Conservation, and Overall Average Score, these metrics failed to reflect its ranking precisely. This finding is confirmed by the prior research19, which has similarly pointed out the inapplicability of these metrics to SAMap. Furthermore, we evaluated the model’s performance in retaining the developmental trajectory by calculating pseudo-time both before and after integration (Methods).

Similar to the integration task, here we applied CAMEX to the datasets and compared its performance against other baseline methods, i.e., SCN18, scANVI34, and CellBLAST35. We conducted ACC (accuracy), Prec (precision), and F1 (F1 score) for the quantitative evaluation of each baseline method.

CAMEX achieves superior integration performance in different multi-species scenarios
We curated three multi-species scRNA-seq datasets, comprising the liver data (involving four species with two technologies), ovary data (encompassing three species with rare populations), and pancreas data (comprising three species with seven batches by six different technologies), to simulate a diverse range of real-world scenarios for benchmarking CAMEX against other baseline methods. In this unsupervised integration setting, CAMEX operates without the need for manually annotated cell type labels.

In the liver dataset, the scRNA-seq data for humans, macaques, and mice originating from the same study were generated using the same sequencing platform36. However, the zebrafish scRNA-seq data were obtained from the zebrafish atlas13, where the detected number of Unique Molecular Identifiers (UMIs) and genes was relatively low (Supplementary Fig. S1a, b). Each of the datasets contains distinct cell types unique to that respective species. While the human, monkey, and mouse datasets exhibit balanced cell types, the zebrafish dataset predominantly consists of hepatocytes (Supplementary Fig. S1c). These two factors significantly increased the complexity of data integration. In this context, CAMEX was the only method capable of accurately integrating scRNA-seq data across different species and achieving the highest overall score (Fig. 2a–c, Supplementary Fig. S1d, Supplementary Data 1). It effectively clustered hepatocytes, fibroblasts, cholangiocytes, endothelial cells, B plasma cells, and Kupffer cells (macrophages). Compared to other methods, CAMEX achieved higher Adjusted Rand Index (ARI) and Normalized Mutual Information (NMI) scores, reflecting its superior ability to cluster cells of the same type. Although scVI also aggregated similar cell types across different species, its large intra-class distances led to lower kBET and ASW batch scores, positioning it below CAMEX for the overall ranking. PyLiger performed well in integrating hepatocytes, but it cannot address other cell types well. Harmony excelled in integrating non-hepatocytes, but failed when dealing with hepatocytes. The integration result of Scanorama was similar to Harmony, but its overall performance lagged behind Harmony. Notably, scGen could integrate hepatocytes from four species while failing to integrate other cell types. In contrast, SATURN, Seurat v4, Harmony, scGPT, Scanorama, and SAMap achieved integration for all cell types except hepatocytes. In addition, all methods successfully integrated liver-specific macrophages (macrophages and Kupffer cells). We examined the expression of typical marker genes in each cell type and verified the identities of macrophages and Kupffer cells (Supplementary Fig. S1e). The genes from each cell type vary considerably across species (Supplementary Fig. S1f). We also performed computational efficiency, hyperparameter comparison, and ablation experiments (Supplementary Figs. S2–S4 and Supplementary Data 2–4). Despite variations in evaluation metrics, most settings successfully isolated distinct cell types in the UMAP. An interesting observation was the clear separation of Kupffer cells and Macrophages when the hidden layer dimension was increased to 512 (Supplementary Fig. S2c).

In the ovary dataset, the scRNA-seq data for humans37–42, macaques43, and mice44 originating from different sequencing platforms resulted in distinct data types (Supplementary Fig. S5a-c). The macaque dataset contains a rare, unique oocyte cell type. In this scenario, CAMEX, scGen, SATURN, Seurat v4, PyLiger, scVI, and scGPT can integrate the scRNA-seq data from three species while preserving this rare cell population, but CAMEX achieved a higher overall score (Supplementary Fig. S5d, e, Supplementary Fig. S6, Supplementary Data 5 and 6).

In the pancreas dataset, the scRNA-seq data came from humans39, macaques45, and mice39, with humans containing five batches (by five different technologies) of data (Supplementary Fig. S7a). We used this to simulate a multi-species, multi-batch scenario. The human batches in this dataset contain roughly the same cell types. In addition, mice and macaques contain cell types that are subsets of humans (Supplementary Fig. S7b, c). In this scenario, Seurat v4, SATURN and CAMEX can preserve the strong biological signal and effectively separate different cell types (Supplementary Fig. S7d, e, Supplementary Fig. S8, Supplementary Data 7, 8). Almost all methods, excluding scGPT, Scanorama, and SAMap, exhibited relatively decent integration performance and can clearly distinguish various cell types. However, Harmony tended to confuse the acinar and ductal cell types. PyLiger confused the acinar and ductal cells and the delta and gamma cells but failed to integrate the alpha and beta cell types. Conversely, scVI cannot distinguish acinar and ductal cells well and faces challenges when integrating the gamma cells.

Overall, CAMEX has demonstrated competitive and robust performance across various multi-species scenarios. It can effectively handle substantial differences and multi-batch data and identify and preserve rare cell populations. Nevertheless, other baseline methods exhibit higher variability in performance, resulting in less stable outcomes.

CAMEX uncovers the conserved differentiation process in the testis across 11 species
The testis is a unique organ with the highest number of tissue-specific genes in its transcriptome, and sperm development represents a continuous and conserved process46. Consequently, comparing these continuous and critical biological processes across different species is essential for evolutionary biology. Here, we collected and utilized the scRNA-seq testis data from 11 species, encompassing five distinct cell types, i.e., somatic, Sertoli, spermatogonia, spermatocytes, and spermatids47. Although this dataset originates from the same source and exhibits consistent cell types (Supplementary Fig. S9a), integrating it across 11 species remains a significant challenge due to the scarcity of one-to-one homologous genes shared among them (Supplementary Fig. S9b). Worth noting is that the original research only merges scRNA-seq data from seven primate datasets via LIGER (v.0.5.0)26 based on primate one-to-one orthologues from Ensembl release 8747.

Firstly, we integrated the testis scRNA-seq data across 11 species using CAMEX and other baseline methods. In this unsupervised integration setting, both CAMEX and the baseline methods are precluded from using any manually annotated cell type labels. We found that SATURN and SAMap could not run normally due to prohibitive memory requirements (GPU memory for SATURN, CPU memory for SAMap), which led us to exclude them from our analysis. From the UMAP, we can see that CAMEX successfully preserved the continuous trajectory of sperm differentiation and separated distinct clusters corresponding to different cell types (Fig. 3a, b and Supplementary Data 9, 10). Furthermore, we calculated the pseudo-times both before and after integration, using their correlation as the trajectory reservation metric. CAMEX achieved the highest value in this regard (Supplementary Fig. S10). Analyzing the results of CAMEX, we observed that the Platypus dataset almost entirely lacked Sertoli cells, consistent with the original findings (Supplementary Fig. S11). PyLiger partially retained the developmental trajectory but failed to integrate spermatids and spermatocytes, erroneously merging somatic cells with spermatids. scVI struggled to effectively integrate Sertoli and spermatocyte cells from chickens with those from other species. Harmony failed to integrate Sertoli cells from opossums and spermatocytes from chickens. Scanorama achieved the worst integration results.

Inspired by previous literature48, to comprehensively evaluate CAMEX’s performance across varying evolutionary distances and increasing numbers of species, we designed two integration scenarios. First, we conducted pairwise integration analyses between humans and each target species, covering a broad phylogenetic spectrum from closely related primates to distant vertebrates. The results clearly show CAMEX excelled in batch correction, while SATURN performed strongly in biological variation conservation. Overall, CAMEX demonstrated an advantage over other methods in both quantitative metrics, rankings, and qualitative UMAP visualization. Second, we created multi-species integration scenarios with progressively expanded species sets, systematically increasing from two to eleven species. CAMEX also consistently maintained high performance across all scenarios (Supplementary Figs. S12–S24 and Supplementary Data 11–21). Overall, most integration algorithms, particularly those relying solely on one-to-one homologous genes, show diminished performance as the evolutionary distance between species increases. Additionally, we conducted an ablation experiment in which human data were integrated separately with each of the other ten species using only one-to-one homologous genes. As shown in Supplementary Fig. S24b, integration performance remained comparable to that achieved with many-to-many homologous genes when closely related primates were involved. However, as evolutionary distance increased—particularly in the human–chicken comparison—a marked decline in performance was observed. These findings suggest that the use of many-to-many homologous genes is a critical factor enabling CAMEX to maintain robust integration performance across large evolutionary spans. Overall, methods like CAMEX, which leverage many-to-many homologous genes, and STAURN, which utilizes sequence similarity, are less impacted by this increasing species divergence.

Secondly, given that the cell subtypes from seven primate species have been manually annotated, we leveraged the cell embeddings obtained from CAMEX to assess the quality of cell embedding in the CAMEX integration task. Using the human dataset as a reference, we applied a simple k-Nearest Neighbors (kNN) classifier to annotate the subtypes for other primate cells. The predicted results are almost consistent with the manual annotations (normalized confusion matrix) (Supplementary Fig. S25). Additionally, we annotated the cell subtypes of four other non-primate species. Cell subtypes exhibit similarities with the annotated major population (Supplementary Fig. S26a). CAMEX identified some cells in the platypus dataset as Sertoli cells. Through further differential analysis, we identified differentially expressed genes (DEGs) (Supplementary Fig. S26b). These genes are highly relevant to spermatogenesis in the testis. For example, the CYP17A1 gene codes for a protein involved in testosterone synthesis and plays a vital role in cholesterol transport49. Similarly, the STAR gene encodes a protein critical for cholesterol transport across mitochondrial membranes, thereby regulating testosterone synthesis50. We speculate that the labeled cells may represent a rare population of potential platypus-specific Sertoli cells.

Thirdly, we obtained the mean embedding for each cell type across all species and conducted hierarchical clustering (Supplementary Fig. S27). The clustering outcomes indirectly confirm the accuracy of cell integration. Notably, opossum and mouse spermatogonia were mistakenly grouped with spermatocytes. We attribute this misclassification to the relatively young age of samples in the original dataset47. Interestingly, spermatogonia appear more similar to somatic cells than spermatocytes and spermatids. This finding suggests that the transcriptional differences during spermatogenesis are crucial47,51. At the cell type level, we established a species expression hierarchy through clustering, although it remains incomplete. This exploration lays the groundwork for further research in this direction.

Finally, we acquired cell and gene embeddings via CAMEX. By clustering the gene embeddings, we identified specific functional modules within clusters 3 and 5 (Fig. 3c, d). Genes in cluster 3 are predominantly expressed in somatic cells. Enrichment analysis revealed associations with the cell surface receptor signaling pathway and extracellular matrix. Somatic cells are all the cells that make up the body of an organism, except for reproductive cells. In the male reproductive organs, somatic cells play multiple roles such as supporting, protecting, secreting and regulating spermatogenesis52. Conversely, genes in cluster 5 are primarily expressed in Spermatocytes, and enrichment analysis linked them to sexual reproduction, meiotic cell cycle process, meiotic cell cycle, meiotic nuclear division, and reproductive process. Meiosis is a well-known process occurring in spermatocytes53 and is consistent with gene function obtained by CAMEX. In summary, CAMEX can integrate conserved developmental processes across multiple and distant species and provide gene embeddings that shed light on cell functions.

CAMEX aligns various development stages of seven organs across seven different species
The evolution of gene expression during organ development of various species remains largely unknown, and different species exhibit distinct patterns. We applied CAMEX and competing methods to integrate and align the RNA-seq data from seven organs across seven species at various developmental stages54. In this unsupervised integration setting, CAMEX does not require any manually annotated cell type labels. We put the bulk data for all developmental stages of all species into CAMEX and other methods to perform the integration task and obtain the embedding together (Fig. 4a). Then, we split and visualize the samples of different organs respectively. Remarkably, CAMEX achieved the highest overall score, and scGen, Harmony, and CAMEX can preserve the continuous development trajectory and separate distinct clusters corresponding to different organs (Fig. 4a, Supplementary Figs. S28, S29, Supplementary Data 22). We noted that the cerebellum and brain were well integrated due to their abundance of neurons and glial cells55, while the testis and ovary, being gonadal cells in early development56, exhibited higher similarity.

While most model organisms display distinct developmental timelines, these are typically mapped to Carnegie stages, spanning from gametogenesis through birth. The Carnegie stages themselves are fundamentally based on observable morphological development during embryogenesis. Crucially, they’ve only been defined by researchers for a handful of model organisms (e.g., human and mouse) and are limited to the period immediately following pregnancy and childbirth. An exciting challenge lies in exploring whether we can expand the Carnegie developmental staging concept beyond this limited set of organisms to encompass any species, or even to predict any point in a species’ lifespan through integrated approaches.

After obtaining organ embeddings from the integration results of CAMEX, we calculated the similarity of these organ embeddings, allowing us to establish correspondence across diverse species. Initially, we selected rats and mice as our subjects and sampled from embryonic day 11 (e11) to postnatal day 112 (p112) for rats and from e10 to p63 for mice. CAMEX aligned the developmental stages of mouse and rat (9/10.5–16/17.5) within the Carnegie stages, consistent with existing literature, and even predicted their stages beyond the specified ranges (Fig. 4b, Supplementary Data 23). Additionally, we adopted CAMEX to analyze the brain transcriptome data from mouse (Carnegie stage 9–16 days) and opossum (Carnegie stage unknown) at multiple developmental stages. While the mouse brain data spanned from e10 to p63, the opossum brain data covered e13.5 to p180. Notably, the opossum central nervous system, including the brain, is markedly premature at birth (p0), resembling an embryonic e11 mouse, as confirmed by the similarity of development time calculated by CAMEX (Supplementary Fig. S30) and a previous study57. However, the developmental stages across different species in the results are not perfectly aligned. We believe there are several contributing factors. Biological Considerations: From a biological perspective, there are two main reasons. First, the developmental speeds across different species might inherently vary. Second, while we align species based on transcriptome expression, the Carnegie developmental stages are primarily defined by morphological development during embryogenesis, which can lead to discrepancies. Model and Computational Considerations: From a modeling and computational standpoint, the observed differences can also stem from the randomness in the neural network optimization process and the various similarity metrics used for embedding calculations. Despite these points, CAMEX still generally achieves roughly alignment of developmental processes across species. For time points that are not well aligned, this can serve as a potential objective for studying the developmental differences across species.

CAMEX could achieve more accurate integration and annotation performance in both relatives and distant species
Integrating and annotating distant species presents a challenging task due to the tendency to share fewer orthologous genes, resulting in more specific cell-type markers. Here, we collected four scRNA-seq datasets from four species, i.e., adult human visual cortex, frontal cortex, cerebellum55, mouse neocortex58, and lizard and turtle pallium59, which were profiled from three technology platforms. Each dataset has a specific cell type that is absent in the datasets of other species (Supplementary Fig. S31a-c).

First, we applied CAMEX and baseline methods to integrate the scRNA-seq data in the absence of cell type annotations. In this setting, all the methods except Scanorama achieved favorable integration performance, with CAMEX receiving the highest overall score (Fig. 5a, Supplementary Fig. S31d, e, Supplementary Data 24 and 25). From the UMAP plot (Supplementary Fig. S32), we can see that all the methods confused the excitatory and inhibitory neurons. CAMEX failed to integrate inhibitory neurons, Harmony and scVI complicated specifically with microglia cells, and scGen failed to mix astrocytes adequately. PyLiger did not properly integrate astrocytes and confused NPCs with excitatory neurons. SATURN misclassified CGCs as NPCs, and scGPT confused ExN with NPCs, in addition to failing to integrate oligodendrocytes.

Traditional cell-type annotation methods heavily depend on clustering and marker genes, posing a significant challenge in non-model species where such prior knowledge is insufficient. To this end, we introduce a GAT classifier after the encoder (Methods). It alleviates the issue by leveraging many-to-many homology relationships, thereby improving annotation accuracy. In this semi-supervised scenario, the human dataset served as the labeled reference. We then applied CAMEX and other methods to annotate the other three unlabeled target datasets, with a focus on the lizard and turtle species. CAMEX and scANVI demonstrated better performance compared to the other two baseline methods, as evaluated by metrics including ACC, Prec, and F1 score (Fig. 5b, c and Supplementary Fig. S33a–f). Unknown is the cell type that does not appear in the reference but appears in the query dataset. NPC does not appear in the reference human dataset but does appear in the query turtle and lizard dataset, so we have defined them as unknown. Due to the absence of neural progenitor cells (NPCs) in the human reference dataset, both CAMEX and scANVI cannot identify them. CAMEX misclassified them as excitatory neurons, whereas scANVI identified them as a mixture of CGCs and other cell types. Interestingly, both predictions offer a degree of reasonableness (Supplementary Fig. S34).

CAMEX facilitates the discovery of new populations and markers in the primate granular dorsolateral prefrontal cortex (DLPFC) data
The discovery of new cell populations and gene markers holds immense importance in understanding evolutionary organisms and identifying novel targets. The granular dorsolateral prefrontal cortex, unique to primates in evolution, plays a pivotal role in cognitive processes. We investigated specific cell types and gene markers within the collected DLPFC dataset across four species, i.e., humans, rhesuses, marmosets, and chimpanzees60. The cell class, subclass, and subtype are manually annotated by clustering and marker genes in the original study60.

Firstly, we integrated all the cells in the DLPFC dataset among the four species by CAMEX in an unsupervised way. Almost all the cell types match accurately, and the major classes (e.g., excitatory neurons, inhibitory neurons, non-neural cells, and glial cells) of cell types are relatively conserved, indicating the capability of CAMEX for cross-species integration (Fig. 6a and Supplementary Fig. S35a).

Then, we simulated the following scenario, assuming we were familiar with one species or model organism. We could adopt its comprehensive reference scRNA-seq atlas to transfer its annotations to other species and reveal several species-specific cell populations relating to species divergence. Considering that during the process of nervous system evolution from a diffuse to a centralized form, glial cells and neurons co-evolved, glial cells, in particular, diversified into a more varied family of cell than the latter61. Here, we use human cells as the reference (training set) to annotate cells in other species (test set) in a semi-supervised way by using a GAT classifier (Methods). Glial cells have four subclasses: astrocytes, oligodendrocytes, oligodendrocyte precursor cells (OPC), and microglia cells. Early studies indicated that OPC and oligodendrocytes among four species shared all subtypes, and the annotation embedding also confirmed that (Supplementary Fig. S35b). As to the astrocytes, CAMEX did not isolate the AQP4 OSMR subtype that exists in humans and chimpanzees (Supplementary Fig. S35b). This might be due to the low proportion of this subtype in the training set (<0.25%).

In addition, microglia are more diverse, and those microglial cells are divided into three distinct cell subtypes in the original study, i.e., APBB1IP, CCL3, and GLDN clusters. Specifically, CCL3 is exclusive to humans, GLDN is found in humans and chimpanzees, and APBB1IP appeared across four species (Fig. 6a left). Leveraging CAMEX with a GAT classifier, we mapped the other three species queries into the human reference dataset, resulting in the distinct isolation of three subtypes (Fig. 6b and Supplementary Fig. S36a, b). At the same time, we investigated the expression of GLDN, CCL3, and APBB1IP in the original manual annotation and the CAMEX prediction (Supplementary Fig. S36c). We found that the prediction of CAMEX has certain rationality, for instance, the expression of CCL3 in the subgroups of chimpanzees. However, there are also false positives, such as CCL3 and GLDN in marmosets. Overall, we therefore wish to clarify that CAMEX should be considered a bioinformatic tool for generating hypotheses about potential cell subpopulations. These computational predictions need further validation through subsequent wet-lab experiments.

Moreover, we conducted differential analysis for each subtype, comparing them with other clusters to identify subtype-specific genes (Supplementary Fig. S37). The differential genes associated with the CCL3 subtype primarily function in immunity. For instance, CD83 encodes a protein known as CD83, which stabilizes MHC II, costimulatory molecules, and CD28 in the cell membrane by counteracting E3 ubiquitin ligases62,63. On the other hand, the differential genes within the GLDN subtype are primarily involved in signal transduction. SPP1 (secreted phosphoprotein 1) is an extracellular protein closely linked to tumor biology, including processes like proliferation, migration, and invasion64,65.

In addition, even within the APBB1IP subtype, we observed abundant differential genes between species (Fig. 6c). NAV2 and NAV3 belong to the neuronal navigator family and are predominantly expressed in the nervous system. They likely play roles in axon guidance, cell migration, and neurodevelopment66. PAM is involved in the amidation process of peptide hormones and neurotransmitters and is crucial for regulating their activity and stability67. Overall, the human and chimpanzee transcriptomes exhibit a high consistency at the single-cell resolution and can serve as model organisms for exploring disease mechanisms or drug development. However, focusing on more species-specific genes is essential when studying the marmoset and the rhesus.

Methods
CAMEX model
Heterogeneous graph of cells and genes. A scRNA-seq dataset comprising N cells and M genes is represented as a matrix X = [x1, x2, ..., xN]T ∈ RN×M. Following a standard preprocessing pipeline (normalization, log-transformation, and selection of highly variable genes, clustering, and obtaining differential genes), we construct a cell-gene bipartite graph A. In this graph, an edge connects cell n and gene m if and only if cell n expresses gene m. Subsequently, cell-cell connections are established by applying the kNN algorithm within the PCA space derived from this processed matrix. In addition, we added a self-loop to each cell and gene, respectively. Then, we connected genes in each A from different species through the many-to-many homologous gene relationships.

In short, there are two types of nodes including cells and genes, and six types of edges (relations) including ‘a cell expresses a gene’, ‘a gene is expressed by a cell’, ‘two cells that are kNN neighbors’, ‘self-loop of a cell’, ‘self-loop of a gene’ and ‘two genes that are orthologous ones’ in the heterogeneous graph. We initialized cell node features with preprocessed expression data (normalized, log-transformed values from the top 2000 highly variable genes (HVGs)), while gene nodes were initialized as zero vectors. The resulting heterogeneous graph incorporated both cell and gene nodes with their respective features.

Heterogeneous graph neural network encoder. CAMEX adopts a heterogeneous graph neural network as the encoder, inspired by the relational graph convolutional network72. We design six independent parameter matrices for each of the above relationships as follows. W (a cell expresses a gene), W (a gene is expressed by a cell), W (two cells that are kNN neighbors), W (two genes that are orthologous ones), W (self-loop of a cell), and W (self-loop of a gene). The encoder transforms the features of the neighbors in a non-linear way and converges to the center node according to the message-passing paradigm. The (l+1)-th layer representation of cell h(l+1) and gene h(l+1) can be defined as follows:

h(l+1) = σ(Σr∈R Σj∈Nr W(l)h(l) + W(l)h(l)) (1)
ni c r j c ni

h(l+1) = σ(Σr∈R Σj∈Nr W(l)h(l) + W(l)h(l)) (2)
mi r j c i,r r j g mi

where h(l) and h(l) represent the l-th layer representation of cells and genes, respectively, W(l) and W(l) represent the l-th layer weight matrices of ‘self-loop of a cell’ and ‘self-loop of a gene’, respectively, Nr denotes the set of neighbor indices of the node (cell or gene) i under relations r ∈ R, c is a normalization factors, W(l) is the weight matrix of the l-th layer transformation of the relationship r, h(l) is the l-th layer representation of neighbor j, and σ is the activation function (i.e., LeakyReLU here).

To reduce the number of parameters to prevent overfitting, we add a hyperparameter ‘share’ that can share weight matrices of the same relationship in different layers, such as W(l) = W(l+1). After the message passes through several layers of neural networks, the cell and gene representations can be mapped into the same embedding space.

Integration task
Edge and feature reconstruction decoder. When integrating different scRNA-seq data in an unsupervised manner, we employ two decoders to reconstruct edges and features, respectively. Specifically, the dot-product decoder aims to reconstruct the edges of the input heterogeneous graph. Edge predictions in the heterogeneous graph are generated using the dot product of cell and gene embeddings:

A0 = H H T (3)
N M

where H and H are the cell and gene embeddings, respectively, and A0 is the reconstructed adjacency matrix.

We measure the quality of each edge in the adjacency matrix using cross-entropy loss:

L = -ΣN ΣM [A logA^ + (1-A )log(1-A^ )] (4)
A i=1 j=1 ij ij ij ij

where A is the edge between the i-th cell and the j-th gene in the original graph, where A^ is obtained by applying a sigmoid function to A0.

The feature reconstruction decoder aims to reconstruct cell features through a fully connected layer:

h0 = W (h(l)) (5)
0 fc ni

where W is the weight of a fully connected neural network and h0 is the reconstructed cell feature.

We utilize the mean squared error to measure the distance between the original h and the predicted h0:

L = 1 ΣN ||h - h0 ||2 (6)
F N i=1 0,i 0,i

Loss function and optimizer in the integration task. The unsupervised integration task loss function can be defined as follows:

L = L + L (7)
integration A F

where L and L are optimized jointly.

The integration task utilizes a graph attention network autoencoder with an encoder-decoder architecture, which is trained end-to-end through backpropagation73. The model optimization employs the Adam optimizer74, a variant of the Stochastic Gradient Descent algorithm, for parameter optimization and is implemented in PyTorch75 using a mini-batch training strategy. We use an NVIDIA GPU laptop RTX 3070 to accelerate the training process. By default, both the encoder and decoder architectures contain two bottleneck layers. We applied a hyperparameter ‘share’ to allow these layers to share parameters, respectively, thereby reducing GPU memory consumption.

Annotation task
Graph attention classifier. Cell labels are not mandatory during the integration task. However, during annotation, it is essential to have one or more reference datasets with pre-defined manual annotations. Consequently, the model could be trained in a semi-supervised manner. The cell annotations provided by the reference dataset can be denoted as: Y = [y1, y2, ..., yNr]T ∈ RC, where Nr represents the number of cells in the reference dataset and C represents the number of cell types. When performing cell annotation tasks, cell embeddings can be input into a GAT classifier. Specifically, the attention coefficients are first computed as follows:

e = a [Wh(l) || Wh(l)] (8)
ij i j

where a is a learnable linear transformation mapping the concatenated features to a scalar attention score, h(l) and h(l) are the l-th layer representations of cells i and j.

Shared parameter W maps cell embedding to a common space.

The concatenation operation (denoted by ||) is applied. Finally, the function a maps the concatenated result to a scalar representing the relevance between nodes i and j. Subsequently, the relevance scores are normalized to obtain attention coefficients:

α = exp(LeakyReLU(e )) / Σ k∈Ni exp(LeakyReLU(e )) (9)
ij ij ik

We use the multi-head self-attention mechanism to predict the final results:

y^ = K 1 ΣK k=1 Σ j∈Nr α(k) W(l) h(l) (10)
i ij r j

where K is the total number of attention-heads, set as eight by default, we employ multi-head attention with K heads (averaging). For heterogeneous relations, the outputs from different relation-specific GAT convolutions are aggregated by mean.

Loss function and optimizer in the annotation task. The semi-supervised annotation task loss function is a cross-entropy loss function that can be defined as follows:

L = - 1 ΣN ΣC y log(y^) (11)
annotation N i=1 c=1 i,c i,c

where N represents sample number, C represents the class number, y represents the true label of sample i on category c, y^ represents the model predicts the probability that sample i belongs to category c.

The optimizer and implementation details are the same as the integration task in ‘Loss function and optimizer in the integration task’.

In a word, integration and annotation tasks share the same encoder architecture. During integration, unsupervised training is achieved without requiring cell type labels, while annotation necessitates a reference dataset with cell type labels for semi-supervised training.

Training details
Graph sampling. For large-scale data with a large number of cells or involving more species, we partition the heterogeneous graph into multiple subgraphs and train the graph neural network using a mini-batch approach. We implement this using the Deep Graph Library (DGL)76.

Checkpoint selection. In deep learning model training, checkpoints are commonly employed to prevent overfitting. However, in the context of unsupervised integration and semi-supervised annotation tasks, these checkpoints can sometimes lead to over-correction or over-fitting. To address this, we introduce a novel metric for evaluating checkpoints. Specifically, we cluster cells based on their embeddings and assess whether to halt training or revert to a checkpoint using the clustering metric ARI.

Graph storage. To efficiently store large-scale graphs, we implemented directed heterogeneous attribute graphs and graph neural networks using the framework DGL library76. The DGL utilizes the coordinate list matrix (COO) matrix for sparse data representation. Its storage logic follows a directed graph structure, mapping relationships from source nodes to target nodes. To maintain consistency with the COO matrix format, we save six kinds of edges independently and opted not to merge the W gc and W cg.

Hyperparameters. CAMEX employs a set of default hyperparameters, including: train_mode: ‘mini_batch’, dim_hidden: 128, batch_size: 1024, gnn_layer_num: 2, epoch_integration: 10, epoch_annotation: 10, and encoder: ‘GCN’. All experiments except for bulk RNA-seq integration used these defaults. For bulk RNA-seq data, which has a small sample size, we recommend setting the train_mode to ‘full_batch’ and setting epoch_integration to 200. Further details are provided in the GitHub tutorials.

Baseline methods
Integration methods
Harmony. Harmony19 has demonstrated superior performance in various benchmarking studies for scRNA-seq data integration. We applied the Python package harmony (v0.0.9) in this study.

scVI. scVI17 is a probabilistic model for scRNA-seq data. We applied the Python package scvi-tools (v0.14.6) in this study.

PyLiger. LIGER26 is a robust online integration model. It has recently been re-implemented as PyLiger25 in Python with improved computational efficiency. We applied the Python package PyLiger (v0.1.3) in this study.

Scanorama. Scanorama27 is a sensitive algorithm that identifies and merges shared cell types between all datasets and accurately integrates heterogeneous scRNA-seq datasets. We applied the Python package Scanorama (v1.7.3) in this study.

SeuratV4. SeuratV429 is one of the most popular single-cell analysis tools based on the R language. We used SeuratV4.2.3 in this study.

scGen. scGen28 is a VAE-based method that can predict gene expression in different perturbations of single-cell data. We use version (2.0.0) in this study.

SAMap. SAMap22 is an algorithm that integrates cross-species single-cell transcriptomic data by utilizing the similarity of gene or protein sequences. We use version (1.0.15) in this study.

SATURN. SATURN30 is a deep learning method for obtaining universal cell embeddings from different species. We download it from GitHub (https://github.com/snap-stanford/saturn) and use it following the tutorials (https://github.com/snap-stanford/SATURN/tree/main/Vignettes/frog_zebrafish_embryogenesis).

scGPT. scGPT19 is a single-cell foundation model trained based on human and mouse scRNA-seq data. We use the latest version (0.2.4) in this study.

Annotation methods
scANVI. scANVI34 is a probabilistic approach to annotate cells in a semi-supervised way. We applied the Python package scvi-tools (v0.14.6).

CellBlast. CellBlast78 is an easy-to-use Python package for cell-type annotation. We applied the Python package CellBlast (v0.5.0) in this study.

SingleCellNet. SingleCellNet18 is a tree model-based R package that can annotate single-cell datasets including cross-species RNA-seq data, and the authors recently migrated this approach into Python (https://github.com/pcahan1/PySingleCellNet/).

It is crucial to note that the benchmarking strategy is task-dependent. For the unsupervised integration and alignment task, CAMEX is compared against label-free methods such as Harmony. For the semi-supervised cell-type annotation task, however, the comparison is made against reference-based methods like SingleCellNet (SCN), which require at least one labeled referenced dataset. In addition, we use default hyperparameters in all experiments for all methods.

Evaluation metrics
Integration metrics. We group the batch correction metrics into two categories: 1) batch correction and 2) conservation of biological variance. Specifically, the batch correction score using metrics including kBET (k-nearest-neighbor Batch Effect Test)31, PCR batch (Principal Component Regression of the batch covariate)31, ASW batch (Average Silhouette Width across batch)31, and Graph connectivity14. The conservation of biological variation score using metrics including ARI (Adjusted Rand Index)32, NMI (Normalized Mutual Information)33, Cell-type ASW14, ISO label F1 (Isolated label F1)14, and ISO label silhouette (Isolated label silhouette)14. We implemented all the above integration metrics by scIB14.

kBET. The kBET algorithm assesses whether the label composition of a cell’s k nearest neighborhood is consistent with the expected label distribution. This evaluation is performed by repeating the test on a random subset of cells, and the outcomes are aggregated as a rejection rate across all tested neighborhoods.

PCR batch. PCR batch quantifies the batch effect removal by comparing the variance contributions of the batch effects of datasets before integration (VC before) and after integration (VC after), respectively.

ASW batch. ASW Batch calculates the silhouette width of cells with respect to batch labels.

Graph connectivity. Graph connectivity assesses whether the graph correctly connects cells of the same cell-type labels among batches.

ARI. The adjusted Rand index (ARI) is a common external evaluation metric for clustering. It assesses the effectiveness of clustering by calculating the number of sample pairs assigned to the same or different clusters in both the true labels and the clustering results.

NMI. NMI computes the normalized mutual information between two clustering results.

Celltype ASW. Cell-type ASW evaluates ASW with respect to cell-type labels, where a higher score means that cells are closer to cells of the same cell type.

Isolated label F1. Isolated label F1 is developed to measure the ability of integration methods to preserve dataset-specific cell types. We adopted the Scanpy pipeline to cluster cells and evaluated the cluster assignment of dataset-specific cell types based on the F1 score.

Isolated label silhouette. Isolated label silhouette measures the conservation of dataset-specific cell types. It evaluates the ASW of dataset-specific cell types.

Annotation metrics. To comprehensively compare the performance of different annotation methods, we employed three kinds of evaluation metrics, including accuracy, precision, and F1 score. Cell annotation is a multi-classification task. Considering the difference in the number of different types of cells in a dataset, all evaluation metrics are weighted. We implemented all the above annotation metrics by Scikit-learn33.

Accuracy. The accuracy can be defined as follows:

accuracy = (TP+TN) / (TP+TN+FP+FN)

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives, and FN is the number of false negatives.

Precision. The precision can be defined as follows:

precision = TP / (TP+FP)

where TP is the number of true positives, and FP is the number of false positives.

F1 score. The F1 score can be interpreted as a harmonic mean of precision and recall, where an F1 score reaches its best value at 1 and worst score at 0. The F1 can be defined as follows:

F1 = 2 * (precision * recall) / (precision + recall)

Batch mean = (kBET + PCR batch + ASW Batch + Graph connectivity) / 4), Biological score (biological mean = (ARI + NMI + Celltype ASW + Isolated label F1 + Isloated label silhouette) / 5), Overall score = (batch mean + biological mean) / 2.

Trajectory conservation analysis
Inspired by previous literature20,48, we applied and adjusted a quantitative metric to evaluate trajectory preservation during integration as Luecken et al.14. This metric is calculated in the following steps:

Pre-integration pseudotime calculation. For each unintegrated (raw) dataset, we first exclude cell types outside our focus. For example, in the testis dataset across 11 species, we exclude somatic and Sertoli cells. We then designate a spermatogonia cell as the trajectory origin (“iroot”) and compute pseudotime values (t_unintegrate) using the scanpy.tl.dpt() function.

For bulk RNA-seq datasets with known developmental time points, we normalize the ground truth temporal values by linearly scaling them to the continuous interval [0, 1], where 0 and 1 represent the start and end points of the developmental process, respectively. This normalization allows for direct comparison with pseudotime estimates while preserving relative temporal relationships.

Post-integration pseudotime analysis. In the integrated dataset (e.g., from CAMEX or Seurat), we select a human spermatogonia cell as the trajectory origin (“iroot”) within the testis dataset. We then compute pseudotime values (t_integrate) using the scanpy.tl.dpt() function.

Trajectory conservation metric. We quantify trajectory preservation using correlation analysis between the pre- and post-integration pseudotime values. Spearman’s rank correlation is used for this quantification.

One-to-one and many-to-many homologous relationships
For the baseline methods, we map the gene symbols to the first or reference species dataset using one-to-one homologous relationships. The many-to-many homology relationships were obtained from the Ensembl database79. The details of how to use it can be found in the GitHub tutorials.

Carnegie stages
Carnegie stages represent a standardized system used to define the 23 stages of embryonic development in vertebrates. The data on ‘the number of days’ in Carnegie stages refer to the days following embryo formation and should be considered an approximate value. These stages track development from the formation of a fertilized egg to an embryo’s progress over approximately 60 days (8 weeks). For instance, a human embryo takes around 60 days to reach Carnegie Stage 23, whereas a chicken embryo completes the same stage in just 10 days. This system enables scientists to compare developmental characteristics across species at equivalent stages, thereby enhancing our understanding of the universal principles governing embryonic development.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "A multi-species single cell integration tool based on graph convolutional framework. For data integration the model consists of\nvariational graph neural network that embeds the data into a shared latent space and decodes, reconstructs, both cellular gene\nexpression profiles and network connectivity between genes and cells. The reconstructed graph can then be used to query by graph\nattention network to annotate the cells against a reference dataset.\nThe method was evaluated by integration of various tissue specific samples from multiple organisms as well as developmental time\ncourse data.\nOverall, this is an interesting method for data integration across evolutionary distinct organisms. The major criticism is that the evaluation\nis not sufficiently thorough. The authors use several metrics of integration for biological conservation and mixing (batch correction).\nHowever, the numerical values are not presented and there are many method comparative statements that are not backed by objective\nvalues. Another limitation of benchmarking was the use of integrated samples in evaluation of cellular annotation rather than using an\nexternal sample.\nMethod:\nDuring the encoder training, what is the optimization function? How are the W matrix updated? What are the initial values of “h” (gene\nand cell)? It’s not clear from the decoder how these values are optimized. The method section is not well organized. It seems that the\nmodel is trained on minimizing L_A and L_F but this is not described sufficiently. Are they optimized in parallel, are there any weights,\ncross terms that are shared?\nIn describing the graph of cell and genes there are six types of edges described. However, “a cell expresses a gene” and a “gene is\nexpressed by a cell” are the same relationship. In effect, there are only five type of edges unless this is a directional graph. If that is the\ncase, what is the meaning of “cell expresses a gene” vs. “gene is expressed in a cell”. Similarly, are there any differences in the Wgc and\nWcg?\nIt seems that the method is restricted to one-to-one gene homology which becomes limiting with more evolutionary distance and\nincreasing number of species added (as the case with sperm). To what extent is the model sensitive to the number of shared\nhomologous genes?\nWhat is the run time performance of CAMEX relative to the other methods?\nEvaluation:\nThe batch correction and biological conservation numerical scores should be included in supplementary material. The bubble plots\nprovide only qualitative views.\nIt will be helpful to provide summary plots of performance as shown in Song et al. Nature Communication, 2023,14:6495 figures 2b and\n7.\nIn the bubble plot it is not clear what constitute a better score. For some measures the score values are very similar to each other. The\nonly measure that shows a range of values is kBET (Suppl. Fig. S3).\nLine 165 - What is the rational for 4:5 weighted score of the batch correction and biological conservation combined score for overall\nmodel performance? Why not give them equal weights?\nIt’s not clear why and when the authors chose to compare (SCN, scANVI Cell BLAST – single cell annotation methods) or when to\ncompare to Harmony, scVI, Scanorama, PyLiger for integration performance?\nLine 249-261, Figure 3b – PyLiger and scVi also seem to preserve the sperm developmental trajectory. The authors state that the\nCAMEX preserved continuous trajectory better than other approaches and that “PyLiger … failed to integrate spermatids and\nspermatocytes erroneously merging somatic cells with spermatids. scVI struggled to effectively integrate Sertoli and spermatocytes cells\nfrom chicken…”. These are all qualitative statements based UMAPs views. There is tendency to over-interpret UMAP plots rather than\nbase evaluation on numerical score metrics (not rank). These statements should be supported by defined evaluation metrics. For\nexample, kBet and PCR batch, ASW batches values for the cell mixing scores or trajectory analysis.\nCross species cell annotation & Supplemental Figure S6&S7- It seems somewhat circular to annotate cells from samples that were used\nin the integration. A more rigorous test would be to annotate new datasets, that was not used in the integration, to evaluate CAMEX’s\nability to transfer labels to new query dataset.\nLine 315: Figure 4 and S9- Similarly to the previous comment, statements about method’s preservation of developmental trajectory or\nlack of (e.g. Harmony) should be supported by an objective metric. It is incorrect to based them on UMAP figures which are strictly\nvisualization.\nFigure 4b and Supple S10b – the developmental time alignments are somewhat inconsistent. While the overall developmental timelines\nare well matched the red line tracing is misleading. For example, mouse brain vs. Opossum brain p6 is most similar to e13 whereas p4\nis most similar to e13&e15. Related inconsistencies are present in Figure 4b where rat e15 matches mouse e12-e14. The authors\nshould comment on this. More importantly, a more rigorous evaluation should be aligning developmental trajectory of datasets that was\nnot included in the original integration.\nMinor:\nWhat is “SOTA” methods ?\nLine 193 – How was “more distinct cell population boundaries” defined and scored?\nLine 195 – What are the kBET and ASW batch score reported?\nLine 223-4 – What does “over-integration” mean in this context?\nLine 261- Figure S6a – What is the color gradient for confusion matrix? Is it similarity scores?\nLine 278 – Figure S8. What is mean embedding? How did they compute the similarity/distances between cell type? Is it based on the\ncells and genes latent space representation or the reconstructed decoder values?\nFigure 3c – Cluster 3 is mostly Somatic cells and cluster 5 is mostly Spermatocytes. This is different from the text in line 293-295.\nLine 310 – Is there a way to quantify the preservation of developmental trajectory? UMAP views are not quantitative. One can draw\nsimilar arrows for scVI, PyLiger and claim that they preserve development.\nFigures S10 and 4b. is the integration of the developmental stages done at once or iteratively, adding each time point at a time?\nLine 350 – “In the scVI results, microglial cells from human, mice, …are not mixed up correctly” – what does that mean? What\ndetermines correct or incorrect mixing? The authors notes deficiencies in PyLiger integration but equally, CAMEX failed to integrate\ninhibitory neurons whereas PyLiger performed better (Figure 5a)\nLine 364 – CAMEX and scANVI are at par. It’s hard to claim CAMEX is better.\nLine 365 – According to Suppl. S12a CAMEX is not better than scANVI. ¬What is the color legend for S12b?\nFigures 5b and S12 – there is no color legends for the pair-wise similarities plots. How are these values calculated?\nFigure 5d – there is a discrepancy between the manual NPC annotation and CAMEX annotation of excitatory neurons (ExN) and\nscANVI. What evidence/corroboration do they show that CAMEX annotation is correct?\nThere are two Figure S13 (one which duplicates S15) which makes it hard to follow the text.\nLine 417-423 This is a hard argument to follow and there are missing references. Also, the claim a correlation between evolutionary\ndistance and proportions of cells is a stretch.\nLine 417 – Looking at Figures 6b and S14 there is very little support to the conclusion that CCL3 and GLDN subtypes are present in\nChimpanzee and Rhesus.", "reviewer_2": "The authors describe CAMEX, an algorithm for cross-species integration of scRNA-seq data based on heterogeneous graph neural\nnetworks. The authors applied the method to diverse datasets and the method showed competitive performance. Here, I outline some\npoints that should be addressed to improve the manuscript:\n1. The authors correctly pointed out in the introduction that some other recent methods, such as SAMap and SATURN, utilise non-one-\nto-one homology between genes. However, there is no comparison between CAMEX and any of these methods, but only with older\nmethods which rely on one-to-one orthologs, which hardly seems fair. Given that including many-to-many relationships between genes\nis essential to the method, the authors should show a benchmark with SAMap or SATURN.\n2. Given that in many cases CAMEX demonstrated comparative performance with other methods, but it costs much more computational\npower as a deep learning method, authors should describe the resource consumption, such as Max memory and running time, in the\nbenchmarks to help future users choose.\n3. Line 130, can the authors please explain what does it mean that the initial cell node features “correspond to its original expression”.\nAren’t gene expression info in the cell-gene edges? What is the added benefit of encoding these into node feature, given the very high\ndimensionality of gene expression, this could be highly resource consuming and generate large objects.\n4. In the small UMAP crops in Fig. 1b, it appears that CAMEX mixed macrophages and Kupffer cells after integration, when scVI kept\nthem separated. Could the authors explain? If this is an issue, how come CAMEX had a higher biological conservation score?\n5. Line 204, does “attention” refer to attention in NNs or simply marker gene expression in Fig. 2d? Please clarify and reformulate the\nsentence because I don’t see that CAMEX HetroGNN has an attention mechanism.\n6. In benchmark UMAPs such as Fig 2b and 3b, there is no clear explanation about what the circled population is referring to in the text\nor in the captions.\n7. In line 264, I wonder what is the reason of using a KNN classifier on embeddings instead of doing graph attention classifier label\ntransfer that is part of CAMEX tool?\n8. In line 308, did the authors indeed use bulk RNA-seq? Is this method also applicable to bulk RNA-seq data, and could the authors\nexplain exactly how it was performed?\n9. In line 327, I am not sure what the authors mean by “predicting their stages beyond the specified ranges”. What is the “prediction”, and\nwhat are the “specified ranges”? Please specify and explain. This mention of “predicting beyond Carnegie stages” appeared a few times in\nthe manuscript as a core result and could greatly benefit from more background and a clearer explanation.\n10. In line 368, it's unclear which population in Fig. 5c-d is the “unknown” population, and I find it hard to understand why these are\n“unknown” but “originally defined as neural progenitor cells”. Please rewrite this paragraph to provide clarity.\n11. The last section on discovering new populations generally seems to lack proof. All results are computational hypotheses, but the\nclaims of discovering new populations are very strong. In Fig. 6b, the GLDN populations in Chimp and Rhesus are highly likely artefacts\ndue to the decision boundary of the GAT classifier. Any small tweak in the decision boundaries can lead to these cells being assigned to\none or another cell type, and there is no proof that these exist in Chimp and Rhesus. Furthermore, there is no visible CCL3 population as\nthere are no purple dots in other species, so I don’t know where the CCL3 population is that the authors claim.\n12. Similarly, the proportion of cell types is likely a computational artefact given the small data size and obscure cell type annotation\ntransfer boundaries. Therefore, the discussions in lines 420-423 cannot be justified. These claims have almost no relation to current\nevolutionary knowledge and cannot be considered reasonable. I suggest removing them.", "reviewer_3": "The authors proposed CAMEX, a heterogeneous Graph Neural Network tool, to address these challenges by leveraging many-to-many\nhomologous relationships for multi-species scRNA-seq data integration, alignment, and annotation. It outperforms other methods on\ncross-species datasets, helps align species at different developmental stages, enables the detection of species-specific cell types and\nmarker genes, and has the potential to offer insights into evolutionary forces across species at single-cell resolution.\n1. What are the noteworthy results?\nCAMEX outperforms state-of-the-art methods in cross-species integration across datasets spanning various species. It enables the\nalignment of developmental stages and the detection of species-specific cell types and marker genes. The tool leverages many-to-many\nhomologous relationships, addressing a key challenge in cross-species integration.\nHowever, there are few methods compared in the results, and further comparison methods need to be added, such as Seurat v4 [1],\nscGen[2], SAMap[3], SATURN [4], etc. At the same time, are the results of the model advantageous compared to the popular large\nlanguage model-based methods that have been compared recently? The authors are requested to supplement the\ncorresponding experiments.\n2. Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not\noriginal, please provide relevant references.\nMulti-species scRNA-seq integration is a critical area for evolutionary biology and comparative genomics. CAMEX’s use of\nheterogeneous GNNs to model many-to-many homology is novel and addresses limitations of existing tools (e.g., Seurat, LIGER, or\nSCOT), which often assume one-to-one orthology.\nHowever, this paper does not compare CAMEX with previous GNN-based methods, especially with the CAME model proposed by the\nsame authors before. The author should clearly state the differences between CAMEX and CAME to prove that it has made an\ninnovative contribution sufficiently. Meanwhile, in the experimental part, the author also does not compare with the results of CAME.\nPlease further explain the reasons for this.\n3. Does the work support the conclusions and claims, or is additional evidence needed?\nThe authors have carried out some experiments on cross-species data sets. However, the distance relationships between different\nspecies were not distinguished in the relevant cross-species experiments. It is recommended that more detailed experiments be\nconducted in reference [5] to verify the effectiveness of the model more comprehensively.\n4. Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit publication or require revision?\n1). When constructing homology relationships, how are homology relationships defined? Are homology databases or sequence\nalignment methods used?\n2). How is the data partitioned during the neural network training process? Is the partition completely random? Or is it based on some\nexperience?\n3). There are relatively few methods for conducting comparative experiments. The authors are requested to further supplement the\ncomparative methods experiments to verify the model performance's superiority.\nBefore the manuscript is published, the authors need to clarify and revise these issues.\n5. Is the methodology sound? Does the work meet the expected standards in your field?\nCAMEX’s use of heterogeneous GNNs to model many-to-many homology is novel and addresses limitations of existing tools, which\noften assume one-to-one orthology. However, specific experiments, such as ablation experiments, are lacking in verifying which\nmodules in the proposed model play a crucial role during the experimental process.\nThe paper's results align with the expectations in the relevant fields. However, it is necessary to supplement the experiments of\ncomparison methods further, especially those proposed in recent years, to verify the advancement of the methods.\n6. Is there enough detail provided in the methods for the work to be reproduced?\nThe author provides the corresponding code and data source and gives relevant hyper-parameter information to replicate the work."}, "author_response": "Rebuttal Summary\nResponse to Review #1\nA multi-species single cell integration tool based on graph convolutional\nframework. For data integration the model consists of variational graph neural\nnetwork that embeds the data into a shared latent space and decodes,\nreconstructs, both cellular gene expression profiles and network connectivity\nbetween genes and cells. The reconstructed graph can then be used to query\nby graph attention network to annotate the cells against a reference dataset.\nThe method was evaluated by integration of various tissue specific samples\nfrom multiple organisms as well as developmental time course data.\nOverall, this is an interesting method for data integration across evolutionary\ndistinct organisms. The major criticism is that the evaluation is not sufficiently\nthorough. The authors use several metrics of integration for biological\nconservation and mixing (batch correction). However, the numerical values are\nnot presented and there are many method comparative statements that are not\nbacked by objective values. Another limitation of benchmarking was the use of\nintegrated samples in evaluation of cellular annotation rather than using an\nexternal sample.\nResponse: We appreciate your kind recognition and detailed comments and\nsuggestions. Following your suggestions, we have extensively revised our\nmanuscript and added more evaluation details.\n1. Method: During the encoder training, what is the optimization function?\nHow are the W matrix updated? What are the initial values of “h” (gene\nand cell)? It’s not clear from the decoder how these values are optimized.\nThe method section is not well organized. It seems that the model is\ntrained on minimizing L_A and L_F but this is not described sufficiently.\nAre they optimized in parallel, are there any weights, cross terms that are\nshared?\nResponse: We have carefully rewritten the loss functions for integration and\nannotation tasks in the revised manuscript. Specifically, we use Adam as the\noptimization function. The parameter matrix 𝑊 is updated using the Adam\noptimizer via backpropagation algorithm. The initial feature of cell is the\nnormalized, log-transformed expression value and the initial feature of gene is\nset to zero. The ℎ is the cell and gene embedding that maps the expression to\nthe latent space through the encoder. The loss function of integration task is\n𝐿 = 𝐿 + 𝐿 seen below and the loss function of annotation task is\n𝑖𝑛𝑡𝑒𝑟𝑔𝑟𝑎𝑡𝑖𝑜𝑛 𝐴 𝐹\n𝐿 seen below. 𝐿 and 𝐿 are optimized in parallel.\n𝑎𝑛𝑛𝑜𝑡𝑎𝑡𝑖𝑜𝑛 𝐴 𝐹\nAccordingly, we added a more detailed description of the data preprocessing\nand model training process in the Methods section (page 15, lines 602-735)\n(marked in red):\n“Heterogeneous graph of cells and genes. Each gene expression matrix with\n𝑁 cells and 𝑀 genes can be defined as 𝑋 = (𝑥 ,𝑥 ,…,𝑥 )𝑇 ∈ 𝑅𝑁×𝑀 . ffter\n1 2 𝑁\nnormalization, log-transformation, selection of highly variable genes, clustering,\nand obtaining differential genes, we convert each expression profile into a cell-\nto-gene bipartite diagram 𝐴, where cell 𝑛 links gene 𝑚 if and only if cell 𝑛\nexpresses gene 𝑚 . hhe cell 𝑛 and cell 𝑛 will be connected if they are\n1 2\nneighbors calculated by the kNN algorithm in the PCf space. In addition, we\nadded a self-loop to each cell and gene, respectively. hhen, we connected\ngenes in each 𝐴 from different species through the many-to-many-\nhomologous gene relationships. In short, there are two types of nodes including\ncells and genes, and six types of edges (relations) including ‘a cell expresses\na gene’, ‘a gene is expressed by a cell’, ‘two cells that are kNN neighbors’, ‘self-\nloop of a cell’, ‘self-loop of a gene’ and ‘two genes that are orthologous ones’ in\nthe heterogeneous graph. We initialized cell node features with preprocessed\nexpression data (normalized, log-transformed values from the top 2,000 highly\nvariable genes (HVGs)), while gene nodes were initialized as zero vectors. hhe\nresulting heterogeneous graph incorporated both cell and gene nodes with their\nrespective features.\nHeterogeneous graph neural network encoder. CfMEX adopts a\nheterogeneous graph neural network [1] as the encoder inspired by the\nrelational graph convolutional network. We design six independent parameter\nmatrices for each of the above relationships as follows. 𝑊 (a cell expresses\n𝑐𝑔\na gene), 𝑊 (a gene is expressed by a cell), 𝑊 (two cells that are kNN\n𝑔𝑐 𝑐𝑐\nneighbors), 𝑊 (two genes that are orthologous ones), 𝑊 (self-loop of a cell),\n𝑔𝑔 𝑐\nand 𝑊 (self-loop of a gene). hhe encoder transforms the features of the\n𝑔\nneighbors in a non-linear way and converges to the center node according to\nthe message-passing paradigm. hhe (𝑙 +1 )-th layer representation of cell\n(𝑙+1) (𝑙+1)\nℎ and gene ℎ can be defined as follows:\n𝑛𝑖 𝑚𝑖\n(𝑙+1) 1 (𝑙) (𝑙) (𝑙) (𝑙)\nℎ = 𝜎(∑ (∑ 𝑊 ℎ +𝑊 ℎ )),\n𝑛𝑖 𝑟∈𝑅 𝑗∈𝑁𝑟𝑐\n𝑖,𝑟\n𝑟 𝑗 𝑐 𝑛𝑖\n𝑖\n(𝑙+1) 1 (𝑙) (𝑙) (𝑙) (𝑙)\nℎ = 𝜎(∑ (∑ 𝑊 ℎ +𝑊 ℎ )),\n𝑚𝑖 𝑟∈𝑅 𝑗∈𝑁𝑟𝑐\n𝑖,𝑟\n𝑟 𝑗 𝑔 𝑚𝑖\n𝑖\n(𝑙) (𝑙)\nwhere ℎ and ℎ represent the 𝑙-th layer representation of cells and genes,\n𝑛𝑖 𝑚𝑖\n(𝑙) (𝑙)\nrespectively, 𝑊 and 𝑊 represent the 𝑙-th layer weight matrices of ‘self-\n𝑐 𝑔\nloop of a cell’ and ‘self-loop of a gene’, respectively, 𝑁𝑟 denotes the set of\n𝑖\nneighbor indices of the node (cell or gene) 𝑖 under relations 𝑟 ∈ 𝑅, 𝑐 is a\n𝑖,𝑟\n(𝑙)\nnormalization factors, 𝑊 is the weight matrix of the 𝑙-th layer transformation\n𝑟\n(𝑙)\nof the relationship 𝑟, ℎ is the 𝑙-th layer representation of neighbor j, and 𝜎\n𝑗\nis the activation function (i.e., Leaky ReLU here).\nho reduce the number of parameters to prevent overfitting, we add a\nhyperparameter ‘share’ that can share weight matrices of the same relationship\nin different layers such as 𝑤𝑙 = 𝑤𝑙+1 . ffter the message passing through\n𝑐𝑐 𝑐𝑐\nseveral layers of neural networks, the cell and gene representation can be\nmapped into the same embedding space.\nIntegration task:\nEdge and feature reconstruction decoder. When integrating different\nscRNf-seq data in an unsupervised manner, we employ two decoders to\nreconstruct edges and features, respectively. Specifically, the dot-product\ndecoder aims to reconstruct the edges of the input heterogeneous graph. Edge\npredictions in the heterogeneous graph are generated using the dot product of\ncell and gene embeddings:\n𝐴′ = 𝐻 ⋅𝐻 ,\n𝑁 𝑀\nwhere 𝐻 and 𝐻 are the cell and gene embeddings, respectively, and 𝐴′ is\n𝑁 𝑀\nthe reconstructed adjacency matrix.\nWe measure the quality of each edge in the adjacency matrix using cross-\nentropy loss:\n𝐿 = −𝛴𝑁 𝛴𝑀 [𝐴 𝑙𝑜𝑔𝐴̂ +(1−𝐴 )𝑙𝑜𝑔 (1−𝐴̂ )],\n𝐴 𝑖=1 𝑗=1 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗\nwhere 𝐴 is the edge between the i-th cell and the j-th gene in the original\n𝑖𝑗\ngraph.\nhhe feature reconstruction decoder aims to reconstruct cell features\nthrough a fully connected layer:\nℎ′ = 𝑊 (ℎ (𝑙) ),\n0 𝑓𝑐 𝑛𝑖\nwhere 𝑤 is the weight of a fully connected neural network and ℎ′ is the\n𝑓𝑐 0\nreconstructed cell feature.\nWe utilize the mean squared error to measure the distance between the\noriginal ℎ and the predicted ℎ′:\n0 0\n𝐿 = 1 𝛴𝑁 (ℎ −ℎ′).\n𝐹 𝑖=1 0 0\n𝑁\nLoss function and optimizer in integration task. hhe unsupervised\nintegration task loss function can be defined as follows:\n𝐿 = 𝐿 + 𝐿\n𝑖𝑛𝑡𝑒𝑟𝑔𝑟𝑎𝑡𝑖𝑜𝑛 𝐴 𝐹\nwhere 𝐿 and 𝐿 are optimized in parallel.\n𝐴 𝐹\nhhe integration task utilizes a graph attention network autoencoder with an\nencoder-decoder architecture, which is trained end-to-end through\nbackpropagation [2]. hhe model optimization uses fdam optimizer [3], a variant\nof the Stochastic Gradient Descent algorithm, for parameter optimization and\nimplemented in Pyhorch [4] using a mini-batch training strategy. We use the\nNVIDIf GPU laptop RhX3070 to accelerate the training process. By default,\nboth the encoder and decoder architectures contain two bottle neck layers. We\napplied a 'share' hyperparameter, which allows these layers to share\nparameters respectively, thereby reducing GPU memory consumption.\nAnnotation task:\nGraph attention classifier. Cell labels are not mandatory during the integration\ntask. However, during annotation, it is essential to have one or more reference\ndatasets with pre-defined manual annotations. Consequently, the model could\nbe trained in a semi-supervised manner. hhe cell annotations provided by the\nreference dataset can be denoted as: 𝑌 = (𝑦 ,𝑦 ,…,𝑦 ) ∈ 𝑅𝐶 , where 𝑁𝑟\n1 2 𝑁𝑟\nrepresents the number of cells in the reference dataset and C denotes the set\nof cell types. When performing cell annotation tasks, cell embeddings can be\ninput into a Gfh classifier. Specifically, the attention coefficients are first\ncomputed as follows:\n(𝑙) (𝑙)\n𝑒 = 𝑎([𝑊ℎ ||𝑊ℎ ]),\n𝑖𝑗 𝑖 𝑗\n(𝑙) (𝑙)\nwhere ℎ and ℎ is the 𝑙-th layer representations of cells i and j.\n𝑖 𝑗\nShared parameter 𝑊 maps cell embedding to a common space. hhe\nconcatenation operation (denoted by ||) is applied. Finally, the function 𝑎 maps\nthe concatenated result to a scalar representing the relevance between nodes\ni and j. Subsequently, the relevance scores are normalized to obtain attention\ncoefficients:\n𝑒𝑥𝑝(𝐿𝑒𝑎𝑘𝑦𝑅𝑒𝐿𝑈(𝑒 ))\n𝑖𝑗\n𝛼 = .\n𝑖𝑗\n𝛴 𝑒𝑥𝑝(𝐿𝑒𝑎𝑘𝑦𝑅𝑒𝐿𝑈(𝑒 ))\n𝑘∈𝑁𝑖 𝑖𝑗\nWe use the multi-head self-attention mechanism to predict the final results:\n𝐾\n𝑦̂ = 1 ∑ 𝐾 𝑦̂𝑘 = 1 ∑ ∑ 𝛼 𝑊 (𝑙) ℎ (𝑙) ,\n𝐾 𝑘=1 𝐾 𝑗∈𝑁𝑟 𝑖𝑗 𝑟 𝑛𝑖\n𝑘=1 𝑖\nwhere K is the total number of attention-heads, set as eight by default.\nLoss function and optimizer in annotation task. hhe semi-supervised\nannotation task loss function is a cross-entropy loss function can be defined as\nfollows:\n𝑁\n𝐶\n1\n𝐿 = − ∑∑𝑦 𝑙𝑜𝑔(𝑦̂ )\n𝑎𝑛𝑛𝑜𝑡𝑎𝑡𝑖𝑜𝑛 𝑁 𝑖,𝑐 𝑖,𝑐\n𝑐=1\n𝑖=1\nwhere 𝑁 represents sample number, 𝐶 represents class number, 𝑦\n𝑖,𝑐\nrepresents the true label of sample 𝑖 on category 𝑐, 𝑦̂ represents the model\n𝑖,𝑐\npredicts the probability that sample 𝑖 belongs to category 𝑐 . hhe optimizer,\nand implementation details is the same in as the integration task in 'Loss\nFunction and Optimizer in Integration hasks'.\nIn a word, integration and annotation tasks share the same encoder\narchitecture. During integration, unsupervised training is achieved without\nrequiring cell type labels, while annotation necessitates a reference dataset with\ncell type labels for semi-supervised training.”\n2. In describing the graph of cell and genes there are six types of edges\ndescribed. However, “a cell expresses a gene” and a “gene is expressed\nby a cell” are the same relationship. In effect, there are only five type of\nedges unless this is a directional graph. If that is the case, what is the\nmeaning of “cell expresses a gene” vs. “gene is expressed in a cell”.\nSimilarly, are there any differences in the Wgc and Wcg?\nResponse: This is a directed heterogeneous attribute graph. Let's take \"Six\ntypes of heterogeneous edges\" in the lower left corner of Figure 1b as an\nexample for illustration. Among them, the circle represents genes and the\ntriangle represents cells. The left half refers to taking cells as the target nodes,\nincluding 𝑊 (two genes that are orthologous ones), 𝑊 (a cell expresses a\n𝑔𝑔 𝑐𝑔\ngene), and 𝑊 (self-loop of a gene); The right half takes genes as the target\n𝑔\nnodes and contains three kinds of relationships including 𝑊 (a gene is\n𝑔𝑐\nexpressed by a cell), 𝑊 (two cells that are kNN neighbors), and 𝑊 (self-\n𝑐𝑐 𝑐\nloop of a cell). There are a total of six kinds of edges. There is no inherent\nbiological difference between 𝑊 and 𝑊 .\n𝑔𝑐 𝑐𝑔\nFurthermore, we have added descriptions in the Training details of the Methods\nsection (page19, lines 729-735):\n“(3) Graph storage: To efficiently store large-scale graphs, we implemented\ndirected heterogeneous attribute graphs and graph neural networks using the\nframework DGL library. hhe DGL utilizes the coordinate list matrix (COO) matrix\nfor sparse data representation. Its storage logic follows a directed graph\nstructure, mapping relationships from source nodes to target nodes. ho\nmaintain consistency with the COO matrix format and align with other edge\ntypes, we save six kinds of edges independently and opted not to merge the\n𝑊 and 𝑊 .”\n𝑔𝑐 𝑐𝑔\n3. It seems that the method is restricted to one-to-one gene homology which\nbecomes limiting with more evolutionary distance and increasing number\nof species added (as the case with sperm). To what extent is the model\nsensitive to the number of shared homologous genes?\nResponse: In fact, we utilized many-to-many orthologous genes when\nconstructing the gene orthologous edges. Additionally, to comprehensively\nevaluate CAMEX's sensitivity with high evolutionary distance and an increasing\nnumber of species, we designed two kinds of integration experiments in the\ntestis dataset across 11 species. Our experimental framework included:\n1. Pairwise integration analyses between humans and each target species,\ncovering a broad phylogenetic spectrum from closely-related primates to distant\nvertebrates.\n2. Multi-species integration scenarios with progressively expanded species sets,\nsystematically increasing from 2 to 11 species.\nWe added Supplementary Figures S12-S24 to clarify this and summarized\nthe results of CAMEX and other baseline methods (page 18, lines 294-311,\nSupplementary Figures S22 and S24). In short, CAMEX is insensitive and\nrobust to the number of shared homologous genes and the distance between\nspecies.\n“Inspired by the previous literature, to comprehensively evaluate CfMEX's\nperformance across varying evolutionary distances and increasing numbers of\nspecies, we designed two integration scenarios. First, we conducted pairwise\nintegration analyses between human and each target species, covering a broad\nphylogenetic spectrum from closely-related primates to distant vertebrates. hhe\nresults clearly show CfMEX excelled in batch correction, while SfhURN\nperformed strongly in biological variation conservation. Overall, CfMEX\ndemonstrated an advantage over other methods in both quantitative metrics,\nrankings, and qualitative UMfP visualization. Second, we created multi-species\nintegration scenarios with progressively expanded species sets, systematically\nincreasing from two to eleven species. CfMEX also consistently maintained\nhigh performance across all scenarios (Supplementary Figs. S12-S24,\nSupplementary Tables S11-S21). Overall, most integration algorithms,\nparticularly those relying solely on one-to-one homologous genes, show\ndiminished performance as the evolutionary distance between species\nincreases. However, methods like CfMEX, which leverage many-to-many\nhomologous genes, and ShfURN, which utilizes sequence similarity, are less\nimpacted by this increasing species divergence.”\nSupplementary Figure S22. Overall integration performance achieved by\nvarious methods when applied to a two-species testis data, allowing for\ncomparative evaluation. The line graph illustrates the integration performance\nof humans with other species. Different colors represent different methods.\nThree distinct lines represent: Batch correction score, Biological variation\nconservation (Bio) score, and Overall score. For the line charts on the left, a\nhigher value indicates better performance. Conversely, for the box plots on the\nright, a lower ranking signifies better performance.\nSupplementary Figure S24. Overall integration performance achieved by\nCAMEX when applied to the testis dataset, ranging from two-species to eleven-\nspecies integrations, allowing for comparative evaluation. Three distinct lines\nrepresent key performance metrics: the Batch correction score, the Biological\nvariation conservation (Bio) score, and the Overall score.\n4. What is the run time performance of CAMEX relative to the other\nmethods?\nResponse: Yes, time consumption is very important. Following your\nsuggestions, we have compared the computational efficiency of different\nmethods by measuring running time on all the datasets. The details can be seen\nin Supplementary Tables S2, S6, S8, S10, and S25. The table indicates that\ndeep learning approaches such as CAMEX and scVI, while slower than\ntraditional methods like Harmony and Scanorama, still operate within an\nacceptable time frame. Furthermore, the CPU memory required by CAMEX\nduring preprocessing and graph building is significant, yet it is lower than or\ncomparable to that of SATURN and SAMap. Notably, due to its implementation\nof graph minibatch training, CAMEX maintains a relatively constant level of\nGPU memory usage, ensuring its scalability.\nThe description of time consumption such as liver dataset result can be seen in\n(page 7, lines 224-226):\n“We also performed computational efficiency, hyperparameter comparison and\nablation experiments (Supplementary Figs. S2-S4, Supplementary hables S2-\nS4).”\nEvaluation:\n5. The batch correction and biological conservation numerical scores should\nbe included in supplementary material. The bubble plots provide only\nqualitative views. It will be helpful to provide summary plots of\nperformance as shown in Song et al. Nature Communication,\n2023,14:6495 figures 2b and 7 [5]. In the bubble plot it is not clear what\nconstitute a better score. For some measures the score values are very\nsimilar to each other. The only measure that shows a range of values is\nkBET (Suppl. Fig. S3).\nResponse: According to your suggestions, we provided the details of the\nscores of different methods for all datasets in the Supplementary Tables S1,\nS5, S7, S9, S11-S22, and S24. Following your suggestions and drawing\ninspiration from Song et al. Nature Communications [5], we created summary\npoint plots of performance where the horizontal and vertical axes represent\nbatch correction score and biological variation conservation score for all\ndatasets. We take the test on the liver dataset as an example (Supplementary\nFigure S1d and Table S1). We provide the scores of 10 methods on the liver\ndataset in the Table S1.\nSupplementary Figure S1. Statistics for the liver dataset across four species.\na, Liver scRNA-seq dataset includes samples from humans, macaques, mice,\nand zebrafishes. b, Number of detected UMIs and genes in different species.\nc, Number of cells in each cell type for each species. Each dataset comprises\ndistinct cell types that are unique to that respective species. d, Summary plots\nillustrate performance, with the horizontal axis representing the average batch\ncorrection score and the vertical axis indicating the average biological variation\nconservation score. e, UMAP visualization of CAMEX integration result, where\ncolor denotes different species. f, Marker genes of Macrophages and Kupffer\ncells in different species.\n6. Line 165 - What is the rational for 4:5 weighted score of the batch\ncorrection and biological conservation combined score for overall model\nperformance? Why not give them equal weights?\nResponse: We calculated the overall score inspired by the previous literature:\nhttps://github.com/jiazhao97/Portal-reproducibility/blob/main/Benchmarking-\nSpleen.ipynb. According to the requirements of the reviewers, we adjusted the\ncalculation method of the overall score. Specifically, we first calculate the\naverage value of the batch score (batch mean = (kBET + PCR batch + ASW\nBatch + Graph connectivity) / 4) and the biological score (biological mean =\n(ARI + NMI + Cell type ASW + Isolated label F1 + Isloated label silhouette) / 5),\nand then the overall score = (batch mean + biological mean) / 2.\nWe recalculated and replotted the bubble charts for all experiments.\nAdditionally, we provide the specific numerical scores for each experiment in\nthe Supplementary Tables (response 5). It is important to note these\nrecalculated results are consistent with, and do not affect the original conclusion.\n7. It’s not clear why and when the authors chose to compare (SCN, scANVI\nCell BLAST – single cell annotation methods) or when to compare to\nHarmony, scVI, Scanorama, PyLiger for integration performance?\nResponse: We designed CAMEX with two distinct modules to handle different\nanalytical scenarios(Figure 1b and Methods, page 15, lines 603-735). For\nintegration tasks, CAMEX operates without needing any cell type labels across\ndatasets, employing an entirely unsupervised training approach. Conversely,\nfor annotation tasks, CAMEX utilizes a semi-supervised approach. It requires\ncell type labels only within the reference dataset and employs a graph attention\nclassifier to transfer these labels from the reference to the query data. We have\ncarefully clarified the differences:\nTherefore, when carrying out the integration task, we compared CAMEX with\nmethods such as Harmony. When performing the annotation task, we\ncompared CAMEX with methods like SCN. We added some descriptions to\nclarify this in the Baseline methods section (page 19, lines 775-777):\n“Notably, when conducting integration tasks, we compare CfMEX with\nintegration methods such as Harmony. When conducting annotation tasks, we\ncompare CfMEX with annotation methods such as SingleCellNet (SCN).”\n8. Line 249-261, Figure 3b – PyLiger and scVI also seem to preserve the\nsperm developmental trajectory. The authors state that the CAMEX\npreserved continuous trajectory better than other approaches and that\n“PyLiger … failed to integrate spermatids and spermatocytes erroneously\nmerging somatic cells with spermatids. scVI struggled to effectively\nintegrate Sertoli and spermatocytes cells from chicken…”. These are all\nqualitative statements based UMAPs views. There is tendency to over-\ninterpret UMAP plots rather than base evaluation on numerical score\nmetrics (not rank). These statements should be supported by defined\nevaluation metrics. For example, kBet and PCR batch, ASW batches\nvalues for the cell mixing scores or trajectory analysis.\nResponse: This is a very valuable comment. Following your advice,\nwe\nproposed a new trajectory conservation metric and it can be calculated as\nfollows (page 21, lines 838-861):\n“Trajectory Conservation Analysis:\nInspired by the previous literature [19, 48], we developed a quantitative metric\nto evaluate trajectory preservation during integration. hhis metric is calculated\nin the following steps:\nPre-integration Pseudotime Calculation: For each unintegrated (raw)\ndataset, we first exclude cell types outside our focus. For example, in the testis\ndataset across 11 species, we excluded somatic and Sertoli cells. We then\ndesignate a spermatogonia cell as the trajectory origin (\"iroot\") and compute\npseudotime values (t_unintegrate) using the scanpy.tl.dpt() function.\nFor bulk RNf-seq datasets with known developmental timepoints, we\nnormalize the ground truth temporal values by linearly scaling them to the\ncontinuous interval [0,1], where 0 and 1 represent the start and end points of\nthe developmental process, respectively. hhis normalization allows for direct\ncomparison with pseudotime estimates while preserving relative temporal\nrelationships.\nPost-integration Pseudotime Analysis: In the integrated dataset (e.g., from\nCfMEX or Seurat), we select a human spermatogonia cell as the trajectory\norigin (\"iroot\") within the testis dataset. We then compute pseudotime values\n(t_integrate) using the scanpy.tl.dpt() function.\nTrajectory Conservation Metric: We quantify trajectory preservation using\ncorrelation analysis between the pre- and post-integration pseudotime values.\nBoth the Pearson correlation coefficient and Spearman rank correlation are\nused for this quantification.”\nWe provide the related result (page 8, lines 283-285):\n“Furthermore, we calculated the pseudo-times both before and after integration,\nusing their correlation as the trajectory reservation metric. CfMEX achieved\nthe highest value in this regard (Supplementary Fig. S10).”\nSupplementary Figure S10. Trajectory conservation analysis for the testis\ndataset across 11 species. This figure presents a UMAP visualization of a joint\nembedding space for the testis dataset, encompassing 11 species. Notably,\nSomatic and Sertoli cell types were excluded from this analysis. The first row\ndisplays cell type annotations by color, while the second row shows pseudo-\ntime after integration. The third row is the correlation dot plot (the vertical axis\nrepresents t_integrate, and the horizontal axis represents t_unintegrate).\n9. Cross species cell annotation & Supplemental Figure S6&S7- It seems\nsomewhat circular to annotate cells from samples that were used in the\nintegration. A more rigorous test would be to annotate new datasets, that\nwas not used in the integration, to evaluate CAMEX’s ability to transfer\nlabels to new query dataset.\nResponse: As described in the Methods, our integration process is\nunsupervised and does not require cell labels. Therefore, annotations assigned\nusing the integrated cell embedding can serve as a method to evaluate the\nquality of the integration. Similar applications have been reported in previous\nliterature. For example, Figure 3c in [6] and Figure 3B in [7].\n[6]:\n[7]:\n10. Line 315: Figure 4 and S9- Similarly to the previous comment, statements\nabout method’s preservation of developmental trajectory or lack of (e.g.\nHarmony) should be supported by an objective metric. It is incorrect to\nbased them on UMAP figures which are strictly visualization.\nResponse: Yes, we proposed a new trajectory conservation metric (page 21,\nlines 838-861). The results of the bulk dataset can be seen in Supplementary\nFigures S28 and S29. Based on this metric, we can see that CAMEX,\nscGen, and Harmony achieved the highest scores. However, some cells cannot be\ncorrectly computed by Harmony. In conclusion, CAMEX and scGen can better\npreserve the developmental trajectory in this scenario.\nSupplementary Figure S28. Integration results of different methods on the\nbulk dataset. a, Overall score evaluation for integration of different methods in\nbulk dataset. b, UMAP visualization of a joint embedding space of different\nmethods. Colors represent species annotations (first row), cell-type annotations\n(second row), stage annotations (third row), pseudo-time (fourth row), and\ncorrelation dot plot (fifth row, the vertical axis represents t_integrate, and the\nhorizontal axis represents t_unintegrate).\nSupplementary Figure S29. Integration results of different methods on the\nbulk dataset. a, UMAP visualization of a joint embedding space of different\nmethods. Colors represent species annotations (first row), cell-type annotations\n(second row), stage annotations (third row), pseudo-time (fourth row), and\ncorrelation dot plot (fifth row, the vertical axis represents t_integrate, and the\nhorizontal axis represents t_unintegrate). b, Summary plots illustrate\nperformance, with the horizontal axis representing the average batch correction\nscore and the vertical axis indicating the average biological variation\nconservation score.\n11. Figure 4b and Supple S10b – the developmental time alignments are\nsomewhat inconsistent. While the overall developmental timelines are well\nmatched the red line tracing is misleading. For example, mouse brain vs.\nOpossum brain p6 is most similar to e13 whereas p4 is most similar to\ne13&e15. Related inconsistencies are present in Figure 4b where rat e15\nmatches mouse e12-e14. The authors should comment on this. More\nimportantly, a more rigorous evaluation should be aligning developmental\ntrajectory of datasets that was not included in the original integration.\nResponse: Yes, you are right. We reconnect the points with the highest\nsimilarity in each column in the figure to avoid misleading (Figure 4b,\nSupplementary Figure S30b). We acknowledge that the developmental\nstages across different species in our results are not perfectly aligned. We\nbelieve this discrepancy stems from both biological factors, such as the\nmismatch between anatomy and gene expression, and computational factors,\nincluding the inherent randomness in the neural network optimization process.\nWe added the following descriptions in the main manuscript. (page 11, lines\n397-411):\n“However, developmental stages across different species in the results are not\nperfectly aligned. We believe there are several contributing factors. Biological\nConsiderations: From a biological perspective, there are two main reasons.\nFirst, the developmental speeds across different species might inherently vary.\nSecond, while we align species based on transcriptome expression, the\nCarnegie developmental stages are primarily defined by morphological\ndevelopment during embryogenesis, which can lead to discrepancies. Model\nand Computational Considerations: From a modeling and computational\nstandpoint, the observed differences can also stem from the randomness in the\nneural network optimization process and the various similarity metrics used for\nembedding calculations. Despite these points, CfMEX still generally achieves\nroughly alignment of developmental processes across species. For time points\nthat are not well aligned, this can serve as a potential objective for studying the\ndevelopmental differences across species.”\nWe recognized that the Methods section might have been insufficiently detailed,\nresulting in a misunderstanding of the integration capabilities of CAMEX. To\nclarify this, its integration does not require any cell type or developmental time\nlabels. Therefore, the inclusion or exclusion of a specific dataset has no impact\non its integration and alignment results.\nMinor:\n12. What is “SOTA” methods ?\nResponse: SOTA is the abbreviation of “State Of The Art”, we added the\ndescription at its first appeared position in the manuscript.\n13. Line 193 – How was “more distinct cell population boundaries” defined and\nscored?\nResponse: The “more distinct cell population boundaries” indicate whether\nthere is only one type of cell present within the same cluster. Following your\ncomment, we provided a more detailed description of this sentence (page 6,\nlines 207-209):\n“CfMEX achieved higher fdjusted Rand Index (fRI) and Normalized Mutual\nInformation (NMI) scores compared to other methods, reflecting its better ability\nto cluster cells of the same type. hhis improved clustering was further\ncorroborated by UMfP visualizations.”\n14. Line 195 – What are the kBET and ASW batch score reported?\nResponse: Previously, we used bubble charts to compare the performance\ndifferences among various methods. According to your suggestions, we have\nincluded point plots where the horizontal and vertical axes represent batch\ncorrection score and biological variation conservation score, respectively, and\nspecific metric values for each metric across different methods in Supplement\nTables S1, S5, S7, S9, S11-S22, and S24.\nFor example, we provide the scores of 10 methods on the liver dataset in the\nSupplementary Table S1.\n15. Line 223-4 – What does “over-integration” mean in this context?\nResponse: Over-integration indicates that some cell types were integrated into\none cluster. In Supplementary Fig. S8, we can see that pyLiger confused two\ntypes of endocrine cells, Delta (yellow) and Gamma (blue) cells. We rewrote\nthis sentence instead of using \"over-integration\" vaguely in the manuscript\n(page 7, lines 249-251):\n“However, Harmony tended to confuse the acinar and ductal cell types. PyLiger\nconfused the acinar and ductal cells and the delta and gamma cells but failed\nto integrate the alpha and beta cell types.”\nSupplementary Figure S8. Benchmarking CAMEX and SOTA methods on the\npancreas dataset across three species. UMAP visualization of a joint\nembedding space on the ovary dataset. Colors represent species annotations\n(first row) and cell-type annotations (second row). Circles represent the cell\npopulation that cannot be integrated very well.\n16. Line 261- Figure S6a – What is the color gradient for confusion matrix? Is\nit similarity scores?\nResponse: It is not similarity scores. Specifically, we first construct a confusion\nmatrix, where the rows represent the true labels and the columns represent the\npredicted labels. Then, we apply row normalization to scale the values in the\nconfusion matrix to a range of [0, 1]. A value closer to 1 is represented by a\ndarker color. An explanation of this and a color bar has been provided in the\nlegend of Supplementary Figure S25.\nSupplementary Figure S25. Integration and kNN annotations of CAMEX on\nthe 11-species testis dataset, focusing on the outcomes for seven primate\nspecies. UMAP visualization of a shared embedding space derived from the\ntestis dataset. The first row of the UMAP uses colors to represent manual\nannotations, while the second row displays colors corresponding to predicted\nannotations. Below the UMAP, a row-normalized confusion matrix is provided\nwhere the rows represent the true labels and the columns correspond to the\npredicted labels. Row normalization scales values within the matrix to a range\nof [0, 1], where a darker color indicates a value closer to 1, signifying higher\naccuracy or agreement.\n17. Line 278 – Figure S8. What is mean embedding? How did they compute\nthe similarity/distances between cell type? Is it based on the cells and\ngenes latent space representation or the reconstructed decoder values?\nResponse: Yes, it is based on the cell and gene latent representation. The\nmean embedding means that we averaged the representations of the same\ntype of cells from the same species. The relationships between cell types can\nbe computed using the scanpy.pl.correlation_matrix function.\nSupplementary Figure S27. Hierarchical clustering of the mean embedding\nfor each cell type across all 11 species. Cell and gene embeddings are from the\nbottleneck layer of CAMEX. Then, we calculated the average embedding for\neach cell type in each species. The Pearson correlation was then applied to\nassess the relationships between these averaged embeddings, with the\nvisualization generated using scanpy.pl.correlation_matrix function.\n18. Figure 3c – Cluster 3 is mostly Somatic cells and cluster 5 is mostly\nSpermatocytes. This is different from the text in line 293-295.\nResponse: Thanks for your kind reminder. Cluster 3 is indeed dominated with somatic\ncells. We have corrected this analysis in the revised manuscript (page 9, lines 342-352).\n“Genes in cluster 3 are predominantly expressed in Somatic cells, and\nenrichment analysis revealed associations with the cell surface receptor\nsignaling pathway and extracellular matrix. Somatic cells are all the cells that\nmake up the body of an organism, except for reproductive cells. In the male\nreproductive organs, somatic cells play multiple roles in supporting, protecting,\nsecreting, and regulating spermatogenesis. Conversely, genes in cluster 5 are\nprimarily expressed in Spermatocytes, and enrichment analysis linked them to\nsexual reproduction, meiotic cell cycle process, meiotic cell cycle, meiotic\nnuclear division, and reproductive process. Meiosis is a well-known process\noccurring in spermatocytes and is consistent with gene function obtained by\nCfMEX.”\n19. Line 310 – Is there a way to quantify the preservation of developmental\ntrajectory? UMAP views are not quantitative. One can draw similar arrows\nfor scVI, PyLiger and claim that they preserve development.\nResponse: Yes. According to your comment, we proposed a new trajectory\nconservation metric and added its definition as follows (page 21, lines 838-\n861). We have provided the results on the bulk dataset (Supplementary\nFigures S28 and S29). Based on this metric, we can see that CAMEX,\nscVI, and Harmony achieved the highest scores. However, some cells cannot be\ncomputed correctly by Harmony. In short, CAMEX and scGen can better\npreserve the developmental trajectory in this scenario. We added the following\nsentences to clarify this analysis (page 10, lines 366-370):\n“Trajectory Conservation Analysis:\nInspired by the previous literature [19, 48], we developed a quantitative metric\nto evaluate trajectory preservation during integration. hhis metric is calculated\nin the following steps:\nPre-integration Pseudotime Calculation: For each unintegrated (raw)\ndataset, we first exclude cell types outside our focus. For example, in the testis\ndataset across 11 species, we excluded somatic and Sertoli cells. We then\ndesignate a spermatogonia cell as the trajectory origin (\"iroot\") and compute\npseudotime values (t_unintegrate) using the scanpy.tl.dpt() function.\nFor bulk RNf-seq datasets with known developmental timepoints, we\nnormalize the ground truth temporal values by linearly scaling them to the\ncontinuous interval [0,1], where 0 and 1 represent the start and end points of\nthe developmental process, respectively. hhis normalization allows for direct\ncomparison with pseudotime estimates while preserving relative temporal\nrelationships.\nPost-integration Pseudotime Analysis: In the integrated dataset (e.g., from\nCfMEX or Seurat), we select a human spermatogonia cell as the trajectory\norigin (\"iroot\") within the testis dataset. We then compute pseudotime values\n(t_integrate) using the scanpy.tl.dpt() function.\nTrajectory Conservation Metric: We quantify trajectory preservation using\ncorrelation analysis between the pre- and post-integration pseudotime values.\nBoth the Pearson correlation coefficient and Spearman rank correlation are\nused for this quantification.”\n“Remarkably, CfMEX achieved the highest overall score, and scGen, Harmony,\nand CfMEX could preserve the continuous development trajectory and\nseparated distinct clusters corresponding to different organs (Fig. 4a,\nSupplementary Figs. S28, S29, and Supplementary Table S22).”\n20. Figures S10 and 4b. is the integration of the developmental stages done at\nonce or iteratively, adding each time point at a time?\nResponse: Yes, the integration of the developmental stages was done in a run.\nWe put the bulk data on all developmental stages of all species into CAMEX to\nperform the integration task without any cell labels or stage labels and obtain\nthe embedding together. The visualization results for all samples across all\nstages and the samples of different organs can be found in Figure 4a and\nSupplementary Figures S28 and S29, respectively.\n21. Line 350 – “In the scVI results, microglial cells from human, mice, …are\nnot mixed up correctly” – what does that mean? What determines correct\nor incorrect mixing? The authors notes deficiencies in PyLiger integration\nbut equally, CAMEX failed to integrate inhibitory neurons whereas PyLiger\nperformed better (Figure 5a)\nResponse: In the scVI results, microglial cells from human, mice, …, are not\nmixed up correctly” means that microglial cells from different species are not\nclustered together in the scVI result (Fig. 5a). According to the UMAP\nvisualization results, the cells do not mix correctly. Yes, the integrated results of\neach method are all imperfect. Following your comment, we have revised it\n(page 11, lines 424-430):\n“fll the methods confused the excitatory and inhibitory neurons. CfMEX failed\nto integrate inhibitory neurons, Harmony and scVI complicated specifically with\nmicroglia cells, and scGen failed to mix astrocytes adequately. PyLiger did not\nproperly integrate astrocytes and confused NPCs with excitatory neurons.\nSfhURN misclassified CGCs as NPCs, and scGPh confused ExN with NPCs,\nin addition to failing to integrate oligodendrocytes.”\n22. Line 364 – CAMEX and scANVI are at par. It’s hard to claim CAMEX is\nbetter. Line 365 – According to Suppl. S12a CAMEX is not better than\nscANVI. ¬What is the color legend for S12b?\nResponse: Yes, CAMEX and scANVI performed comparably. In the original\ntext, we clarified the competitive performance of CAMEX and scANVI compared\nto other methods (page 28, lines 439-441).\n“CfMEX and scfNVI demonstrated better performance compared to the other\ntwo baseline methods, as evaluated by metrics including fCC, Prec, and F1\nscore (Fig. 5b-c, Supplementary Fig. S33a-f).”\nFigure S12b (S33) illustrates the confusion matrix, showing the predicted\nresults versus the actual results. We construct a confusion matrix, where the\nrows represent the true labels and the columns correspond to the predicted\nlabels. Then, we apply row normalization to scale the values in the confusion\nmatrix to a range of 0 to 1. A value closer to 1 is represented by a darker color.\nSupplementary Figure S33. Predicted annotation by CAMEX and baseline\nmethods. The cell types are abbreviated as follows: Astrocyte (Ast), Brain\nPericyte (BP), Endothelial (End), Excitatory Neuron (ExN), Inhibitory Neuron\n(InN), Microglial cell (Micro), Neural Progenitor Cell (NPC), Oligodendrocyte\n(Oligo), Oligodendrocyte Precursor Cell (OPC), Purkinje cell (PC), Cerebellar\nGranule Cell (CGC), Ependymal cell (EC), Macrophage (mø), UnKnown (UK).\na, ACC, PRC, and F1 scores achieved by CAMEX and baseline methods on\nthe mouse dataset. b, Row-normalized confusion matrix of CAMEX and\nbaseline methods on the mouse dataset. A row-normalized confusion matrix is\nprovided where the rows represent the true labels and the columns correspond\nto the predicted labels. Row normalization scales values within the matrix to a\nrange of [0, 1], where a darker color indicates a value closer to 1, signifying\nhigher accuracy or agreement. c, UMAP visualization of the mouse dataset.\nColors denote manual and predicted annotations. Circles represent the End\ncells. d, ACC, PRC, and F1 scores achieved by CAMEX and baseline methods\non the lizard dataset. e, Row-normalized confusion matrix of CAMEX and\nbaseline methods on the lizard dataset. f, UMAP visualization of the lizard\ndataset. Colors denote manual and predicted annotations. Circles represent the\nNPC cells.\n23. Figures 5b and S12 – there is no color legends for the pair-wise similarities\nplots. How are these values calculated?\nResponse: The different color legends in Fig 5b represent different methods.\nI\nbelieve what you are referring to is the color bar in Fig 5c. Our description of\nthe \"color bar\" can be referred to in the previous comment (comment 22).\nFigure 5. CAMEX integrates and annotates relatives and distant species\ncortex scRNA-seq dataset. a, UMAP visualization of joint embedding space.\nColors denote species annotations (first row) and cell type annotations (second\nrow). Circles represent the cell population that can not be integrated very well.\nThe cell types are abbreviated as follows: Astrocyte (Ast), Brain Pericyte (BP),\nEndothelial (End), Excitatory Neuron (ExN), Inhibitory Neuron (InN),\nMicroglial cell (Micro), Neural Progenitor Cell (NPC), Oligodendrocyte (Oligo),\nOligodendrocyte Precursor Cell (OPC ), Purkinje cell (PC), Cerebellar Granule\nCell (CGC), Ependymal cell (EC), Macrophage (mø), UnKnown (UK). b, ACC,\nPRC, and F1 scores achieved by CAMEX and baseline methods on the turtle\ndataset. Unknown is the cell type dost not appear in reference but appear in\nquery dataset. NPC dost not appear in reference human dataset but appear in\nquery turtle and lizard dataset, so we defined them as unknown. c, Row-\nnormalized confusion matrix of CAMEX and baseline methods on turtle dataset.\nRow-normalized confusion matrix is provided where the rows represent the true\nlabels and the columns correspond to the predicted labels. Row normalization\nscales values within the matrix to a range of [0, 1], where a darker color\nindicates a value closer to 1, signifying higher accuracy or agreement. d, UMAP\nvisualization of turtle dataset. Colors denote manual and predicted annotations.\nCircles represent the cell population that cannot be integrated very well.\n24. Figure 5d – there is a discrepancy between the manual NPC annotation\nand CAMEX annotation of excitatory neurons (ExN) and scANVI. What\nevidence/corroboration do they show that CAMEX annotation is correct?\nResponse: Indeed, this is an interesting question. Specifically, since the\ntraining set does not include any NPCs (humans), the model cannot accurately\npredict the NPCs in the test set. However, the prediction results of CAMEX and\nscANVI are both reasonably valid. We have added the following descriptions\n(page 12, lines 443-450):\n“Unknown is the cell type that does not appear in the reference but appears in\nthe query dataset. NPC does not appear in the reference human dataset but\ndoes appear in the query turtle and lizard dataset, so we have defined them as\nunknown. Due to the absence of neural progenitor cells (NPCs) in the human\nreference dataset, both CfMEX and scfNVI struggled to identify them\naccurately. CfMEX misclassified them as excitatory neurons, whereas scfNVI\nidentified them as a mixture of CGCs and other cell types. Interestingly, both\npredictions offer a degree of reasonableness (Supplementary Fig. S34).”\nSupplementary Figure S34. The manual annotation, CAMEX, and scANVI\nprediction for turtle. UMAP visualization of cells in the turtle dataset.\nColors represent species annotations (first row) and marker genes (second row).\nAmong them, SOX4 and SOX11 are the markers of NPC, SLC17A6 is the\nmarker of ExN, and SLC17A6 is the marker of CGC. Circles represent the NPC\ncell population.\n25. There are two Figure S13 (one which duplicates S15) which makes it hard\nto follow the text.\nResponse: Thanks for your kind reminder. We have deleted the redundant one\nand carefully proofread the legends in the Supplementary File.\n26. Line 417-423 This is a hard argument to follow and there are missing\nreferences. Also, the claim a correlation between evolutionary distance\nand proportions of cells is a stretch.\nResponse: According to the suggestion of you and Reviewer 2, we have\nremoved these sentences (line 417-423) in the manuscript.\n27. Line 417 – Looking at Figures 6b and S14 there is very little support to the\nconclusion that CCL3 and GLDN subtypes are present in Chimpanzee and\nRhesus.\nResponse: We agree with your comment. To address this, we have provided a\nmore detailed analysis of the expression patterns for GLDN, CCL3, and\nfPBB1IP. We compared our CAMEX predictions with the original\nliterature's manual annotations (see Supplementary Figure S36 and below).\nOur results show that CAMEX obtained some interesting findings, such as the\nexpression of CCL3 in some cell populations in chimpanzees. However, we also\nacknowledge it can produce some false positives, such as the predicted\nexpression of CCL3 and GLDN in marmosets.\nOverall, we clarify that CAMEX should be considered a bioinformatics tool\nfor generating hypotheses about potential cell subpopulations. These\ncomputational predictions must be validated through subsequent wet-lab\nexperiments. Accordingly, we toned down our conclusion and modified related\nsentences (page13, lines 486-494):\n“ft the same time, we investigated the expression of GLDN, CCL3, and\nfPBB1IP in the original manual annotation and the CfMEX prediction\n(Supplementary Fig. S36c). We found that the prediction of CfMEX has\ncertain rationality, for instance, the expression of CCL3 in the subgroups of\nchimpanzees. However, there are also false positives, such as CCL3 and\nGLDN, in marmosets. Overall, we therefore wish to clarify that CfMEX should\nbe considered a bioinformatics tool for generating hypotheses about potential\ncell subpopulations. hhese computational predictions must be validated\nthrough subsequent wet-lab experiments.”\nSupplementary Figure S36. Integrating the microglia cell datasets of four\nspecies by CAMEX. c, Violin plot of GLDN, CCL3, and fPBB1IP expression\nin original manual annotation and CAMEX prediction clusters.\nReviewer #1 (Remarks on code availability):\nI only had a brief view of the code. It seems well organized and they provide\nseveral tutorial examples. I did not have a chance to run the code independently.\nResponse: The code includes installation and user manuals, making the\ninstallation and use process convenient. Meanwhile, the author also provides\ncorresponding tutorials and hyperparameters, facilitating the use and repetition\nof experiments.\nReference:\n[1] M. Schlichtkrull, T. N. Kipf, P. Bloem, R. Van Den Berg, I. Titov, and M.\nWelling, \"Modeling relational data with graph convolutional networks,\" in The\nsemantic web: 15th international conference, ESWC 2018, Heraklion, Crete,\nGreece, June 3–7, 2018, proceedings 15, 2018: Springer, pp. 593-607.\n[2] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning representations\nby back-propagating errors, Nature, vol. 323, no. 6088, pp. 533-536, 1986.\n[3] K. Diederik, \"Adam: A method for stochastic optimization,\" (No Title), 2014.\n[4] A. Paszke et al., \"Automatic differentiation in pytorch,\" 2017.\n[5] Song Y, Miao Z, Brazma A, et al. Benchmarking strategies for cross-species\nintegration of single-cell RNA sequencing data. Nature Communications, 2023,\n14(1): 6495.\n[6] Xiong L, Tian K, Li Y, et al. Online single-cell data integration through\nprojecting heterogeneous datasets into a common cell-embedding space[J].\nNature Communications, 2022, 13(1): 6118.\n[7] Stuart T, Butler A, Hoffman P, et al. Comprehensive integration of single-cell\ndata[J]. cell, 2019, 177(7): 1888-1902. e21.\nResponse to Review #2\nThe authors describe CAMEX, an algorithm for cross-species integration of\nscRNA-seq data based on heterogeneous graph neural networks. The authors\napplied the method to diverse datasets and the method showed competitive\nperformance. Here, I outline some points that should be addressed to improve\nthe manuscript:\n1. The authors correctly pointed out in the introduction that some other recent\nmethods, such as SAMap and SATURN, utilise non-one-to-one homology\nbetween genes. However, there is no comparison between CAMEX and any of\nthese methods, but only with older methods which rely on one-to-one orthologs,\nwhich hardly seems fair. Given that including many-to-many relationships\nbetween genes is essential to the method, the authors should show a\nbenchmark with SAMap or SATURN.\nResponse: We sincerely thank the reviewer for this valuable comment.\nFollowing your suggestion, we compared CAMEX with several typical and\nrecently proposed single-cell integration algorithms. These include widely used\nmethods such as Seurat v4 [1] and scGen [2], and cross-species integration\nalgorithms like SAMap [4] and SATURN [5], and large language model-based\nmethods scGPT [5].\nSpecifically, we have added the results in the four sections (page 6, lines 186-\n258; page 7, lines 260-354; page 10, lines 356-411; page 11, lines 413-450)\nand Supplementary Tables S1, S5, S7, S9, S11-S22, and S24. Here, we take\nthe result of the liver dataset as an example for demonstration.\n\"In the integration task, we applied CfMEX to the datasets and compared its\nperformance against other SOhf methods, i.e., Harmony [16], PyLiger [24, 25],\nscVI [17], Scanorama [26], scGEN [27], scGPh [28], Seurat [29], SfMap [21],\nand SfhURN [30]. For SfMap and SfhURN, we conducted our experiments\nfollowing their respective tutorials (Methods).”\n“Compared to other methods, CfMEX achieved higher fdjusted Rand Index\n(fRI) and Normalized Mutual Information (NMI) scores compared to other\nalgorithms, reflecting its superior ability to cluster cells of the same type.”\n“Notably, scGen could integrate hepatocytes cells from four species while failing\nto integrate other cell types. In contrast, SfhURN, Seurat v4, Harmony, scGPh,\nScanorama, and SfMap achieved integration for all cell types except\nhepatocytes cells.”\n“hhe genes from each cell type vary considerably across species, and they\nexhibit high cell-type specificity. We also performed computational efficiency,\nhyperparameter comparison, and ablation experiments (Supplementary Figs.\nS2-S4, Supplementary Tables S2-S4). Despite variations in evaluation\nmetrics, most settings successfully isolated distinct cell types in the UMfP. fn\ninteresting observation was the clear separation of Kupffer cells and\nMacrophages when the hidden layer dimension was increased to 512\n(Supplementary Fig. S2c).”\nFigure 2. Benchmarking CAMEX and SOTA methods on the liver dataset\nacross four species. a, Overall score evaluation on the liver. The overall score\nconsiders both batch correction and biological variation conservation\nperformance. Rankings are depicted using a color gradient, where lighter colors\nindicate better performance. b, UMAP visualization of joint embedding space\non the liver dataset. Circles represent the Hepatocytes cell population from\ndifferent species. Colors represent species annotations. c, UMAP visualization\nof joint embedding space on the liver dataset. Colors represent cell type\nannotations. Squares represent the Kupffer cells and Macrophages as well as\ncholangio cell population from different species.\n2. Given that in many cases CAMEX demonstrated comparative performance\nwith other methods, but it costs much more computational power as a deep\nlearning method, authors should describe the resource consumption, such as\nMax memory and running time, in the benchmarks to help future users choose.\nResponse: We sincerely thank the reviewer for this valuable comment. We\nhave supplemented the comparison for experiments in the first four sections by\nincluding peak memory usage (CPU and GPU, if available) and running time.\nSpecifically, we compared the maximum CPU memory, maximum GPU memory\n(when applicable), and running time of the different methods.\nWe have provided the details in Supplementary Tables S2, S6, S8, S10, and\nS25. These results indicate that deep learning approaches such as CAMEX\nand scVI, while slower than traditional methods like Harmony and Scanorama,\nstill operate within an acceptable time frame. Furthermore, the CPU memory\nrequired by CAMEX during preprocessing and graph building is significant, yet\nit is lower than or comparable to that of SATURN and SAMap. Notably, due to\nits implementation of graph minibatch training, CAMEX maintains a relatively\nconstant level of GPU memory usage, ensuring its scalability.\nThe description of time consumption such as liver dataset result can be seen in\n(page 7, lines 224-226):\n“We also performed computational efficiency, hyperparameter comparison and\nablation experiments (Supplementary Figs. S2-S4, Supplementary hables S2-\nS4).”\n3. Line 130, can the authors please explain what does it mean that the initial\ncell node features “correspond to its original expression”. Aren’t gene\nexpression info in the cell-gene edges? What is the added benefit of encoding\nthese into node feature, given the very high dimensionality of gene expression,\nthis could be highly resource consuming and generate large objects.\nResponse: The expression matrix serves two key functions in the data\npreprocessing stage: constructing heterogeneous graphs and initializing node\nfeatures. The details can be seen in Methods (page 15, lines 603-620). If the\nexpression matrix is transformed solely into a graph structure (e.g., an\nadjacency matrix), the continuous expression values will be reduced to binary\ninteractions. This results in a loss of quantitative transcriptional information.\nTherefore, we use gene expression as the nodes’ initial feature.\nAs the reviewer pointed out, the expression matrix is indeed high-\ndimensional. Therefore, in the first layer of the graph neural network, we use\nthe embedding layer to reduce the input dimension (usually selected highly\nvariable genes, such as 2000 dimensions) to the hidden layer dimension\n(default is 128 dimensions), which significantly reduces the computational load\nof the subsequent graph neural network. In addition, we use graph mini-batch\nstrategy to reduce the GPU memory consumption. The number of nodes in\neach minibatch is 1024, which occupies approximately 1300MB of GPU\nmemory.\n4. In the small UMAP crops in Fig. 1b, it appears that CAMEX mixed\nmacrophages and Kupffer cells after integration, when scVI kept them\nseparated. Could the authors explain? If this is an issue, how come CAMEX\nhad a higher biological conservation score?\nResponse: This is indeed as the reviewers said. scVI has more advantages\nover CAMEX in distinguishing macrophages from kupffer cells in this dataset.\nIn addition, CAMEX can improve this result by adjusting the number of\ndimensions of the hidden layer in Supplementary Figure S2c.\nSupplementary Figure S2. Hyperparameter comparison of CAMEX on the\nliver dataset. a, Full-batch versus mini-batch comparison. b, Hidden dimension\ncomparison. c, Three experimental repeats with a hidden dimension of 512.\nIn fact, the biological conservation score of scVI is higher than that of CAMEX\n(Figure 2a). We used the scIB [6] Python package to calculate these metrics.\nFor ARI and NMI, scIB employs the Louvain algorithm to derive clustering labels\nfrom the data, using a resolution range of 0.1 to 2 (with a step of 0.1). These\ngenerated labels are then compared against the true labels to compute ARI and\nNMI. The relatively low ARI and NMI scores for scVI are attributable to its less\neffective integration of certain cell types, such as Hepatocytes and Endothelial\ncells (Figure 2b), which led to suboptimal clustering. Despite this, scVI\noutperformed CAMEX in the other three metrics related to biological variation\nconservation score. Conversely, CAMEX excelled in batch correction and\nultimately achieved a better overall score.\nThe calculation process of the metrics is as follows: Specifically, we first\ncalculate the average value of the batch score (batch mean = (kBET + PCR\nbatch + ASW Batch + Graph connectivity) / 4) and the biological score\n(biological mean = (ARI + NMI + Cell type ASW + Isolated label F1 + Isloated\nlabel silhouette) / 5), and then the overall score = (batch mean + biological mean)\n/ 2.\nFigure 2. Benchmarking CAMEX and SOTA methods on the liver dataset\nacross four species. a, Overall score evaluation on the liver. The overall score\nconsiders both batch correction and biological variation conservation\nperformance. Rankings are depicted using a color gradient, where lighter colors\nindicate better performance. b, UMAP visualization of joint embedding space\non the liver dataset. Circles represent the Hepatocytes cell population from\ndifferent species. Colors represent species annotations. c, UMAP visualization\nof joint embedding space on the liver dataset. Colors represent cell type\nannotations. Squares represent the Kupffer cells and Macrophages as well as\ncholangio cell population from different species.\nWe provide the scores of 10 methods on the liver dataset in the Supplementary\nTable S1.\n5. Line 204, does “attention” refer to attention in NNs or simply marker gene\nexpression in Fig. 2d? Please clarify and reformulate the sentence because I\ndon’t see that CAMEX HetroGNN has an attention mechanism.\nResponse: We acknowledge that this sentence may lead to misunderstanding.\nWhile CAMEX indeed utilizes a heterogeneous GAT network with attention, in\nthis context, it specifically refers to differential expression or marker genes. To\navoid ambiguity, we have reorganized the sentences accordingly (page 7, line\n222-224):\n“hhe genes from each cell type vary considerably across species, and they\nexhibit high cell-type specificity.”\n6. In benchmark UMAPs such as Fig 2b and 3b, there is no clear explanation\nabout what the circled population is referring to in the text or in the captions.\nResponse: We have added descriptions of these circles in the legend to\nfacilitate readers' clearer and easier understanding of these pictures.\n“Figure 2. Benchmarking CAMEX and SOTA methods on the liver dataset\nacross four species. a, Overall score evaluation on the liver. hhe overall score\nconsiders both batch correction and biological variation conservation\nperformance. Rankings are depicted using a color gradient, where lighter colors\nindicate better performance. b, UMfP visualization of joint embedding space\non the liver dataset. Circles represent the hepatocytes cell population from\ndifferent species. Colors represent species annotations.”\n“Figure 3. CAMEX achieved superior integration performance and\npreserved the differentiation process in the testis dataset across 11\nspecies. a, Overall score evaluated on testis and phylogenetic tree generated\nby TimeTree (v.5) (http://www.timetree.org/). b, UMAP visualization of a joint\nembedding space on the testis dataset. Colors denote species annotations (first\nrow) and cell type annotations (second row).”\n7. In line 264, I wonder what is the reason of using a KNN classifier on\nembeddings instead of doing graph attention classifier label transfer that is part\nof CAMEX tool?\nResponse: KNN classifier is not the part of CAMEX tool. Through unsupervised\nintegration, we can obtain the embeddings of cells. In addition to the\nconventional batch correction and biological variation conversation score. We\nhope to prove through other methods that CAMEX can generate high-quality\nembeddings. The sperm development dataset spanning 11 species that we\nused contains manually annotated cell types and sub cell types. We took\nhuman dataset as a reference, used the KNN-based classifier as the prediction\nmodel, annotated the cell types of other species based on the embeddings\nunsupervised by CAMEX, and achieved a relatively high accuracy rate. This\nreflects from the side that CAMEX can generate high-quality embeddings.\n8. In line 308, did the authors indeed use bulk RNA-seq? Is this method also\napplicable to bulk RNA-seq data, and could the authors explain exactly how it\nwas performed?\nResponse: Yes, we did use the bulk transcriptome data as described in the\noriginal text and reference [7]. CAMEX can be applied to bulk data. For bulk\ndata, we treat each sample as used like a single cell and use the data\npreprocessing pipeline of the CAMEX, and graph neural network training. The\nobtained embeddings are analyzed for downstream tasks.\n9. In line 327, I am not sure what the authors mean by “predicting their stages\nbeyond the specified ranges”. What is the “prediction”, and what are the\n“specified ranges”? Please specify and explain. This mention of “predicting\nbeyond Carnegie stages” appeared a few times in the manuscript as a core\nresult and could greatly benefit from more background and a clearer\nexplanation.\nResponse: Thanks for your comment. After integration, we can obtain the\nembedding of each sample. “Prediction” means that we align the developmental\nstages of different samples through the similarity of embedding and transfer the\ndevelopment stage label. Model organisms are often extensively studied and\ntheir developmental stages are relatively well-defined and are named after the\nCarnegie Stages. However, the developmental stages of some non-model\norganisms are uncertain. By aligning, we hope to be able to transfer the\ndevelopmental stage labels from one well studied organism to unfamiliar\norganisms. This is called “predictions”.\nThe ambiguity in these statements arose from our insufficient description of the\nCarnegie stages. To clarify, we have further introduced the motivation of the\ndevelopment stage label prediction (page 10, lines 374-383):\n“While most model organisms display distinct developmental timelines, these\nare typically mapped to Carnegie stages, spanning from gametogenesis\nthrough birth. hhe Carnegie stages themselves are fundamentally based on\nobservable morphological development during embryogenesis. Crucially,\nthey've only been defined by researchers for a handful of model organisms (e.g.,\nhuman and mouse) and are limited to the period immediately following\npregnancy and childbirth. fn exciting challenge lies in exploring whether we\ncan expand the Carnegie developmental staging concept beyond this limited\nset of organisms to encompass any species, or even to predict any point in a\nspecies' lifespan through integrated approaches.”\nand the definition of Carnegie stages in the main text (page 22, lines 870-878):\n“Carnegie stages represent a standardized system used to define the 23 stages\nof embryonic development in vertebrates. hhese stages track development\nfrom the formation of a fertilized egg to an embryo's progress over\napproximately 60 days (8 weeks). For instance, a human embryo takes around\n60 days to reach Carnegie Stage 23, whereas a chicken embryo completes the\nsame stage in just 10 days. hhis system enables scientists to compare\ndevelopmental characteristics across species at equivalent stages, thereby\nenhancing our understanding of the universal principles governing embryonic\ndevelopment.”\n10. In line 368, it's unclear which population in Fig. 5c-d is the “unknown”\npopulation, and I find it hard to understand why these are “unknown” but\n“originally defined as neural progenitor cells”. Please rewrite this paragraph to\nprovide clarity.\nResponse: This description of manuscript is indeed not clear enough.\nSpecifically, since the reference dataset (training set) contains a limited number\nof cell types, the cell types sometimes cannot include all the cell types in the\nquery dataset. For example, here humans are taken as a reference. However,\nthe human dataset does not contain neural progenitor cells, while the lizard\ndataset used as the query dataset contains neural progenitor cells\n(Supplementary Figure S31c and below). At this point, we define the neural\nprogenitor cell as unknown. As requested by the reviewers, we rewrote this\nparagraph to facilitate readers' understanding (page 12, line 443-445):\n“Unknown is the cell type dost not appear in reference but appear in query\ndataset. NPC dost not appear in reference human dataset but appear in query\nturtle and lizard dataset, so we defined them as unknown.”\nSupplementary Figure S31. Benchmarking CAMEX and baseline methods on\nthe cortex dataset across four species. a, Cortex scRNA-seq dataset includes\nsamples from humans, mice, lizards, and turtles. b, Number of detected UMIs\nand genes in different species. c, Each dataset comprises distinct cell types\nthat are unique to that respective species.\n11. The last section on discovering new populations generally seems to lack\nproof. All results are computational hypotheses, but the claims of discovering\nnew populations are very strong. In Fig. 6b, the GLDN populations in Chimp\nand Rhesus are highly likely artefacts due to the decision boundary of the GAT\nclassifier. Any small tweak in the decision boundaries can lead to these cells\nbeing assigned to one or another cell type, and there is no proof that these exist\nin Chimp and Rhesus. Furthermore, there is no visible CCL3 population as\nthere are no purple dots in other species, so I don’t know where the CCL3\npopulation is that the authors claim.\nResponse: We completely agree with the reviewer that our initial claims of\n‘discovering’ new cell populations were too strong, and identifying novel\nsubtypes from UMAP plots alone is insufficient. To address this, we have\nprovided a more detailed analysis of the expression patterns for GLDN, CCL3,\nand fPBB1IP. We compared our CAMEX predictions with the original\nliterature's manual annotations (see Supplementary Figure S36 and below).\nOur results show that CAMEX obtained some interesting findings, such as the\nexpression of CCL3 in some cell populations in chimpanzees. However, we also\nacknowledge it can produce some false positives, such as the predicted\nexpression of CCL3 and GLDN in marmosets.\nOverall, we clarify that CAMEX should be considered a bioinformatics tool\nfor generating hypotheses about potential cell subpopulations. These\ncomputational predictions must be validated through subsequent wet-lab\nexperiments. Accordingly, we toned down our conclusion and modified related\nsentences (page13, lines 486-494):\n“ft the same time, we investigated the expression of GLDN, CCL3, and\nfPBB1IP in the original manual annotation and the CfMEX prediction\n(Supplementary Fig. S36c). We found that the prediction of CfMEX has\ncertain rationality, for instance, the expression of CCL3 in the subgroups of\nchimpanzees. However, there are also false positives, such as CCL3 and\nGLDN, in marmosets. Overall, we therefore wish to clarify that CfMEX should\nbe considered a bioinformatics tool for generating hypotheses about potential\ncell subpopulations. hhese computational predictions must be validated\nthrough subsequent wet-lab experiments.”\nSupplementary Figure S36. Integrating the microglia cell datasets of four\nspecies by CAMEX. c, Violin plot of GLDN, CCL3, and fPBB1IP expression\nin original manual annotation and CAMEX prediction clusters.\n12. Similarly, the proportion of cell types is likely a computational artefact given\nthe small data size and obscure cell type annotation transfer boundaries. Therefore,\nthe discussions in lines 420-423 cannot be justified. These claims have almost no\nrelation to current evolutionary knowledge and cannot be considered reasonable.\nI suggest removing them.\nResponse: We sincerely thank the reviewer for this valuable comment.\nFollowing your suggestion, we have removed them.\nReviewer #2 (Remarks on code availability):\nCodes are available with limited documentation. There are some notebooks,\nhowever, there are no annotation transfer examples, as all notebooks end at\ndata integration. Given that annotation transfer is a key part of the method, it is\nessential to include at least one example.\nResponse: According to the requirements of the reviewers, we improved the\ntutorial and provided detailed steps for annotation transfer with two examples.\nThe URLs are:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/4cortex_annotati\non/integration_annotation_in_relatives_distant_species.ipynb\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/5micro_mapping\n/discovery_new_populations_markers.ipynb\nReference:\n[1] Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W.M., Zheng, S., Butler, A.,\nLee, M.J., Wilk, A.J., Darby, C. and Zager, M. (2021) Integrated analysis of\nmultimodal single-cell data. Cell, 184, 3573–3587.\n[2] Lotfollahi, M., Wolf, F.A. and Theis, F.J. (2019) scGen predicts single-cell\nperturbation responses. Nat. Methods, 16, 715–721.\n[3] Wang J, Ma A, Chang Y, et al. (2021) scGNN is a novel graph neural network\nframework for single-cell RNA-Seq analyses. Nature Communications, 12(1):\n1882.\n[4] Tarashansky, A.J., Musser, J.M., Khariton, M., Li, P., Arendt, D., Quake, S.R.\nand Wang, B. (2021) Mapping single-cell atlases throughout Metazoa unravels\ncell type evolution. eLife, 10, e66747.\n[5] Rosen, Y., Brbic, M., Roohani, Y., Swanson, K., Li, Z. and Leskovec, J. (2024)\nToward universal cell embeddings: integrating single-cell RNA-seq datasets\nacross species with SATURN. Nat. Methods, 21, 1492–1500.\n[6] Luecken, M.D., Büttner, M., Chaichoompu, K., et al. (2022) Benchmarking\natlas-level data integration in single-cell genomics. Nature methods, 19(1): 41-\n50.\n[7] Cardoso-Moreira, M., Halbert, J., Valloton, D., et al. (2019) Gene expression\nacross mammalian organ development. Nature, 571(7766): 505-509.\nResponse to Review #3\nThe authors proposed CAMEX, a heterogeneous Graph Neural Network tool,\nto address these challenges by leveraging many-to-many homologous\nrelationships for multi-species scRNA-seq data integration, alignment, and\nannotation. It outperforms other methods on cross-species datasets, helps\nalign species at different developmental stages, enables the detection of\nspecies-specific cell types and marker genes, and has the potential to offer\ninsights into evolutionary forces across species at single-cell resolution.\n1. What are the noteworthy results?\nCAMEX outperforms state-of-the-art methods in cross-species integration\nacross datasets spanning various species. It enables the alignment of\ndevelopmental stages and the detection of species-specific cell types and\nmarker genes. The tool leverages many-to-many homologous relationships,\naddressing a key challenge in cross-species integration.\nHowever, there are few methods compared in the results, and further\ncomparison methods need to be added, such as Seurat v4 [1], scGen [2],\nSAMap [3], SATURN [4], etc. At the same time, are the results of the model\nadvantageous compared to the popular large language model-based methods\nthat have been compared recently? The authors are requested to supplement\nthe corresponding experiments.\nResponse: We appreciate your kind recognition on this study and thanks for\nthe valuate suggestion. Following your suggestion, we have compared CAMEX\nwith several typical and recently proposed single-cell integration algorithms.\nThese include widely used methods such as Seurat v4 [1] and scGen [2], and\ncross-species integration algorithms like SAMap [3] and SATURN [4], and large\nlanguage model-based methods scGPT [5]. Specifically, we have added the\nresults in the four sections and related Supplement Figures (page 6, lines\n186-258; page 7, lines 260-354; page 10, line 356-411; page 11, lines 413-\n450) and Supplementary Tables S1, S5, S7, S9, S11-S22, and S24. Here, we\ntake the result of liver dataset as an example for demonstration.\nFigure 2. Benchmarking CAMEX and SOTA methods on the liver dataset\nacross four species. a, Overall score evaluation on the liver. The overall score\nconsiders both batch correction and biological variation conservation\nperformance. Rankings are depicted using a color gradient, where lighter colors\nindicate better performance. b, UMAP visualization of joint embedding space\non the liver dataset. Circles represent the hepatocyte cell population from\ndifferent species. Colors represent species annotations. c, UMAP visualization\nof joint embedding space on the liver dataset. Colors represent cell type\nannotations. Squares represent the Kupffer cells and Macrophages, as well as\nthe cholangiocyte cell population from different species.\nIn addition, if the large language models referenced by the reviewers pertain\nto tools like ChatGPT or DeepSeek for annotating cell type by differentially\nexpressed genes, we believe this is a highly promising direction, with the\ncommunity actively exploring this area. However, we don't have sufficient\ncomputational resources, especially GPUs with large memory, to run these\nlarge language models. Accordingly, we have expanded the discussion section\nto include the potential of large language models in single-cell data annotation,\nparticularly the advantages and disadvantages in cross-species scenarios\n(page 14, lines 583-599):\n“In addition, Large Language Models (LLMs) are quickly being adopted in\nsingle-cell data analysis and will bring novel insight for cross-species annotation.\nhhe methods can be roughly divided into several categories: Marker-based\nmethods leverage the natural language understanding of LLMs to interpret lists\nof differentially expressed (marker) genes from single-cell cluster analysis,\ninferring cell types. hext sequence conversion methods transform high-\ndimensional single-cell gene expression data into text sequences. hools like\nscBERh exemplify this for annotation. hhe LLM fgent framework (e.g., scfgent)\ncreates autonomous agents with planning and memory modules to\nautomatically annotate scRNf-seq data through multi-round interactions.\nDespite their promise, data availability in cross-species studies remains limited,\nposing challenges for robust computational analysis. Key concerns include:\nEnsuring the biological accuracy and consistency of LLM-generated\nannotations and preventing hallucinations. Managing the computational costs\nand resource demands of deploying and running large LLMs. Verifying the\nbiological significance and reliability of LLM explanations to build researcher\ntrust in the results.”\n2. Will the work be of significance to the field and related fields? How does it\ncompare to the established literature? If the work is not original, please provide\nrelevant references.\nMulti-species scRNA-seq integration is a critical area for evolutionary biology\nand comparative genomics. CAMEX’s use of heterogeneous GNNs to model\nmany-to-many homology is novel and addresses limitations of existing tools\n(e.g., Seurat, LIGER, or SCOT), which often assume one-to-one orthology.\nHowever, this paper does not compare CAMEX with previous GNN-based\nmethods, especially with the CAME model proposed by the same authors\nbefore. The author should clearly state the differences between CAMEX and\nCAME to prove that it has made an innovative contribution sufficiently.\nMeanwhile, in the experimental part, the author also does not compare with the\nresults of CAME. Please further explain the reasons for this.\nResponse: Thanks for kind comment. Our previous method CAME is a semi-\nsupervised annotation algorithm specifically designed for cross-species\nanalysis between two species. CAME necessitates labeled reference data and\ndoes not support unsupervised integration of single-cell data from multiple\nspecies. In contrast, CAMEX enables unsupervised integration across any\nnumber of species without requiring labels. Furthermore, when labeled\nreference data is available, CAMEX can also transfer these labels to the query\ndata.\nWe have further highlighted and emphasized the distinctions between CAME\nand CAMEX in the Introduction (page 3, lines 80-89) and Discussion\nsections (page 14, lines 538-546):\nIntroduction:\n“SfMfP relies on sequence alignment, a computationally intensive process,\nand is primarily used for visualization purposes rather than generating\nembeddings for downstream tasks. SfhURN measures inter-species gene\nsimilarity using protein-language models. But SfhURN faces challenges from\nproteome-transcriptome mismatches (notably from alternative splicing events)\nand requires cell annotation inputs that might limit its application. Our previous\nannotation method, CfME, utilizes many-to-many homologous gene\nrelationships, restricting to pairwise species alignment (two species) within a\nshared embedding space. CfME cannot perform integration without cell labels\nacross multiple species datasets.”\nDiscussion:\n“ho our knowledge, SfMap pioneered this area, but its reliance on sequence\nalignment leads to significant computational time consumption. Furthermore,\nSfMap's CPU memory usage increases sharply with a large number of species,\nposing a practical limitation. SfhURN, while an advanced algorithm that\ncombines protein language models with transcriptional expression, requires cell\ntype annotation as input. hhis prerequisite significantly limits its broad\napplication. Our previous method, CfME, is a semi-supervised model that\nrelies on reference cell labels for annotation and lacks the capability to perform\nintegration tasks.”\n3. Does the work support the conclusions and claims, or is additional evidence\nneeded?\nThe authors have carried out some experiments on cross-species data sets.\nHowever, the distance relationships between different species were not\ndistinguished in the relevant cross-species experiments. It is recommended\nthat more detailed experiments be conducted in reference [7] to verify the\neffectiveness of the model more comprehensively.\nResponse: We carefully read the method in reference [7]. To comprehensively\nevaluate CAMEX's performance with a large evolutionary distance and an\nincreasing number of species, we designed two types of integration\nexperiments in the testis dataset across 11 species. Our experimental\nframework included:\n(1) Pairwise integration analyses between humans and each target species,\ncovering a broad phylogenetic spectrum from closely related primates to distant\nvertebrates.\n(2) Multi-species integration scenarios with progressively expanded species\nsets, systematically increasing from 2 to 11 species.\nWe added the Supplementary Figures S12-S24 to clarify this. Furthermore, the\nresults of experiments conducted by CAMEX and other baseline methods are\nsummarized in (page18 line 294-311 and below) and in the figure\n(Supplementary Figures S22 and S24 and below). In conclusion, CAMEX is\ninsensitive and robust to the number of shared homologous genes and the\ndistance between species.\n“Inspired by the previous literature, to comprehensively evaluate CfMEX's\nperformance across varying evolutionary distances and increasing numbers of\nspecies, we designed two integration scenarios. First, we conducted pairwise\nintegration analyses between human and each target species, covering a broad\nphylogenetic spectrum from closely related primates to distant vertebrates. hhe\nresults clearly show CfMEX excelled in batch correction, while SfhURN\nperformed strongly in biological variation conservation. Overall, CfMEX\ndemonstrated an advantage over other methods in both quantitative metrics,\nrankings, and qualitative UMfP visualization. Second, we created multi-species\nintegration scenarios with progressively expanded species sets, systematically\nincreasing from two to eleven species. CfMEX also consistently maintained\nhigh performance across all scenarios (Supplementary Figs. S12-S24,\nSupplementary Tables S11-S21). Overall, most integration algorithms,\nparticularly those relying solely on one-to-one homologous genes, show\ndiminished performance as the evolutionary distance between species\nincreases. However, methods like CfMEX, which leverage many-to-many\nhomologous genes, and ShfURN, which utilizes sequence similarity, are less\nimpacted by this increasing species divergence.”\nSupplementary Figure S22. Overall integration performance achieved by\nvarious methods when applied to a two-species testis data, allowing for\ncomparative evaluation. The line graph illustrates the integration performance\nof humans with other species. Different colors represent different methods.\nThree distinct lines represent: Batch correction score, Biological variation\nconservation (Bio) score, and Overall score. For the line charts on the left, a\nhigher value indicates better performance. Conversely, for the box plots on the\nright, a lower ranking signifies better performance.\nSupplementary Figure S24. Overall integration performance achieved by\nCAMEX when applied to the testis dataset, ranging from two-species to eleven-\nspecies integrations, allowing for comparative evaluation. Three distinct lines\nrepresent key performance metrics: the Batch correction score, the Biological\nvariation conservation (Bio) score, and the Overall score.\nAccording to reference [7], we discussed the relationship between distance and\nintegration effect and compared the differences in results among different\nalgorithms. Overall, CAMEX is insensitive and robust to the number of shared\nhomologous genes and the distance between species. We added the details to\nthe supplement section (page 8, lines 294-311):\n“Inspired by the previous literature [7], to comprehensively evaluate CfMEX's\nperformance across varying evolutionary distances and increasing numbers of\nspecies, we designed two integration scenarios. First, we conducted pairwise\nintegration analyses between humans and each target species, covering a\nbroad phylogenetic spectrum from closely related primates to distant\nvertebrates. hhe results clearly show that CfMEX excelled in batch correction,\nwhile SfhURN performed better in biological variation conservation. Overall,\nCfMEX demonstrated an advantage over other methods in quantitative metrics,\nrankings, and qualitative UMfP visualization. Second, we created multi-species\nintegration scenarios with progressively expanded species sets, systematically\nincreasing from two to eleven species. CfMEX also consistently maintained\nhigh performance across all scenarios (Supplementary Figs. S12-S24,\nSupplementary Tables S11-S21). Overall, most integration algorithms,\nparticularly those relying solely on one-to-one homologous genes, show\ndiminished performance as the evolutionary distance between species\nincreases. However, methods like CfMEX, which leverage many-to-many\nhomologous genes, and ShfURN, which utilizes sequence similarity, are less\nimpacted by this increasing species divergence.”\n4. Are there any flaws in the data analysis, interpretation and conclusions? Do\nthese prohibit publication or require revision?\n1). When constructing homology relationships, how are homology relationships\ndefined? Are homology databases or sequence alignment methods used?\n2). How is the data partitioned during the neural network training process? Is\nthe partition completely random? Or is it based on some experience?\n3). There are relatively few methods for conducting comparative experiments.\nThe authors are requested to further supplement the comparative methods\nexperiments to verify the model performance's superiority.\nBefore the manuscript is published, the authors need to clarify and revise these\nissues.\nResponse: 1). The homology relationships were obtained from the Ensembl\ndatabase [9]. We have enhanced the data section to include details on how to\nobtain one-to-one and many-to-many homologous genes (page 22, lines 868-\n878). Additionally, examples and usage instructions have been provided on\nGitHub for reference.\n“hhe many-to-many homology relationships were obtained from the Ensembl\ndatabase [7]. hhe details of how to use it can be found in the GitHub tutorials.”\n2). In the integration task, none of the datasets have labels. This is an\nunsupervised (self-supervised) learning task, so there is no need to divide the\ndata set. The main results can be seen in the Results sections 2, 3, and 4. For\nthe annotation task, a reference dataset with labels (comprising data from a\nspecific species) was employed as the training set to annotate query datasets\nwithout labels (from other species). The results can be seen in section “CAMEX\ncould achieve more accurate integration and annotation performance in both\nrelatives and distant species” (page 11, lines 436-438):\n“In this scenario, we used the human dataset as the reference (with cell labels)\nand applied CfMEX to annotate the other three datasets (without cell labels),\nparticularly those from lizard and turtle samples.”\n3). As mentioned previously in Response to Comment 1, we have expanded\nour study with new comparative experiments. In addition to the existing baseline\nmethods (scVI, pyliger, Harmony, and Scanorama), we have now included\nSeurat v4, scGen, SAMap, SATURN, and scGPT in the method comparison.\n5. Is the methodology sound? Does the work meet the expected standards in\nyour field?\nCAMEX’s use of heterogeneous GNNs to model many-to-many homology is\nnovel and addresses limitations of existing tools, which often assume one-to-\none orthology. However, specific experiments, such as ablation experiments,\nare lacking in verifying which modules in the proposed model play a crucial role\nduring the experimental process.\nThe paper's results align with the expectations in the relevant fields. However,\nit is necessary to supplement the experiments of comparison methods further,\nespecially those proposed in recent years, to verify the advancement of the\nmethods.\nResponse: (1) Thanks for your valuable comment. According to your\nsuggestion, we have added new ablation experiments to demonstrate the\neffectiveness of key modules. These include replacing different kinds of edges\nduring network construction. Furthermore, we conducted multiple\nhyperparameter comparison tests, such as varying the initial number of features,\nthe number of GNN layers, and the dimensionality. We added the revisions in\nthe main/supplementary text (page 7, lines 224-226, Supplementary Figs. S2,\nS3, S4, and Table S3):\n“We also performed computational efficiency, hyperparameter comparison, and\nablation experiments (Supplementary Figs. S2-S4, Supplementary Tables\nS2-S4). Despite variations in evaluation metrics, most settings successfully\nisolated distinct cell types in the UMfP. fn interesting observation was the clear\nseparation of Kupffer cells and Macrophages when the hidden layer dimension\nwas increased to 512 (Supplementary Fig. S2c).”\n(2) For baseline methods, we have provided a more detailed experiment\ncompared with the latest cross-species integration algorithm. The details can\nbe seen in Response to comment 1.\nSupplementary Figure S2. Hyperparameter comparison of CAMEX on the\nliver dataset. a, Full-batch versus mini-batch comparison. b, Hidden dimension\ncomparison. c, Three experimental repeats with a hidden dimension of 512.\nSupplementary Figure S3. Hyperparameter comparison of CAMEX on the\nliver dataset. a, Illustration of the effect of batch size during mini-batch training.\nb, Evaluation of the influence of GNN layer depth within the encoder. c,\nComparison of the impact of different training epoch numbers.\nSupplementary Figure S4. Ablation study of edge contributions in CAMEX.\nPerformance impact of removing individual edge types from the original six-\nedge configuration.\n6. Is there enough detail provided in the methods for the work to be reproduced?\nThe author provides the corresponding code and data source and gives\nrelevant hyper-parameter information to replicate the work.\nResponse: We appreciate your careful review and inspection.\nReviewer #3 (Remarks on code availability):\nThe code provides installation and user manuals, which make the installation\nand use convenient. Meanwhile, the author also gives corresponding tutorials\nand hyperparameters, facilitating users to use or repeat experiments.\nResponse: We sincerely appreciate your kind help and insightful comments.\nReference:\n[1] Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W.M., Zheng, S., Butler, A.,\nLee, M.J., Wilk, A.J., Darby, C. and Zager, M. (2021) Integrated analysis of\nmultimodal single-cell data. Cell, 184, 3573–3587.\n[2] Lotfollahi, M., Wolf, F.A. and Theis, F.J. (2019) scGen predicts single-cell\nperturbation responses. Nat. Methods, 16, 715–721.\n[3] Tarashansky, A.J., Musser, J.M., Khariton, M., Li, P., Arendt, D., Quake, S.R.\nand Wang, B. (2021) Mapping single-cell atlases throughout Metazoa unravels\ncell type evolution. eLife, 10, e66747.\n[4] Rosen, Y., Brbic, M., Roohani, Y., Swanson, K., Li, Z. and Leskovec, J. (2024)\nToward universal cell embeddings: integrating single-cell RNA-seq datasets\nacross species with SATURN. Nat. Methods, 21, 1492–1500.\n[5] Cui H, Wang C, Maan H, et al. scGPT: toward building a foundation model\nfor single-cell multi-omics using generative AI. Nature Methods, 2024, 21(8):\n1470-1480.\n[6] Zhong, Huawen, et al. \"Benchmarking cross-species single-cell RNA-seq\ndata integration methods: towards a cell type tree of life.\" Nucleic Acids\nResearch 53.1 (2025): gkae1316.\n[7] Harrison P W, Amode M R, Austine-Orimoloye O, et al. Ensembl 2024.\nNucleic acids research, 2024, 52(D1): D891-D899.\nRebuttal Summary\nResponse to Review #2\nI appreciate the authors’ response to my comments and their efforts to improve\nthe manuscript. The manuscript has been greatly improved. I only have a few\nadditional comments:\nResponse: We appreciate your kind recognition and detailed comments.\nFollowing your suggestions, we have extensively revised our manuscript to\nimprove its quality.\n1. Line 227 - “The genes from each cell type vary considerably across species,\nand they exhibit high cell-type specificity.” Please explain what genes the\nauthors are referring to here?\nResponse: Here, we refer to the genes shown in Supplementary Figure S1f,\nwhich have been identified as marker genes of Macrophage and Kupffer cells.\nSpecifically, we searched CellMarker 2.0 [1] and verified the results in the\noriginal literature. For example, APOC1 and APOE are identified as\nMacrophage cell marker genes in humans [2], and CD5L is a marker gene of\nKupffer cells in mice [3]. Following your comments, we have cited the\nsupplementary subfigure in the related sentences to clarify this (page 7, lines\n235-236):\n“The genes from each cell type vary considerably across species\n(Supplementary Fig. S1f).”\nSupplementary Figure S1. f, Marker genes of Macrophages and Kupffer cells\nin different species, which are searched in CellMarker2.0 [1] and verified in the\noriginal literature. For example, APOC1 and APOE are identified as the\nMacrophage cell marker genes in humans [2], and CD5L is a marker gene of\nthe Kupffer cell in mice [3].\n[1] Hu C, Li T, Xu Y, et al. CellMarker 2.0: an updated database of manually\ncurated cell markers in human/mouse and web tools based on scRNA-seq\ndata[J]. Nucleic acids research, 2023, 51(D1): D870-D876.\n[2] Zhao J, Zhang S, Liu Y, et al. Single-cell RNA sequencing reveals the\nheterogeneity of liver-resident immune cells in human[J]. Cell discovery, 2020,\n6(1): 22.\n[3] Wang J, Hu W, Shen Z, et al. Dissecting the single-cell transcriptome\nunderlying chronic liver injury[J]. Molecular Therapy Nucleic Acids, 2021, 26:\n1364-1373.\n2. It would be good to explicitly mention WHEN the semi-supervised GAT-based\nannotation transfer functionality of CAMEX is used for various work shown in\nthis manuscript. Line 436 seems to be an application of the GAT, but it came\nafter the KNN annotation in line 320. The author’s rationale for using KNN to\nshow the embedding quality is reasonable, which can be explicitly mentioned\nfor line 320. Along the same line, please add code annotation to explain line In\n[58] of the notebook:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/4cortex_annotati\non/integration_annotation_in_relatives_distant_species.ipynb.\nResponse: Yes, as the reviewers pointed out, we employed the CAMEX\nintegration in the corresponding sections of Figures 2, 3, and 4. We used both\nthe CAMEX integration and annotations in the corresponding sections of\nFigures 5 and 6. When we use the terms “integration” or “alignment”, it indicates\nthat we are conducting an unsupervised integration task. When using the term\n“annotation”, it suggests that we are conducting a semi-supervised label\ntransfer task. We clarified this point in the \"Baseline methods\" section (page 20,\nlines 811-816):\n“It is crucial to note that the benchmarking strategy is task-dependent. For the\nunsupervised integration and alignment task, CAMEX is compared against\nlabel-free methods such as Harmony. For the semi-supervised cell-type\nannotation task, however, the comparison is made against reference-based\nmethods like SingleCellNet (SCN), which require at least one labeled reference\ndataset.”\nFurthermore, we explicitly clarify the task in each section. For example, in the\nResult section 2 “CAMEX achieves superior integration performance in different\nmulti-species scenarios”, we added the sentence (page 6, lines 204-205): “In\nthis unsupervised integration setting, CAMEX operates without the need for\nmanually annotated cell type labels.\nIn the Result section 3 “CAMEX uncovers the conserved differentiation process\nin the testis across 11 species”, we added the sentence (page 8, lines 290-\n292): “In this unsupervised integration setting, both CAMEX and the baseline\nmethods are precluded from using any manually annotated cell type labels.”\nIn the Result section 4 “CAMEX aligns various development stages of seven\norgans across seven different species”, we added the sentence (page 10, lines\n387-388): In this unsupervised integration setting, CAMEX does not require any\nmanually annotated cell type labels.\nIn the Result section 5 “CAMEX could achieve more accurate integration and\nannotation performance in both relatives and distant species”, we added the\nsentences (page 12, lines 449-450): “First, we applied CAMEX and baseline methods to\nintegrate the scRNA-seq data in the absence of cell type annotations”.\n(page 12, lines 465-467): “In this semi-supervised scenario, the human dataset\nserved as the labeled reference. We then applied CAMEX to annotate the three\nunlabeled target datasets, with a focus on the lizard and turtle species.”\nIn the Result section 6 “CAMEX facilitates the discovery of new populations and\nmarkers in the primate granular dorsolateral prefrontal cortex (DLPFC) data”,\nwe added the sentence (page 13, lines 491-492): “Firstly, we integrated all the\ncells in the DLPFC dataset among the four species by CAMEX in an\nunsupervised method.”\n(page 13, lines 503-505): “Here, we use the cells in humans as the reference\n(training set) to annotate cells in other species (test set) in a semi-supervised\nway by using a GAT classifier (Method).”\nIn addition, we have added a sentence to clarify KNN (page 9, lines 336-337):\n“Secondly, … to assess the quality of cell embedding in the CAMEX integration\ntask.”\nIn response to your request, we have provided code annotations on each line\nof code in the corresponding sections of the Jupyter notebook:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/4cortex_annotati\non/integration_annotation_in_relatives_distant_species.ipynb\n3. Please specify which dates are included in the Carnegie stages for the\nmouse and rat (9/10.5 - 16/17.5) in line 394, so readers can understand what\n‘beyond the specified ranges’ means. It will also help to mark the range within\nCarnegie stages on Fig. 4b.\nResponse: Following your suggestion, we added \"mouse and rat (9/10.5 -\n16/17.5)\" to the sentence (page 11, lines 414-415):\n“CAMEX aligned the developmental stages of mouse and rat (9/10.5 – 16/17.5)\nwithin the Carnegie stages.”\nSecondly, in Fig. 4b and Supplementary Figure S30, we have marked the\nexisting range of the Carnegie stages. In the legend of Fig. 4b and\nSupplementary Figure S30, we explain that \"ek\" and \"pk\" denote the k-th\nembryonic day and postnatal day, respectively.\nFigure 4b. Similarities among all organs in mice and rats, with the mouse as a\nreference, ensuring their developmental times are well aligned. The Carnegie\nstages in mouse and rat correspond to embryonic days 9-16 and 10.5-17.5,\nrespectively. Here, \"ek\" and \"pk\" denote the k-th embryonic day and postnatal\nday, respectively.\nSupplementary Figure S30b. Similarities among the brains of mice and\nopossums, with the mouse as a reference, ensuring their developmental times\nare well aligned. The opossum central nervous system, including the brain, is\nmarkedly premature at birth (p0), resembling the embryonic e11 mice. The\nCarnegie stages in the mouse correspond to gestational days 9-16, while the\ncorresponding stages for the opossum have not been established. Here, \"ek\"\nand \"pk\" denote the k-th embryonic day and postnatal day, respectively.\nFinally, in the methods section, we clarify that \"Carnegie stages\" refers to the\nkth embryonic day (page 23, lines 912-914):\n“The data on ‘the number of days’ in Carnegie stages refers to the days\nfollowing embryo formation and should be considered an approximate value.”\n4. Line 437 - this description of the current limitation of cell type annotation\nmethods is confusing - what single group?\nResponse: Thanks for your reminder. Here, we mean to conduct cell\nannotation based on clustering and marker genes. We have revised this\nsentence (page 12, lines 460-462):\n“Traditional cell-type annotation methods are heavily dependent on clustering\nand marker genes, posing a significant challenge in non-model species where\nsuch prior knowledge is insufficient.”\n5. Line 429: I wonder what “All the methods confused the excitatory and\ninhibitory neurons” means, and which figures show this confusion? Line 442 -\nIf all methods confused excitatory and inhibitory neurons, how can they be\nannotated?\nResponse: According to the second comment from the reviewer. We have once\nagain clarified that in this section, two tasks - integration and annotation - have\nbeen carried out separately.\n(page 12, lines 449-450): “First, we applied CAMEX and baseline methods to\nintegrate the scRNA-seq data in the absence of cell type annotations”.\n(page 12, lines 465-467): “In this semi-supervised scenario, the human dataset\nserved as the labeled reference. We then applied CAMEX to annotate the three\nunlabeled target datasets, with a focus on the lizard and turtle species.”\nIn addition, as the reviewer said, Line 429 “All the methods confused the\nexcitatory and inhibitory neurons” means the results of the integration task\n(Supplementary Fig. S32). In the UMAP plot, we believe that the excitatory\nand inhibitory neurons resulting from all methods are not perfect. Line 442\nrefers to the annotation task. We have revised this sentence (page 12, lines\n453-455):\n“From the UMAP plot (Supplementary Fig. S32), we can see that all the\nmethods confused the excitatory and inhibitory neurons.”\n6. I understand the author’s addition of discussion about LLM models. However,\nit seems quite descriptive for a discussion section and disconnected from the\nmain idea of this work; therefore, it is not a good ending to the discussion.\nPerhaps move the current text to the introduction and shorten it, and consider\ncomparing CAMEX with scGPT-based approaches in discussions.\nResponse: We agreed with your comment and shortened it. And we moved it\nto the introduction section and streamlined it (page 3, lines 80-83):\n“Furthermore, inspired by the rapid advancement of large models in natural\nlanguage processing, analogous foundation models have emerged in the\nsingle-cell field, such as scGPT [19], which facilitates downstream tasks like\nbatch effect removal.”\nIn addition, we compared CAMEX and scGPT in each section. For example,\n7. In line 616, it would be good to clarify whether the PCA space used for KNN\nis from a standard scRNA-seq dimensionality reduction analysis.\nResponse: Yes, as you noted, we followed the common scRNA-seq analysis\npipeline, which includes normalization, log-transformation, and selection of\nhighly variable genes, followed by PCA and the application of the KNN\nalgorithm. According to the reviewers' suggestions, we have provided a more\ndetailed description of this section (page 16, lines 631-638):\n“A scRNA-seq dataset comprising 𝑁 cells and 𝑀 genes is represented as a\nmatrix 𝑿 = (𝒙 ,𝒙 ,…,𝒙 )𝑇 ∈ 𝑅𝑁×𝑀. Following a standard preprocessing\n1 2 𝑁\npipeline (normalization, log-transformation, and selection of highly variable\ngenes, clustering, and obtaining differential genes), we construct a cell-gene\nbipartite graph 𝑨. In this graph, an edge connects cell 𝑛 and gene 𝑚 if and\nonly if cell 𝑛 expresses gene 𝑚. Subsequently, cell-cell connections are\nestablished by applying the kNN algorithm within the PCA space derived from\nthis processed matrix.”\nResponse to Review #1\nThe editor also asked me to comment on behalf of Reviewer #1 regarding\ntheir comments.\nResponse: We appreciate your time and effort in reviewing our manuscript.\n1. How many layers are there in the heterogeneous graph neural network\nencoder? Please state this explicitly.\nResponse: We employed a two-layer heterogeneous GNN network by default\nand used this default parameter in all our experiments. We have clarified this in\nthe methods section (page 19, lines 764-771):\n“(4) Hyperparameters:\nCAMEX employs a set of default hyperparameters, including: train_mode:\n'mini_batch', dim_hidden: 128, batch_size: 1024, gnn_layer_num: 2,\nepoch_integration: 10, epoch_annotation: 10, and encoder: 'GCN'. All\nexperiments except for bulk RNA-seq integration used these defaults. For bulk\nRNA-seq data with a small sample size, we recommend setting the\ntrain_mode to ‘full_batch’ and setting epoch_integration to 200. Further details\nare provided in the GitHub tutorials.”\n8. In Luecken et al. (Nature Methods, 2021), pseudotime was calculated before\nand after integration using DPT, and Spearman’s rank correlation between the\ntwo pseudotime vectors was used as a trajectory conservation score. This\napproach is already established and very similar to the author’s approach;\nhowever, it is important to note that this correlation-based measure is valid\nbecause pseudotime is a relative ordering, not an absolute metric. Spearman’s\ncorrelation compares ranks, effectively checking whether the ordering of cells\nalong the trajectory remains consistent after integration. Therefore, the author\nis advised to remove Pearson’s method for calculating pseudotime\nconservation and credit Luecken et al., if the methods are the same.\nResponse: Thanks for your reminder. We have cited this paper and removed\nPearson’s correlation (page 22, lines 881-883).\n“We applied and adjusted a quantitative metric to evaluate trajectory\npreservation during integration from Luecken et al. [1]. This metric is calculated\nin the following steps:\nSpearman’s rank correlation is used for this quantification.”\n[1] Luecken M D, Büttner M, Chaichoompu K, et al. Benchmarking atlas-level\ndata integration in single-cell genomics[J]. Nature methods, 2022, 19(1): 41-50.\nSupplementary Figure S10. Trajectory conservation analysis for the testis\ndataset across 11 species.\nSupplementary Figure S28. Integration results of different methods on the\nbulk dataset.\nSupplementary Figure S29. Integration results of different methods on the\nbulk dataset.\n24. I would rephrase “both CAMEX and scANVI struggled to identify them\naccurately” because the methods can't annotate a cell type that does not exist\nin the training set, and this method doesn’t infer. I would highlight the nature\nthat these methods rely on existing annotations so they simply cannot annotate.\nResponse: Yes, as you pointed out, this is indeed a limitation of the existing\nannotation methods. We have clarified the shortcomings of the existing\nmachine learning-based annotation methods in the manuscript (page 12, lines\n475-477).\n“Due to the absence of neural progenitor cells (NPCs) in the human reference\ndataset, both CAMEX and scANVI cannot identify them.”\n27. The authors mentioned CCL3 is exclusive to humans, but later said that the\nexpression of CCL3 in the subgroups of chimpanzees is of certain rationality\nfrom CAMEX. These seem to contradict. Meanwhile, please also cite\nreferences to support the claim from lines 487-489 about which cell types exist\nin which species.\nResponse: The cell types and subtypes are manually annotated in the original\nliterature [1] through clustering and marker genes. In the original literature, the\nauthor believes that CCL3 is unique to humans. Here, we used the human\ndataset as a reference (training set) and employed semi-supervised learning to\nperform cell annotation with CAMEX. From the results, we speculate that CCL3\nalso exists in Chimpanzees, and we can confirm this through the UMAP plot\nand verification of gene expression levels in Supplementary Figure S36.\nSupplementary Figure S36. Integrating the microglia cell dataset of four\nspecies by CAMEX. a, UMAP visualization of the microglial dataset across four\nspecies. b, Visualizations of CCL3 and GLDN expression in the microglia\nacross four species. c, Violin plot of GLDN, CCL3, and APBB1IP expression in\noriginal manual annotation and CAMEX prediction clusters.\nFollowing your comments, we rewrote this sentence (page 12-13, line 488-\n490): “The cell class, subclass, and subtype are manually annotated by\nclustering and marker genes in the original study [1].”\n[1] Ma S, Skarica M, Li Q, et al. Molecular and cellular evolution of the primate\ndorsolateral prefrontal cortex[J]. Science, 2022, 377(6614): eabo7257.\nReviewer #2 (Remarks on code availability):\nThe authors have provided examples of the annotation capability using CAMEX\naccordingly.\nResponse: We sincerely thank you for this valuable comment.\nResponse to Review #3\nThe revised manuscript demonstrates significant improvements in\nmethodological description, dataset construction, model architecture, and\nexperimental validation. The overall logic is now more rigorous and coherent,\nand the conclusions are more convincing. Below are several minor revision\nsuggestions aimed at further enhancing the work:\nResponse: We appreciate your kind recognition and detailed comments.\nFollowing your suggestions, we have carefully revised our manuscript and\nadded more evaluation details.\n1. Although the authors have supplemented their work with comparative\nexperiments involving methods such as Seurat v4, scGen, SAMap, SATURN,\nand scGPT, it is still recommended to enhance the transparency of\nexperimental design consistency further and ensure the fairness of results.\nSpecifically, please clarify whether each baseline method underwent\nsystematic hyperparameter optimization or was run with default configurations.\nIn addition, it is advisable to include comparisons of computational resource\nconsumption and runtime efficiency to enable a more comprehensive\nassessment of each method’s feasibility in practical application scenarios.\nResponse: Thanks for your suggestion. As you said, for all experiments, we\nused the default hyperparameters for each of the baseline algorithms and\nCAMEX. Specifically, for baseline methods (page 21, lines 817):\n“In addition, we use default hyperparameters in all experiments for all methods\nexcept the bulk RNA-seq dataset.”\nFor CAMEX (page 19, lines 764-771):\n“(4) Hyperparameters:\nCAMEX employs a set of default hyperparameters, including: train_mode:\n'mini_batch', dim_hidden: 128, batch_size: 1024, gnn_layer_num: 2,\nepoch_integration: 10, epoch_annotation: 10, and encoder: 'GCN'. All\nexperiments except for bulk RNA-seq integration used these defaults. For bulk\nRNA-seq data, which has a small sample size, we recommend setting the\ntrain_mode to ‘full_batch’ and setting epoch_integration to 200. Further details\nare provided in the GitHub tutorials.”\nIn addition, the training hyperparameters of the 5 main experiments have been\nmade public on the GitHub repository, such as:\nFigure 2 liver dataset:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/1liver/params.py\nFigure 3 testis dataset:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/2testis/params.p\ny\nFigure 4 bulk dataset:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/3bulk/params.py\nFigure 5 cortex dataset:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/4cortex_annotati\non/params.py\nFigure 6 whole brain dataset:\nhttps://github.com/zhanglabtools/CAMEX/blob/main/analysis/5micro_mapping\n/params.py\nRegarding the comparisons of computational resource consumption and\nruntime efficiency, Supplementary Tables 2, 6, 8, 10, and 25 present the\nresource consumption of all algorithms for various tasks.\n2. In the multi-species integration experiments, the authors demonstrate that\nCAMEX maintains stable performance as evolutionary distance increases. The\nunderlying causes or biological mechanisms driving performance variation with\nincreasing evolutionary distance should be further explored in order to\nstrengthen the interpretability of the results.\nResponse: Thanks for your kind suggestion. Among the 11 species included in\nthe testis dataset, seven (i.e., Bonobo, Chimpanzee, Human, Gorilla, Gibbon,\nMacaque, and Marmoset) are primates, while the remaining three (Mouse,\nOpossum, and Platypus) are mammals. Overall, apart from chickens, the\nevolutionary distances between species are not significant. As shown in\nSupplementary Figure S22, the performance of most integration methods\ndeclines markedly at the phylogenetic position of Mouse, consistent with the\nexpectation that integration becomes more challenging with increasing\nevolutionary distance.\nTo further investigate this, we performed an ablation experiment in which\nthe human data were integrated separately with each of the other ten species\nusing only one-to-one homologous genes. Supplementary Figure S24b\nillustrates that integration performance remained comparable to that achieved\nwith many-to-many homologous genes when closely related primates were\ninvolved. However, as evolutionary distance expanded, particularly in the case\nof human–chicken integration, a significant decrease in performance was\nobserved.\nThese results indicate that the use of many-to-many homologous genes is\na key factor enabling CAMEX to maintain robust integration performance\nacross larger evolutionary spans. A detailed description of this ablation study\nhas been included in the manuscript (page 9, lines 324-333).\n“Additionally, we conducted an ablation experiment in which human data were\nintegrated separately with each of the other ten species using only one-to-one\nhomologous genes. As shown in Supplementary Figure S24b, integration\nperformance remained comparable to that achieved with many-to-many\nhomologous genes when closely related primates were involved. However, as\nevolutionary distance increased—particularly in the human–chicken\ncomparison—a marked decline in performance was observed. These findings\nsuggest that the use of many-to-many homologous genes is a critical factor\nenabling CAMEX to maintain robust integration performance across large\nevolutionary spans.”\nSupplementary Figure S22. Illustration of the overall integration\nperformance achieved by various methods when applied to a two-species\ntestis data, allowing for comparative evaluation. The line graph illustrates\nthe single-cell data integration performance of humans with other species.\nDifferent colors represent different methods. Three distinct lines represent:\nBatch correction score, Biological variation conservation (Bio) score, and\nOverall score. For the line charts on the left, a higher value indicates better\nperformance. Conversely, for the box plots on the right, a lower ranking signifies\nbetter performance.\nSupplementary Figure S24. Summary of the overall integration\nperformance achieved by CAMEX when applied to the testis dataset,\nranging from two-species to eleven-species integrations, allowing for\ncomparative evaluation. a, Three distinct lines represent key performance\nmetrics: the Batch correction score, the Biological variation conservation (Bio)\nscore, and the Overall score. b, On the testis development dataset, the result\nobtained by integrating the single-cell data of humans and ten other species\nseparately using only one-to-one homologous genes.\n3. It is recommended that the newly added experimental results be organized\nand integrated in a more structured manner to improve the overall coherence\nand clarity of the manuscript. At the same time, please consider including\nspecific examples of application scenarios for cross-species analysis in the\nIntroduction and Discussion sections, to further highlight the research value and\npotential applications of the work.\nResponse: We appreciate your detailed comment. First, in the method section,\nwe emphasized that integration and annotation are two distinct tasks. The main\nfeature is that the integration is unsupervised and does not require any cell\nlabels. At the same time, the annotation is semi-supervised and requires at least\none labeled reference dataset (page 20, lines 811-817):\n“It is crucial to note that the benchmarking strategy is task-dependent. For the\nunsupervised integration and alignment task, CAMEX is compared against\nlabel-free methods such as Harmony. For the semi-supervised cell-type\nannotation task, however, the comparison is made against reference-based\nmethods like SingleCellNet (SCN), which require at least one labeled reference\ndataset.”\nSecond, we explicitly clarify the task of each section. In the result section 2\n“CAMEX achieves superior integration performance in different multi-species\nscenarios”, we added the description (page 6, lines 204-205): “In this\nunsupervised integration setting, CAMEX operates without the need for\nmanually annotated cell type labels.\nIn section 3 “CAMEX uncovers the conserved differentiation process in the\ntestis across 11 species”, we added the description (page 8, lines 290-292):\n“In this unsupervised integration setting, both CAMEX and the baseline\nmethods are precluded from using any manually annotated cell type labels.”\nIn section 4 “CAMEX aligns various development stages of seven organs\nacross seven different species”, we added the description (page 10, lines 387-\n388): “In this unsupervised integration setting, CAMEX does not require any\nmanually annotated cell type labels.”\nIn section 5 “CAMEX could achieve more accurate integration and annotation\nperformance in both relatives and distant species”, we added the description\n(page 12, lines 449-450): “First, we applied CAMEX and baseline methods to\nintegrate the scRNA-seq data in the absence of cell type annotations”.\n(page 12, lines 465-467): “In this semi-supervised scenario, the human dataset\nserved as the labeled reference. We then applied CAMEX to annotate the three\nunlabeled target datasets, with a focus on the lizard and turtle species.”\nIn section 6 “CAMEX facilitates the discovery of new populations and markers\nin the primate granular dorsolateral prefrontal cortex (DLPFC) data”, we added\nthe description (page 13, lines 491-492): “Firstly, we integrated all the cells in\nthe DLPFC dataset among the four species by CAMEX in an unsupervised\nmethod.”\n(page 13, lines 503-505): “Here, we use the cells in humans as the reference\n(training set) to annotate cells in other species (test set) in a semi-supervised\nway by using a GAT classifier (Method).”\nAs to the specific examples of application scenarios, CAMEX mainly addresses\nthe issues of batch effect removal, alignment of the Carnegie stages, and the\ndiscovery of new cell types or subtypes. We added some descriptions.\nIntroduction (page 3, lines 73-75): “For instance, batch effects make it\nchallenging to align the developmental stages of different species, hinder cross-\nspecies cell annotation, and the identification of novel cell types and subtypes.”\nDiscussion (page 15, lines 615-627): “Fourthly, Single-cell RNA sequencing\nstudies in plants face greater challenges than those in animals due to gene\nfamily expansions from whole-genome duplications, which complicate cross-\nspecies integration [1]. In addition, recent studies have extended cross-species\nintegration models to encompass the construction of developmental\nevolutionary trees and perturbation predictions [2], representing highly\npromising applications in both biological and pharmaceutical research. It is\nanticipated that CAMEX will incorporate additional components in the future to\nbetter adapt to emerging tasks. Last, as technology evolves, we envision\nCAMEX broadening its scope to cover more cross-species data [3] and to\nencompass multi-omics and spatial transcriptomics. This approach would\nprovide a more comprehensive map for species origin and evolution by\nleveraging multi-dimensional information.”\n[1] Passalacqua M J, Gillis J. Coexpression enhances cross-species integration\nof single-cell RNA sequencing across diverse plant species[J]. Nature Plants,\n2024, 10(7): 1075-1080.\n[2] Zhong H, Han W, Cui G, et al. Unify: Learning Cellular Evolution with\nUniversal Multimodal Embeddings[J]. bioRxiv, 2025: 2025.09. 07.674681.\n[3] Wang P, Liu W, Wang J, et al. scCompass: An integrated cross-species\nscRNA-seq database for AI-ready[J]. bioRxiv, 2024: 2024.11. 12.623138.\n4. The formulas should be reformatted according to standard conventions,\nensuring consistent notation throughout the manuscript.\nResponse: Following your suggestion, we standardized the mathematical\nnotation by the bold italic uppercase letters for matrices, bold italic lowercase\nletters for vectors, italic lowercase letters for scalars, and calligraphic\nuppercase letters for sets. Superscripts were enclosed in parentheses to\ndistinguish layer indices from exponents. For example, we revised the notations\nin the following sentences as (page 16, lines 631-638):\n“A scRNA-seq dataset comprising 𝑁 cells and 𝑀 genes is represented as a\nmatrix 𝑿 = (𝒙 ,𝒙 ,…,𝒙 )𝑇 ∈ 𝑅𝑁×𝑀 . Following a standard preprocessing\n1 2 𝑁\npipeline (normalization, log-transformation, and selection of highly variable\ngenes, clustering, and obtaining differential genes), we construct a cell-gene\nbipartite graph 𝑨. In this graph, an edge connects cell 𝑛 and gene 𝑚 if and\nonly if cell 𝑛 expresses gene 𝑚 . Subsequently, cell-cell connections are\nestablished by applying the kNN algorithm within the PCA space derived from\nthis processed matrix.”"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The author's revision is well done.", "reviewer_3": "The authors have provided clear and comprehensive responses to all the questions raised in the review. The revision appropriately\naddresses the suggested modifications, and the manuscript has been improved accordingly. I find the changes satisfactory, and I have\nno further comments or concerns.\n(Remarks on code availability)\nThe code accompanying the paper is well-prepared and serves as a valuable resource for the research community. It includes a clear\nand comprehensive README file with detailed instructions for installation, dependencies, and execution, which made the setup\nprocess straightforward. The documentation also highlights key parameters and configurations, enabling users to adapt the code to their\nown datasets and experiments."}, "author_response": ""}]

Decision letter