ScATAC-seq-based ecDNA structure prediction method, method for identifying cells carrying ecDNA and medium

Abstract

The invention provides an ecDNA structure prediction method based on scATAC-seq, a method for identifying cells carrying ecDNA and a medium, and relates to the technical field of biological information. According to the ecDNA structure prediction method based on scATAC-seq, inconsistent readings of ecDNA breakpoint sources are extracted through comparison with a reference genome, enriched peak regions are called in sequence and converted into a graph, a breakpoint graph is combined with a Bayesian model, and accurate prediction of structural information such as ecDNA breakpoints, ecDNA regions and connection modes is achieved. The invention also provides a method for identifying cells carrying the ecDNA, and the method is used for accurately acquiring the cells carrying the ecDNA through standardization treatment, formation of a gene activity matrix and non-single-peak test on the ecDNA structure based on the scATAC-seq. According to the present invention, the ecDNA structure can be predicted, the cells carrying the ecDNA can be accurately identified, and the method can be used for the accurate analysis and research of the tumor heterogeneity, evolution and drug resistance process of the ecDNA.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and more specifically, to a method for predicting ecDNA structure based on scATAC-seq and a method and medium for identifying cells carrying ecDNA. Background Art

[0002] Extrachromosomal circular DNAs (eccDNAs) are a type of circular DNA molecule that originates from genomic DNA but exists independently of linear chromosomes. Among these, large extrachromosomal circular DNA (ranging from hundreds of kilobases to several megabases) is commonly found in cancer cells. Studies have shown that ecDNA plays an important role in tumor heterogeneity, tumor evolution, and drug resistance, and has the potential to become a target for tumor therapy.

[0003] The earliest discovery of extrachromosomal DNA can be traced back to 1965, when Cox et al. observed small double-chromosome structures without visible centromeres in tumor cells and named them "double minutes" (DMs). As research deepened, the term "DMs" was gradually replaced by the more inclusive "ecDNA". For decades, researchers have widely used experimental techniques such as fluorescence in situ hybridization (FISH) and electrophoresis to label, separate and analyze ecDNA. However, experimental methods are usually applicable to a small number of samples or cell lines, and there are certain limitations in the study of patient-derived ecDNA, mainly reflected in cumbersome operations and high costs.

[0004] In recent years, the rapid development of genomic sequencing technology has made it possible to identify ecDNA in large patient cohorts. Due to the high accessibility of ecDNA, identifying potential ecDNA based on assay for transposase-accessible chromatin using sequencing (ATAC-seq) data has become a viable approach. However, conventional ATAC-seq samples often contain millions of cells, which results in the averaging of ecDNA signals, making it difficult to accurately reflect whether individual cells carry ecDNA.

[0005] The emergence of single-cell sequencing assay for transposase-accessible chromatin (scATAC-seq) has enabled accessibility signals to be traced back to each specific cell, providing new opportunities for ecDNA research at the single-cell level. However, existing technologies still face challenges in determining whether a single cell carries ecDNA. Therefore, there is an urgent need for a method that can predict the presence of ecDNA at the single-cell level to more accurately analyze the role of ecDNA in tumor heterogeneity, evolution, and drug resistance. Summary of the Invention

[0006] The purpose of the present invention is to solve the above-mentioned technical problems and provide a method for predicting ecDNA structure based on scATAC-seq. This method can predict ecDNA structure based on scATAC-seq. Another purpose of the present invention is to provide a method for further identifying cells carrying ecDNA based on this prediction method; another purpose of the present invention is to provide a medium for the above two methods.

[0007] To achieve the above object, the present invention adopts the following technical solutions:

[0008] The present invention provides a method for predicting ecDNA structure based on scATAC-seq, comprising the following steps:

[0009] S1. Align and quantify the sample scATAC-seq sequencing results with the reference genome to generate an alignment file. Extract discordant reads derived from ecDNA breakpoints from the alignment file.

[0010] S2. calling enriched peak regions from discordant reads, wherein the peak regions are peak regions that appear in both the forward and reverse directions;

[0011] S3. Convert the peak regions into a graph, where the nodes of the graph are the enriched peak regions; the discordant edges of the graph are the connections between the peak regions enriched by the discordant reads;

[0012] S4. Estimate the locations of ecDNA breakpoints based on the enriched regions of discordant reads using a Bayesian model. Convert the graph to a breakpoint graph based on the ecDNA breakpoints, where the nodes of the breakpoint graph represent base-accurate breakpoints. Add continuous edges to the discordant edges of the breakpoint graph, where the continuous edges are defined as the sequencing depth and coverage of all reads between the breakpoints.

[0013] S5. Save the reasonable nodes in the graph, where the reasonable nodes need at least one inconsistent edge and one continuous edge; and output the complete structure of ecDNA.

[0014] Furthermore, the inconsistent reads in step S1 are reads that are aligned to the reference genome by both paired-end sequencing and meet one of the following conditions:

[0015] (1) Aligned to the same chromosome and in the same direction;

[0016] (2) aligned to the same chromosome, and the insert length was greater than 5 kb;

[0017] (3) Aligned to different chromosomes.

[0018] Furthermore, the threshold k of the enriched peak region in step S2 is calculated as the data points greater than 95%, and the calculation formula is shown in Formula I:

[0019] Formula I:

[0020] In Formula I, X is the sequencing depth observed in a certain region; λ is the average sequencing depth of the whole genome, that is, the expectation of the Poisson distribution; and k is a threshold to be determined so that the probability P of a depth less than k in the context of the whole genome is greater than or equal to 95%.

[0021] The present invention provides a method for identifying cells carrying ecDNA, comprising the following steps:

[0022] P1. Align and quantify the sample's scATAC-seq sequencing results with the reference genome, obtaining a result file that contains the alignment information, including the source cell, chromosome, chromosome start and end points, cell label, the number of reads supporting the region, and whether the cell is qualified.

[0023] P2. Based on the comparison result file obtained in step P1, perform ecDNA structure prediction, and obtain the coordinates of ecDNA in the genome according to the reference genome; the ecDNA structure prediction method is the ecDNA structure prediction method according to any one of claims 1 to 3;

[0024] P3. Select the ecDNA predicted and identified in step P2 and normalize the alignment result file in P1 to obtain the accessibility signal of the ecDNA region of the normalized cell;

[0025] P4. Generate a gene activity matrix for each cell based on the ecDNA selected in step P3, and extract the matrix subset containing genes in the ecDNA region;

[0026] P5. Based on the ecDNA selected in step P3, extract the inconsistent reads of the ecDNA breakpoints in the alignment result file obtained in step P1, and map the inconsistent reads to cell labels;

[0027] P6. Perform a non-singlet test on the accessibility signal of the ecDNA region of the cell after normalization in step P3 to obtain a non-singlet test result for each ecDNA region;

[0028] P7. Based on the non-singlet test results for each ecDNA region obtained in step P6, fit the accessibility signals of the normalized ecDNA regions of the cells obtained in step P3 and the gene activity matrix of each cell obtained in step P4 to distinguish potential ecDNA+ and ecDNA- cell subpopulations.

[0029] P8. Integrate the results obtained in steps P5 and P7 to obtain cells carrying ecDNA in the sample; the integration is to perform Fisher's exact test on the components of cells with inconsistent readings at ecDNA breakpoints obtained in step P5 and ecDNA+ cells in step P7 to verify the statistical significance of the results.

[0030] Furthermore, the step P1 includes the following steps:

[0031] P11. Using cell labeling and next-generation sequencing, perform 10× scATAC-seq sequencing on the Illumina platform to generate fastq files.

[0032] P12. Use cellranger-atac to align and quantify the fastq file obtained in P11 with the human hg38 genome to obtain the aligned bam file, fragment.tsv.gz file, and singlecell.csv file. The bam file marks the cell from which each read originated. The fragment.tsv.gz file records the chromosome, chromosome start, chromosome end, cell label, and the number of reads supporting the region. The singlecell.csv file records the cell label and whether the cell is qualified.

[0033] Preferably, the standardization processing method in step P3 is: according to the singlecell.csv file in step P1, qualified cells are retained; using the fragment.tsv.gz file obtained in step P1 as input, the CPM values ​​(ecCPM values) of these qualified cells in the ecDNA region are calculated and standardized to obtain the accessibility signal of the ecDNA region of the standardized cells.

[0034] More preferably, the normalization process is as shown in Formula II:

[0035] Formula II:

[0036] Normalized ecCPMi is the normalized ecCPM value of the i-th cell; the number of reads in the ecDNA region is the number of reads in the ecDNA region in the cell; the number of reads in all alignments is the total number of reads in the cell.

[0037] Furthermore, the step P4 includes the following steps:

[0038] P41. Use the fragment.tsv.gz file from step P1 as input to ArchR software and use the addGeneScoreMatrix function to construct a matrix of all gene activity for each cell.

[0039] P42. Extract the genes contained in the ecDNA selected in step P3 within the range of genomic coordinates;

[0040] P43. From the complete gene activity matrix of each cell obtained in P41, the activity matrix of genes in the ecDNA region of qualified cells is retained to obtain a matrix subset;

[0041] The genes in the ecDNA region are extracted through step P42.

[0042] Furthermore, in the non-unimodal test process described in step P6, if the P value of the accessibility signal of the ecDNA region of the normalized cell is less than 0.05, the alternative hypothesis is accepted and the data obeys a multimodal distribution; if the P value is greater than or equal to 0.05, the null hypothesis is accepted and the data obeys a unimodal distribution.

[0043] Furthermore, if the results of step P6 show that the data follows a multimodal distribution, a Gaussian mixture model is used to decompose the accessibility signal of the ecDNA region of the normalized cells obtained in step P3 into a component number of 2. Otherwise, k-means clustering is used to cluster the ecDNA gene activity matrix of step P4 into a component number of 2.

[0044] The formula of the Gaussian mixture model is shown in Formula III:

[0045] Formula III:

[0046] Among them, μ1 and μ2 are the means of the two Gaussian models, σ1 and σ2 are the estimated standard deviations of the two Gaussian models, and w is the estimated weight parameter of the first Gaussian model.

[0047] The present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program implements any of the above methods when executed by a processor.

[0048] The beneficial effects of the present invention are as follows:

[0049] (1) The present invention provides an ecDNA structure prediction method based on scATAC-seq, which extracts inconsistent reads derived from ecDNA breakpoints by comparing with a reference genome, and sequentially calls the enriched peak regions, which are converted into graphs and breakpoint graphs combined with a Bayesian model to achieve accurate prediction of structural information such as ecDNA breakpoints, ecDNA regions, and connection modes.

[0050] (2) The present invention also provides a method for identifying cells carrying ecDNA. This method utilizes standardized processing and the formation of a gene activity matrix for the above-mentioned scATAC-seq-based ecDNA structure, and accurately obtains cells carrying ecDNA through non-single-peak testing.

[0051] (3) The present invention can predict the structure of ecDNA and accurately identify cells carrying ecDNA. It can be used to accurately analyze and study the role of ecDNA in tumor heterogeneity, evolution, and drug resistance, and has important medical detection application prospects for non-therapeutic purposes. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 This is a flowchart of the ecDNA structure prediction method based on scATAC-seq in Example 1 of the present invention.

[0053] Figure 2 This is a brief flow chart of the method for identifying cells carrying ecDNA according to Example 2 of the present invention.

[0054] Figure 3 This is a detailed flow chart of the method for identifying cells carrying ecDNA according to Example 2 of the present invention.

[0055] Figure 4 This is a schematic diagram of part of the results of executing step P1 in Example 4 of the present invention.

[0056] Figure 5A This is a schematic diagram of the ecEGFR results of executing step P2 in Example 4 of the present invention.

[0057] Figure 5B This is a schematic diagram of the ecMDM4 result of executing step P2 in Example 4 of the present invention.

[0058] Figure 5C This is a schematic diagram of the ecPDGFRA results of executing step P2 in Example 4 of the present invention.

[0059] Figure 6 This is a schematic diagram of part of the results of executing step P5 in Example 4 of the present invention.

[0060] Figure 7 This is a schematic diagram of part of the results of executing steps P6 and P7 in Example 4 of the present invention.

[0061] Figure 8 This is a schematic diagram of part of the results of executing step P8 in Example 4 of the present invention. DETAILED DESCRIPTION

[0062] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0063] Example 1

[0064] This example provides a method for predicting ecDNA structure based on scATAC-seq. Figure 1 As shown, the following steps are included:

[0065] S1. Align and quantify the sample scATAC-seq sequencing results with the reference genome to generate an alignment file. Extract discordant reads derived from ecDNA breakpoints from the alignment file.

[0066] Alternatively, the sample scATAC-seq sequencing results may be processed by aligning and quantifying the sample scATAC-seq sequencing results with a reference genome, generating a comparison result file that includes the alignment information, including the source cell, chromosome, chromosome start, chromosome end, cell label, the number of sequencing results supporting the region, and whether the cell is qualified.

[0067] As a preferred embodiment, the sample scATAC-seq sequencing results are obtained by the following method: using cell labeling and second-generation sequencing, the sample to be analyzed is sequenced at 10×scATAC-seq on the Illumina platform to obtain a fastq file; the fastq file is aligned and quantified with the human hg38 genome using cellranger-atac to obtain a bam file after alignment, in which each label marks the cell from which the read originated; and also includes a fragment..tsv.gz file containing information about chromosomes, chromosome start, chromosome end, cell labels, and the number of reads supporting the region; as well as a singlecell.csv file recording cell labels and whether the cell is qualified, and a bam index file.

[0068] S2. calling enriched peak regions from discordant reads, wherein the peak regions are peak regions that appear in both the forward and reverse directions;

[0069] As an optional solution, the discordant reads in step S2 need to be paired and meet one of the following conditions:

[0070] (1) Aligned to the same chromosome and in the same direction;

[0071] (2) aligned to the same chromosome, and the insert length was greater than 5 kb;

[0072] (3) Aligned to different chromosomes.

[0073] As a preferred solution, in step S2, the enriched peak region is called from the inconsistent reads, and the threshold k of the enriched peak region is calculated to be greater than 95% of the data points, and the calculation formula is shown in Formula I:

[0074] Formula I:

[0075] In Formula I, X is the observed sequencing depth in a certain region, which is assumed to follow a Poisson distribution; λ is the average sequencing depth of the entire genome, that is, the expectation (mean) of the Poisson distribution; k is a threshold we want to determine so that the probability P of a depth less than k in the context of the entire genome is greater than or equal to 95%; i is an intermediate variable in the cumulative sum, representing different values ​​in the Poisson distribution (from 0 to k-1).

[0076] S3. Convert the peak regions into a graph, where the nodes of the graph are the enriched peak regions; the discordant edges of the graph are the connections between the peak regions enriched by the discordant reads;

[0077] S4. Estimate the locations of ecDNA breakpoints based on the enriched regions of discordant reads using a Bayesian model. Convert the graph to a breakpoint graph based on the ecDNA breakpoints, where the nodes of the breakpoint graph represent base-accurate breakpoints. Add continuous edges to the discordant edges of the breakpoint graph, where the continuous edges are defined as the sequencing depth and coverage of all reads between the breakpoints.

[0078] S5. Save the reasonable nodes in the graph, where the reasonable nodes need at least one inconsistent edge and one continuous edge; and output the ecDNA connection method.

[0079] Example 2

[0080] This embodiment provides a method for identifying cells carrying ecDNA. Figure 2 and attached Figure 3 As shown, the following steps are included:

[0081] P1. Align and quantify the sample scATAC-seq sequencing results with the reference genome, generating a result file that records the source cell, chromosome, chromosome start and end points, cell label, the number of sequencing results supporting the region, and whether the cell is qualified.

[0082] As an optional solution, step P1 includes the following steps:

[0083] P11. Using cell labeling and next-generation sequencing, perform 10× scATAC-seq sequencing on the Illumina platform to generate fastq files.

[0084] P12. Use cellranger-atac to align and quantify the fastq file with the human hg38 genome, generating a bam file. Each tag identifies the cell from which the read originated. The file also includes a fragment.tsv.gz file, which contains information about chromosomes, chromosome start and end points, cell labels, and the number of reads supporting that region. The singlecell.csv file, which records cell labels and cell eligibility, and the bam index file, bai.

[0085] P2. Perform ecDNA structure prediction on the alignment results in step P1, and obtain the ecDNA coordinates in the genome based on the reference genome; the ecDNA structure prediction is performed as described in Example 1 using the scATAC-seq-based ecDNA structure prediction method;

[0086] P3. Select the ecDNA predicted and identified in step P2 and normalize the accessibility signals of the cells in S1 at the ecDNA region;

[0087] As an optional solution, the standardization processing method of step P3 is: retain qualified cells based on the singlecell.csv file in step P1; use the fragment.tsv.gz file of P1 as input to calculate the CPM values ​​of these qualified cells standardized in the ecDNA region.

[0088] As a preferred solution, the CPM value of the cell in the ecDNA region is normalized using Normalized ecCPM i It is defined as the number of raw reads in the cellular ecDNA region divided by the total number of reads, multiplied by 1e 6 Then take the logarithm and standardize it as shown in Formula II:

[0089] Formula II:

[0090] In formula II, ecCPMi is the ecCPM value of the i-th cell; the number of reads in the ecDNA region is the number of reads in the ecDNA region in the cell; the number of reads in all alignments is the total number of reads in the cell.

[0091] P4. Generate a gene activity matrix for each cell based on the ecDNA selected in step P3, and extract the matrix subset containing genes in the ecDNA region;

[0092] As an optional solution, step P4 includes the following steps:

[0093] P41. Use the fragment.tsv.gz file from step P1 as input to the ArchR software and use the addGeneScoreMatrix function to construct a matrix of all gene activity for each cell.

[0094] P42. Extract the genes included in the coordinate axis range of the ecDNA selected in step P3;

[0095] P43. From the complete gene activity matrix of each cell obtained in P41, the activity matrix of genes in the ecDNA region of qualified cells is retained to obtain a matrix subset;

[0096] The genes in the ecDNA region are extracted through step P42.

[0097] P5. Based on the alignment result file of step P1 and the ecDNA selected in step P3, extract the inconsistent reads of the ecDNA breakpoints and map these reads to cell labels;

[0098] P6. Perform a non-singlet test on the normalized signals of the ecDNA regions of all cells in step P3;

[0099] As an optional solution, if the P value of the non-unimodal test in step P6 is less than 0.05, the alternative hypothesis is accepted, that is, the data follow a multimodal distribution; if the P value of the non-unimodal test is greater than or equal to 0.05, the null hypothesis is accepted, that is, the data are considered to follow a unimodal distribution.

[0100] P7. Based on the results of step P6, fit the data from steps P3 and P4 to distinguish potential ecDNA+ and ecDNA- cell subpopulations;

[0101] As a preferred solution, if the result of step P6 shows that the data follows a multimodal distribution, a Gaussian mixture model is used to decompose the ecDNA regional signals of all cells normalized in step P3 with a component number of 2; otherwise, k-means clustering is used to fit the ecDNA gene activity matrix of step P4 with a component number of 2;

[0102] The formula of the Gaussian mixture model is shown in Formula III:

[0103] Formula III:

[0104] Where μ1 and μ2 are the means of the two Gaussian models, σ1 and σ2 are the estimated standard deviations of the two Gaussian models, and w is the estimated weight parameter of the first Gaussian model.

[0105] P8. Integrate the results from steps P5 and P7, evaluate the results using a Fisher's exact test, verify the statistical significance of the results, and identify cells carrying extrachromosomal circular DNA. Preferably, a Fisher's exact test is performed on the cells with discordant ecDNA breakpoint reads from step P5 and the ecDNA+ cells from step P7. The Fisher's exact test determines the significance of the ecDNA+ cells. If the significance P value is < 0.05, the ecDNA+ cells are considered ecDNA-carrying cells.

[0106] Example 3

[0107] Based on the above-mentioned ecDNA structure prediction method and method for identifying cells carrying ecDNA, this embodiment provides a computer-readable storage medium for executing the above-mentioned method.

[0108] Example 4

[0109] This example uses the scATAC-seq sequencing data of glioblastoma patient GBM4349 (SRA database accession number: SRR10315838) as input, and the method for predicting cells carrying ecDNA is as follows:

[0110] Step P1: Using cell labeling (DNA barcode) and next-generation sequencing (NGS), specifically, performing 10× scATAC-seq sequencing on the sample to be analyzed using the Illumina platform to obtain a fastq file.

[0111] The specific method is to use cellranger-atac (10×Genomics official alignment software) to align and quantify the fastq file with the hg38 genome (the download address of the 10×genomics hg38 reference genome is https: / / cf.10xgenomics.com / supp / cell-atac / refdata-cellranger-arc-GRCh38-2020-A-2.0.0.tar.gz) (the basic parameters when executing cellranger-atac are as follows: The command "count --id = sampleID --reference = / path / to / reference_dir --fastqs = / path / to / fastq_dir --sample = fastq_prefix" returns the default values ​​for all other parameters. The user is required to provide four parameters: --id is the unique run ID and output folder name; --reference is the path to the folder containing the cellranger-atac reference; --fastqs is the path to the input fastq data; and --sample is the prefix of the selected fastq filenames. This generates an aligned bam file (possorted_bam.bam), where each read is tagged with a barcode (CB) identifying the cell from which it originated. Output from this process also includes the file fragment.tsv.gz (a compressed spreadsheet file similar to the BED format, with five columns after decompression: one for each chromosome, chromosome start, chromosome end, cell label, and the number of reads supporting that region); singlecell.csv (a spreadsheet file containing information such as cell labels and cell eligibility); and a bam index file.

[0112] Figure 4 This is a schematic diagram of some of the comparison and quantification results obtained by executing step P1 in this embodiment, as shown in FIG. Figure 4 As shown, the quality of alignment and quantification is assessed. The quality report obtained by running cellranger-atac on the scATAC-seq sample from GBM4349 is as follows: Figure 4 A is the quality analysis report of cellranger-atac, including information such as the number of cells that passed quality control and the average number of fragments per cell. Figure 4 B is the insert size distribution of the fragments; Figure 4 C is the distribution of fragments from the transcription start site.

[0113] Step P2: Identify the ecDNA in the sample using the possorted_bam.bam file obtained in Step P1 and the blacklisted regions of the hg38 genome. Perform structural prediction of the ecDNA based on the alignment results in Step P1 and obtain the genomic coordinates of the ecDNA based on the hg38 reference genome. The specific method is as follows:

[0114] S1. Align and quantify the sample's scATAC-seq sequencing results (fastq file) with the reference genome using cellranger-atac to generate an alignment file. Extract discordant reads derived from ecDNA breakpoints from the alignment file.

[0115] (1) Obtain the bam file obtained by comparison in step P1 and the blacklist region of hg38 (hg38.blacklist.bed).

[0116] (2) The bam file obtained in step P1 is read jointly by multiprocess (a Python module that supports multithreaded processing) and pysam (a Python module for processing genomic data). Each thread processes a fixed-size genomic window (64 Mb by default).

[0117] (3) Extract the reads that may be the source of ecDNA breakpoints (discordant reads) from the bam file to obtain the discordant reads; the discordant reads are reads that are aligned to the reference genome by both paired-end sequencing and meet one of the following conditions:

[0118] 1) Aligned to the same chromosome and in the same direction;

[0119] 2) Aligned to the same chromosome, with the insert length greater than 5 kb;

[0120] 3) Aligned to different chromosomes.

[0121] (4) Among the inconsistent reads obtained in step (3), the reads that are aligned to the autosomes (chr1-chr22) and sex chromosomes (chrX, chrY) of the hg38 reference genome are retained, and the genomic coordinates of the reads cannot exceed the range of the chromosome; the reads with cell barcode labels are retained; the reads that overlap with the blacklist regions of the hg38 reference genome are excluded, and the filtered inconsistent reads are obtained;

[0122] S2. Call the enriched peak regions from the filtered discordant reads. The peak regions are the peak regions that appear in both the forward and reverse directions. The specific steps are as follows:

[0123] (1) Calculate the threshold for enrichment of inconsistent read regions. Use the average depth of the entire genome as λ and calculate the threshold for data points greater than 95%. The formula is shown in Equation I:

[0124] Formula I:

[0125] In Formula I, X is the observed sequencing depth at any region in the entire genome, assumed to follow a Poisson distribution. λ is the average sequencing depth across the entire genome, i.e., the expected (mean) value of the Poisson distribution. k is a threshold we wish to determine such that, given the entire genome, the probability of a depth less than k is greater than or equal to 95%. i is an intermediate variable in the cumulative sum, representing different values ​​in the Poisson distribution (ranging from 0 to k-1).

[0126] (2) Obtaining the enriched peak region from the discordant reads obtained in step (1); the enrichment threshold is obtained by step (1), and the extension length of the peak region to both sides is defaulted to 500 bp.

[0127] (3) Retain the peak regions that appear in both the forward and reverse directions in the results of step (2).

[0128] S3. Convert the peak regions into a graph, where the nodes of the graph are the enriched peak regions; the discordant edges of the graph are the connections between the peak regions enriched by inconsistent reads; specifically, convert the peak regions of step S2 into a graph, where each enriched peak region is regarded as a node of the graph; and the connections between the peak regions enriched by inconsistent reads are regarded as discordant edges of the graph.

[0129] S4. Estimate the precise locations of ecDNA breakpoints based on the enriched regions of discordant reads using a Bayesian model. Based on the precise breakpoint locations calculated in step S2, convert the graph into a breakpoint map, where nodes represent base-accurate breakpoints. Add continuous edges to the discordant edges of the breakpoint map. A continuous edge is defined as the sequencing depth and coverage of all reads (including discordant and normal reads) between the breakpoint nodes.

[0130] S5. Delete invalid nodes in the breakpoint graph. A valid node must have at least one discordant edge and one continuous edge. Output all valid ecDNA connection methods from the graph.

[0131] The algorithm of steps S1-S5 is named circlehunter2.

[0132] The predicted ecDNAs in the samples are shown in Table 1; Figure 5 shows the structures and regions of some ecDNAs identified in step P2. The ecDNA structures detected from the scATAC-seq data of GBM4349 include ecEGFR ( Figure 5A ), ecMDM4( Figure 5B ) and ecPDGFRA( Figure 5C ).

[0133] Table 1 Predicted ecDNA in samples

[0134]

[0135]

[0136] Step P3: Select ecDNAs numbered ecDNA_4_8, ecDNA_6_1, and ecDNA_7_1 in turn, and standardize the reads of all cells in the ecDNA regions (ecDNA_4_8, ecDNA_6_1, and ecDNA_7_1) (as shown in Formula II). That is, according to the singlecell.csv file in step P1, qualified cells are retained; the fragment.tsv.gz file in step P1 is used as input to calculate the CPM value (count per million) of the standardized ecDNA region in these qualified cells.

[0137] Normalized ecCPM is used to normalize the CPM values ​​of cells in the ecDNA region. i It is defined as the number of raw reads in the cellular ecDNA region divided by the total number of reads, multiplied by 1e 6 Then take the logarithm and calculate the formula as shown in Formula II:

[0138] Formula II:

[0139] Normalized ecCPM i is the normalized ecCPM value of the i-th cell. The reads of the ecDNA region are the reads that overlap with the ecDNA region in the cell; the reads of all alignments are all the reads in the cell.

[0140] Step P4: Using the fragment.tsv.gz file obtained in Step P1 as input, use ArchR software to generate a gene activity matrix for GBM4349 patients and extract the matrix subsets containing genes in the ecDNA regions selected in Step P3: ecDNA_4_8, ecDNA_6_1, and ecDNA_7_1. The specific steps are as follows:

[0141] P41. Use the fragment.tsv.gz file from step P1 as input to the ArchR software and use the addGeneScoreMatrix function to construct a matrix of all gene activity in cells from GBM4349 patients.

[0142] P42. Extract the genes contained within the genomic coordinate range of the ecDNA regions selected in step P3 (ecDNA_4_8, ecDNA_6_1, and ecDNA_7_1). Next, from the complete gene activity matrix for each cell obtained in step P41, obtain the subset of the matrix containing the genes contained within each ecDNA region for each cell.

[0143] P43. From the total gene activity matrix of each cell obtained in P41, the activity matrix of genes in the ecDNA region of qualified cells is retained to obtain a matrix subset; the genes in the ecDNA region are extracted through step P42.

[0144] Table 2 shows the genes carried by some ecDNAs identified in step P2. Among them, ecDNA_6_1 (ecMDM4) carries 25 genes, including MDM4, SOX13, ZBED6, and REN; ecDNA_4_8 (ecEGFR) carries genes including EGFR; and ecDNA_7_1 (ecPDGFRA) carries 44 genes, including PDGFRA.

[0145] Table 2 Genes carried by some ecDNA identified in step P2

[0146]

[0147]

[0148] Step P5: Based on the bam file of the alignment result in step P1 and the ecDNA selected in step P3, use samtools software to extract inconsistent reads from ecDNA breakpoints and map these inconsistent reads to cell labels.

[0149] Take ecDNA_6_1_1 as an example, Figure 6 As shown, the inconsistent reads enriched regions of ecDNA_6_1_1, namely the inconsistent reads in the start_peak (chr1: 203, 462, 575-203, 463, 809) and end_peak (chr1: 205, 164, 836-205, 166, 125) regions, were extracted, where Figure 6The positions of breakpoints representing ecMDM4, discordant reads, and ecDNA breakpoints were marked using red boxes and arrows, respectively, and the cells corresponding to these reads were extracted.

[0150] Step P6: Perform a non-unimodal test on the normalized reads from step P3. If the P value for this non-unimodal test is less than 0.05, the alternative hypothesis is accepted, indicating a multimodal distribution. If the P value for this non-unimodal test is greater than or equal to 0.05, the null hypothesis is accepted, indicating a unimodal distribution. The results show a significant P value < 0.05, supporting the alternative hypothesis, indicating a multimodal distribution of signals in this ecDNA region.

[0151] Step P7: If the data follows a multimodal distribution according to the results of step P6, a Gaussian mixture model (GMM) is used to decompose the ecDNA regional signals of all cells normalized in step P3 into a component of 2. Otherwise, k-means clustering is used to cluster the ecDNA gene activity matrix in step P4 into a component of 2.

[0152] The formula of the Gaussian mixture model is shown in Formula III:

[0153] Formula III:

[0154] Among them, μ1 and μ2 are the means of the two Gaussian models, σ1 and σ2 are the estimated standard deviations of the two Gaussian models, and w is the estimated weight parameter of the first Gaussian model.

[0155] Based on the results of step P6 above, this example uses a Gaussian mixture model to decompose the normalized readings of all cells, setting the number of components to 2, resulting in two components. We consider the component with the higher mean value to be the signal from ecDNA+ cells, and the component with the lower mean value to be the signal from ecDNA- cells. Gaussian decomposition is performed on the normalized signals of ecEGFR, ecPDGFRA, and ecMDM4, respectively. Figure 7 The results show that a Gaussian mixture model was used to fit the normalized reads of scATAC-seq cells from GBM4349 patients, respectively for ecEGFR ( Figure 7 A), ecMDM4( Figure 7 B) and ecPDGFRA ( Figure 7 C) were fitted, and the P values ​​for ecEGFR and ecPDGFRA were 0, 0, and 0.04, respectively ( Figure 7 A. Figure 7 B. Figure 7 C).

[0156] The dip.test was used to test for non-unimodality, with a P value < 0.05 indicating statistical significance. The red curve represents the distribution of reads from cells carrying ecDNA, and the black curve represents the distribution of reads from cells without ecDNA. The red dots indicate cells with discordant reads at the ecDNA breakpoint region.

[0157] Step P8. Perform Fisher's exact test on the cells with inconsistent reads at the ecDNA breakpoints obtained in Step P5 and the ecDNA+ cells obtained in Step P7 to verify the statistical significance of the results and obtain cells carrying extrachromosomal circular DNA;

[0158] Specifically: Fisher's exact test was used to test whether the ecDNA+ cell population obtained in step P7 was enriched with inconsistent reads, and the significance P value of ecDNA+ cells was obtained. At this time, the significance P value was equal to 0 (P < 0.05), accepting the alternative hypothesis, supporting that inconsistent reads tended to the ecDNA+ cell population, and proving the reliability of the classification results.

[0159] Figure 8 The comprehensive results of steps P3-P7 are shown, including the decomposition of Gaussian mixture models, activity heatmaps of genes carried by ecDNA, cells carrying inconsistent reads, and the final prediction results, showing the GBM4349 patient carrying ecEGFR ( Figure 8 A), ecMDM4( Figure 8 B) and ecPDGFRA ( Figure 8 C) Predicted status of cells.

[0160] The present invention discloses a method for identifying cells carrying ecDNA based on scATAC-seq data. By analyzing inconsistent alignment reads related to ecDNA in the scATAC-seq library, candidate ecDNA is identified. Then, based on the high copy characteristics of ecDNA, the accessibility signals of all cells in the ecDNA region are fitted and classified. Finally, statistical tests are performed to ensure the accuracy of the prediction results, thereby accurately identifying cells carrying ecDNA.

[0161] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit the solution. Although the present invention has been described in detail with reference to the preferred embodiments, ordinary technicians in this field should understand and can modify or replace the technical solution of the present invention based on the understanding of this solution without departing from the purpose and scope of the technical solution of the present invention, which should be covered by the scope of the claims of the present invention.

Claims

1. A method for predicting ecDNA structure based on scATAC-seq, characterized in that: The steps include: S1. Align and quantify the sample scATAC-seq sequencing results with the reference genome to generate an alignment file. Extract discordant reads derived from ecDNA breakpoints from the alignment file. S2. calling enriched peak regions from discordant reads, wherein the peak regions are peak regions that appear in both the forward and reverse directions; S3. Convert the peak regions into a graph, where the nodes of the graph are the enriched peak regions; the discordant edges of the graph are the connections between the peak regions enriched by the discordant reads; S4. Estimate the locations of ecDNA breakpoints based on the enriched regions of discordant reads using a Bayesian model. Convert the graph to a breakpoint graph based on the ecDNA breakpoints, where the nodes of the breakpoint graph represent base-accurate breakpoints. Add continuous edges to the discordant edges of the breakpoint graph, where the continuous edges are defined as the sequencing depth and coverage of all reads between the breakpoints. S5. Save the reasonable nodes in the graph, where the reasonable nodes need at least one inconsistent edge and one continuous edge; and output the complete structure of ecDNA.

2. The ecDNA structure prediction method based on scATAC-seq according to claim 1, characterized in that The inconsistent reads in step S1 are reads that are aligned to the reference genome for both paired-end sequencing and meet one of the following conditions: (1) Aligned to the same chromosome and in the same direction; (2) aligned to the same chromosome, and the insert length was greater than 5 kb; (3) Aligned to different chromosomes.

3. The ecDNA structure prediction method based on scATAC-seq according to claim 1, characterized in that The threshold k of the enriched peak region in step S2 is calculated as the data points greater than 95%, and the calculation formula is shown in Formula I: Formula I: In Formula I, X is the sequencing depth observed in a certain region; λ is the average sequencing depth of the whole genome, that is, the expectation of the Poisson distribution; and k is a threshold to be determined so that the probability P of a depth less than k in the context of the whole genome is greater than or equal to 95%.

4. A method for identifying cells carrying ecDNA, characterized in that: The steps include: P1. Align and quantify the sample's scATAC-seq sequencing results with the reference genome, obtaining a result file that contains the alignment information, including the source cell, chromosome, chromosome start and end points, cell label, the number of reads supporting the region, and whether the cell is qualified. P2. Based on the comparison result file obtained in step P1, perform ecDNA structure prediction, and obtain the coordinates of ecDNA in the genome according to the reference genome; the ecDNA structure prediction method is the ecDNA structure prediction method according to any one of claims 1 to 3; P3. Select the ecDNA predicted and identified in step P2 and normalize the alignment result file in P1 to obtain the accessibility signal of the ecDNA region of the normalized cell; P4. Generate a gene activity matrix for each cell based on the ecDNA selected in step P3, and extract the matrix subset containing genes in the ecDNA region; P5. Based on the ecDNA selected in step P3, extract the inconsistent reads of the ecDNA breakpoints in the alignment result file obtained in step P1, and map the inconsistent reads to cell labels; P6. Perform a non-singlet test on the accessibility signal of the ecDNA region of the cell after normalization in step P3 to obtain a non-singlet test result for each ecDNA region; P7. Based on the non-singlet test results for each ecDNA region obtained in step P6, fit the accessibility signals of the normalized ecDNA regions of the cells obtained in step P3 and the gene activity matrix of each cell obtained in step P4 to distinguish potential ecDNA+ and ecDNA- cell subpopulations. P8. Integrate the results from steps P5 and P7 to identify cells carrying ecDNA in the sample; The integration is to perform Fisher's exact test on the cells with inconsistent readings at the ecDNA breakpoints obtained in step P5 and the components of the ecDNA+ cells in step P7 to verify the statistical significance of the results.

5. The method for identifying cells carrying ecDNA according to claim 4, characterized in that: The step P1 comprises the following steps: P11. Using cell labeling and next-generation sequencing, perform 10× scATAC-seq sequencing on the Illumina platform to generate fastq files. P12. Use cellranger-atac to align and quantify the fastq file obtained in P11 with the human hg38 genome to obtain the aligned bam file, fragment.tsv.gz file, and singlecell.csv file. The bam file marks the cell from which each read originated. The fragment.tsv.gz file records the chromosome, chromosome start, chromosome end, cell label, and the number of reads supporting the region. The singlecell.csv file records the cell label and whether the cell is qualified.

6. The method for identifying cells carrying ecDNA according to claim 4, characterized in that: The standardization processing method described in step P3 is: according to the singlecell.csv file in step P1, qualified cells are retained; using the fragment.tsv.gz file obtained in step P1 as input, the CPM values ​​of these qualified cells in the ecDNA region are calculated and standardized to obtain the accessibility signal of the ecDNA region of the standardized cells.

7. The method for identifying cells carrying ecDNA according to claim 4, characterized in that: The step P4 comprises the following steps: P41. Use the fragment.tsv.gz file from step P1 as input to ArchR software and use the addGeneScoreMatrix function to construct a matrix of all gene activity for each cell. P42. Extract the genes contained in the ecDNA selected in step P3 within the range of genomic coordinates; P43. From the complete gene activity matrix of each cell obtained in P41, the activity matrix of genes in the ecDNA region of qualified cells is retained to obtain a matrix subset; The genes in the ecDNA region are extracted through step P42.

8. The method for identifying cells carrying ecDNA according to claim 4, characterized in that: During the non-unimodal test described in step P6, if the P value of the accessibility signal of the ecDNA region of the normalized cells is less than 0.05, the alternative hypothesis is accepted and the data obeys a multimodal distribution; if the P value is greater than or equal to 0.05, the null hypothesis is accepted and the data obeys a unimodal distribution.

9. The method for identifying cells carrying ecDNA according to claim 4, characterized in that: The results of step P6 showed that the data followed a multimodal distribution. A Gaussian mixture model was used to decompose the accessibility signal of the ecDNA region of the normalized cells obtained in step P3 into a component number of 2. Otherwise, k-means clustering was used to cluster the ecDNA gene activity matrix in step P4 into a component number of 2. The formula of the Gaussian mixture model is shown in Formula III: Formula III: Among them, μ1 and μ2 are the means of the two Gaussian models, σ1 and σ2 are the estimated standard deviations of the two Gaussian models, and w is the estimated weight parameter of the first Gaussian model.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 3 or any one of claims 4 to 9 is implemented.

Images