
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae187
bbae187
Problem Solving Protocol
AcademicSubjects/SCI01060
CoT: a transformer-based method for inferring tumor clonal copy number substructure from scDNA-seq data
Liu Furui School of Information Engineering, Ningxia University, 750021, Ningxia, China

Shi Fangyuan School of Information Engineering, Ningxia University, 750021, Ningxia, China
Collaborative Innovation Center for Ningxia Big Data and Artificial Intelligence Co-founded by Ningxia Municipality and Ministry of Education, Ningxia University, 750021, Ningxia, China

Du Fang School of Information Engineering, Ningxia University, 750021, Ningxia, China
Collaborative Innovation Center for Ningxia Big Data and Artificial Intelligence Co-founded by Ningxia Municipality and Ministry of Education, Ningxia University, 750021, Ningxia, China

Cao Xiangmei Basic Medical School, Ningxia Medical University, 750001, Ningxia, China

https://orcid.org/0000-0001-6526-6991
Yu Zhenhua School of Information Engineering, Ningxia University, 750021, Ningxia, China
Collaborative Innovation Center for Ningxia Big Data and Artificial Intelligence Co-founded by Ningxia Municipality and Ministry of Education, Ningxia University, 750021, Ningxia, China

Corresponding author. Zhenhua Yu, School of Information Engineering, Ningxia University, 750021, No. 489 Helan West Road, Xixia District, Yinchuan City, Ningxia Province, China. Tel: +86-15121975933. Email: zhyu@nxu.edu.cn
5 2024
25 4 2024
25 4 2024
25 3 bbae18706 12 2023
08 3 2024
16 4 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Single-cell DNA sequencing (scDNA-seq) has been an effective means to unscramble intra-tumor heterogeneity, while joint inference of tumor clones and their respective copy number profiles remains a challenging task due to the noisy nature of scDNA-seq data. We introduce a new bioinformatics method called CoT for deciphering clonal copy number substructure. The backbone of CoT is a Copy number Transformer autoencoder that leverages multi-head attention mechanism to explore correlations between different genomic regions, and thus capture global features to create latent embeddings for the cells. CoT makes it convenient to first infer cell subpopulations based on the learned embeddings, and then estimate single-cell copy numbers through joint analysis of read counts data for the cells belonging to the same cluster. This exploitation of clonal substructure information in copy number analysis helps to alleviate the effect of read counts non-uniformity, and yield robust estimations of the tumor copy numbers. Performance evaluation on synthetic and real datasets showcases that CoT outperforms the state of the arts, and is highly useful for deciphering clonal copy number substructure.

copy number alteration
intra-tumor heterogeneity
single-cell sequencing
deep learning
Natural Science Foundation of Ningxia Province 10.13039/501100004772 2023AAC05006 Key Research and Development Program of Ningxia 10.13039/100016692 2023BEG02009 National Natural Science Foundation of China 10.13039/501100001809 62062058
==== Body
pmcINTRODUCTION

Tumor genomes are often featured with various types of mutations, such as copy number alteration (CNA), single nucleotide variation and indel [1, 2]. CNA refers to the amplification or deletion of the chromosomal regions, and its size can scale from several kilobases to the entire chromosome. CNAs in cancers can be divided into large-scale CNAs and focal CNAs, and both types of CNAs frequently occur across multiple cancer types [3, 4]. Krijgsman et al. [5] define a cutoff of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\le $\end{document}3Mb for the size of focal CNAs. CNA and other genomic mutations are continuously accumulated throughout tumor evolution, supporting tumor persistence and growth [1]. As a result, tumor clones that show different mutational profiles, known as the intra-tumor heterogeneity (ITH), emerge in the tumor microenvironment [6]. ITH is a critical factor contributing to drug resistance and cancer treatment failure [7]. To better understand the clonal architecture, one can infer the underlying tumor phylogenetic tree, which depicts the composition of tumor clones (clonal substructure), their respective genomic mutations and their lineage relationship [8, 9]. Therefore, precisely inferring the clonal substructure is essential for the reconstruction of tumor phylogeny, understanding disease state and finding effective drug targets, thus helping design personalized treatment plans for cancers. Specifically, how to accurately identify clonal copy number substructure has been a long-standing interest in the cancer research community.

Single-cell DNA sequencing (scDNA-seq) [10] has provided an unprecedented view of CNAs at single-cell resolution, and shows advantages over bulk sequencing in unscrambling underlying tumor clonal substructure. However, technical issues such as amplification bias and low sequencing coverage make it highly challenging to precisely infer mutations from noisy sequencing data, thus impeding reliable elucidation of ITH [11]. For instance, amplification bias causes the distribution of read counts to be over-dispersed across the genome, and probably introduces systematic errors when measuring copy numbers; ultralow sequencing coverage leads to high sparsity of the read counts, and makes it difficult to distinguish between different tumor clones. Due to these specific complex features of scDNA-seq data, computational methods originally designed for subclonal copy number analysis from bulk sequencing data are not suitable for processing scDNA-seq data. This promotes the emergence of a batch of new methods [12–20] for specifically calling single-cell CNAs.

Single-cell CNA calling methods can be grouped into two classes according to whether negative controls are required as part of the input. The negative controls are typically created by pooling normal cells, whose identities are inferred from or directly specified in the input data. To relieve the effects of local signal fluctuation and high sparsity of the read counts, single-cell based methods usually measure read counts in large-sized bins (e.g. 500kb). In addition, some common preprocessing steps, such as removing outlier bins and cells, correction for GC-content bias and library size normalization, are adopted to enhance the quality of read counts before CNA calling [11]. Nevertheless, existing single-cell CNA calling methods differ significantly in computational models for calling copy numbers. For instance, SCOPE [14] uses a modified Bayesian information criterion to select the optimal number of breakpoints, which divide the genome into distinct segments; SeCNV [18] deduces copy number segments by minimizing structural entropy of a depth congruent map, which quantifies the similarity between any pair of genomic bins. Despite that these statistical methods perform acceptably well in most cases, they have limited generalization ability on complex datasets such as hyperploid tumors [21]. In addition, most of the aforementioned methods do not consider clonal substructure information when calling CNAs, and accurate inference of tumor clonal copy number substructure from scDNA-seq data remains challenging. Identifying clonal substructure can be formulated as a clustering problem, that is, how to aggregate cells having similar copy number profiles into the same subpopulation.

As deep representation learning has superior ability of extracting informative features through nonlinear transformation functions, it has been widely used to learn latent representations and cluster cells from single-cell sequencing data [22–26]. For instance, Seurat [27] is a representative method for analyzing single-cell data and the authors have released several improved versions of the tool [28, 29]; scTAG [30] models count data with the zero-inflated negative binomial (ZINB) distribution and uses a graph embedding autoencoder to learn representations for scRNA-seq cells. Despite that a plenty of methods are available for clustering scRNA-seq data, these methods may not be suitable for identifying tumor clones from scDNA-seq data due to two reasons: (1) the statistical characteristics of scRNA-seq and scDNA-seq data are largely different. scRNA-seq can only quantify copy numbers of expressed genes, and thus often yields highly sparse counts data, while scDNA-seq can profile copy numbers of the whole genome, and the sparsity of counts data can be eschewed using large bins to measure the read counts; (2) clustering methods for scRNA-seq data often implicitly assume independency between genes, and select a predefined number of highly variable genes for downstream analysis, while this strategy is not suitable for preprocessing scDNA-seq data as copy numbers of adjacent genomic regions are closely related [31]. To model read counts from scDNA-seq data, rcCAE [21] leverages a convolutional autoencoder to find a low-dimensional latent space where tumor clones can be easily distinguished between each other, and simultaneously denoise the read counts for better identification of single-cell copy numbers. Despite that rcCAE gains advantages over conventional statistical models in disentangling distinct subpopulations from noisy read counts, the convolutional kernels used in rcCAE tend to capture local copy number patterns, and may be less effective in learning global representations that contribute to better clustering of the cells.

Here, we introduce a novel deep learning method for reasoning tumor clonal copy number substructure from scDNA-seq data. The backbone of our method is a Copy number Transformer (CoT) autoencoder that takes single-cell read counts as input and learns latent global representations of the cells (as shown in Figure 1). Compared with the most relevant competitor rcCAE, our method shows two advantages: (1) CoT is better at learning global representation of the cells. As CoT treats read counts of each cell as sequential data and uses multi-head attention module [32] to model dependencies between genomic regions, copy number patterns of each tumor clone can be more easily captured, which contributes to better representation of the cells. By comparison, rcCAE employs convolutional layers to reason local copy number patterns and then aggregates these features into a latent representation, which makes it less effective in representing the relationship between genomic regions that are far away from each other but closely related in copy number state; (2) CoT performs more robustly against read counts fluctuation when calling single-cell CNAs. rcCAE does not exploit clonal substructure information when calling CNAs, and calculates copy numbers under a cell-by-cell manner. By comparison, CoT jointly analyzes read counts data of the cells from the same subpopulation, thus being able to get more robust estimation of the copy numbers. We showcase the superior performance of CoT on various simulated datasets with different tumor ploidy and data heterogeneity. In addition, clonal copy number substructure inference results on three 10X Genomics datasets demonstrate that our method yields more accurate identification of tumor clones than the state-of-the-art methods.

Figure 1 The model architecture of CoT. The pipeline of CoT consists of three subprocesses including learning cell embeddings with a Transformer autoencoder, identifying cell subpopulations with a GMM and inferring single-cell copy numbers with clone-specific HMMs. The encoder of the Transformer autoencoder consists of a stack of six identical layers, each of which is formed by a multi-head self-attention sub-layer and a feed forward sub-layer. The decoder has the same structure as the encoder. The encoder maps the read counts \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} of a cell into a latent representation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $z$\end{document}, and the decoder reconstructs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $z$\end{document}. The reshaping operation aggregates read counts of every \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document} consecutive bins into a vector, and thus reduces the length of read count sequence to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lceil \frac{M}{p} \rceil $\end{document} (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} denotes the number of bins). Given the inferred cell labels by the GMM, CoT employs an HMM for each of the subclones to estimate copy number segments.

MATERIALS AND METHODS

Simulated and real scDNA-seq datasets

As evolutionary process of tumor can be characterized by a phylogenetic tree, where genomic mutations such as CNAs accumulate along the branches of the tree, we emulate tumor phylogeny by generating a tumor tree followed by attachment of the cells. Each node except root of the tree represents a tumor clone, and each edge <\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document}> denotes newly acquired CNAs by the child node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document} compared with the parent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document}. The tree is constructed by following a finite-site evolutionary model of tumor, and copy number gains/losses of predefined segments are randomly generated for each edge of the tree. Once the tree is constructed, the cells are placed on the tree to emulate different-sized tumor clones. Given copy number profiles of each cell, we use SCSsim [33] to generate scDNA-seq data under 0.02 sequencing coverage. The produced reads are further processed with BWA [34] and SAMtools [35] to generate BAM files. More details about the tree construction and simulation can be found in Supplementary Methods.

We generate various datasets to emulate different tumor ploidy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P$\end{document} (diploidy, triploidy and tetraploidy) and data heterogeneity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $K$\end{document} (i.e. the number of tumor clones). Concretely, the number of cells is fixed to 200, the number of clones is selected from {4, 9, 14, 19} and the simulation is repeated five times for each parameter combination of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $K$\end{document}. This results in 60 scDNA-seq datasets for performance evaluation.

We also get three 10X Genomics datasets from a breast tumor tissue [19] to further evaluate the performance of CoT.

Preprocessing read counts data

We get read count data from fixed-size windows (or bins) along the genome, and employ same approaches as used in [21] to normalize and correct the read counts. After preprocessing, a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N\times M$\end{document} data matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R$\end{document} is obtained to represent the normalized read counts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document} cells in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} bins. Due to the low sequencing coverage of scDNA-seq data, the read counts are measured in large-sized bins (e.g. 500kb), resulting in high-dimensional read counts sequences with more than 5000 bins on human genome. To alleviate the difficulty of handling sequences with high dimensionality, we reshape the size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N\times p\times T$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $T=\lceil \frac{M}{p} \rceil $\end{document} denotes the length of resized sequence and extended elements of the matrix are set to 0. This reshaping operation aggregates read counts of every \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document} consecutive windows to form a feature vector, which makes it appropriate to use current off-shelf Transformer-based approaches to learn representations from sequences. We denote the resulted read counts data by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $X$\end{document} for downstream analysis.

Reasoning clonal copy number substructure

The Transformer models [32, 36–38] leverage attention mechanism to capture the context dependencies between different positions of the input sequence. Inspired by this, we use a Transformer autoencoder to learn representations from single-cell read counts data. Given the reshaped read counts sequence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} (a matrix with size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p\times T$\end{document}) of a cell, the encoder maps it to a latent representation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $z$\end{document}, then the decoder reconstructs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $z$\end{document}. With the self-attention mechanism, the Transformer autoencoder model is able to capture the global dependencies between different genomic regions, and thus deliver better representations of the cells.

An illustration of the model architecture of the Transformer autoencoder is given in Figure 1. Specifically, the encoder consists of a stack of six identical layers, each of which is formed by a multi-head self-attention sub-layer and a feed-forward sub-layer. In addition, residual connection and layer normalization are applied to each sub-layer. The latent layer is a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $d$\end{document}-dimensional fully connected layer, and the decoder has a same structure as the encoder. We add a fully connected layer before the decoder network to restore the data dimension. ReLU activation function is applied to all intermediate layers except the latent layer, and the exponential function is applied to the output of the decoder network. Notably, our model architecture does not contain embedding layer for the input, and all sub-layers produce outputs of dimension equal to the input dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document}.

The network weights are optimized with the ’Adam’ algorithm by minimizing MSE loss between the input and output. After training of the Transformer autoencoder is finished, we get latent representation of each cell with the encoder network, and employ a Gaussian mixture model (GMM) based clustering approach [21] to identify cell subpopulations. We run the GMM 100 times with randomly initialized parameters to find the best result on each dataset.

Inferring single-cell copy numbers

We employ a hidden Markov model (HMM) to analyze the log normalized read counts (LRC) data for each subpopulation, and hidden states of the HMM are defined as copy numbers up to 10. The original LRC data are divided into \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $K$\end{document} subgroups and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document}-th subgroup contains LRC data of the cells belonging to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document}-th subpopulation. Without loss of generality, a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $V\times M$\end{document} matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $Y$\end{document} is used to represent the LRC data of an inferred subpopulation that consists of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $V$\end{document} cells. Given a hidden state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $s$\end{document}, we assume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $Y_{vj}$\end{document} follows a normal distribution with standard deviation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\sigma $\end{document} and mean of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mu _{s}=\log \left ( \frac{c_{s}}{2} \right ) +o$\end{document}; here different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $o$\end{document} correspond to different tumor ploidy and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c_{s}$\end{document} is the copy number in hidden state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $s$\end{document}. Assuming the initial hidden state distribution and state transition probability matrix of the HMM are denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\pi $\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document}, respectively, we aim to infer the maximum likelihood estimations of the model parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\theta =(\pi , A, \sigma , o)$\end{document}, which is resolved by the widely used expectation-maximization (EM) algorithm. Specifically, we update the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\sigma $\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $o$\end{document} in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $n$\end{document}-th iteration of the algorithm as follows:

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{align*} & o^{(n)}=\frac{\sum_{v=1}^{V}\sum_{j=1}^{M}\sum_{s=1}^{S} \gamma_{vjs}^{(n)} \left( Y_{vj}-\log\left(\frac{c_{s}}{2}\right) \right)}{\sum_{v=1}^{V}\sum_{j=1}^{M}\sum_{s=1}^{S} \gamma_{vjs}^{(n)}}; \end{align*}\end{document}

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{align*} & \sigma^{(n)}=\sqrt{\frac{\sum_{v=1}^{V}\sum_{j=1}^{M}\sum_{s=1}^{S} \gamma_{vjs}^{(n)} \left( Y_{vj}-\log\left(\frac{c_{s}}{2}\right)-o^{(n)} \right)^{2}}{\sum_{v=1}^{V}\sum_{j=1}^{M}\sum_{s=1}^{S} \gamma_{vjs}^{(n)}}}, \end{align*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma _{vjs}$\end{document} denotes the posterior probability that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document}-th bin of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $v$\end{document}-th cell is in state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $s$\end{document}. After the model converges, we infer the copy number of each bin for each cell from the hidden state that gains the largest posterior probability, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $z_{vj}=\mathop{\arg \max }_{s} \gamma _{vjs}$\end{document}.

Performance comparison

We compare CoT with several recently published state-of-the-art methods for calling single-cell CNAs, including SCOPE [14], SCICoNE [15], SeCNV [18] and rcCAE [21]. These methods are selected for performance comparison as they differ in methodological basis, and perform well on different datasets as previously reported [21]. CHISEL [19] is excluded from performance evaluation as it additionally calculates allele frequency data under a large bin size (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\sim $\end{document}5Mb). In addition, we also select two methods including Seurat [29] and scTAG [30] that can be used to cluster scRNA-seq data, to assess their performance of identifying tumor clones from scDNA-seq data. For clustering scRNA-seq data, Seurat is a representative method and scTAG is an advanced approach that surpasses the state of the arts. We use these two scRNA-seq-based methods to check if they will still be highly effective when applied to scDNA-seq data, and examine their performance limits regardless of the difference between the two omics data. When evaluating the clonal substructure inference performance, we calculate adjusted rand index (ARI) and normalized mutual information (NMI) by comparing the inferred cell labels with the ground truth, and compute the difference (denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document}) between the ground truth and predicted number of clusters. Similar to the performance metrics adopted in rcCAE, we use accuracy, balanced accuracy and F-measure to indicate copy number accuracy. All experimental results are obtained from a computational server having 1 NVIDIA RTX 3090 GPU, 256GB RAM and 32 CPU cores. More details regarding the calculation of performance metrics and usage of each method are provided in the Supplementary Methods.

RESULTS

Results on simulated datasets

We first apply CoT to the simulated datasets and compare with other methods. Figure 2 shows the clustering accuracy of rcCAE, Seurat, scTAG and CoT. The results indicate that rcCAE and CoT perform similarly on datasets with four tumor clones, while CoT consistently surpasses rcCAE when the data heterogeneity increases. For instance, the median \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document} on triploid datasets with 19 clones is -3 for rcCAE and 0 for CoT. Seurat can yield perfect clustering results on simple datasets with only four or nine tumor clones, while it suffers significant performance degradation on complex datasets with larger number of clones, and tends to underestimate the number of clusters (e.g. mean ARI is 0.75, mean NMI is 0.87 and median \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document} is 5 on datasets containing 19 clones). Despite that the true cell labels are given as part of the input for scTAG (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document} = 0), its ARI and NMI scores change significantly, and tend to decrease when data heterogeneity increases on hyperploid samples. It is also observed that rcCAE’ clustering accuracy tends to decrease with the number of clones, especially on hyperploid samples. By comparison, CoT better maintains the clustering performance, and accurately aggregates the cells into distinct subpopulations on all triploid samples. For instance, the mean ARIs of rcCAE, Seurat, scTAG and CoT on tetraploid samples containing 19 tumor clones are 0.97, 0.76, 0.72 and 0.99, respectively. These results suggest that the self-attention module adopted in CoT is more effective than the convolutional kernels used by rcCAE, and scRNA-seq-based methods are less effective in clustering scDNA-seq data with complex clonal substructure.

Figure 2 Cell clustering performance comparison between CoT and other methods on simulated scDNA-seq datasets. The results on the simulated five diploid datasets, five triploid datasets and five tetraploid datasets are analyzed. Three performance metrics, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document} (difference between the ground truth and predicted number of clusters), ARI and NMI, are compared when the simulated number of clones increases from 4 to 19. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Delta K$\end{document} of scTAG is 0 on all datasets as the ground truth number of clusters is given as part of the input of scTAG.

We further compare the copy number estimation performance of CoT and other methods, and the results are shown in Figure 3. In general, our method performs comparably with or better than existing methods. On the diploid datasets, SeCNV, rcCAE and CoT achieve high accuracies across different numbers of tumor clones. By comparison, SCOPE gets high overall accuracies but yields lower balanced accuracies and F-measures, which may result from misestimation of some specific copy number states; SCICoNE performs poorly as it overestimates the ploidy of most of the cells (Supplementary Fig. 1). Similarly, CoT exhibits high accuracies on the triploid datasets, and generally performs better than the most competitive existing methods including SeCNV and rcCAE. For instance, the mean F-measures of SCOPE, SCICoNE, SeCNV, rcCAE and CoT on datasets having nine tumor clones are 0.57, 0.27, 0.89, 0.92 and 0.97, respectively. For the complex tetraploid samples, CoT also consistently yields high balanced accuracies and F-measures. We also observe that the performance of SCOPE and SeCNV improve when the number of tumor clones increases. The increased data heterogeneity makes it easier to identify genomic breakpoints shared by multiple tumor clones, which may probably contribute to the elevated performance of SCOPE and SeCNV, as they calculate global statistics from the read counts of all cells to segment the genome.

Figure 3 Copy number estimation accuracy of the methods on simulated datasets. Given the estimated and ground truth copy numbers of a cell, three metrics including accuracy, balanced accuracy and macro-averaged F-measure are calculated to make performance comparison. The boxplots showcase how the performance of each method changes with respect to the number of clones.

Additional benchmark testing is performed for the methods on small-sized datasets, checking if they can still deliver reliable results with limited number of cells. The datasets are generated by following the same simulation procedure, with the number of cells and number of clones being set to 100 and 4, respectively. The results in Supplementary Fig. 2 show that SeCNV, rcCAE and CoT achieve high accuracies on diploid samples, while SeCNV suffers from degraded performance on hyperploid datasets. SCOPE and SCICoNE tend to misestimate the cell ploidy (Supplementary Fig. 3), thus yielding lower accuracies in recognizing different copy number states. By comparison, CoT outperforms other methods by providing more robust estimations of the copy number segments and cell ploidy. Taken together, these simulation results demonstrate CoT’s superiority of identifying clonal copy number substructure from complex scDNA-seq data.

Finally, we investigate the effects of hyper-parameters such as bin size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $s$\end{document}, latent dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $d$\end{document} and input dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document}. Specifically, values in {100, 200, 300, 400, 500 kb}, {10, 15, 20} and {32, 64, 128, 256, 512} are tested for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $b$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $d$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document}, respectively. The evaluation results shown in Supplementary Figures 4–12 suggest that our method yields consistently good results across different test conditions, which demonstrates CoT’s robustness against the changes of hyper-parameters. More details about the ablation study can be found in Supplementary Results.

Results on three 10X Genomics datasets

We further apply CoT to three 10X Genomics datasets (denoted as C, D and E, respectively) generated from a breast tumor tissue [19]. These datasets are obtained by sequencing single cells from different tumor sections with sequencing coverages in 0.02X-0.05X. The numbers of cells contained in each dataset are 1254, 1384 and 1446, respectively. The merged BAM file containing read alignments of all cells can be downloaded from https://support.10xgenomics.com/single-cell-dna/datasets/1.0.0/breast_tissue_{C,D,E}_2k. A previous study has employed CHISEL algorithm [19] to jointly analyze all cells from the five sequenced sections of the tumor, and found that there are eight tumor clones and one normal subpopulation. Specifically, there are one normal diploidy and five tumor clones in dataset C, one normal subpopulation accompanied with four tumor clones in dataset D and six subpopulations formed by one normal diploidy and five tumor clones in dataset E. On dataset C, the normal diploidy (463 cells) and two tumor clones (529 and 212 cells, respectively) comprise a large proportion of the cells, while the remaining three clones contain 19, 16 and 15 cells, respectively. Similarly, the dominant subpopulations in dataset D are the normal diploidy (542 cells) and two tumor clones (520 and 270 cells, respectively). For the dataset E, two of the clones have similar copy number profiles but differ in allele-specific copy numbers. Considering the fact that read counts data can only be used to quantify the total copy number of a segment, these two clones could be treated as the same cluster when clustering the cells based on read counts, and therefore we regroup the cells into five clusters, each of which consists of 388, 782, 168, 78 and 30 cells. Cluster 1 is predominantly comprised of diploid normal cells, while clusters 2–5 are aneuploid tumor clones that mainly differ in copy number profiles on chromosomes 2–4. As the results of CHISEL are obtained by joint analysis of read counts and allele frequency data, we think it is appropriate to exploit the results of CHISEL as the ground truth to assess the performance of CoT.

The cell clustering results of CoT are shown in Figure 4. CoT aggregates the cells into five clusters on datasets C and D, and finds eight clusters on dataset E. We also obtain the clustering results of rcCAE, Seurat and scTAG on the same datasets to make a performance comparison (Figure 4 and Supplementary Table 1). The results show that our method outperforms other methods on all three datasets by delivering higher ARI and NMI scores. For instance, ARIs of CoT, rcCAE, Seurat and scTAG on dataset E are 0.945, 0.687, 0.885 and 0.939, respectively. In addition, cells belonging to different clusters are more clearly separated in the latent space learned by CoT. For instance, CoT achieves almost perfect clustering results on dataset D by correctly distinguishing between distinct cell subpopulations, rcCAE wrongly assigns the cells from two clones to the cluster 5 and Seurat fails to identify a minor clone. Despite that the ground truth cell labels are given as part of the input of scTAG, it is less effective than CoT in identifying clonal substructure. These results suggest performance bottleneck occurs when employing scRNA-seq-based methods to cluster scDNA-seq data due to intrinsic differences between the two omics data. To check if the attention model used in CoT does improve the clustering performance, we replace each multi-head self-attention module with a convolutional layer with similar number of parameters, and compare the clustering accuracy. The results in Supplementary Table 1 imply the attention model enables more accurate clustering of the cells, suggesting that our model architecture with the Transformer component is more effective.

Figure 4 Cell clustering results of CoT and other methods on three 10X Genomics datasets. The predicted cell labels of CHISEL algorithm [19] are used as the ground truth to compare the clustering performance. Based on the learned cell embeddings, the first two principal components are used to visualize the clustering results. The results show that CoT outperforms other methods on all three datasets.

We further analyze the label distribution of CoT, rcCAE, Seurat and scTAG by comparing the predicted and ground truth cell labels, and the results are shown in Supplementary Fig. 13. All methods are able to identify the most abundant subpopulations on each dataset, and CoT is more powerful in recognizing minor clusters. In general, CoT yields more pure clusters of which each is dominantly comprised of the cells from a same subpopulation, while other methods have higher probability to aggregate cells from different subpopulations into a same cluster. For instance, the predicted cluster 5 of rcCAE on dataset C completely contains CHISEL cluster 3 and partially encompasses the cells from CHISEL cluster 6; inferred cluster 3 by Seurat on dataset D entirely comprises CHISEL clusters 2 and 3; cluster 2 of scTAG on dataset E fully covers CHISEL cluster 5 and majority of the cells from CHISEL cluster 4. The elevated clustering performance of CoT benefits from the utilization of multi-head attention mechanism, which is profit for better learning of global representations of the cells.

We then evaluate the copy number estimation accuracy of rcCAE and CoT (other methods are not evaluated as the required per-cell BAMs are not available for the 10X Genomics datasets). As the ground truth single-cell copy number profiles are unknown, we take CHISEL’s results as the golden standard to calculate accuracy, balanced accuracy and F-measure. As shown in Supplementary Fig. 14, CoT generally obtains higher values than rcCAE on three performance metrics. For instance, the mean F-measure of CoT on three datasets are 0.84, 0.83 and 0.81, respectively, while the corresponding values of rcCAE are 0.75, 0.73 and 0.78, respectively. The results imply that our joint modeling of the cells belonging to the same cluster enables improved inference of copy number segments.

We further investigate the difference of copy number profiles between distinct clusters. Specifically, for each cluster we select the cell that has the shortest distance to the center of the cluster as a representative cell, and make a comparison between these representative cells. The copy number estimation results of CoT on dataset E are depicted in Figure 5, and the results on datasets C and D are given in Supplementary Figures 15–16. As shown in Figure 5 and Supplementary Fig. 13, the cells belonging to CHISEL cluster 1 are divided into two subgroups (clusters 1 and 2) by CoT. Cluster 1 is the normal diploid subpopulation that is featured with neutral copy number across the genome, while cluster 2 contains only four cells that have hyperploid genomes. These cells of cluster 2 show large distances to cluster 1 in the latent space (as shown in Figure 4), and thus are assigned to a separate cluster by CoT. Similarly, our method also divides the cells of CHISEL cluster 4 into two subgroups (clusters 3 and 4) with ploidy of 3.55 and 3.50, respectively. The main difference is observed on chromosome 11, where cluster 3 has copy number of 4 on p15.5-q23.1 and copy number of 2 on q23.1-q25, while the corresponding copy numbers of cluster 4 are 3 and 1, respectively. It is noted that these differences are not revealed by rcCAE’ results (clusters 6 and 7 in figure 6 of the rcCAE paper). This suggests that our method can identify more refined copy number difference between tumor clones. Cluster 5 aligns well with the CHISEL cluster 2, and this tumor shows copy number of 3 on chromosomes 2–3, copy number of 4 on chromosomes 4–5 and copy number of 3 on chromosome 6p. There are similar copy number differences between clusters 6 and 7 as observed between clusters 3 and 4, and these two clusters align well with the CHISEL cluster 3. The main difference between clusters 3–4 and 6–7 is shown on chromosome 2, where the former has copy number of 3, while the latter has copy number of 4. Clusters 4 and 7 identified by CoT do not match any clusters of rcCAE, and they all show copy number deletion on q23.1–q25. To give a more detailed illustration of the difference, we select three cells from cluster 4 to compare copy number estimation results of CoT and rcCAE (Supplementary Fig.17). As LRC data of the genomic regions with low-prevalence copy number states cannot be well reconstructed, rcCAE tends to misestimate copy numbers of these regions (Supplementary Fig. 17). For instance, reconstructed LRC data of the chromosome 7q of the cell ‘AGAGTGGGTGTGATCT” heavily deviate from the expected values, and this makes rcCAE predict the copy number as 4, while CoT still recognizes the correct copy number of 3. Finally, cluster 8 completely encompasses CHISEL cluster 5, which suggests they represent a same tumor clone. The distinctive feature of this clone is that it has a copy number of 3 on chromosomes 2–4. By making a comparison between the predicted clusters on datasets C and E, we find the following mapping relationships: (C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{1}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{1}$\end{document}), (C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{3}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{3}$\end{document}), (C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{4}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{5}$\end{document}) and (C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{5}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{8}$\end{document}); here C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{i}$\end{document} denotes cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} of C and E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{j}$\end{document} denotes cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} of E. Similarly, mapping relationships (D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{1}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{1}$\end{document}), (D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{3}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{3}$\end{document}), (D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{4}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{5}$\end{document}) and (D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{5}$\end{document}:E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $_{6}$\end{document}) are found between the clusters of D and E. In addition, cluster 2 of C and cluster 2 of D seem to represent a same tumor clone as they present highly similar copy number profiles.

Figure 5 Copy number estimation results of CoT on real dataset E. CoT identifies eight clusters including one normal subpopulation and seven tumor clones. For each cluster, the cell that has the shortest distance to the center of the cluster is selected as the representative cell. The copy number profiles of the representative cells are compared. Cluster 1 represents the normal diploid subpopulation, cluster 2 denotes a tumor clone with hyperploid genome and other tumor clones show similar copy number profiles on most of the chromosomes. The main differences among clusters 3–8 are observed on chromosomes 2, 3, 4 and 11.

DISCUSSION AND CONCLUSIONS

As evolution of tumor can be seen as a reiterative clonal expansion process, tumor clones exhibiting genetic diversity emerge in the tumor during this process [2]. ITH is one of the major determinants of drug resistance and cancer treatment failure, and is associated with poor prognosis of the patient [39]. In addition, different cancer types can exhibit high discrepancy in the degree of ITH [6]. Therefore, accurately identifying clonal substructure is important to reconstruct the tumor phylogeny, and reveal the key genomic mutations that drive tumor progression [40].

Single-cell DNA sequencing technology makes it convenient to identify CNAs at single-cell resolution, and unmask tumor clonal substructure without a complicated deconvolution procedure as required in bulk sequencing. However, specific technical issues such as amplification bias and low sequencing coverage make it difficult to accurately elucidate ITH from scDNA-seq data. To decipher clonal copy number substructure, we introduce a new method CoT that utilizes a Transformer autoencoder as the backbone, to learn cell embeddings with strong differentiation from scDNA-seq data, and thus enable accurate identification of distinct subpopulations by clustering. Furthermore, through joint analysis of read counts data of the cells from the same subpopulation, CoT is able to alleviate the effect of read counts non-uniformity, thus delivering reliable and robust copy number estimation results. Through comprehensive evaluations of CoT on simulated datasets and the 10X Genomics breast cancer datasets, we demonstrate that it surpasses the other methods in calling tumor clones and single-cell CNAs. We would like to emphasize that our method can also be applied to other types of tumors.

Despite that the adopted Transformer model shows better capability in distinguishing between tumor clones, there are still two deficiencies that may limit the practicality of CoT. First, CoT uses read counts in genomic bins with size of hundreds of kilobases to calculate copy numbers, while accurately measuring copy numbers of small-sized bins (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\le $\end{document}10 kb) is essential for finding focal CNAs and is still an unresolved problem. Due to the ultra-low sequencing coverage and large size of the genome, the resulting read counts data will be of high sparsity and high dimensionality, which complicates the downstream ITH and CNA analysis. The ZINB-based model may be a feasible solution as it has shown advantages in modeling the distribution of sparse scRNA-seq data [26, 41]. Secondly, the multi-head attention mechanism employed in CoT is conducive for learning global features of cells, while copy numbers are highly relevant between adjacent genomic regions and often appear as local patterns, which cannot be fully captured by the Transformer submodules. As convolutional neural network (CNN) is highly effective in recognizing local patterns, hybrid architectures that combine CNN with Transformer have presented higher capability in fusing local features with the global representations of the data [42, 43]. Therefore, building a hybrid CNN-Transformer architecture is another feasible direction to improve CoT.

In summary, we present a new bioinformatics method CoT, to identify cell subpopulations as well as detect single-cell CNAs from scDNA-seq data. We think it will be useful for gaining deeper understanding of the ITH as well as the evolutionary mechanism of cancer, and aid in the prediction of prognosis and outcomes of cancer patients.

Key Points

We introduce a transformer-based method CoT for dissecting tumor clonal copy number substructure from scDNA-seq data.

CoT is able to capture the copy number relationship between different genomic regions and thus more effective in yielding global representations for single cells.

CoT accurately reasons tumor clonal substructure and estimates single-cell copy numbers on various simulated datasets and three 10X Genomics datasets.

Supplementary Material

supp_bbae187

FUNDING

This work was supported in part by the Natural Science Foundation of Ningxia Province (2023AAC05006), Key Research and Development Program of Ningxia (2023BEG02009) and National Natural Science Foundation of China (62062058).

AUTHORS’ CONTRIBUTIONS

Z.Y. and F.D. conceived the study. F.L. implemented the method, analyzed the data and drafted the manuscript. F.S., X.C. and Z.Y. improved the manuscript. All authors read and approved the final manuscript for publication.

DATA AVAILABILITY

The sequencing data of three 10X Genomics datasets are available at https://support.10xgenomics.com/single-cell-dna/datasets/1.0.0/breast_tissue_{C,D,E}_2k.

SOFTWARE AVAILABILITY

An implementation of the CoT method is available at https://github.com/zhyu-lab/cot.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE

Not applicable.

Author Biographies

Furui Liu is a Master’s student at School of Information Engineering, Ningxia University. His research interests include deep learning and bioinformatics.

Fangyuan Shi is an assistant professor at School of Information Engineering, Ningxia University. Her research interests include bioinformatics and cancer genomics.

Fang Du is a professor at School of Information Engineering, Ningxia University. Her research interests include biomedical data mining and knowledge graph.

Xiangmei Cao is an associate professor at Basic Medical School, Ningxia Medical University. Her research interests include bioinformatics and pathogenic mechanism of glioma.

Zhenhua Yu is an associate professor at School of Information Engineering, Ningxia University. His research interests include bioinformatics, computational biology and deep learning.
==== Refs
References

1. Stratton MR , CampbellPJ, AndrewP, Futreal. The cancer genome. Nature 2009;458 (7239 ): 719–24.19360079
2. Turajlic S , SottorivaA, GrahamT, SwantonC. Resolving genetic heterogeneity in cancer. Nat Rev Genet 2019;20 (7 ):404–16.30918367
3. Beroukhim R , MermelCH, PorterD, et al. The landscape of somatic copy-number alteration across human cancers. Nature 2010;463 (7283 ):899–905.20164920
4. Boutros PC , FraserM, HardingNJ, et al. Spatial genomic heterogeneity within localized, multifocal prostate cancer. Nat Genet 2015;47 (7 ):736–45.26005866
5. Krijgsman O , CarvalhoB, MeijerGA, et al. Focal chromosomal copy number aberrations in cancer—needles in a genome haystack. Biochim. Biophys. Acta, Mol. Cell. Res. 2014;1843 (11 ):2698–704.
6. McGranahan N , SwantonC. Clonal heterogeneity and tumor evolution: past, present, and the future. Cell 2017;168 (4 ):613–28.28187284
7. Vasan N , BaselgaJ, HymanDM. A view on drug resistance in cancer. Nature 2019;575 (7782 ):299–309.31723286
8. Kuipers J , JahnK, BeerenwinkelN. Advances in understanding tumour evolution through single-cell sequencing. Biochim. Biophys. Acta, Rev. Cancer 2017;1867 (2 ):127–38.28193548
9. Jiaqi W , El-KebirM. Clonarch: visualizing the spatial clonal architecture of tumors. Bioinformatics 2020;36 (Supplement_1 ):i161–8.32657368
10. Gawad C , KohW, QuakeSR. Single-cell genome sequencing: current state of the science. Nat Rev Genet 2016;17 (3 ):175–88.26806412
11. Mallory XF , EdrisiM, NavinN, NakhlehL. Methods for copy number aberration detection from single-cell dna-sequencing data. Genome Biol 2020;21 (1):208.
12. Garvin T , AboukhalilR, KendallJ, et al. Interactive analysis and assessment of single-cell copy-number variations. Nat Methods 2015;12 (11 ):1058–60.26344043
13. Wang X , ChenH, ZhangNR. Dna copy number profiling using single-cell sequencing. Brief Bioinform 2018;19 (5 ):731–6.28159966
14. Wang R , LinD-Y, JiangY. Scope: a normalization and copy-number estimation method for single-cell dna sequencing. Cell Syst 2020;10 (5 ):445–452.e6.32437686
15. Kuipers J , TuncelMA, FerreiraP, et al. Single-cell copy number calling and event history reconstruction. bioRxiv 2020;2020.04.28.065755. 10.1101/2020.04.28.065755.
16. Feng X , ChenL, QingY, et al. Scyn: single cell cnv profiling method using dynamic programming. BMC Genomics 2021;22 (5 ):651.34789142
17. Hui S , NielsenR. Sconce: a method for profiling copy number alterations in cancer evolution using single-cell whole genome sequencing. Bioinformatics 2022;38 (7 ):1801–8.35080614
18. Ruohan W , YuweiZ, MengboW, et al. Resolving single-cell copy number profiling for large datasets. Brief Bioinform 2022;23 (4):bbac264.
19. Zaccaria S , RaphaelBJ. Characterizing allele-and haplotype-specific copy numbers in single cells with chisel. Nat Biotechnol 2021;39 (2 ):207–14.32879467
20. Hui S , NielsenR. SCONCE2: jointly inferring single cell copy number profiles and tumor evolutionary distances. BMC Bioinformatics 2022;23 (1 ):348.35986254
21. Zhenhua Y , LiuF. Fangyuan Shi, and Fang Du. Rccae: a convolutional autoencoder method for detecting intra-tumor heterogeneity and single-cell copy number alterations. Brief Bioinform 2023;24 (3 ):bbad108.36961311
22. Grønbech CH , VordingMF, TimshelPN, et al. Scvae: variational auto-encoders for single-cell gene expression data. Bioinformatics 2020;36 (16 ):4415–22.32415966
23. Ciortan M , DefranceM. GNN-based embedding for clustering scRNA-seq data. Bioinformatics 2021;38 (4 ):1037–44.
24. Li G , ShaliuF, WangS, et al. A deep generative model for multi-view profiling of single-cell rna-seq and atac-seq data. Genome Biol 2022;23 (1 ):20.35022082
25. Yan J , MaM, ZhenhuaY. Bmvae: a variational autoencoder method for clustering single-cell mutation data. Bioinformatics 2023;39 (1 ):btac790.36478203
26. Ren L , WangJ, LiZ, et al. Scmcs: a framework for single-cell multi-omics data integration and multiple clusterings. Bioinformatics 2023;39 (4 ):btad133.36929930
27. Satija R , FarrellJA, GennertD, et al. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol 2015;33 (5 ):495–502.25867923
28. Hao Y , HaoS, Andersen-NissenE, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184 (13 ):3573–3587.e29.34062119
29. Hao Y , StuartT, KowalskiMH, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42 (2):293–304.
30. Zhuohan Y , YifuL, WangY, et al. Zinb-based graph embedding autoencoder for single-cell rna-seq interpretations. In Proceedings of the AAAI Conference on Artificial Intelligence 2022;36 (4):4671–4679, published by AAAI Press, Palo Alto, California USA.
31. Markowska M , CakałaT, MiasojedowB, et al. Conet: copy number event tree model of evolutionary tumor history for single-cell data. Genome Biol 2022;23 (1 ):1–35.34980209
32. Vaswani A , ShazeerN, ParmarN, et al. Attention is all you need. In Proceedings of the 31st International Conference on Neural Information Processing Systems 2017;30 :6000–6010, published by Curran Associates Inc, Morehouse Lane, Red Hook, NY, USA.
33. Zhenhua Y , FangD, SunX, LiA. Scssim: an integrated tool for simulating single-cell genome sequencing data. Bioinformatics 2020;36 (4 ):1281–2.31584615
34. Li H , DurbinR. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics 2009;25 (14 ):1754–60.19451168
35. Li H , HandsakerB, WysokerA, et al. The sequence alignment/map format and samtools. Bioinformatics 2009;25 (16 ):2078–9.19505943
36. Chen J , HaoX, TaoW, et al. Transformer for one stop interpretable cell type annotation. Nat Commun 2023;14 (1):223.
37. Yang F , WangW, WangF, et al. scBERT as a large-scale pretrained deep language model for cell type annotation of single-cell RNA-seq data. Nat Mach Intell 2023;4 (10 ):852–66.
38. Jing X , ZhangA, LiuF, ZhangX. Stgrns: an interpretable transformer-based method for inferring gene regulatory networks from single-cell transcriptomic data. Bioinformatics 2023;39 (4 ):btad165.37004161
39. Patel AP , TiroshI, TrombettaJJ, et al. Single-cell rna-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 2014;344 (6190 ):1396–401.24925914
40. Schwartz R , SchäfferAA. The evolution of tumour phylogenetics: principles and practice. Nat Rev Genet 2017;18 (4 ):213–29.28190876
41. Wan H , ChenL, DengM. scNAME: neighborhood contrastive clustering with ancillary mask estimation for scRNA-seq data. Bioinformatics 2022;38 (6 ):1575–83.34999761
42. Zhang J , LiC, LiuG, et al. A cnn-transformer hybrid approach for decoding visual neural activity into text. Comput Methods Programs Biomed 2022;214 :106586.
43. He Q , YangQ, XieM. Hctnet: a hybrid cnn-transformer network for breast ultrasound image segmentation. Comput Biol Med 2023;155 :106629.
