Method for identifying single cell ecDNA cloning trajectory

Abstract

The invention discloses an identification method of a single cell ecDNA cloning trajectory. Potential tumor markers in single cell ATAC data are identified based on ecDNA cloning characteristics. The cancer marker is recognized through the ecDNA, tumor cells can be sensitively recognized, and early diagnosis can be carried out; the copy number, the structure and the distribution of the ecDNA often have great heterogeneity in tumor cells. Different types of ecDNA clones may be carried among different tumor cell populations, so that the evolution trajectory of tumors and the diversity of the cell populations can be disclosed; according to the present invention, the gene structure change of the tumor cell population in different stages can be revealed by researching the clone trajectory of the ecDNA so as to provide the potential marker for the early screening and monitoring of cancers.

Description

Technical Field

[0001] The present invention relates to the field of tumor marker identification, and in particular to a method for identifying single-cell ecDNA clone trajectories. Background Art

[0002] Recent studies have shown that ecDNA often contains specific oncogenes, drug resistance genes, immune regulatory genes, and other regulatory elements, making it an important factor in cancer pathogenesis and poor prognosis in multiple cancer types. However, many aspects of ecDNA remain unclear, including its role in aberrant gene expression, cancer evolution, biogenesis, maintenance, and clearance. Given the high oncogene copy number and enhanced chromatin accessibility of ecDNA, ATAC-seq is a powerful tool for studying the functional role of ecDNA in enhancer hijacking, cross-domain regulation, tumor heterogeneity, and uneven distribution.

[0003] A key feature of ecDNA is its high copy number and enriched promoter regions, which can significantly increase the expression levels of related oncogenes. For example, oncogenes such as KLF7, MYC, and EGFR are often found on ecDNA, and their expression is closely related to tumor aggressiveness and drug resistance (Warburton, D., et al. (2020). "Extrachromosomal DNA in cancer." Nature Reviews Cancer, 20(4), 201-213.; Mertens, F., et al. (2021). "Extrachromosomal DNA and its role in cancer evolution and drug resistance." Nature Communications, 12(1), 2713.). The application of scATAC-seq and single-cell genomics technologies has enabled researchers to deeply explore the structure and function of ecDNA, revealing its complex role in cancer. In particular, in various cancer types, ecDNA has been found to regulate gene expression through "enhancer hijacking", thereby supporting tumor evolution (Lee J., et al. (2019). "EcDNA-mediated oncogene expression and cancer evolution." Nature Genetics, 51(3), 359-366).

[0004] Despite the increasing importance of ecDNA in cancer research, many questions about its biological functions, formation mechanisms, and how ecDNA is maintained and eliminated remain unclear. Existing research on ecDNA heterogeneity and clonal dynamics at the single-cell level is still relatively limited. This invention identifies ecDNA based on single-cell ATAC-seq data and identifies tumor target molecules based on their clonal characteristics. By dynamically tracking the clonal characteristics of ecDNA, we can reveal the expression regulation patterns of specific oncogenes and their role in tumor heterogeneity. Summary of the Invention

[0005] The purpose of this application is to provide a tumor target molecule identification method based on ecDNA cloning characteristics, and to identify tumor target molecules based on ecDNA and its cloning characteristics identified by single-cell ATAC-seq data.

[0006] In order to solve the above technical problems, the present application provides the following technical solution: a method for identifying single-cell ecDNA clone trajectories, the method operates as follows: S101. Collect paired normal tissue samples, precancerous lesion tissue samples, and cancer tissue samples from multiple patients for single-cell sequencing; S102, performing cell annotation on the pre-processed single-cell sequencing data; S103, identifying cells containing ecDNA and characteristic genes of ecDNA in the single-cell sequencing data, constructing an ecDNA feature matrix, performing dimensionality reduction clustering on the feature matrix, and identifying the central cells of the clones; S104: ecDNA clone trajectory analysis and visualization, and selected diffuse genes for trajectory analysis to identify single-cell ecDNA clone trajectories.

[0007] As a further description of the above scheme: The detailed operations of single-cell data preprocessing and cell annotation in S102 are as follows: S201, Cell Ranger (cell sequencing) ATAC data upstream processing pipeline; scATAC-seq data were processed using the 10× Genomics Cell Ranger ATAC pipeline, which includes aligning FASTQ files to the reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin regions (peaks), and outputting fragment files. The resulting fragments.tsv.gz and peaks.bed files were used for downstream analysis. S202, construct Seurat object; The SignacR package was used to process the output of Cell Ranger. The peaks.bed file for each sample was read and converted into a GenomicRanges object. The singlecell.csv metadata file containing cell filter information was loaded, and high-quality cells with passed_filters > 500 were retained to ensure data reliability. Peak data from all samples were then merged into a unified combined.peaks set. The reduce function was used to eliminate overlapping regions and filter out atypical peak regions with lengths < 20 bp or > 10,000 bp. For each sample, a FragmentObject was created to store cell-specific fragment information. A feature matrix (FeatureMatrix) was generated based on the fragment information, with rows representing peaks and columns representing cells. This matrix was encapsulated as a ChromatinAssay object and integrated into a Seurat object, along with sample metadata. Finally, the Seurat objects of all samples were merged to form an integrated dataset. A dataset field was added to identify the sample source and supplemented with human genome annotation information. S203, perform quality control on Seurat objects; First, the nucleosome signal was calculated to assess the characteristics of DNA open regions, and cells with a nucleosome signal of nucleosome_signal>4 were filtered out to remove nucleosome contamination. Second, the TSS enrichment score was calculated to assess the degree of chromatin openness near the transcription start site, and high-quality cells with TSS.enrichment>3 were retained. At the same time, the peak fragment ratio (pct_reads_in_peaks>15%) was analyzed to ensure that the cells had sufficient chromatin open region signal. Finally, the blacklist region ratio was calculated to exclude low-quality cells with a blackbox region ratio of blacklist_ratio<0.05. This series of rigorous quality control steps was verified using multiple visualization methods such as DensityScatter, TSSPlot, FragmentHistogram, and VlnPlot. S204, cell type annotation; The preprocessed peaks data were normalized using the signacR package for TF-IDF, and highly expressed characteristic peaks with a total count > 10 were selected. Singular value decomposition (SVD) was then used to reduce the data dimensionality. A K-nearest neighbor graph was constructed based on the dimensionality reduction results, and the Louvain algorithm was used to cluster cells into multiple characteristic clusters. The gene activity matrix was then calculated by quantifying the number of open chromatin fragments in gene bodies and promoter regions, and then normalized by LogNormalization. Finally, the differential activity characteristics of marker genes of various cell types in the tumor microenvironment were analyzed in each cluster, and each cluster was annotated as a specific cell type based on the characteristic expression profile.

[0008] As a further description of the above scheme, the process of S103 identifying cells containing ecDNA and the characteristic genes of ecDNA in single-cell data is as follows: S301, screening of amplified fragments and copy number variation analysis; To identify cells carrying ecDNA in single-cell ATAC-seq data, it is first necessary to detect abnormally amplified regions on the genome. Using the AMULET tool, CNV analysis is performed on the sequencing data of each single cell to identify genomic intervals with abnormally high copy numbers and CNVs greater than 5. To ensure data reliability, stringent filtering criteria were set, including alignment quality MAPQ > 5, reasonable insert size ≤ 300 bp, complete alignment of paired-end reads to the reference genome hg38, and removal of potential interference factors such as PCR duplicates. Based on these criteria, a binary matrix was constructed, in which rows represent detected amplification sites and columns represent different single cells, thereby systematically characterizing the distribution pattern of amplified fragments in the cell population. S302, detection and characterization of genomic breakpoints; The formation of ecDNA is often accompanied by genomic rearrangement events, so detecting breakpoints is a key step in identifying ecDNA. Potential breakpoints are identified by analyzing the characteristics of aligned reads, mainly in the following three situations: (1) Large inserts >500 bp, indicating that there may be large-scale rearrangements in the genomic sequence; (2) Cross-chromosome alignment, that is, paired reads are aligned to different chromosomes; (3) Co-directional alignment, that is, the alignment direction of the two reads is the same, which suggests the possible presence of inversions or circular DNA structures; all detected breakpoints retain the cell barcode information, ensuring that subsequent analysis can be performed at single-cell resolution; S303, identification and structural reconstruction of ecDNA; To accurately identify ecDNA, amplified regions with CNVs (gene copy number variations) greater than 2 are expanded upstream and downstream, with a default extension of 1000-3000 bp to cover the possible breakpoint junction regions. Subsequently, the bedtools tool is used to calculate the genomic overlap between the expanded amplified fragments and the breakpoints. Cell barcode matching is then used to verify that the amplified fragments and breakpoints originate from the same cell. Upon successful verification, the circular structure of the ecDNA is reconstructed, and the complete genomic composition is inferred. S304, characteristic gene annotation and functional analysis; ecDNA often carries genes that are closely related to the occurrence and development of cancer. To analyze the biological function of ecDNA, the ChIPseeker tool is used to annotate the expanded amplified fragments and identify the important characteristics of the coding genes and regulatory elements contained therein. In addition, the ecDNAscope tool outputs detailed structural information of ecDNA in each cell, presenting it in the form of an ordered list of fragments, and sorting and visualizing ecDNA according to copy number, number of breakpoints, and carrier cell frequency parameters.

[0009] As a further description of the above scheme: Constructing the ecDNA feature matrix specifically includes the following steps: After obtaining the ecDNA-related gene list, the ecDNA feature matrix is ​​constructed. The ecDNA feature matrix is ​​a numerical matrix with rows representing genes and columns representing samples. It is obtained by traversing the target gene list, extracting the rows corresponding to each gene from the copy number data frame of multiple samples in turn, and calculating the total copy number of the gene in all samples.

[0010] As a further description of the above scheme: Dimensionality reduction clustering of the feature matrix specifically includes the following steps: Based on the Seurat R package, the ecDNA matrix was used to create a Seurat object. The NormalizeData function was used to normalize the data to eliminate the effects of sequencing depth or library size, making the expression of ecDNA signature genes between different cells comparable. The FindVariableFeatures function was used to calculate signature genes with high coefficients of variation. After normalizing the data using the ScaleData function, principal component analysis was performed on the signature matrix, and the first 20 principal components were selected as the basis for subsequent analysis. In the cluster analysis stage, the FindClusters function was used to set the resolution parameter to 2 for cell clustering.

[0011] As a further description of the above scheme: the identification of clonal central cells specifically includes the following steps: S01, identifying marker genes of ecDNA based on threshold values; S02. The cell with the highest expression level of marker genes in each cluster is defined as the clone center cell; The threshold values ​​are: the gene must be expressed in at least 1% of cells; and the fold change value after log2 conversion must be greater than or equal to 0.1.

[0012] As a further description of the above protocol, the detailed steps for ecDNA clone trajectory analysis and visualization are as follows: For each cluster, the sum of the expression values ​​of each ecDNA feature in all cells in the cluster was calculated; Pearson correlation coefficients between different clusters were calculated based on the sum; Two clusters with a Pearson correlation coefficient greater than 0.3 were defined as continuous clonal trajectories; Use UMAP (non-linear dimensionality reduction algorithm) to visualize the clonal evolution of clusters with low ecDNA copy numbers to clusters with higher counts; The detailed steps for selecting diffuse genes for trajectory analysis are as follows: After obtaining the ecDNA cloning trajectory, the clonal relationships between all clusters are collected. After selecting a specific cloning trajectory, ecDNA characteristic genes with significant differences in copy number between clusters are calculated as diffuse genes. Then, a heat map of diffuse genes between clusters is drawn to visualize the changes in diffuse genes between clusters. The significant difference is defined as the cumulative change in copy number of a certain ecDNA feature between different clusters in the cloning trajectory with a step size of n, where n can be customized.

[0013] The characteristics or innovations of the present invention are as follows: 1) Highly sensitive identification at the single-cell level. By integrating single-cell ATAC-seq data with ecDNA clonal profiling, tumor-specific ecDNA and the oncogenes it carries (such as KLF7 and MYC) can be accurately identified at single-cell resolution, significantly improving the sensitivity of early cancer diagnosis. Compared to traditional methods, this approach can capture subtle heterogeneity in ecDNA within tumor cell populations, providing molecular markers for early screening.

[0014] 2) Dynamically analyze the evolutionary trajectory of tumor clones. Based on the copy number, structural heterogeneity, and clonal distribution characteristics of ecDNA, the diversity of tumor cell subpopulations and the dynamic evolution of their clones are revealed. By tracking the spatial distribution and temporal evolution of different ecDNA clones, the molecular mechanisms of key processes such as tumor resistance and metastasis can be analyzed, providing new strategies for targeted therapy.

[0015] Compared with existing technologies, the present invention has the following beneficial technical effects: By identifying cancer markers through ecDNA, the present invention can sensitively identify tumor cells and perform early diagnosis; the copy number, structure, and distribution of ecDNA often have great heterogeneity in tumor cells. Different tumor cell populations may carry different types of ecDNA clones, and the present invention helps to reveal the evolutionary trajectory of tumors and the diversity of their cell populations. By studying the cloning trajectory of ecDNA, the present invention can reveal the changes in the genetic structure of tumor cell populations at different stages, thereby providing potential markers for early screening and monitoring of cancer. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 Figure 3 is the ecDNA clonal trajectory; a shows the ecDNA clonal trajectory of 27 clusters identified based on single-cell ATAC-seq dimensionality reduction clustering; b shows the correlation of ecDNA clones between different clusters; c shows the ecDNA genes shared between different clusters; d shows the correlation between ecDNA and cell cycle markers; e shows the diffuse gene expression along the ecDNA clonal trajectory; f shows the distribution of all sample types containing ecDNA cells; g shows the key ecDNA features in the evolution of skin cancer along the clonal trajectory. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all of the embodiments. In the following description, specific details such as specific configurations and components are provided only to help fully understand the embodiments of the present application. Therefore, it should be clear to those skilled in the art that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of the present application. In addition, for clarity and brevity, the description of known functions and structures has been omitted in the embodiments.

[0018] It should be understood that references throughout this specification to "one embodiment" or "this embodiment" mean that a particular feature, structure, or characteristic associated with the embodiment is included in at least one embodiment of the present application. Therefore, the appearance of "one embodiment" or "this embodiment" throughout this specification does not necessarily refer to the same embodiment. Furthermore, these particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0019] In addition, the present application may repeat reference numerals and / or letters in different examples. This repetition is for the purpose of simplicity and clarity and does not in itself indicate the relationship between the various embodiments and / or settings discussed.

[0020] The term "and / or" in this article is only a description of the association relationship of associated objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, B exists alone, and A and B exist at the same time. The term " / and" in this article describes another type of association object relationship, indicating that two relationships can exist. For example, A / and B can mean: A exists alone, and A and B exist alone. In addition, the character " / " in this article generally indicates that the previous and subsequent associated objects are in an "or" relationship.

[0021] The term "at least one" in this article is merely a description of the association relationship between associated objects, indicating that three relationships may exist. For example, at least one of A and B can mean: A exists alone, A and B exist at the same time, and B exists alone.

[0022] It should also be noted that, in this document, relational terms such as first and second are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "include," "comprises," or any other variations thereof are intended to cover non-exclusive inclusion.

[0023] Example 1 A tumor target molecular identification method based on ecDNA clonal characteristics, including ecDNA identification, ecDNA clonal trajectory identification, and diffuse gene identification, including: S101: Collect paired normal skin tissue samples, precancerous actinic keratosis tissue samples, and cutaneous squamous cell carcinoma tissue samples from three patients and perform single-cell ATAC sequencing on them; S102, preprocessing the single-cell data; performing cell annotation on the preprocessed single-cell sequencing data; S103, identifying cells containing ecDNA in single-cell sequencing data; identifying characteristic genes of ecDNA; and constructing an ecDNA characteristic matrix; S104. Perform dimensionality reduction and clustering on the feature matrix; identify clonal central cells; analyze and visualize ecDNA cloning trajectories; and select diffuse genes for trajectory analysis.

[0024] The detailed operations for single-cell data preprocessing and cell annotation in S102 are as follows: S201, Cell Ranger ATAC data upstream processing flow; Freshly sequenced scATAC-seq data were processed using the 10× Genomics Cell Ranger ATAC pipeline. This pipeline includes aligning FASTQ files to a reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin regions (peaks), and outputting files such as fragments. The resulting files, such as fragments.tsv.gz and peaks.bed, were used for downstream analysis.

[0025] S202, construct Seurat object; The SignacR package was used to process the Cell Ranger output. The peaks.bed file for each sample was read and converted into a GenomicRanges object. The singlecell.csv metadata file containing cell filter information was also loaded. High-quality cells with a passed_filters count greater than 500 were retained to ensure data integrity. Peak data from all samples were then merged into a unified combined.peaks collection. The reduce function was used to eliminate overlapping regions and filter out atypical peak regions with lengths <20 bp or >10,000 bp. For each sample, a FragmentObject object was created to store cell-specific fragment information. From this, a feature matrix (FeatureMatrix) was generated, with rows representing peaks and columns representing cells. These matrices were encapsulated as ChromatinAssay objects and integrated into Seurat objects, along with sample metadata. Finally, the Seurat objects for all samples were merged to form a combined dataset. The dataset field was added to identify the sample source and supplemented with human genome annotation information.

[0026] S203, perform quality control on Seurat objects; First, the nucleosome signal was calculated to assess the characteristics of DNA open regions, and cells with excessively high nucleosome signals (nucleosome_signal>4) were filtered out to remove nucleosome contamination. Next, the TSS enrichment score was calculated to assess the degree of chromatin openness near the transcription start site, retaining high-quality cells with TSS.enrichment>3. The ratio of reads within peaks (pct_reads_in_peaks>15%) was analyzed to ensure that cells had sufficient open chromatin signal. Finally, the blacklist ratio (blacklist_ratio<0.05) was calculated to exclude low-quality cells with excessively high reads in difficult-to-map regions such as repetitive sequences. This stringent quality control process was validated using various visualization methods, including DensityScatter, TSSPlot, FragmentHistogram, and VlnPlot.

[0027] Specifically, preprocessing of single-cell data includes the following requirements: Number of open chromatin fragments (nCount_ATAC): retain 1,000 to 30,000 valid fragments per cell and exclude low-quality cells (too few fragments) or potential multiple cells (too many fragments).

[0028] Peak region read ratio (pct_reads_in_peaks): >15%, ensuring that most sequencing fragments are located in the open chromatin region and reducing background noise.

[0029] Blacklist region ratio (blacklist_ratio): <5%, filtering out the impact of highly repetitive sequences or technical interference regions.

[0030] Nucleosome signal (nucleosome_signal): <4, excluding abnormal fragment length caused by undigested nucleosomes.

[0031] TSS enrichment score (TSS.enrichment): >3, ensuring that the data are significantly enriched in the gene promoter region, reflecting the reliability of the experimental quality.

[0032] S204, cell type annotation; The preprocessed peaks data were normalized using the signacR package. Highly expressed peaks with a total count greater than 10 were selected, and data dimensionality reduction was achieved using singular value decomposition (SVD). A K-nearest neighbor graph was constructed based on the dimensionality reduction results, and the Louvain algorithm was used to cluster cells into multiple characteristic clusters. A gene activity matrix was then calculated by quantifying the number of open chromatin fragments in gene bodies and promoter regions, and then normalized using LogNormalize. Finally, the differential activity profiles of marker genes within various cell types within the tumor microenvironment (such as immune cells, stromal cells, and tumor cells) were analyzed within each cluster. Each cluster was annotated as a specific cell type based on its characteristic expression profiles.

[0033] The process of S103 identifying cells containing ecDNA and the characteristic genes of ecDNA in single-cell data is as follows: S301, screening of amplified fragments and copy number variation analysis; Identifying cells harboring ecDNA in single-cell ATAC-seq data first requires detecting aberrantly amplified regions across the genome. Because the human genome is diploid, ecDNA typically manifests as a significant local increase in copy number (CNV>2). CNV analysis was performed on the sequencing data for each single cell using the AMULET tool to identify genomic regions with abnormally high copy number. To ensure data reliability, stringent filtering criteria were established, including alignment quality (MAPQ>5), reasonable insert size (≤300 bp), full alignment of paired-end reads to the reference genome (hg38), and removal of potential confounding factors such as PCR duplicates. Based on these criteria, a binary matrix was constructed, with rows representing detected amplification sites and columns representing distinct single cells, thereby systematically characterizing the distribution pattern of amplified fragments within the cell population.

[0034] S302, detection and characterization of genomic breakpoints; The formation of ecDNA is often accompanied by genomic rearrangement events, so detecting breakpoints is a key step in identifying ecDNA. This study identifies potential breakpoints by analyzing the characteristics of aligned reads, primarily in the following three situations: (1) large inserts (>500 bp, this threshold can be adjusted based on research needs), indicating the presence of large-scale genomic sequence rearrangements; (2) cross-chromosomal alignments, where paired reads align to different chromosomes; and (3) co-directional alignments, where two reads align in the same direction, which may indicate the presence of inversions or circular DNA structures. All detected breakpoints retain cellular barcode information, ensuring that subsequent analysis can be performed at single-cell resolution.

[0035] S303, identification and structural reconstruction of ecDNA; To accurately identify ecDNA, the present invention performs upstream and downstream extensions (default 2000bp, adjustable) of amplified regions with CNV > 2 to cover the possible breakpoint junctions. Subsequently, the bedtools tool is used to calculate the genomic overlap between the extended amplified fragments and the breakpoints. A key verification step is to confirm whether these amplified fragments and breakpoints originate from the same cell, which is achieved by matching their cellular barcodes. Once verification is successful, the circular structure of the ecDNA can be reconstructed and its complete genomic composition can be inferred.

[0036] S304, characteristic gene annotation and functional analysis; ecDNA often carries genes closely related to the occurrence and development of cancer (such as oncogenes or drug-resistance genes). To analyze the biological functions of ecDNA, the present invention uses the ChIPseeker tool to perform gene annotation on the expanded amplified fragments, identifying important features such as coding genes and regulatory elements contained therein. In addition, the ecDNAscope tool outputs detailed structural information of ecDNA in each cell (presented in the form of an ordered list of fragments) and sorts and visualizes ecDNA based on parameters such as copy number, number of breakpoints, and frequency of carrying cells. These data not only reveal the distribution pattern of ecDNA on the genome (such as whether it is enriched in specific oncogene sites), but also provide important clues for studying the role of ecDNA in tumor heterogeneity and evolution.

[0037] Constructing the ecDNA feature matrix specifically includes the following steps: The ecDNA gene list was generated in the previous step. The ecDNA feature matrix is ​​a numerical matrix with rows representing genes and columns representing samples. It is generated by iterating through the target gene list, extracting the row corresponding to each gene from the copy number data frame of multiple samples, and calculating the sum of the copy number of each gene across all samples.

[0038] The dimensionality reduction clustering of the feature matrix specifically includes the following steps: Seurat objects were created using the ecDNA matrix based on the Seurat R package. The data were normalized using the NormalizeData function to eliminate the effects of sequencing depth or library size and ensure comparable expression of ecDNA signature genes across different cells. Signature genes with high coefficients of variation were calculated using the FindVariableFeatures function. After normalizing the data using the ScaleData function, principal component analysis (PCA) was performed on the signature matrix, and the top 20 principal components were selected as the basis for subsequent analyses. For cluster analysis, the FindClusters function was used with a resolution parameter of 2 to perform cell clustering.

[0039] The identification of the cloning center cell specifically comprises the following steps: S01, identifying ecDNA marker genes based on thresholds; S02. The cell with the highest expression level of marker genes in each cluster is defined as the clone center cell.

[0040] The thresholds are: the gene must be expressed in at least 1% of cells; the fold change value after log2 conversion is greater than or equal to 0.1.

[0041] As a further description of the above scheme, the detailed steps for ecDNA clone trajectory analysis and visualization are as follows: For each cluster, the sum of the expression values ​​of each ecDNA feature in all cells in the cluster was calculated; Pearson correlation coefficients between different clusters were calculated based on the sum; Two clusters with a Pearson correlation coefficient greater than 0.3 were defined as continuous clonal trajectories; Use UMAP to visualize clone trajectories; The detailed steps for selecting diffuse genes for trajectory analysis are as follows: After obtaining the ecDNA clone track, the clonal relationships between all clusters are collected. After selecting a specific clone track (e.g., cluster3 -> cluster7 -> cluster12), the ecDNA characteristic genes with significant copy number differences between clusters are calculated as diffuse genes. Then, a heat map of diffuse genes between clusters is drawn to visualize the changes in diffuse genes between clusters. A significant difference is defined as a cumulative change in the copy number of a certain ecDNA feature between different clusters in the cloning trajectory, with a step size of n, where n can be customized.

[0042] Example 2 This embodiment is carried out on the basis of the above-mentioned embodiment 1, and the similarities with the above-mentioned embodiment 1 are not repeated here.

[0043] Constructing the ecDNA feature matrix specifically includes the following steps: (1) A total of 8713 ecDNA signatures were identified in 8453 cells containing ecDNA; (2) Based on the CNV values ​​of 8713 identified features, behavioral signatures were constructed and listed as the ecDNA feature matrix of cells.

[0044] The dimensionality reduction clustering of the feature matrix specifically includes the following steps: performing principal component analysis (PCA) on the feature matrix, and selecting the first 20 principal components to perform dimensionality reduction clustering analysis on the matrix.

[0045] Identifying the central cell of the clone specifically includes the following steps: (1) Identify ecDNA marker genes based on a threshold value; the threshold value is: the gene must be expressed in at least 1% of cells; the fold change value after log2 conversion is greater than or equal to 0.1.

[0046] (2) The cell with the highest expression of marker genes in each cluster is defined as the clone center cell.

[0047] ecDNA clone trajectory analysis and visualization specifically includes the following steps: (1) For each cluster, the sum of the expression values ​​of each ecDNA feature in all cells in the cluster is calculated; (2) Calculate the Pearson correlation coefficient between different clusters based on the sum; (3) Two clusters with a Pearson correlation coefficient greater than 0.3 were defined as continuous clone trajectories; (4) Use UMAP to visualize the clone trajectory. The results are shown in the attached figure. Figure 1 shown.

[0048] A. The figure shows the ecDNA clone trajectories of 27 clusters identified by dimensionality reduction clustering based on single-cell ATAC-seq. B. The figure shows the correlation between ecDNA clones in different clusters. The Pearson correlation coefficients between clusters range from 0.3 to 0.7, indicating that ecDNA clones can be tracked well between cell clusters. C. The figure shows the ecDNA genes shared between different clusters. Among them, genes such as NCL, CCNYL1, METTL21A, and MIR033 are well shared between different clusters. D. The figure shows the correlation between ecDNA and cell cycle markers, indicating that the amount of ecDNA is positively correlated with the cell's proliferation capacity. The stronger the cell's proliferation capacity, the greater the amount and characteristics of ecDNA. E. The figure shows diffuse gene expression along ecDNA clone tracks. Clusters 13, 7, and 3 share a large amount of ecDNA. The transmission of these shared ecDNA tracks illustrates the clonal relationships between cells. F. The figure shows the distribution of all sample types containing ecDNA cells, ranging from paracancerous lesions to precancerous lesions to skin squamous cell carcinomas, illustrating the evolution of ecDNA clonal trajectories throughout the disease process. G. The figure shows key ecDNA signatures during skin cancer progression along clonal trajectories. Different ecDNAs drive the progression from peritumoral cancer to precancerous lesions and finally to squamous cell carcinoma. For example, CTDSP1 and MCM7 ecDNAs specifically drive the development of squamous cell carcinoma, while BLC2L11 and GTDC1 ecDNAs drive the development of precancerous lesions.

[0049] The foregoing description is merely a preferred embodiment of the present invention and does not limit the scope of protection of the present invention. Those skilled in the art will readily appreciate that the present invention is susceptible to various modifications and variations. Any variation, modification, replacement, integration, or parameter change to these embodiments, which is within the spirit and principles of the present invention and which achieves the same functionality through conventional substitutions, without departing from the principles and spirit of the present invention, falls within the scope of protection of the present invention.

Claims

1. A method for identifying single-cell ecDNA clone trajectories, characterized in that: The method operates as follows: S101. Collect paired normal tissue samples, precancerous lesion tissue samples, and cancer tissue samples from multiple patients for single-cell sequencing; S102, performing cell annotation on the pre-processed single-cell sequencing data; S103, identifying cells containing ecDNA and characteristic genes of ecDNA in the single-cell sequencing data, constructing an ecDNA feature matrix, performing dimensionality reduction clustering on the feature matrix, and identifying the central cells of the clones; S104: ecDNA clone trajectory analysis and visualization, and selected diffuse genes for trajectory analysis to identify single-cell ecDNA clone trajectories.

2. The method for identifying single-cell ecDNA clone trajectories according to claim 1, characterized in that: The detailed operations for single-cell data preprocessing and cell annotation in S102 are as follows: S201, Cell Ranger ATAC data upstream processing flow; scATAC-seq data were processed using the 10× Genomics Cell Ranger ATAC pipeline, which includes aligning FASTQ files to the reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin peaks, and outputting fragments files. The resulting fragments.tsv.gz and peaks.bed files were used for downstream analysis. S202, construct Seurat object; The SignacR package was used to process the output of Cell Ranger. The peaks.bed file of each sample was read and converted into a GenomicRanges object. The singlecell.csv metadata file containing cell filtering information was loaded. High-quality cells with passed_filters > 500 were retained to ensure data reliability. The peak data from all samples are then merged into a unified combined.peaks set. The reduce function is used to eliminate overlapping regions and filter out atypical peak regions with lengths <20bp or >10,000bp. For each sample, a FragmentObject is created to store cell-specific fragment information. Based on this fragment information, a count matrix FeatureMatrix is ​​generated, with rows representing peaks and columns representing cells. This matrix is ​​encapsulated as a ChromatinAssay object and integrated into a Seurat object, along with sample metadata. Finally, an integrated dataset is formed by merging the Seurat objects of all samples, adding a dataset field to identify the sample source and supplementing with human genome annotation information. S203, perform quality control on Seurat objects; First, the characteristics of DNA open regions were evaluated by calculating the nucleosome signal, and cells with a nucleosome signal of nucleosome_signal>4 were filtered out to remove nucleosome contamination. Second, the TSS enrichment score was calculated to assess the degree of chromatin openness near the transcription start site, and high-quality cells with TSS.enrichment>3 were retained. At the same time, the peak fragment ratio (pct_reads_in_peaks>15%) was analyzed to ensure that the cells had sufficient chromatin open region signals. Finally, the blacklist region ratio was calculated to exclude low-quality cells with a blackbox region ratio of blacklist_ratio<0.

05. The quality control steps were verified by multiple visualization methods including DensityScatter, TSSPlot, FragmentHistogram and VlnPlot; S204, cell type annotation; The preprocessed peaks data were normalized using the signacR package for TF-IDF, and highly expressed characteristic peaks with a total count > 10 were selected. Data dimensionality reduction was achieved using singular value decomposition (SVD). A K-nearest neighbor graph was constructed based on the dimensionality reduction results, and the Louvain algorithm was used to cluster cells into multiple characteristic clusters. The gene activity matrix was then calculated by quantifying the number of open chromatin fragments in gene bodies and promoter regions, and then normalized using LogNormalize. Finally, the differential activity characteristics of marker genes of various cell types in the tumor microenvironment were analyzed in each cluster, and each cluster was annotated as a specific cell type based on the characteristic expression profile.

3. The method for identifying single-cell ecDNA clone trajectories according to claim 1, characterized in that: The process of S103 identifying cells containing ecDNA and the characteristic genes of ecDNA in single-cell data is as follows: S301, screening of amplified fragments and copy number variation analysis; To identify cells carrying ecDNA in single-cell ATAC-seq data, it is first necessary to detect aberrantly amplified regions on the genome; The AMULET tool was used to perform CNV analysis on the sequencing data of each single cell to identify genomic intervals with abnormally high copy numbers and CNVs > 5. To ensure data reliability, stringent filtering criteria were set, including alignment quality MAPQ > 5, reasonable insert size ≤ 300 bp, complete alignment of paired-end reads to the reference genome hg38, and removal of potential interference factors such as PCR duplicates. Based on these criteria, a binary matrix was constructed, in which rows represent detected amplification sites and columns represent different single cells, thereby systematically characterizing the distribution pattern of amplified fragments in the cell population. S302, detection and characterization of genomic breakpoints; The formation of ecDNA is often accompanied by genomic rearrangement events, so detecting breakpoints is a key step in identifying ecDNA. Potential breakpoints are identified by analyzing the characteristics of aligned reads, mainly in the following three situations: (1) Large inserts >500 bp, indicating that there may be large-scale rearrangements in the genomic sequence; (2) Cross-chromosome alignment, that is, paired reads are aligned to different chromosomes; (3) Co-directional alignment, that is, the alignment direction of the two reads is the same, which suggests the possible presence of inversions or circular DNA structures; all detected breakpoints retain the cell barcode information, ensuring that subsequent analysis can be performed at single-cell resolution; S303, identification and structural reconstruction of ecDNA; To accurately identify ecDNA, amplified regions with CNV > 2 were extended 1000-3000 bp upstream and downstream to cover the possible breakpoint junction regions. The bedtools tool was then used to calculate the genomic overlap between the extended amplified fragment and the breakpoint. Cell barcode matching was then used to verify that the amplified fragment and the breakpoint were from the same cell. Upon successful verification, the circular structure of the ecDNA was reconstructed and the complete genomic composition was inferred. S304, characteristic gene annotation and functional analysis; To understand the biological functions of ecDNA, the ChIPseeker tool was used to annotate the expanded amplified fragments and identify the important features of the coding genes and regulatory elements contained therein. In addition, the ecDNAscope tool outputs detailed structural information of the ecDNA in each cell, presenting it in the form of an ordered list of fragments, and sorting and visualizing the ecDNA according to copy number, number of breakpoints, and carrier cell frequency parameters.

4. The method for identifying single-cell ecDNA clone trajectories according to claim 3, characterized in that: Constructing the ecDNA feature matrix specifically includes the following steps: After obtaining the ecDNA-related gene list, the ecDNA feature matrix is ​​constructed. The ecDNA feature matrix is ​​a numerical matrix with rows representing genes and columns representing samples. It is obtained by traversing the target gene list, extracting the rows corresponding to each gene from the copy number data frame of multiple samples in turn, and calculating the total copy number of the gene in all samples.

5. The method for identifying single-cell ecDNA clone trajectories according to claim 3, characterized in that: The dimensionality reduction clustering of the feature matrix specifically includes the following steps: Based on the Seurat R package, the ecDNA matrix was used to create a Seurat object. The NormalizeData function was used to normalize the data to eliminate the effects of sequencing depth or library size, making the expression of ecDNA signature genes between different cells comparable. The FindVariableFeatures function was used to calculate signature genes with high coefficients of variation. After normalizing the data using the ScaleData function, principal component analysis was performed on the signature matrix, and the first 20 principal components were selected as the basis for subsequent analysis. In the cluster analysis stage, the FindClusters function was used to set the resolution parameter to 2 for cell clustering.

6. The method for identifying single-cell ecDNA clone trajectories according to claim 3, characterized in that: The identification of the cloning center cell specifically comprises the following steps: S01, identifying marker genes of ecDNA based on threshold values; S02. The cell with the highest expression level of marker genes in each cluster is defined as the clone center cell; The threshold values ​​are: the gene must be expressed in at least 1% of cells; and the fold change value after log2 conversion must be greater than or equal to 0.

1.

7. The method for identifying single-cell ecDNA clone trajectories according to claim 1, characterized in that: The detailed steps for ecDNA clone trajectory analysis and visualization are as follows: For each cluster, the sum of the expression values ​​of each ecDNA feature in all cells in the cluster was calculated; Pearson correlation coefficients between different clusters were calculated based on the sum; Two clusters with a Pearson correlation coefficient greater than 0.3 were defined as continuous clonal tracks; The nonlinear dimensionality reduction algorithm UMAP was used to visualize the clonal evolution process from clusters with low ecDNA copy numbers to clusters with higher counts; The detailed steps for selecting diffuse genes for trajectory analysis are as follows: After obtaining the ecDNA cloning trajectory, the clonal relationships between all clusters are collected. After selecting a specific cloning trajectory, ecDNA characteristic genes with significant differences in copy number between clusters are calculated as diffuse genes. Then, a heat map of diffuse genes between clusters is drawn to visualize the changes in diffuse genes between clusters. The significant difference is defined as the cumulative change in copy number of a certain ecDNA feature between different clusters in the cloning trajectory with a step size of n, where n can be customized.

Images