Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks

Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks View ORCID Profile Sishir Subedi , View ORCID Profile Yongjin P. Park doi: https://doi.org/10.1101/2025.06.04.25328900 Sishir Subedi 1 Bioinformatics Graduate Program, University of British Columbia , BC, Canada 2 BC Cancer Research, Part of Provincial Health Care Authority , BC, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sishir Subedi Yongjin P. Park 2 BC Cancer Research, Part of Provincial Health Care Authority , BC, Canada 3 Department of Pathology and Laboratory Medicine, University of British Columbia , BC, Canada 4 Department of Statistics, University of British Columbia , BC, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yongjin P. Park For correspondence: ypark{at}bccrc.ca Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Gene expression variation in cancer cells is attributed to many inherited and environmental factors, including genetic variants and cellular landscapes. Decomposing different sources of information is intractable with single-cell RNA-seq alone. However, we show that our new approach can split them with the help of multiple patients, assuming that cell types are widely shared and genetic effects are specifically present in a particular patient. Our approach based on a cross-attention neural network was applied to three different cancer types to identify cell types and patient-specific genetic effects in transcriptomic data. Residual expressions, excluding cell types, can implicate patient-specific disease mechanisms. Background Cancer cells emerge as a result of long-time accumulation of mutations on several driver genes. As cancer cells progress and adapt to tissue of origin and metastasize to new tumour microenvironments, the molecular mechanisms of primary mutations are shaped by cellular landscapes [ 1 ] and germline mutations [ 2 ], showing a much broader degree of polymorphism than the initial genetic variation. Conversely, different mutations may converge to a similar cancer type due to strong influences from tumour microenvironments and the evolutionary selection process [ 3 , 4 ]. Such interactions between genetic and tumour microenvironmental factors are often omitted in observational cancer studies. Therefore, identifying the underpinning mechanisms of cancer development and progression still stands a substantial challenge in basic and translational research. Single-cell genomic techniques emerge as an instrumental tool in our efforts for unbiased profiling of tumour cells and their microenvironments, providing a necessary basis to disentangle cancer heterogeneity [ 5 ] and gene regulatory networks and gene-gene interaction circuitry beyond direct mutational targets [ 6 ]. Traditional single-cell analysis methods often assume that the data generation process of cancer data is not very different from many other single-cell data measured on reference tissue samples. However, single-cell transcriptomics derived from multiple cancer patients shows different characteristics from the corresponding reference tissues and cell types of non-cancer samples. Here, we recognize that extracting cellular commonality shared across many patient-derived samples in contrast to patient-specific (potentially genetic) patterns is a key question to be addressed in high-dimensional single-cell data analysis. Related work Although this work does not entirely hinge on batch adjustment, a single-cell omics data integration is often considered a batch correction problem. Large-scale single-cell datasets provide unbiased resources to study intrinsic biological mechanisms and patient-specific environmental/genetic factors. However, integrating datasets over tens of patients poses significant computational challenges, making existing batch adjustment methods less ideal. A large body of existing work focuses on reconciling low-dimensional latent representations (cellular embedding) across different batches and data modalities [ 7 ] in order to result in unbiased results in downstream clustering and cell-type-specific differential gene expression (DGE) analysis. Several notable methods build on graph-based data integration framework, where cell-cell interaction graphs come from k-nearest neighbour search, including the first version of a mutual nearest neighbour method [ 8 ], Scanorama [ 9 ], BBKNN (batch-balancing k-nearest neighbour) [ 10 ], Harmony [ 11 ], and the weighted nearest neighbour method tailored for multimodal integration [ 12 ]. The performance of these distance graph-based methods arguably relies more on quality latent representation learning than the algorithmic matching and adjustment procedures, as they take substantially similar steps. Yet, by default, almost all of them learn latent embedding by principal component analysis (PCA), or equivalently singular value decomposition (SVD) on log1p-transformed genes. Non-negative matrix factorization (NMF) identifies more interpretable and robust latent feature space [ 13 , 14 ]. An atlas-scale benchmark study focusing on data integration emphasizes that model-based approaches, such as a variational autoencoder model [ 15 ], yield generally richer and robust latent feature space [ 16 ]. The notable examples include SCVI [ 17 ] and scANVI [ 18 ], where we can directly incorporate batch information and covariates in the generative model to remove potential biases present in single-cell data. However, as the research community has gradually become cautious about the risks of directly manipulating latent space to produce visually appealing two-dimensional scatter plots [ 19 ], demands for objective and systematic evaluation criteria also grew. A recent study warns potentially harmful artifacts could be seen among several model-based methods [ 20 ]. Model parameters can be easily unidentifiable if a deep learning algorithm has to rely on noisy observations while not encouraging the models to pick up invariant features. More recently, several approaches have been proposed for variance decomposition in multi-donor single-cell data integration problems. CellANOVA [ 21 ] employs a linear mixed effect modelling approach to separating out condition-specific and/or batch effects from common (fixed) effects estimated by existing batch adjustment methods, such as Harmony [ 11 ] trained on SVD features. Factorization methods, GEDI [ 22 ] and Lemur [ 23 ], seek to isolate donor-specific or condition-specific signals, which are often marginalized in traditional data integration methods. The idea is to include two a priori independent components in the matrix factorization models: gene-to-latent dictionary parameters and latent-to-cell latent variables. A deep “fair” autodencoder model [ 24 ], such as ScANVI [ 18 ] and Biolord [ 25 ], also modify latent-to-cell variables to disentangle multiple sources of variation in latent space. We find a similarity with our previous work, deltaTopic [ 26 ], modelling distribution shifts. However, environmental or conditional effects considered in the previous work are intrinsically different from patient-specific effects that we normally observe in cancer data. We found that a careful experimental design with paired control-condition samples is much desired to estimate batch and condition effects. ConstrastiveVI [ 27 ] is designed to address similar questions, hoping that shared and condition-specific latent dimensions are statistically independent with the knowledge of control/background cells, similar to RUV (removing unwanted variation) methods based on control genes and samples [ 28 ]. Again, we normally have little knowledge of control cells in heterogeneous cancer data; hence, the method is not always applicable. Our contributions We propose a deep learning approach, PICASA, short for Partitioning Inter-patient Cellular Attributes by Shared Attention networks, to systematically dissect multiple sources of variation embedded in cancer transcriptomic data. Our novelty lies in a tailored neural network architecture and training algorithm that can effectively tackle practical challenges practitioners frequently encounter in multi-donor cancer single-cell data analysis. We hypothesize that intrinsic cellular effects are invariantly manifested in the majority of patients/donors/samples. Built on the success of foundation models in single-cell data [ 29 – 32 ], we show that ubiquitous cell type signatures can be robustly redeemed by training cell pairs as if clinicians understand tumour heterogeneity in light of matched normal samples. Having the common cellular contexts established, we show patient or donor-specific factors can be better understood. For each of the three different cancer examples, including breast [ 33 ], ovarian [ 34 ], and lung cancer [ 35 ], we show that a cell type commonality and the residual patient-specific signature models are effectively decomposed by PICASA while being markedly independent of each other. Unlike many existing distance-based batch adjustment methods (unable to recover cell-type-specific generative models), PICASA learns transferrable gene/feature embedding coordinates and cell-type-specific gene-gene interaction patterns as an attention layer (to be discussed later). We also report that many cancer patient-specific signatures are deemed somatic and genetic, such as copy number variants (CNV). However, we regretfully expose the fundamental limitations of an existing CNV calling strategy based on transcriptomic data, as gene expression vectors normalized along genome axes may contain mixed information, obfuscating between genetic and cellular factors. Methods Study design Our primary interest lies in a systematic approach to disentangling gene expression variation of single-cell RNA-seq data measured in multiple patients/donors according to their potential sources of origin. As a single cell in cancer studies may be affected by genetic and cellular environmental factors, we assume that cellular effects are to be shared across multiple patients, whereas genetic factors can be more isolated within each patient and a specific treatment condition. Each cell’s patient-of-origin and its treatment assignments are largely considered labelled information. However, cell type labels are held out from our algorithm, but later used to confirm that the signatures shared across patients indeed correspond to tumour microenvironments. Data preprocessing All the datasets used in this study are publicly available. We obtained two groups of datasets–one group contains simulated and real data for integration benchmark experiments, and the other group contains cancer patient single-cell data from different cancer types. We conducted a series of standard single-cell data preprocessing and quality-control steps. First, we filtered out mitochondrial and spike genes along with genes detected in less than three cells. Then, we normalized cells by a fixed sequencing depth, i.e., 10 4 , followed by log normalization and high variable gene selection. Throughout all the experiments, we used 2,000 highly variable genes for the initial unsupervised learning for the commonality model and brought all the features back for the subsequent analysis with quality control procedures recommended in Scanpy library [ 36 ]. Lung cancer The lung cancer dataset [ 35 ] was retrieved from Gene Expression Omnibus (GEO) with the accession number GSE148071 . The dataset consists of 42 tissue biopsy samples, ranging from stage III to IV non-small cell lung cancer (NSCLC). We selected 23 patient samples that contain at least 1000 cells and constructed a dataset with 72,514 total cells, including 26,077 genes. Ovarian cancer The high-grade serous ovarian cancer (HGSOC) dataset, collected from eleven patients before and after chemotherapy [ 34 ], was retrieved from Gene Expression Omnibus GSE165897 . All the patients were homogeneously treated with post-neoadjuvant chemotherapy (NACT) as recommended for inoperable patients with poor prognosis. Breast cancer The breast cancer single-cell dataset of multiple subtypes [ 33 ], including ER +, HER2 +, and Triple Negative Breast Cancer (TNBC), is publicly available at GEO with the accession number GSE176078 . Of the 26 patients, we retained cells collected from 21 patients, of which we found at least 1,000 cells. Normal pancreas single-cell data The total dataset consists of 14,890 cells, each of which contains 2,000 genes, across eight different sequencing batches. We downloaded them from Seuret data integration tutorial, https://satijalab.org/seurat/archive/v3.2/integration.html , but each of them is also publicly available from NCBI GEO GSE81076 , GSE85241 , GSE86469 , and GSE84133 , and EBI Array Express: E-MTAB-5061 . AI model for partitioning inter-patient cellular attributes Overview of PICASA We design PICASA in order to decompose high-dimensional single-cell gene expression vectors. PICASA first seeks to identify inherent cell type patterns robustly shared across multiple patients and then characterize patient-specific gene expression patterns by contrasting them with the common signatures. Our premise is that gene expression patterns derived from cell type signatures are widely present across multiple cancer patients and matched normal samples, as they reflect inherent tissue-of-origin characteristics. However, we argue that somatic mutations and environmental factors are manifested in a patient-specific manner. PICASA consists of two types of neural network models with the following purposes: A commonality network shared across patients : Briefly, we designed a specialized neural network architecture that recognizes the commonality across nearest neighbour pairs of cells between matched patients and trained it with a contrastive learning method [ 37 , 38 ] (see below for details). The model learns patient-invariant latent representations ( Fig. 1 ) and gene-gene interaction networks of each cell type or cell state (e.g., Fig. 2B ). A patient-specific network : Provided that this recognition network can faithfully encode latent states of cell states robustly shared across patients, we capture the remaining patient-specific (unique) gene expression signals with the second neural network built on a transitional variational autoencoder model [ 39 ]. Here, we used negative binomial autoencoder model (scvi) [ 17 ] while conditioning on a part of latent representation learned from the commonality network. Download figure Open in new tab Fig. 1. Overview of PICASA PICASA model consists of two networks: A commonality network for identifying shared biological cell state representation across patients and a patient-specific network for capturing patient-specific effects. (A) The approximate neighbor search method estimates pairwise single-cell data from two patients. (B) The cell-pair data is embedded in high-dimension at the gene level to share gene interactions and projected onto common latent space, capturing shared biological cell states using contrastive learning. (C) A normalized gene expression matrix and pre-trained common latent space by the commonality network are used in the patient-specific network consisting of encoder-decoder architecture. The model learns patient-specific representation from the original expression data not present in the shared common representation. Download figure Open in new tab Fig. 2. PICASA disentangle sample-specific and sample-invariant effects in pancreas and breast cancer data PICASA captures shared cell type and sequencing platform-specific batch effects in normal pancreas data. UMAP of (A) mixed and decomposed (B) common and (C) unique representations (upper panel is cell type label and lower panel is sample label). PICASA captures shared cell type and hormone receptor-driven patient-specific effects in breast cancer. UMAP of (D) mixed and decomposed (E) common representations. UMAP of (F) cancer and (G) non-cancer cells from unique representation space learned by the patient-specific network (upper panel is cell type label and lower panel is sample label). Single-cell data The model takes a gene expression count matrix, X , with n rows (cells) and m columns (genes); each element takes non-negative value, namely X ig ∈ ℝ ≥0 , for a cell i and a gene g . Some studies share a count matrix of normalized values, so the values can be a rational number. Since such a count matrix was collected on each patient p , we specifically mark each count matrix with its patient of origin, namely, . Likewise, we denote a row vector of a cell i by x i with the patient of origin p and use for each element of a gene g in a cell i observed in a patient p . Learning commonality with cross-attention mechanisms between patients The commonality network architectures consists of three components ( Fig. 1B ): (1) gene embedding, (2) cross-attention between a pair of cells derived from different patients, and (3) cell encoding. All the parameters of these components are estimated by a contrastive learning algorithm, which puts generative training to a classification problem [ 37 , 38 ]. Gene embedding space The goal of gene embedding is to locate a gene expression value in semantically identifiable locations. Roughly speaking, a traditional gene expression mapping locates values in the same real line, but our embedding scheme systematically allocates each gene and its expression values into fine-grained positions. We may consider gene embedding as an automatic quantization step. For each cell, we first take log-transformation, namely followed by the depth normalization to 10,000, so that . Without loss of too much generality, we partition the transformed and normalized expression values into 3,000 bins, meaning that our expression vocabulary size is kept at that level. An embedding layer allocates a cell vector to enlarged feature space ℝ m × d , where d is 25 in all the models. We used torch.nn.Embedding implemented in pyTorch library. We may introduce an additional embedding scheme based on gene symbols, as scGPT [ 30 ] does, if we want our embedding layer to be universally applicable across different data sets and studies. Since we assume commonalities across cancer cells is captured by top 2,000 highly expressed genes, and the downstream attention layers always see the same order of genes, such a gene symbol-based embedding does not play a significant role. Each gene’s identity is implicitly identifiable by its position. Cross-attention layers to recognize the commonality How can we represent common cellular mechanisms? Gene regulatory mechanisms are arguably best represented by multiple snapshots of gene-gene interaction networks. In a large language model (LLM), attention mechanisms have been shown to be powerful enough to generate realistic and meaningful sentences by simply tracking which pairs of words (tokens) are associated within a sentence and a paragraph [ 40 ], and so can be the meaning of a cell (gene expression tokens). In order to capture invariant gene interaction networks, we train a pair of gene token vectors (cells) sampled across different patients. For instance, our cross-attention layer seeks to associate expression tokens of a cell i derived from a patient p with another, yet similar, cell j of a different patient q , and vice versa. More precisely, we have the embedded genes/tokens, , for the cell i , where m = 2000, the number of genes (fixed order), and d =25, the embedding dimension, and E ( j ) for the counterpart cell j ( Fig. 1A -B). In order to feed the embedding vector for each gene into the scaled-dot product attention layer, we construct query, key, and value matrices with the dense model parameter matrix W : Query: , Key: , Value: , Having two sets of query, key, and value tensors on both side of cells, i and j , we can synthesize a cross-attention vector for a gene g from the cell i to j as if we are mapping a source cell i to a target cell j ’s context ( Fig. 1B ): where an unscaled cross-attention matrix A takes each element, with the normalizer matrix, These cross-attention intermediate representations, C ( i → j ) ∈ ℝ m × d , are then pooled across the embedding dimension, and constitute a cell-level score vector ∈ ℝ m : Using two fully-connected (FC) layers (the first with 100 and the second with 25 units), we can then map the pooled gene vector into low-dimensional latent space: just as in a typical single-cell variational autoencoder model, such as scvi [ 17 ]. For donors p and q matched up, say that we have two sets of cells, 𝒢 p and 𝒢 q , respectively. For a pair of cells, i ∈ 𝒢 p and j ∈ 𝒢 q , found by k -nearest neighbour (kNN) search, the cross-attention layers bring about two latent representation vectors, z i and z j , respectively. Interestingly, the cross-attention layer force genes in the cell i to be associated with the genes in the matched cell j , and vice versa. Contrastive learning of common cellular representations Unlike generative model training, contrastive learning involves little modelling assumption. Since the latent representation in the bottleneck layer, z ∈ ℝ 25 , is somewhat restrictive, we linearly project them a onto large space: where we use two FC layers–the first and second with 50 units. For each layer, we linearly project the input data and passed through a series of 1-dimensional batch normalization [ 41 ], rectified linear unit (ReLu) activation [ 42 ], and dropout regularization [ 43 ] with the rate 0.1. The similarity of two projected vectors, for instance, h i and h j , can be measured by the following cosine similarity: In order to train the model parameters via self-supervised learning framework [ 38 , 44 ], or noise contrastive estimation (NCE) [ 37 ], we need to define positive and negative examples. Without loss of generality, suppose we pair cells in the sample p , hence a cell with a cell , sampled from the matched sample q . We give positive labels, i and j , to the originally paired cells and find noisy (background) pairs by simply sampling other cell k within the same cell group, , except for itself. Then, we have similarity scores, for the positive pairs and the other scores, for the negative ones. By sampling negative/background cells within the cell group derived from the same donor/patient, we achieve sharp contrasts between the positive and negative pairs. For each mini-batch, mutual information loss [ 44 ] desired by NCE can be bounded by the following: where is a set of negative examples for a cell i found within the same patient group. Learning patient-specific unique representations Our next question is: Had we transported a cell into cell-type-less contexts, what would be the patient-specific signatures? We will address this question by designing a conditional model. Built on a typical variational autoencoder architecture, more precisely, a sufficiently deep encoder and zero-inflated negative binomial (ZINB) decoder proposed by scvi [ 17 ], we include all the sources of gene expression variation. The technical details of ZINB models can be found in the original [ 17 ] and related papers [ 18 ]. Model definition Letting y i be a latent representation vector, we can design an autoencoder model with the following information flow: Encoder: x i → y i Decoder: where the reconstructed signals follow ZINB. Following the idea of the fair autoencoder model [ 24 ], we can introduce the latent features z i , each of which can be evaluated (without gradient) by a forward pass of a gene expression vector x i from the bottleneck layer of the previous network ( Fig. 1C ). Encoder: x i → y i Decoder: Now, our goal boils down to training a ZINB autoencoder model, while conditioning on the “observed” commonality features, so that the new latent representation vector, y , can pick up new signatures independent of common cell type signals. Once we have this conditional decoder model, we can predict patient-specific gene expression patterns in a counterfactual manner. For the estimation of cell-type-free, thus patient-specific, perhaps environmental, factors, we can simply switch off z , such as Decoder( y i , do( z i = 0 )). Multi-task learning to delineate the sources of information In addition to the generative loss of the ZINB model, say L generative , we enforce strong conditional independence between the two types of latent variables Y and Z by adding two additional types of loss functions. All pairwise cosine similarity within a minibatch: Donor/patient label prediction (cross entropy): where if and only if a cell i belongs to a patient p , i.e., , otherwise . We may consider putting different weights on different types of losses. Here, we simply equally weighted, constituting the following total loss function: For each loss, we take into account of different number of data points by taking the average using torch.mean function. Results PICASA disentangle sample-specific and sample-invariant effects Unsupervised clustering on sample-ubiquitous factors aligns with cell type annotations As a proof of concept, we first tested PICASA on a well-known batch correction problem, where single-cell expressions of similar pancreatic tissue samples were surveyed by multiple sequencing platforms and technologies (see Methods). Here, the first goal is to learn robust latent representations pertinent to known cell types and states in normal pancreatic samples in batch and platform independent manners ( Fig 2A-C ). We reused cell type annotations provided by Seurat v3 package [ 45 ]. Since cell type factors are considered one of the largest variance components, a powerful autoencoder model [ 25 ] is generally expected to yield good data integration results. However, we discovered that cells are influenced by both known cell types and batch membership, resulting in many small clusters consisting of cells of the same cell type stemming from one batch. Although two-dimensional visualization, such as UMAP (Uniform Manifold Approximation and Projection [ 46 ]), often exaggerates potential algorithmic artifacts [ 19 ], many disconnected components scattered across multiple batches and cell types ( Fig. 2A ) suggest cells of the same cell types are parent across multiple batches and not connected by 15-nearest neighbourhood graph (UMAP’s default). We found similar results with ScVi [ 17 ]. Even though we designed PICASA with more general utilities in mind, PICASA clearly separates out two types of latent representations–the sample-specific ( Fig. 2B ) and commonality (thus cell types) latent space ( Fig. 2C ). We not only see that cells in the commonality space located in proximity with other cells from the same cell types ( Fig. 2C ), while staying away from the different cell types, but we also find that the sample-specific latent features encourage the cells of the same types and batches are grouped together ( Fig. 2B ). As another proof of concept, we trained deep autoencoder models on the breast cancer dataset [ 33 ] (see Methods) derived from 26 tumour samples. We constructed multi-donor data consisting of 87,012 cells after removing patient samples with small cell counts, retaining samples with at least 1,000 cells. Here, we reused cell-type annotation labels provided by the original authors. The authors used Seurat v3 to perform Leiden clustering with default parameters followed by the supervised classification method Garnett [ 47 ] to annotate the resolved clusters based on canonical cell type markers. Of the total 21,327 genes/features, we selected the top 2,000 highly variable genes to train the first model of PICASA and brought back the full feature set for the second stage. An autoencoder model trained without donor information generally results in cell clusters of homogeneous cell types and of a single patient ( Fig. 2D ). Both cell-type variation and patient heterogeneity seem equally dominant, although cell types are generally considered a stronger binding force in a typical batch correction algorithm [ 8 , 10 ]. PICASA, however, resulted in strong cell-type-specific clustering patterns in the UMAP space ( Fig. 2E , top), mixing cells of multiple patients together within each cell-type cluster/neighbour ( Fig. 2E , bottom). On the other hand, the unique and patient-specific model of PICASA induces tight clustering of cancer cells of the same cancer patient ( Fig. 2F , top) and the same subtype ( Fig. 2F , bottom). However, interestingly, non-cancer cells (tumour microenvironment) are less patient-specific, forming larger clusters including cells across multiple patients. PICASA captures therapy-induced changes in high-grade serous ovarian cancer (HGSOC) Robust commonality model Next, we applied PICASA to cancer data with a more sophisticated experimental design. We obtained 51,786 cells from the publicly available ovarian cancer study [ 34 ]. This study traced changes over 32,847 genes for eleven patients between pre- and post-neoadjuvant chemotherapy (NACT). We used major cell type annotations provided by the original authors, which had been resolved by Seurat v3 [ 45 ] shared nearest neighbourhood clustering pipeline. While not explicitly handling and adjusting patient heterogeneity, autoencoder model fitting resulted in fragmented clustering patterns. Cells generally form well-separated groups by the major cell type categories ( Fig. 3A , left), but we also noticed granular structures primarily determined by the patient-of-origin labels ( Fig. 3A , right), confirming a large degree of patient heterogeneity, even for non-cancer cells. Download figure Open in new tab Fig. 3. PICASA captures therapy-induced changes in high-grade serous ovarian cancer data PICASA captures shared cell type and treatment based patient-specific effects in ovarian cancer. UMAP of (A) mixed and decomposed (B) common and (C) unique representations. (C) Heatmap of gene interaction scores of known marker genes for different cell types in the ovarian cancer dataset learned by the cross-attention module. UMAP of (D) cancer and (E) non-cancer cells from unique representation space learned by the patient-specific network (left panel is patient label and right panel is treatment label). (F) Patient-wise contour plots of CNV profiles generated using copyKAT from patient-specific and mixed representations of ovarian cancer data. (G) Cell type-wise contour plots of CNV profiles for selected patients generated using copyKAT from patient-specific and mixed representations. (H) Heatmap of correlation between CNV profiles from patient-specific and mixed representations for different cell types at treatment-naive and post-NACT stages. In the space shared across multiple patient samples, however, we found cells annotated as the same cell type form bigger clusters, one cluster for each cell type, and we also note that relatively more minor cell types generally form independent groups ( Fig. 3B , left). No batch-or sample-specific clusters are found on the same latent space ( Fig. 3B , right, coloured by the patient labels), contrary to the results of the previous autoencoder model ( Fig. 3A , right). In order to better understand how cross-attention neural architecture can achieve robust cell-type annotations, we examined the patterns of the attention layer within each cell type ( Fig. 3C ), ascertaining the most frequent gene-gene interactions captured during the model training ( Fig. 1 ). For epithelial ovarian cancer (EOC) cells, two marker genes, EPCAM and WFDC2 tend to strongly co-expressed/interacting with each other. For HGSOC, EpCAM is a well-known marker gene, particularly highly expressed in chemotherapy-resistant tumour samples [ 48 ]. WFDC2 (or HE4 ) is also another well-established marker gene of ovarian cancer [ 49 ], exhibiting albeit subtype specificity [ 50 ]. For cancer-associated fibroblasts (CAF), two marker genes, PDPN and DCN , tend to be co-expressed in our model, strongly aligned with CAF studies in pancreatic ductal adenocarcinoma [ 51 ]. For the immune cell types, we generally observe model the model parameters corresponding to marker genes are highly activated– CD14 for monocytes, CD8A and PTPRC /CD45 for T-cells, and MS4A1 /CD20 along with CD79A (B-cell receptor complex) for B-cells. The cross-attention mechanisms link these marker genes with other marker genes. For monocytes, we found strong connection between CD14 and THBD ( Fig. 3C ). As THBD are specifically expressed in hypoxic tumour regions [ 52 ], the high activity of CD14 together with THBD suggests a dominant role M2 macrophages for this HGSOC cohort. Moreover, the co-occurrence of FCER1G corroborates that the monocytes/macrophages are highly tumour-infiltrating [ 53 ], showing clear pathological characteristics of tumour microenvironments. Patient-specific genetic changes across multiple cell types While conditioning on these cell-type latent states, PICASA estimated latent representations of the residual expressions. For clarity, we projected cancer ( Fig. 3D ) and non-cancer cells ( Fig. 3E ) separately. The visual inspection of UMAP embedding suggests cancer cells generally form clusters in a highly patient-specific manner, whereas we find non-cancer cells tend to mingle across different patients. Interestingly, for the cancer cells ( Fig. 3D , right), we can spot markedly distinctive cell groups within EOC733, EOC87, and EOC372, according to the treatment regimen. We then further investigated how these patient-specific signals had been generated in light of somatic mutations. Not having access to the genome sequencing data, we estimated genomic copy number variant (CNV) profiles for all the cells using CopyKAT [ 54 ], which consistently outperformed in the recent benchmark study [ 55 ]. We first estimated CNV profiles on the mixed (raw) single-cell RNA-seq data and then on the patient-specific transcriptomic signals after eliminating the effect of cell type effects (depicted in Fig. 1C ; see Methods for details). Briefly, we trained the full model parameters, including the commonality signatures and later set the commonality to zero values to estimate the residual expression values in forward generation. As recommended, we used CopyKAT patient by patient within the same cell type to avoid potential batch effects. Using the two types of cell-level CNV profiles, we first compared them within each patient ( Fig. 3E ), asking whether the residual signals indeed corresponded to patient-specific genetic effects estimated by CopyKAT or whether there were any possibilities that cellular effects (which we removed) being read as “somatic” copy number changes. For six out of the eleven patients, we found a strong genome-wide correlation between the CNV expressions based on the raw and residual data (R > 0.4); for the other four cases, we found weakly correlated or anti-correlated patterns. In order to better understand what was left after suppressing the cell type effect (the commonality latent states), we aggregated the CNV profiles within each cell type and treatment regimen to compare the correlation between the raw and residual data. A majority of copy number changes could be attributed to occur in EOC cells (R = median: 0.89; mean: 0.48 ±0.84), rather than other cell types in the tumour microenvironments, such as the macrophages (R = median: 0.63; mean: 0.18 ±0.83), the T-cells (R = median: 0.52; mean: 0.26 ±0.6), and the B-cells (R = median: 0.65; mean: 0.57 ±0.41). Especially for the cells in the samples showing weaker CNV correlations between the raw and residual (EOC733, EOC136, EOC540, EOC443, EOC1005), we found that the estimated CNV profiles more widely fluctuate depending on cellular environments and treatment effects ( Fig. 3G ). In EOC773, we measured substantially strong Spearman’s correlation between the raw and residual CNV values before treatment (R=0.99) and the majority of cell populations were considered the EOC type. After the neoadjuvant chemotherapy (NACT), we still found a strong correlation for the EOC (R=0.84). However, we also found that non-cancer cells report high correlation values within each own cell type, such as the cases of T-cells (R=0.83) and the natural killer (NK) cells (R=0.83). Interestingly, we found an inverse correlation for the macrophages (R=-0.86 and very weak correlation for the B-cells (R=-0.06). For some of the cancer-associated fibroblasts (CAF) samples, we also found strong correlations, such as EOC349 (R=0.97) before treatment, EOC136 (R=0.86), EOC153 (R=0.84), EOC3 (R=0.95), EOC443 (R=0.94), and EOC87 (R=0.95) after treatment. Arguably, CNV results relying on gene expression profiles are not the best way to track cancer clonal status. Here, our results also suggest that a substantial fraction of the CNV variance can be affected by intrinsic and ubiquitous cell type effects rather than individualized somatic changes. Nonetheless, it is reassuring to see generally stronger correlations among the malignant cell types. The previous copy number signature analysis [ 56 ] demonstrates the importance of CNV profiles in tumour samples, implicating markedly different prognoses and relapse depending on varying levels of genetic exposures. While looking at copy number changes in non-cancer single-cell data could be controversial, our results indicate that the CNV profiles in the CAFs for several patients seem more persistent than the original cancer cells. It is possible that these CAF cells were somehow misclassified in the original study. Alternatively, it is possible that cell-cell communications between the tumour and CAF cells could have resulted in genomic instability, potentially contributing to relapse and poor prognosis after treatment. We also note that differentially gene expression patterns often highly overlap with copy number changes in many cancer studies [ 57 , 58 ]. PICASA characterizes subtypes of non-small cell lung cancer (NSCLC) cells Robust cell type clusters shared across patients and different subtypes As we have seen in the previous cancer studies, PICASA consistently map 72,514 cells derived from 23 NSCLC patients, encompassing both lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) [ 35 ], and strongly induce distinctive cell type clusters while not showing any patient-specific effects ( Fig. 4A ). Each cell type’s commonality is characterized by interactions between marker genes ( Fig. 4B ). The immune cell types enrich well-known cell surface and other marker genes/proteins: a surface marker gene CD14 and lysosome ( LYZ ) gene for monocytes, immunoglobulin genes ( IGHG3 and JCHAIN ) for plasma cells, and CD2 and IGHG3 for T-cells. Moreover, for endothelial cells, a key autocrine regulator gene, ANGPT2 (angiopoietin-2), is highly expressed [ 59 ] and interacts with a wide spectrum of genes, including extracellular matrix protein genes, COL1A2 and DCN2 , which are more pertinent to epithelial and fibroblasts. Download figure Open in new tab Fig. 4. PICASA characterizes subtypes of non-small cell lung cancer PICASA captures shared cell type and histologically defined cancer subtype-specific effects in lung cancer. (A) UMAP of decomposed common representation with cell type label. (B) Heatmap of gene interaction scores of known marker genes for different cell types in the lung cancer dataset learned by the cross-attention module. (C) Gene set enrichment analysis of lung cancer data based on ranked genes in the cross-attention matrix with the known cell-type marker genes in the PanglaoDB database. (D) Heatmap of correlation between CNV profiles from patient-specific and mixed representations for different cell types. NAPSA (napsin A) is highly coexpressed with KRT8A and LYZ in the epithelial cells’ patterns, as implicated in independent experimentally validated study [ 60 ]. We notice a switching pattern of NAPSA and KRT6A expressions between the epithelial and malignant cells. Considering that KRT6A promotes epithelial-to-mesenchymal transition (EMT) in many cancer types [ 61 ], including NSCLC, we can speculate systematic rearrangement of gene-gene interaction networks during cancer progression. On the other hand, DCN (decorin), an extracellular matrix protein/gene, suppresses the EMT process and functions as a tumour suppressor [ 62 , 63 ]. We found DCN is specifically interacting with the other extracellular matrix protein/gene, COL1A2 , and coexpressed with KRT6A , perhaps inhibiting KRT6A ’s function in EMT, maintaining the fibroblasts non-cancerous. Cell type-specific cross-attention activities also recapitulate known cell-type marker genes ( Fig. 4C ). We conducted a gene set enrichment analysis of highly variable genes (2,000 genes) within each common cell type cluster against the marker genes in PanglaoDB [ 64 ] using rank-based Fast Gene Set Enrichment Analysis (GSEA) [ 65 ] implemented in R’s fgsea package. High-activity genes in the endothelial cluster are significantly over-represented in endothelial cells (p=0.001), stromal cells (p=0.008), and precursor mesothelial cells (p=0.008). We found that the top genes in the fibroblast cluster enrich mesangial cells (p=0.005) and myofibroblast cells (p=0.006). For the malignant cells, our GSEA results highlight the true cell type of origin, which are epithelial cell types, such as alveolar goblet cells (p<0.001), epithelial cells (p=0.001), and basal cells (p=0.005). We observe similar enrichment patterns for the non-cancer epithelial cells, whose gene activity ranks are associated with goblet cells (p=0.015), epithelial cells (p=0.02), alveolar cells (p=0.03), and basal cells (p=0.04). As for the plasma cells and other immune cells, we found the plasma cell clusters match with plasma (p=0.001)/B cells (p=0.006), the monocyte/macrophage cell cluster with alveolar macrophages (p=0.009), and T-cells with T-cells (p=0.002), cytotoxic T-cells (p=0.002), helper T-cells (p=0.006) and naive T-cells (p=0.04). The residual (patient-specific) expression patterns recapitulate genetic effects As we did in the ovarian cancer study, we made CNV calls for all the cells using CopyKAT [ 54 ]. We summarized Spearman’s correlation between the two types of CNV profiles estimated on the raw (mixed) and residual (patient-specific) expression, respectively. For the Alveolar cells derived from the four LUAD (adenocarcinoma) patients, we found strong correlations (median: 0.9; mean: 0.81 ±0.26), suggesting that the genetic effects of the cells in these patients had not been “adjusted” by cell type effects (thus, independent). Likewise, we found strong positive correlations for cancer cells derived from LUSC (squamous carcinoma) patients (eight out of thirteen), yielding R=median: 0.89; mean: 0.47 ±0.74. In summary, much of the genetic effects were seen in alveolar cell clusters for LUAD patients, whereas genetic effects of cancer cells were quite robustly sustained in LUSC patients. Considering substantially higher mutation frequencies for TP53 (genomic integrity) in LUSC patients (63%) than LUAD (34%) [ 66 ], we can expect more prevalent structural variation in chromosomal architectures, spanning over a long stretch of many gene bodies. However, interestingly, we found negative correlation values for P39 and P8 (adenocarcinoma) for the monocytes/macrophages, respectively, R=-0.93 and -0.56. We also found negative correlations for the malignant cells in P16 and P28, R=-0.79 and -0.69, respectively. Such negative correlations between the raw and residual data leave room for further discussion. Single-cell data alone measured within each patient and cell type seem to lack a necessary basis to discern transcriptomic signals derived from genetic/somatic effects and developmental/differentiation processes. Especially for cancer cells, genetic structural variants can easily manifest in gene expression activities [ 57 ]. If cell-type-specific signatures stretch over long consecutive blocks of the genome, normal transcriptomic patterns may as well be misinterpreted as genetic/somatic. Benchmark analysis Simulation studies for data integration tasks Here, we focused on the performance of the first step of PICASA ( Fig. 1B ). To the best of our knowledge, no simulation method can specifically simulate heterogeneous patient-specific genetic effects along with multiple cell type contexts. We used benchmark data designed for atlas-level data integration tasks [ 16 ]. For the total of 19,318 cells, we first defined equally-sized (technical) batches and independently grouped them into four cell-type groups. The batch membership was provided to unsupervised learning algorithms, but the cell-type labels were held out for benchmark evaluation. For simplicity, we focused on 2,000 genes that are highly expressed and thus consistently available for all the cells. We compared the quality of latent variables induced by multiple stages of our PICASA framework with other methods both qualitatively ( Fig. 5A-B ) and quantitatively ( Fig. 5C-D ), including PCA: PCA-based data integration implemented in Scanpy [ 36 ] BBKNN: Batch-balancing K-nearest neighbour [ 10 ] Combat [ 67 ]: factorization-based batch adjustment method implemented in Scanpy Harmony [ 11 ]: cluster-based batch adjustment Scanorama [ 9 ]: mutual nearest neighbour method Liger [ 13 , 14 ]: a non-negative matrix factorization method SCVI [ 17 ]: variational autoencoder model with batch labels PICASA-C (this work): latent commonality space PICASA-U (this work): residual, sample-unique space PICASA-UC (this work): the commonality and residual states combined Download figure Open in new tab Fig. 5. PICASA accurately delineates between common cell type and patient-specific unique patterns UMAP representation of Simulation-2 (large-scale, balanced) dataset capturing shared cell type effects and integrating cells across samples. (A) Cell clusters are labelled by cell type. (B) Cell clusters are labelled by sample label. (C) Comparative analysis of PICASA capturing shared cell type effects across patients in different datasets. The plot shows cell type Local Inverse Simpson’s Index (LISI) (lower value suggests better mixing of cells from the same cell type across samples). (D) Comparative analysis of PICASA integrating cells across patients in different datasets. The plot shows batch label Local Inverse Simpson’s Index (LISI) (higher value suggests better mixing of cells across samples). As a representative example, we projected the latent states induced by all these methods in two-dimensional UMAP space and coloured each cell according to cell types ( Fig. 5A ) and samples ( Fig. 5B ). BBKNN, Liger, and PICASA-C form a single (nearly perfect) cluster for each cell type group. The other methods, such as SCVI and COMBAT, generally put one cell type in one cluster but often result in disjoint and scattered tiny clusters unable to reconcile the differences between samples ( Fig. 5B ). Clusters generated by PICASA-C are tightly close to each other within each cell type while mixing multiple samples. On the other hand, PICASA-U works in the opposite way, mixing multiple cell types while separating out samples forming sample-specific clusters. PICASA induces cell-type-specific clusters mixing multiple samples/batches We then quantified how well PICASA and other methods perform in data integration tasks by calculating LISI (Local Inverse Simpson’s Index) scores [ 11 ], namely, where p is the probability of cell i ’s neighbours belong to the same cell type or batch. A low LISI value means a homogeneous mixture of neighbouring cells with the minimum value=1; a high LISI implies the opposite, heterogeneity within the neighbours. Our desiderata for PICASA-C is to have low LISI values with respect to cell type labels and high LISI regarding batch/sample labels. For PICASA-U, we expect the opposite–high cell type LISI and low batch LISI. All the methods we considered here show strong performance in terms of grouping cells into clusters of a homogeneous cell type, notably except for PICASA-U ( Fig. 5C ). PICASA-U was expected to find clusters independent of cell type labels, so it did for the simulation, pancreas, and ovarian cancer data. However, in the breast and lung cancer studies, the LISI values of PICASA-U were not significantly different from the other methods. For the batch LISI, we find that PICASA-U and PICASA-UC clearly pick up single-patient (or batch) neighbours as intended by design ( Fig. 1C ). PICASA-C, along with LIGER, identified cell type-specific clusters on the basis of cells from diverse batches. Taken together, benchmark studies suggest that two types of latent variables that PICASA generates are empirically identifiable and statistically independent from each other. Discussion We note several limitations of this study. One of our key ideas in neural network architecture builds on a simplified version of a transformer network [ 40 ]. For simplicity, we only implemented expression embedding schemes restricting the 2,000 most highly variable genes/features in modelling gene-gene interaction attention matrices. In many cases, unless our goal is to build a foundation model, a thousand genes often suffice to capture general cell type signatures. However, we can increase the number of features in our cross-attention layers by adopting gene symbol-based embedding and tokenization as in other single-cell foundation models [ 30 , 32 ], where we restrict the number of tokens with different genes for different cells. Moreover, the attention mechanism layer itself can be more efficiently implemented, reducing memory footprint and computation cost [ 68 ]. More efficient transformer architectures will enable us to include ten times more features in the model [ 69 ]. We intentionally designed our approach to take two steps in getting final patient-specific expression values, later used to make genome-wide CNV calls. In our preliminary analysis, resolving both shared and patient-specific latent representations poses challenges in model identifiability. In fact, without enforcing strong prior probability, finding latent variables mutually independent of one another is a non-trivial problem [ 70 – 72 ], perhaps an open question. In a sense, estimating the parameters of a complete joint model would necessarily involve somewhat arbitrary and laborious hyperparameter optimization. Instead, joint analysis with genomic data (DLP+) [ 73 ] and probabilistic mapping between DNA and RNA samples [ 58 ] is expected to pave a new avenue to more principled approaches. Alternatively, we can preprocess to have single nucleotide variant calls and read depth profiles from RNA-seq data to make sure that different data modalities capture genetic effects. Conclusions We present a new AI approach to delineating patient-specific and patient-invariant signals from single-cell RNA-seq data derived from tens of patients. We demonstrated our approach, PICASA, which can first identify the commonality latent variables shared across patients/donors using specialized cross-attention neural network layers, and the residual signals, not explained by the commonality, generally correspond to CNV calling results. Here, we also challenged the notion of patient specificity, albeit revealing that copy number calling purely based on the raw count data can be based on cell type signatures shared across multiple patients. Nonetheless, much of cancer cells’ patient-specificity can be explained by genetic factors independent of cell type variation. Independent genetic screening will provide complementary data sources to corroborate our observations. Data Availability All the datasets used in this study are publicly available in Gene Expression Omnibus (GEO) database. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE148071 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE165897 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE176078 Footnotes Lead contact: Yongjin P. Park, ypp{at}stat.ubc.ca , yongjin.park{at}ubc.ca Reference 1. ↵ Patel AS , Yanai I. A developmental constraint model of cancer cell states and tumor heterogeneity . Cell . 2024 ; 187 : 2907 – 18 . OpenUrl CrossRef PubMed 2. ↵ Vali-Pour M , Park S , Espinosa-Carrasco J , Ortiz-Martínez D , Lehner B , Supek F. The impact of rare germline variants on human somatic mutation processes . Nat Commun . 2022 ; 13 : 3724 . OpenUrl PubMed 3. ↵ Binnewies M , Roberts EW , Kersten K , Chan V , Fearon DF , Merad M , et al. Understanding the tumor immune microenvironment (TIME) for effective therapy . Nat Med . 2018 ; 24 : 541 – 50 . OpenUrl CrossRef PubMed 4. ↵ Gavish A , Tyler M , Greenwald AC , Hoefflin R , Simkin D , Tschernichovsky R , et al. Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours . Nature . 2023 ; 618 : 598 – 606 . OpenUrl CrossRef PubMed 5. ↵ Sellars MC , Wu CJ , Fritsch EF . Cancer vaccines: Building a bridge over troubled waters . Cell . 2022 ; 6. ↵ Hahn WC , Bader JS , Braun TP , Califano A , Clemons PA , Druker BJ , et al. An expanded universe of cancer targets . Cell . 2021 ; 184 : 1142 – 55 . OpenUrl CrossRef PubMed 7. ↵ Argelaguet R , Cuomo ASE , Stegle O , Marioni JC . Computational principles and challenges in single-cell data integration . Nat Biotechnol . 2021 ; 8. ↵ Haghverdi L , Lun ATL , Morgan MD , Marioni JC . Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors . Nat Biotechnol . 2018 ; 36 : 421 – 7 . OpenUrl CrossRef PubMed 9. ↵ Hie B , Bryson B , Berger B. Efficient integration of heterogeneous single-cell transcriptomes using scanorama . Nat Biotechnol . 2019 ; 37 : 685 – 91 . OpenUrl CrossRef PubMed 10. ↵ Polański K , Park J-E , Young MD , Miao Z , Meyer KB , Teichmann SA . BBKNN: Fast batch alignment of single cell transcriptomes . Bioinformatics . 2019 ; 11. ↵ Korsunsky I , Millard N , Fan J , Slowikowski K , Zhang F , Wei K , et al. Fast, sensitive and accurate integration of single-cell data with harmony . Nat Methods . 2019 ; 12. ↵ Hao Y , Hao S , Andersen-Nissen E , Mauck WM 3rd . , Zheng S , Butler A , et al. Integrated analysis of multimodal single-cell data . Cell . 2021 ; 184 : 3573 – 3587 .e29. OpenUrl CrossRef PubMed 13. ↵ Welch JD , Kozareva V , Ferreira A , Vanderburg C , Martin C , Macosko EZ . Single-Cell multi-omic integration compares and contrasts features of brain cell identity . Cell . 2019 ; 177 : 1873 – 1887 .e17. OpenUrl CrossRef PubMed 14. ↵ Kriebel AR , Welch JD . UINMF performs mosaic integration of single-cell multi-omic datasets using nonnegative matrix factorization . Nat Commun . 2022 ; 13 : 780 . OpenUrl PubMed 15. ↵ Kingma DP , Welling M. Auto-Encoding variational bayes . arXiv [stat ML] . 2013 ; 16. ↵ Luecken MD , Büttner M , Chaichoompu K , Danese A , Interlandi M , Mueller MF , et al. Benchmarking atlas-level data integration in single-cell genomics . Nat Methods . 2022 ; 19 : 41 – 50 . OpenUrl CrossRef PubMed 17. ↵ Lopez R , Regier J , Cole MB , Jordan MI , Yosef N. Deep generative modeling for single-cell transcriptomics . Nat Methods . 2018 ; 15 : 1053 – 8 . OpenUrl CrossRef PubMed 18. ↵ Xu C , Lopez R , Mehlman E , Regier J , Jordan MI , Yosef N. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models . Mol Syst Biol . 2021 ; 17 : e9620 . OpenUrl CrossRef PubMed 19. ↵ Kobak D , Linderman GC . Initialization is critical for preserving global data structure in both t -SNE and UMAP . Nat Biotechnol . 2021 ; 1 – 2 . 20. ↵ Antonsson SE , Melsted P. Batch correction methods used in single cell RNA-sequencing analyses are often poorly calibrated . bioRxiv . 2024 ;2024.03.19.585562. 21. ↵ Zhang Z , Mathew D , Lim TL , Mason K , Martinez CM , Huang S , et al. Recovery of biological signals lost in single-cell batch integration with CellANOVA . Nat Biotechnol . 2024 ; 1 – 17 . 22. ↵ Madrigal A , Lu T , Soto LM , Najafabadi HS . A unified model for interpretable latent embedding of multi-sample, multi-condition single-cell data . Nat Commun . 2024 ; 15 : 6573 . OpenUrl PubMed 23. ↵ Ahlmann-Eltze C , Huber W. Analysis of multi-condition single-cell data with latent embedding multivariate regression . Nat Genet . 2025 ; 1 – 9 . 24. ↵ Louizos C , Swersky K , Li Y , Welling M , Zemel R. The variational fair autoencoder . arXiv [statML] . 2015 ; 25. ↵ Piran Z , Cohen N , Hoshen Y , Nitzan M. Disentanglement of single-cell data with biolord . Nat Biotechnol . 2024 ; 42 : 1678 – 83 . OpenUrl PubMed 26. ↵ Zhang Y , Khalilitousi M(sam) , Park YP . Unraveling dynamically encoded latent transcriptomic patterns in pancreatic cancer cells by topic modeling . Cell Genomics . 2023 ; 0 . 27. ↵ Weinberger E , Lin C , Lee S-I. Isolating salient variations of interest in single-cell data with contrastiveVI . Nat Methods . 2023 ; 20 : 1336 – 45 . OpenUrl CrossRef PubMed 28. ↵ Risso D , Ngai J , Speed TP , Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples . Nat Biotechnol . 2014 ; 32 : 896 – 902 . OpenUrl CrossRef PubMed 29. ↵ Szałata A , Hrovatin K , Becker S , Tejada-Lapuerta A , Cui H , Wang B , et al. Transformers in single-cell omics: A review and new perspectives . Nat Methods . 2024 ; 21 : 1430 – 43 . OpenUrl CrossRef PubMed 30. ↵ Cui H , Wang C , Maan H , Pang K , Luo F , Duan N , et al. scGPT: Toward building a foundation model for single-cell multi-omics using generative AI . Nat Methods . 2024 ; 1 – 11 . 31. Theodoris CV , Xiao L , Chopra A , Chaffin MD , Al Sayed ZR , Hill MC , et al. Transfer learning enables predictions in network biology . Nature . 2023 ; 618 : 616 – 24 . OpenUrl CrossRef PubMed 32. ↵ Yang F , Wang W , Wang F , Fang Y , Tang D , Huang J , et al. scBERT as a large-scale pretrained deep language model for cell type annotation of single-cell RNA-seq data . Nature Machine Intelligence . 2022 ; 1 – 15 . 33. ↵ Wu SZ , Al-Eryani G , Roden DL , Junankar S , Harvey K , Andersson A , et al. A single-cell and spatially resolved atlas of human breast cancers . Nat Genet . 2021 ; 53 : 1334 – 47 . OpenUrl CrossRef PubMed 34. ↵ Zhang K , Erkan EP , Jamalzadeh S , Dai J , Andersson N , Kaipio K , et al. Longitudinal single-cell RNA-seq analysis reveals stress-promoted chemoresistance in metastatic ovarian cancer . Sci Adv . 2022 ; 8 : eabm1831 . OpenUrl CrossRef PubMed 35. ↵ Wu F , Fan J , He Y , Xiong A , Yu J , Li Y , et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer . Nat Commun . 2021 ; 12 : 2540 . OpenUrl CrossRef PubMed 36. ↵ Wolf FA , Angerer P , Theis FJ . SCANPY: Large-scale single-cell gene expression data analysis . Genome Biol . 2018 ; 19 : 15 . OpenUrl CrossRef PubMed 37. ↵ Gutmann MU , Hyvärinen A. Noise-contrastive estimation of unnormalized statistical models, with applications to natural image statistics . J Mach Learn Res . 2012 ; 13 : 307 – 61 . OpenUrl 38. ↵ Chen T , Kornblith S , Norouzi M , Hinton G. A simple framework for contrastive learning of visual representations . International conference on machine learning . PMLR ; 2020 . p. 1597 – 607 . 39. ↵ Kingma DP , Welling M. Auto-Encoding variational bayes . arXiv [statML] . 2013 ; 40. ↵ Guyon I , Luxburg UV , Bengio S , Wallach H , Fergus R , Vishwanathan S , et al. , editors Vaswani A , Shazeer N , Parmar N , Uszkoreit J , Jones L , Gomez AN , et al. Attention is all you need . In: Guyon I , Luxburg UV , Bengio S , Wallach H , Fergus R , Vishwanathan S , et al. , editors. Advances in neural information processing systems . Curran Associates, Inc .; 2017 . 41. ↵ Ioffe S , Szegedy C. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift . arXiv [csLG] . 2015 ; 42. ↵ Nair V , Hinton GE . Rectified linear units improve restricted boltzmann machines . Proceedings of the 27th international conference on international conference on machine learning . Madison, WI, USA : Omnipress ; 2010 . p. 807 – 14 . 43. ↵ Srivastava N , Hinton GE , Krizhevsky A , Sutskever I , Salakhutdinov R. Dropout: A simple way to prevent neural networks from overfitting . J Mach Learn Res . 2014 ; 15 : 1929 – 58 . OpenUrl CrossRef 44. ↵ Oord A van den , Li Y , Vinyals O. Representation learning with contrastive predictive coding . arXiv [csLG] . 2018 ; 45. ↵ Stuart T , Butler A , Hoffman P , Hafemeister C , Papalexi E , Mauck WM 3rd, et al. Comprehensive integration of single-cell data . Cell . 2019 ; 177 : 1888 – 1902 .e21. OpenUrl CrossRef PubMed 46. ↵ Becht E , McInnes L , Healy J , Dutertre C-A , Kwok IWH , Ng LG , et al. Dimensionality reduction for visualizing single-cell data using UMAP . Nat Biotechnol . 2018 ; 47. ↵ Pliner HA , Shendure J , Trapnell C. Supervised classification enables rapid annotation of cell atlases . Nat Methods . 2019 ; 16 : 983 – 6 . OpenUrl CrossRef PubMed 48. ↵ Tayama S , Motohara T , Narantuya D , Li C , Fujimoto K , Sakaguchi I , et al. The impact of EpCAM expression on response to chemotherapy and clinical outcomes in patients with epithelial ovarian cancer . Oncotarget . 2017 ; 8 : 44312 – 25 . OpenUrl CrossRef PubMed 49. ↵ Schummer M , Ng WV , Bumgarner RE , Nelson PS , Schummer B , Bednarski DW , et al. Comparative hybridization of an array of 21,500 ovarian cDNAs for the discovery of genes overexpressed in ovarian carcinomas . Gene . 1999 ; 238 : 375 – 85 . OpenUrl CrossRef PubMed Web of Science 50. ↵ Schaner ME , Ross DT , Ciaravino G , Sorlie T , Troyanskaya O , Diehn M , et al. Gene expression patterns in ovarian carcinomas . Mol Biol Cell . 2003 ; 14 : 4376 – 86 . OpenUrl Abstract / FREE Full Text 51. ↵ Elyada E , Bolisetty M , Laise P , Flynn WF , Courtois ET , Burkhart RA , et al. Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts . Cancer Discov . 2019 ; 9 : 1102 – 23 . OpenUrl Abstract / FREE Full Text 52. ↵ Tang W , Du J , Li L , Hu S , Ma S , Xue M , et al. Hypoxia-related THBD+ macrophages as a prognostic factor in glioma: Construction of a powerful risk model . J Cell Mol Med . 2024 ; 28 : e18393 . OpenUrl CrossRef PubMed 53. ↵ Yang R , Chen Z , Liang L , Ao S , Zhang J , Chang Z , et al. Fc fragment of IgE receptor ig (FCER1G) acts as a key gene involved in cancer immune infiltration and tumour microenvironment . Immunology . 2023 ; 168 : 302 – 19 . OpenUrl 54. ↵ Gao R , Bai S , Henderson YC , Lin Y , Schalck A , Yan Y , et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes . Nat Biotechnol . 2021 ; 39 : 599 – 608 . OpenUrl CrossRef PubMed 55. ↵ Song M , Ma S , Wang G , Wang Y , Yang Z , Xie B , et al. Benchmarking copy number aberrations inference tools using single-cell multi-omics datasets . Brief Bioinform . 2025 ; 26 : bbaf076 . OpenUrl PubMed 56. ↵ Macintyre G , Goranova TE , De Silva D , Ennis D , Piskorz AM , Eldridge M , et al. Copy number signatures and mutational processes in ovarian carcinoma . Nat Genet . 2018 ; 50 : 1262 – 70 . OpenUrl CrossRef PubMed 57. ↵ Shao X , Lv N , Liao J , Long J , Xue R , Ai N , et al. Copy number variation is highly correlated with differential gene expression: A pan-cancer study . BMC Med Genet . 2019 ; 20 : 175 . OpenUrl CrossRef PubMed 58. ↵ Campbell KR , Steif A , Laks E , Zahn H , Lai D , McPherson A , et al. Clonealign: Statistical integration of independent single-cell RNA and DNA sequencing data from human cancers . Genome Biol . 2019 ; 20 : 54 . OpenUrl CrossRef PubMed 59. ↵ Hegen A , Koidl S , Weindel K , Marmé D , Augustin HG , Fiedler U. Expression of angiopoietin-2 in endothelial cells is controlled by positive and negative regulatory promoter elements . Arterioscler Thromb Vasc Biol . 2004 ; 24 : 1803 – 9 . OpenUrl Abstract / FREE Full Text 60. ↵ Ueno T , Linder S , Elmberger G. Aspartic proteinase napsin is a useful marker for diagnosis of primary lung adenocarcinoma . Br J Cancer . 2003 ; 88 : 1229 – 33 . OpenUrl CrossRef PubMed 61. ↵ Yang B , Zhang W , Zhang M , Wang X , Peng S , Zhang R. KRT6A promotes EMT and cancer stem cell transformation in lung adenocarcinoma . Technol Cancer Res Treat . 2020 ; 19 : 1533033820921248 . OpenUrl PubMed 62. ↵ Kolb M , Margetts PJ , Sime PJ , Gauldie J. Proteoglycans decorin and biglycan differentially modulate TGF-beta-mediated fibrotic responses in the lung . Am J Physiol Lung Cell Mol Physiol . 2001 ; 280 : L1327 – 34 . OpenUrl PubMed Web of Science 63. ↵ Mao L , Yang J , Yue J , Chen Y , Zhou H , Fan D , et al. Decorin deficiency promotes epithelial-mesenchymal transition and colon cancer metastasis . Matrix Biol . 2021 ; 95 : 1 – 14 . OpenUrl PubMed 64. ↵ Franzén O , Gan L-M , Björkegren JLM. PanglaoDB: A web server for exploration of mouse and human single-cell RNA sequencing data . Database . 2019; 2019 . 65. ↵ Korotkevich G , Sukhov V , Budin N , Shpak B , Artyomov MN , Sergushichev A. Fast gene set enrichment analysis . Cold Spring Harbor Laboratory . 2021 . p. 060012 . 66. ↵ Mendiratta G , Ke E , Aziz M , Liarakos D , Tong M , Stites EC . Cancer gene mutation frequencies for the U.S. population . Nat Commun . 2021 ; 12 : 5961 . OpenUrl CrossRef PubMed 67. ↵ Johnson WE , Li C , Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods . Biostatistics . 2007 ; 8 : 118 – 27 . OpenUrl CrossRef PubMed Web of Science 68. ↵ Dao T , Fu DY , Ermon S , Rudra A , Ré C. FlashAttention: Fast and memory-efficient exact attention with IO-awareness . arXiv [csLG] . 2022 ; 69. ↵ Hao M , Gong J , Zeng X , Liu C , Guo Y , Cheng X , et al. Large-scale foundation model on single-cell transcriptomics . Nat Methods . 2024 ; 21 : 1481 – 91 . OpenUrl CrossRef PubMed 70. ↵ Higgins I , Matthey L , Pal A , Burgess CP , Glorot X , Botvinick M , et al. Beta-VAE: Learning basic visual concepts with a constrained variational framework . Int Conf Learn Represent . 2016 ; 71. Chen RTQ , Li X , Grosse R , Duvenaud D. Isolating sources of disentanglement in variational autoencoders . arXiv [csLG] . 2018 ; 72. ↵ Kim H , Mnih A. Disentangling by factorising. In: Dy J, Krause A, editors. Proceedings of the 35th international conference on machine learning . Stockholmsmässan, Stockholm Sweden : PMLR ; 2018 . p. 2649 – 58 . 73. ↵ Laks E , McPherson A , Zahn H , Lai D , Steif A , Brimhall J , et al. Clonal decomposition and DNA replication states defined by scaled single-cell genome sequencing . Cell . 2019 ; 179 : 1207 – 1221 .e22. OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted June 06, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks Sishir Subedi , Yongjin P. Park medRxiv 2025.06.04.25328900; doi: https://doi.org/10.1101/2025.06.04.25328900 Share This Article: Copy Citation Tools Decomposing patient heterogeneity of single-cell cancer data by cross-attention neural networks Sishir Subedi , Yongjin P. Park medRxiv 2025.06.04.25328900; doi: https://doi.org/10.1101/2025.06.04.25328900 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Oncology Subject Areas All Articles Addiction Medicine (568) Allergy and Immunology (863) Anesthesia (300) Cardiovascular Medicine (4436) Dentistry and Oral Medicine (444) Dermatology (382) Emergency Medicine (608) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1509) Epidemiology (15229) Forensic Medicine (30) Gastroenterology (1124) Genetic and Genomic Medicine (6600) Geriatric Medicine (668) Health Economics (997) Health Informatics (4538) Health Policy (1368) Health Systems and Quality Improvement (1613) Hematology (542) HIV/AIDS (1264) Infectious Diseases (except HIV/AIDS) (15916) Intensive Care and Critical Care Medicine (1103) Medical Education (623) Medical Ethics (146) Nephrology (667) Neurology (6599) Nursing (346) Nutrition (998) Obstetrics and Gynecology (1144) Occupational and Environmental Health (957) Oncology (3333) Ophthalmology (974) Orthopedics (369) Otolaryngology (420) Pain Medicine (436) Palliative Medicine (130) Pathology (663) Pediatrics (1693) Pharmacology and Therapeutics (691) Primary Care Research (711) Psychiatry and Clinical Psychology (5447) Public and Global Health (9232) Radiology and Imaging (2198) Rehabilitation Medicine and Physical Therapy (1370) Respiratory Medicine (1196) Rheumatology (593) Sexual and Reproductive Health (712) Sports Medicine (530) Surgery (712) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a00e9053d8661b23',t:'MTc3OTY0OTQ0MQ=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();