A method for identifying single-cell ecDNA clonal trajectories

Abstract

The application discloses a single-cell ecDNA clone track recognition method, and recognizes potential tumor markers in single-cell ATAC data based on ecDNA clone characteristics. The application can sensitively recognize tumor cells and perform early diagnosis by recognizing cancer markers through ecDNA; 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 application helps to reveal the evolution track of tumors and the diversity of cell populations thereof; by studying the clone track of ecDNA, the application can reveal the gene structure changes of tumor cell populations in different stages, thereby providing potential markers for early screening and monitoring of cancers.

Description

Technical Field

[0001] This invention relates to the field of tumor marker identification, and more particularly to a method for identifying the cloning trajectory of single-cell ecDNA. Background Technology

[0002] Recent studies have shown that ecDNA typically contains specific oncogenes, drug resistance genes, immune regulatory genes, and other regulatory elements, making it a crucial factor in cancer pathogenesis and poor prognosis across various 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 copy number of tumor genes and enhanced chromatin accessibility of ecDNA, ATAC-seq is a powerful tool for investigating the functional roles of ecDNA in enhancer hijacking, trans-functional regulation, tumor heterogeneity, and uneven distribution.

[0003] A key characteristic of ecDNA is its high copy number and enriched promoter regions, which can significantly enhance the expression levels of related oncogenes. For example, oncogenes such as KLF7, MYC, and EGFR are commonly found on ecDNA, and their expression is closely related to tumor invasiveness 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 explore the structure and function of ecDNA in depth, revealing its complex role in cancer. In particular, ecDNA has been found to regulate gene expression through "enhancer hijacking" in various cancer types, 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 regarding its biological functions, formation mechanisms, and how to maintain and eliminate it 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, the expression regulatory patterns of specific oncogenes and their roles in tumor heterogeneity can be revealed. Summary of the Invention

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

[0006] To address the aforementioned technical problems, this application provides the following technical solution: a method for identifying the cloning trajectory of single-cell ecDNA, the method comprising the following steps:

[0007] S101. Collect paired normal tissue samples, precancerous lesion tissue samples, and cancer tissue samples from multiple patients for single-cell sequencing;

[0008] S102. Perform cell annotation on the preprocessed single-cell sequencing data;

[0009] S103. Identify cells containing ecDNA and characteristic genes of ecDNA in single-cell sequencing data, construct ecDNA feature matrix, perform dimensionality reduction and clustering on the feature matrix, and identify cloning center cells.

[0010] S104, ecDNA cloning trajectory analysis and visualization, and selected diffuse genes for trajectory analysis to identify single-cell ecDNA cloning trajectories.

[0011] As a further description of the above scheme: the detailed operations for preprocessing single-cell data and cell annotation in S102 are as follows:

[0012] S201, Cell Ranger (cell sequencing) ATAC data upstream processing workflow;

[0013] The scATAC-seq data from sequencing was processed using the 10× Genomics Cell Ranger ATAC workflow. The workflow included aligning the FASTQ file to a reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin regions (peaks), and outputting fragments files. The resulting fragments.tsv.gz and peaks.bed files were used for downstream analysis.

[0014] S202, Construct the Seurat object;

[0015] The output of Cell Ranger was processed using the SignacR package. This involved reading the peaks.bed file for each sample and converting it into GenomicRanges objects, while simultaneously loading the singlecell.csv metadata file containing cell filtering information. High-quality cells with passed_filters > 500 were retained to ensure data reliability. The 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 > 10000 bp. For each sample, a FragmentObject was created to store cell-specific fragment information. Based on this fragment information, a FeatureMatrix was generated, 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 from all samples were merged to form an integrated dataset. A dataset field was added to identify the sample origin, and human genome annotation information was supplemented.

[0016] S203. Perform quality control on Seurat objects;

[0017] First, the characteristics of open DNA regions were assessed by calculating nucleosome signals, and cells with nucleosome signal > 4 were filtered out to remove nucleosome contamination. Second, TSS enrichment scores were calculated to assess the degree of chromatin openness near transcription start sites, retaining high-quality cells with TSS enrichment > 3. Simultaneously, the proportion of intra-peak fragments (pct_reads_in_peaks) > 15% was analyzed to ensure that cells had sufficient open chromatin region signals. Finally, the proportion of blacklisted regions was calculated, and low-quality cells with a blacklist ratio < 0.05 were excluded. This series of rigorous quality control steps was validated using multiple visualization methods, including DensityScatter, TSS Plot, Fragment Histogram, and VlnPlot.

[0018] S204, Cell type annotation;

[0019] The preprocessed peaks data were normalized using TF-IDF using the signacR package, and high-expression characteristic peaks with a total count >10 were selected. Dimensionality reduction was achieved through singular value decomposition (SVD). Based on the dimensionality reduction results, a K-nearest neighbor graph was constructed, and cell clustering was performed using the Louvain algorithm to divide the cells into multiple characteristic clusters. Subsequently, the gene activity matrix was calculated by quantifying the number of open chromatin fragments in the gene body and promoter regions, and LogNormalize normalization was performed. Finally, by analyzing the differential activity characteristics of marker genes of various cell types in each cluster within the tumor microenvironment, each cluster was annotated as a specific cell type based on the characteristic expression profile.

[0020] As a further description of the above scheme: The process by which S103 identifies cells containing ecDNA in single-cell data and the characteristic genes of ecDNA is as follows:

[0021] S301. Screening and copy number variation analysis of amplified fragments;

[0022] To identify cells carrying ecDNA in single-cell ATAC-seq data, it is first necessary to detect anomalously amplified regions on the genome; the AMULET tool is used to perform CNV analysis on the sequencing data of each single cell to identify genomic regions with abnormally high copy numbers of CNV > 5.

[0023] To ensure data reliability, strict 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 from PCR duplication. Based on these criteria, a binary matrix was constructed, where rows represent detected amplification sites and columns represent different single cells, thereby systematically characterizing the distribution pattern of amplified fragments in the cell population.

[0024] S302. Detection and feature analysis of genomic breakpoints;

[0025] The formation of ecDNA is usually accompanied by genomic rearrangement events, therefore, detecting breakpoints is a key step in identifying ecDNA. Potential breakpoints are identified by analyzing the characteristics of aligned reads, mainly including the following three scenarios:

[0026] (1) Large insertion fragments > 500 bp indicate that there may be large-scale rearrangement of the genome sequence;

[0027] (2) Cross-chromosome alignment, that is, the paired reads are aligned to different chromosomes;

[0028] (3) Same direction alignment, that is, the alignment direction of the two reads is the same, which suggests that there may be inverted or circular DNA structures; all detected breakpoints retain cell barcode information, ensuring that subsequent analysis can be performed at single-cell resolution;

[0029] Identification and structural reconstruction of S303 and ecDNA;

[0030] To accurately identify ecDNA, amplified regions with CNV (copy number variation) > 2 are extended upstream and downstream, with a default extension range of 1000bp-3000bp, to cover possible breakpoint connection regions. Subsequently, the bedtools tool is used to calculate the genomic overlap between the amplified fragment and the breakpoint. Then, cell barcode matching is used to verify whether the amplified fragment and the breakpoint originate from the same cell. After successful verification, the circular structure of the ecDNA is reconstructed, and the complete genomic composition is inferred.

[0031] S304, Feature gene annotation and functional analysis;

[0032] ecDNA typically carries genes closely related to cancer development and progression. To analyze the biological functions of ecDNA, the ChIPseeker tool is used to annotate the expanded fragments, identifying the coding genes and regulatory elements contained within. In addition, the ecDNAscope tool outputs detailed structural information of ecDNA in each cell, presented as an ordered list of fragments, and sorts and visualizes the ecDNA based on parameters such as copy number, number of breakpoints, and frequency of carrying cells.

[0033] As a further description of the above scheme: Constructing the ecDNA feature matrix specifically includes the following steps:

[0034] After obtaining the list of genes related to ecDNA, an ecDNA feature matrix is ​​constructed. The ecDNA feature matrix is ​​a numerical matrix where rows represent genes and columns represent samples. It is obtained by traversing the list of target genes, 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 the gene in all samples.

[0035] As a further description of the above scheme: the dimensionality reduction and clustering of the feature matrix specifically includes the following steps:

[0036] Based on the Seurat R package, Seurat objects are created using ecDNA matrices. The NormalizeData function is used to standardize the data to eliminate the influence of sequencing depth or library size, making the expression of ecDNA characteristic genes between different cells comparable. The FindVariableFeatures function is used to calculate characteristic genes with high coefficients of variation. After normalizing the data using the ScaleData function, principal component analysis is performed on the feature matrix, and the first 20 principal components are selected as the basis for subsequent analysis. In the cluster analysis stage, the FindClusters function is used to set the resolution parameter to 2 for cell clustering.

[0037] As a further description of the above scheme: the identification of clonal center cells specifically includes the following steps:

[0038] S01. Identify marker genes for ecDNA based on threshold values;

[0039] S02. The cell with the highest expression level of marker genes in each cluster is defined as the cloning center cell;

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

[0041] As a further description of the above scheme, the detailed operation of ecDNA cloning trajectory analysis and visualization is as follows:

[0042] For each cluster, the sum of the expression values ​​of each ecDNA feature across all cells in that cluster is calculated.

[0043] Pearson correlation coefficients between different clusters are calculated based on summation.

[0044] Two clusters with a Pearson correlation coefficient greater than 0.3 are defined as consecutive clonal trajectories;

[0045] The clonal evolution process from clusters with low ecDNA copy numbers to clusters with high copy numbers was visualized using UMAP (a non-linear dimensionality reduction algorithm).

[0046] The detailed steps for selecting diffuse genes for trajectory analysis are as follows:

[0047] After obtaining the ecDNA cloning trajectory, the cloning relationships between all clusters are collected. After selecting a specific cloning trajectory, the ecDNA feature genes with significant differences in copy number between clusters are calculated as diffuse genes. Then, a heatmap of diffuse genes between clusters is plotted to visualize the changes of diffuse genes between clusters. 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.

[0048] The features or innovations of this invention are as follows:

[0049] 1) Highly sensitive identification at the single-cell level. By integrating single-cell ATAC-seq data with ecDNA clonal feature analysis, tumor-specific ecDNA and its oncogenes (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 method can capture subtle heterogeneity of ecDNA in tumor cell populations, providing molecular markers for early screening.

[0050] 2) Dynamic analysis of tumor clonal evolution trajectories. Based on the copy number, structural heterogeneity, and clonal distribution characteristics of ecDNA, the diversity of tumor cell subpopulations and their dynamic clonal evolution patterns can be revealed. By tracking the spatial distribution and temporal evolution of different ecDNA clones, the molecular mechanisms of key processes such as tumor drug resistance and metastasis can be elucidated, providing new strategies for targeted therapy.

[0051] Compared with existing technologies, this invention offers the following advantages: By identifying cancer biomarkers through ecDNA, this invention can sensitively identify tumor cells and enable early diagnosis. The copy number, structure, and distribution of ecDNA often exhibit significant heterogeneity within tumor cells. Different tumor cell populations may carry different types of ecDNA clones, and this invention helps to reveal the evolutionary trajectory of tumors and the diversity of their cell populations. Furthermore, by studying the clonal trajectory of ecDNA, this invention can reveal changes in the gene structure of tumor cell populations at different stages, thereby providing potential biomarkers for early cancer screening and monitoring. Attached Figure Description

[0052] Figure 1 The following sections illustrate ecDNA cloning trajectories: a) shows the ecDNA cloning trajectories of 27 clusters identified by single-cell ATAC-seq dimensionality reduction clustering; b) shows the correlation between ecDNA clones across different clusters; c) shows the shared ecDNA genes among different clusters; d) shows the correlation between ecDNA and cell cycle markers; e) shows the diffuse gene expression along the ecDNA cloning trajectories; f) shows the sample type distribution of all cells containing ecDNA; and g) shows key ecDNA features in the evolution of skin cancer along the cloning trajectories. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. In the following description, specific details such as specific configurations and components are provided merely to help fully understand the embodiments of this application. Therefore, those skilled in the art should understand that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of this application. In addition, for clarity and brevity, descriptions of known functions and structures are omitted in the embodiments.

[0054] It should be understood that the phrase "an embodiment" or "this embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "an embodiment" or "this embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments.

[0055] Furthermore, reference numerals and / or letters may be repeated in different examples within this application. Such repetition is for the purpose of simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or settings discussed.

[0056] In this article, the term "and / or" is merely a description of the relationship between related 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 simultaneously. The term " / and" in this article describes another type of relationship between related objects, 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 related objects before and after it are in an "or" relationship.

[0057] In this article, the term "at least one" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, "at least one of A and B" can mean: A exists alone, A and B exist simultaneously, or B exists alone.

[0058] 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, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion.

[0059] Example 1

[0060] A method for identifying tumor target molecules based on ecDNA clonal characteristics, including ecDNA identification, ecDNA clonal trajectory identification, and diffuse gene identification, comprising:

[0061] S101. Collect 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.

[0062] S102. Preprocess the single-cell data; perform cell annotation on the preprocessed single-cell sequencing data;

[0063] S103. Identify cells containing ecDNA in single-cell sequencing data; identify characteristic genes of ecDNA; construct an ecDNA feature matrix;

[0064] S104. Perform dimensionality reduction and clustering on the feature matrix; identify clonal center cells; analyze and visualize ecDNA clonal trajectory; select diffuse genes for trajectory analysis.

[0065] The detailed steps for preprocessing and annotating single-cell data in S102 are as follows:

[0066] S201, Cell Ranger ATAC data upstream processing flow;

[0067] The scATAC-seq data immediately after sequencing was processed using the 10×Genomics Cell Ranger ATAC workflow. This workflow includes aligning the FASTQ file to a reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin regions (peaks), and outputting fragments files. The final generated fragments.tsv.gz and peaks.bed files are used for downstream analysis.

[0068] S202, Construct the Seurat object;

[0069] The output of Cell Ranger was processed using the SignacR package. This involved reading the peaks.bed file for each sample and converting it into GenomicRanges objects, while simultaneously loading the singlecell.csv metadata file containing cell filtering information. High-quality cells with more than 500 passed filters were retained to ensure data reliability. The 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 shorter than 20 bp or longer than 10000 bp. For each sample, a FragmentObject was created to store cell-specific fragment information. Based on this, a 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 from all samples were merged to form an integrated dataset, with a dataset field added to identify the sample origin and supplemented with human genome annotation information.

[0070] S203. Perform quality control on Seurat objects;

[0071] First, the characteristics of open DNA regions were assessed by calculating nucleosome signal, and cells with excessively high nucleosome signals (nucleosome_signal>4) were filtered out to remove nucleosome contamination. Second, TSS enrichment scores were calculated to assess the degree of chromatin openness near transcription start sites, retaining high-quality cells with TSS.enrichment>3. Simultaneously, the proportion of fragments within peaks (pct_reads_in_peaks>15%) was analyzed to ensure that cells had sufficient open chromatin region signals. Finally, the proportion of blacklisted regions (blacklist_ratio<0.05) was calculated to exclude low-quality cells with excessively high proportions of fragments in difficult-to-map regions such as repetitive sequences. This series of rigorous quality control steps was validated using multiple visualization methods, including DensityScatter, TSSSPlot, FragmentHistogram, and VlnPlot.

[0072] Specifically, preprocessing single-cell data includes the following requirements:

[0073] Cellular chromatin open fragment number (nCount_ATAC): retain 1,000 to 30,000 valid fragments per cell, excluding low-quality cells (too few fragments) or potentially multicellular cells (too many fragments).

[0074] Peak reads percentage (pct_reads_in_peaks): >15%, ensuring that most sequencing fragments are located in open chromatin regions, reducing background noise.

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

[0076] Nucleosome signal: <4, excluding fragment length abnormalities caused by undigested nucleosomes.

[0077] 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.

[0078] S204, Cell type annotation;

[0079] The preprocessed peaks data were normalized using TF-IDF using the signacR package, and high-expression characteristic peaks with a total count >10 were selected. Dimensionality reduction was achieved through singular value decomposition (SVD). A K-nearest neighbor graph was constructed based on the dimensionality reduction results, and the Louvain algorithm was used for cell clustering, dividing the cells into multiple characteristic clusters. Subsequently, the gene activity matrix was calculated by quantifying the number of open chromatin fragments in gene bodies and promoter regions, and then normalized using LogNormalize. Finally, by analyzing the differential activity characteristics of marker genes of various cell types (such as immune cells, stromal cells, and tumor cells) in each cluster within the tumor microenvironment, each cluster was annotated as a specific cell type based on its characteristic expression profile.

[0080] The S103 identification process for cells containing ecDNA and the characteristic genes of ecDNA in single-cell data is as follows:

[0081] S301. Screening and copy number variation analysis of amplified fragments;

[0082] Identifying cells carrying ecDNA in single-cell ATAC-seq data first requires detecting anomalous amplification regions on the genome. Since the human genome is diploid, ecDNA typically exhibits significantly increased local copy numbers (CNV>2). The AMULET tool was used to perform CNV analysis on the sequencing data of each single cell to identify genomic regions with abnormally high copy numbers. To ensure data reliability, strict filtering criteria were established, 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 duplications. Based on these criteria, a binary matrix was constructed, where rows represent detected amplification sites and columns represent different single cells, thus systematically characterizing the distribution pattern of amplified fragments in the cell population.

[0083] S302. Detection and feature analysis of genomic breakpoints;

[0084] The formation of ecDNA is often accompanied by genomic rearrangement events, therefore, detecting breakpoints is a crucial step in identifying ecDNA. This invention identifies potential breakpoints by analyzing the characteristics of aligned reads, mainly including the following three situations: (1) large insertion fragments (>500 bp, this threshold can be adjusted according to research needs), indicating that there may be large-scale rearrangement of the genomic sequence; (2) cross-chromosome alignment, that is, paired reads are aligned to different chromosomes; (3) same-direction alignment, that is, the alignment directions of two reads are the same, which may indicate the presence of inversions or circular DNA structures. All detected breakpoints retain cell barcode information, ensuring that subsequent analysis can be performed at single-cell resolution.

[0085] Identification and structural reconstruction of S303 and ecDNA;

[0086] To accurately identify ecDNA, this invention extends the amplified region with CNV>2 upstream and downstream (default extension 2000bp, this parameter is adjustable) to cover possible breakpoint junction regions. Subsequently, the genomic overlap between the extended amplified fragments and breakpoints is calculated using the bedtools tool. A crucial validation step is confirming whether these amplified fragments and breakpoints originate from the same cell, which is achieved by matching their cell barcodes. Once validation is successful, the circular structure of the ecDNA can be reconstructed, and its complete genomic composition can be inferred.

[0087] S304, Feature gene annotation and functional analysis;

[0088] ecDNA typically carries genes closely related to cancer development and progression (such as oncogenes or drug resistance genes). To elucidate the biological functions of ecDNA, this invention uses the ChIPseeker tool to annotate expanded fragments, identifying important features such as coding genes and regulatory elements. Furthermore, the ecDNAscope tool outputs detailed structural information of ecDNA in each cell (presented as 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 patterns of ecDNA across the genome (e.g., whether it is enriched at specific oncogene sites) but also provide important clues for studying the role of ecDNA in tumor heterogeneity and evolution.

[0089] Constructing the ecDNA feature matrix specifically includes the following steps:

[0090] The ecDNA-related gene list has already been obtained in the previous process. The ecDNA feature matrix is ​​a numerical matrix where rows represent genes and columns represent samples. It is obtained by traversing the target gene list, extracting the row corresponding to each gene from the copy number data frames of multiple samples, and calculating the sum of the copy numbers of that gene in all samples.

[0091] The specific steps for dimensionality reduction and clustering of the feature matrix include:

[0092] Seurat objects are created using the ecDNA matrix based on the Seurat R package. The NormalizeData function is used to standardize the data, eliminating the influence of sequencing depth or library size and ensuring comparability of ecDNA characteristic gene expression across different cells. The FindVariableFeatures function is used to calculate characteristic genes with high coefficients of variation. After normalizing the data using the ScaleData function, principal component analysis (PCA) is performed on the feature matrix, and the top 20 principal components are selected as the basis for subsequent analyses. In the clustering analysis phase, the FindClusters function is used to cluster cells with a resolution parameter set to 2.

[0093] The identification of clonal center cells specifically includes the following steps:

[0094] S01. Identify ecDNA marker genes based on thresholds;

[0095] S02. The cell with the highest expression level of marker genes in each cluster is defined as the cloning center cell.

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

[0097] As a further description of the above scheme, the detailed operation of ecDNA cloning trajectory analysis and visualization is as follows:

[0098] For each cluster, the sum of the expression values ​​of each ecDNA feature across all cells in that cluster is calculated.

[0099] Pearson correlation coefficients between different clusters are calculated based on summation.

[0100] Two clusters with a Pearson correlation coefficient greater than 0.3 are defined as consecutive clonal trajectories;

[0101] Use UMAP to visualize the clone trajectory;

[0102] The detailed steps for selecting diffuse genes for trajectory analysis are as follows:

[0103] After obtaining the ecDNA clonal trajectory, clonal relationships between all clusters were collected. After selecting a specific clonal trajectory (e.g., cluster 3 -> cluster 7 -> cluster 12), ecDNA characteristic genes with significant differences in copy number between clusters were calculated as diffuse genes. A heatmap of diffuse genes across clusters was then plotted to visualize the changes in diffuse genes between clusters.

[0104] Significant differences are defined as the cumulative change in the copy number of a certain ecDNA feature across different clusters in the cloning trajectory with a step size of n, where n can be customized.

[0105] Example 2

[0106] This embodiment is based on the above embodiment 1, and the similarities with the above embodiment 1 will not be repeated.

[0107] Constructing the ecDNA feature matrix specifically includes the following steps:

[0108] (1) A total of 8,713 ecDNA features were identified in 8,453 cells containing ecDNA;

[0109] (2) Based on the CNV values ​​of the 8713 identified features, behavioral features were constructed and listed as the ecDNA feature matrix of cells.

[0110] 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.

[0111] Identifying clonal center cells specifically includes the following steps:

[0112] (1) Identify ecDNA marker genes based on thresholds; the thresholds are: the gene must be expressed in at least 1% of cells; the fold change value after log2 transformation is greater than or equal to 0.1.

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

[0114] ecDNA cloning trajectory analysis and visualization specifically includes the following steps:

[0115] (1) For each cluster, calculate the sum of the expression values ​​of each ecDNA feature across all cells in that cluster;

[0116] (2) Calculate the Pearson correlation coefficient between different clusters based on summation;

[0117] (3) Two clusters with a Pearson correlation coefficient greater than 0.3 are defined as consecutive clonal trajectories;

[0118] (4) The cloning trajectory was visualized using UMAP, and the results are shown in the appendix. Figure 1 As shown.

[0119] A. The figure shows the ecDNA cloning trajectories of 27 clusters identified based on single-cell ATAC-seq dimensionality reduction clustering;

[0120] B. The figure shows the correlation of ecDNA clones between different clusters. The Pearson correlation coefficients between various clusters range from 0.3 to 0.7, indicating that ecDNA has good clonal tracking ability among cell clusters.

[0121] C. The figure shows the ecDNA genes shared between different clusters, among which NCL, CCNYL1, METTL21A, and MIR033 have good sharing among different clusters;

[0122] 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 proliferative capacity. The stronger the cell's proliferative capacity, the greater the amount and characteristics of ecDNA.

[0123] E. The figure shows the diffuse gene expression along the ecDNA clonal trajectory, from cluster 13 to 7 and then to cluster 3, sharing a large amount of ecDNA. The trajectory transfer process of these shared ecDNAs illustrates the clonal relationship between cells.

[0124] F. The figure shows the distribution of all sample types containing ecDNA, including adjacent normal, precancerous lesions and squamous cell carcinoma of the skin, demonstrating the evolution of the clonal trajectory of ecDNA throughout the entire pathogenesis process;

[0125] G. The figure illustrates key ecDNA features in the evolution of skin cancer along the clonal trajectory. Different ecDNAs drive the process from adjacent normal tissue to precancerous lesions and then 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.

[0126] The above description is merely a preferred embodiment of the present invention and does not limit the scope of protection of the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any changes, modifications, substitutions, integrations, and parameter alterations to these embodiments within the spirit and principles of the present invention, achieved through conventional substitutions or by achieving the same function without departing from the principles and spirit of the present invention, fall within the scope of protection of the present invention.

Claims

1. A method for identifying the cloning trajectory of single-cell ecDNA, characterized in that, The method operates as follows: S101, collect paired normal tissue samples, precancerous lesion tissue samples and cancer tissue samples of a plurality of patients for single cell sequencing; S102, cell annotation is performed on the pretreated single cell sequencing data; S103, cells containing ecDNA and ecDNA characteristic genes in the single cell sequencing data are identified, an ecDNA characteristic matrix is constructed, dimensionality reduction clustering is performed on the characteristic matrix, and clonal center cells are identified; The process of identifying cells containing ecDNA and ecDNA characteristic genes in single cell data is as follows: S301, screening of amplified fragments and copy number variation analysis; In order to identify cells carrying ecDNA in single cell ATAC-seq data, it is necessary to first detect abnormal amplification regions on the genome; CNV analysis is performed on the sequencing data of each single cell using the AMULET tool, and genomic intervals with abnormally high copy number CNV>5 are identified; In order to ensure data reliability, strict filtering criteria are set, including alignment quality MAPQ>5, reasonable insert size ≤300 bp, double-end read alignment to the reference genome hg38, and removal of potential interference factors such as PCR duplication; based on the criteria, a binary matrix is constructed, in which the rows represent the detected amplification sites and the columns represent different single cells, thereby systematically characterizing the distribution pattern of amplified fragments in the cell population; S302, detection and feature analysis of genomic breakpoints; The formation of ecDNA is 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, including the following three cases: (1) Large insert size >500 bp, indicating large-scale rearrangement of genomic sequences; (2) Cross-chromosome alignment, i.e., paired reads align to different chromosomes; (3) Same direction alignment, i.e., the alignment directions of the two reads are the same, which suggests the presence of an inversion or circular DNA structure; all detected breakpoints retain cell barcode information, ensuring that subsequent analysis can be performed at single cell resolution; S303, identification and structure reconstruction of ecDNA; In order to accurately identify ecDNA, the amplification region with CNV>2 is extended 1000-3000 bp upstream and downstream to cover the breakpoint connection region; then the bedtools tool is used to calculate the genomic overlap between the extended amplified fragments and the breakpoints; then it is verified whether the amplified fragments and breakpoints come from the same cell by matching the cell barcode; after successful verification, the circular structure of ecDNA is reconstructed, and the complete genomic composition is inferred; S304, characteristic gene annotation and functional analysis; To analyze the biological function of ecDNA, ChIPseeker tool was used to annotate the extended amplified fragments, identify the important features of coding genes and regulatory elements contained therein; In addition, ecDNAscope tool outputs detailed structural information of ecDNA in each cell in the form of an ordered fragment list, and sorts and visualizes ecDNA according to copy number, breakpoint number and cell frequency parameters; S104, ecDNA clone trajectory analysis and visualization, and selected diffuse genes for trajectory analysis to identify single-cell ecDNA clone trajectories; The detailed operation of selecting diffuse genes for trajectory analysis is as follows: After obtaining the ecDNA clone trajectory, the cloning relationship between all clusters cluster is collected; After selecting a specific clone trajectory, the ecDNA feature genes with significant differences in copy number between clusters cluster are calculated as diffuse genes; Then draw the heatmap of diffuse genes between clusters cluster, visualize the change of diffuse genes between clusters cluster, and the significant difference is defined as the cumulative change of the copy number of a certain ecDNA feature between different clusters cluster in the clone trajectory with n as the step, and the value of n is defined.

2. The method of claim 1, wherein, The detailed operation of preprocessing and cell annotation of single-cell data in S102 is as follows: S201, upstream processing flow of Cell Ranger ATAC data; The scATAC-seq data just sequenced off-line is processed using the Cell Ranger ATAC process of 10x Genomics, which includes aligning FASTQ files to the reference genome, removing PCR duplicates, identifying valid cells, generating open chromatin region peaks, and outputting fragments files; The final fragments.tsv.gz, peaks.bed files are used for downstream analysis; S202, construct Seurat object; The output results of Cell Ranger are processed using the SignacR package, by reading the peaks.bed file of each sample and converting it to a GenomicRanges object, while loading the singlecell.csv metadata file containing cell filtering information, and retaining high-quality cells with passed_filters>500 to ensure data reliability; Subsequently, the peak data of all samples are merged into a unified combined.peaks set, the reduce function is used to eliminate overlapping regions, and non-typical peak regions with a length <20 bp or >10000 bp are filtered out; for each sample, a FragmentObject fragment object is created to store cell-specific fragment information, a count matrix FeatureMatrix is generated based on the fragment information, which represents peaks in rows and cells in columns, the matrix is encapsulated as a ChromatinAssay object and integrated into a Seurat object, while sample metadata is attached; finally, an integrated dataset is formed by merging Seurat objects of all samples, dataset field is added to identify sample sources, and human genome annotation information is supplemented; S203, quality control of the Seurat object; Firstly, the DNA open region feature is evaluated by calculating the nucleosome signal, and cells with nucleosome_signal>4 are filtered out to remove nucleosome contamination; secondly, the TSS enrichment score is calculated to evaluate the chromatin openness near the transcription start site, and high-quality cells with TSS.enrichment>3 are retained; at the same time, the proportion of peak reads pct_reads_in_peaks>15% is analyzed to ensure that the cells have enough chromatin open region signal; finally, the blacklist ratio is calculated to exclude low-quality cells with blacklist_ratio<0.05; The quality control step is verified by DensityScatter, TSSPlot, FragmentHistogram and VlnPlot visualization methods; S204, cell type annotation; The peaks data after preprocessing is processed by TF-IDF standardization using the signacR package, and high-expression feature peaks with total count>10 are selected, and data dimensionality reduction is realized by singular value decomposition SVD; based on the dimensionality reduction result, a K-nearest neighbor graph is constructed, and Louvain algorithm is used for cell clustering, and the cells are divided into multiple characteristic clusters; then the gene activity matrix is calculated by quantifying the number of open chromatin fragments in the gene body and promoter region, and LogNormalize standardization processing is performed; finally, the differential activity characteristics of marker genes of various cells in the tumor microenvironment in each cluster are analyzed, and each cluster is annotated as a specific cell type according to the characteristic expression profile.

3. The method of claim 1, wherein, Constructing the ecDNA feature matrix specifically includes the following steps: After obtaining the ecDNA-related gene list, the ecDNA feature matrix is constructed, which is a numerical matrix with rows representing genes and columns representing samples; it is obtained by traversing the target gene list, extracting the corresponding row of 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.

4. The method of claim 1, wherein, The dimensionality reduction and clustering of the feature matrix specifically include the following steps: Seurat objects were created based on the ecDNA matrix using the Seurat R package. The data was normalized to eliminate the effects of sequencing depth or library size using the NormalizeData function, so that the expression of ecDNA signature genes between different cells is comparable. The FindVariableFeatures function was used to calculate the signature genes with high coefficient of variation. After the data was normalized using the ScaleData function, principal component analysis was performed on the feature matrix, and the first 20 principal components were selected as the basis for subsequent analysis. In the clustering analysis stage, the FindClusters function was used to set the resolution parameter to 2 for cell clustering.

5. The method of claim 1, wherein, The identification of the clonal center cell specifically includes the following steps: S01, identifying the marker genes of ecDNA according to the threshold value; S02, defining the cell with the highest expression of marker genes in each cluster as the clonal center cell; The threshold value is that the gene must be expressed in at least 1% of the cells, and the log2 transformed fold change value is greater than or equal to 0.

1.

6. The method of claim 1, wherein, The detailed operation of ecDNA clonal trajectory analysis and visualization is as follows: For each cluster, the expression values of each ecDNA feature in all cells in the cluster are added up; Based on the addition calculation, the Pearson correlation coefficient between different clusters is calculated; Two clusters with a Pearson correlation coefficient greater than 0.3 are defined as a continuous clonal trajectory; The nonlinear dimension reduction algorithm UMAP is used to visualize the clonal evolution process from the cluster with lower ecDNA copy number to the cluster with higher count.

Images