A method for analyzing the epigenetic characteristics of chromosomal instability (CIN) in tumor cells

Abstract

The present invention relates to the technical field of single-cell multi-omics sequencing, and discloses a method for analyzing the epigenetic characteristics of chromosomal instability (CIN) in tumor cells. The steps include: preparing single-cell RNA libraries and DNA libraries for sequencing; respectively preprocessing the DNA and RNA sequencing data, aligning them to the human genome reference sequence, identifying differentially methylated regions, performing karyotype analysis, dividing tumor cells into different subgroups, and quantifying gene expression levels; analyzing the correlations among DNA methylation differences, gene expression differences, and karyotype heterogeneity between different subgroups. The method of the present invention can simultaneously measure DNA methylation, gene expression, and karyotype of single cells, and on this basis, analyze the correlations among DNA methylation, gene expression, and karyotype heterogeneity of tumor cells, which helps to more deeply understand the mechanisms of tumorigenesis and evolution.

Description

Technical Field

[0001] The present invention relates to the technical field of single-cell multi-omics sequencing, and particularly relates to a method for analyzing the epigenetic characteristics of chromosomal instability (CIN) in tumor cells. Background Art

[0002] Chromosomal instability (CIN) is characterized by an increased frequency of chromosomal number gain and loss among cells, manifested as karyotypic heterogeneity. As is well known, CIN plays a key role in cancer initiation and evolution and has been widely documented in various cancers. Despite the strong correlation between CIN and cancer, the precise relationship between them remains a topic of intense debate.

[0003] The limitations of existing single-cell analysis tools have restricted the research and understanding of karyotypic heterogeneity. For example, the HeLa cell line has attracted much attention due to its CIN and has been widely studied. Previous studies have revealed significant heterogeneity in karyotype and chromosomal stability among different HeLa subtypes, different generations, and within the same clone. Although many researchers have documented the chromosomal heterogeneity of HeLa, the underlying mechanisms have not been explored in depth. Notably, in the HeLa-CCL2 subtype, the chromosomal number differences among individual cells are particularly obvious. Although these differences have been observed in multiple studies, few studies have focused on identifying and characterizing different subpopulations within the HeLa-CCL2 cell line and exploring the potential mechanisms of their differences. This is due to the limitations of existing single-cell analysis tools, and advanced single-cell multi-omics methods are needed to fill the key gaps in the current understanding of karyotypic heterogeneity.

[0004] Recently, the advancement of multi-omics sequencing technologies has made it possible to simultaneously analyze the DNA methylome, transcriptome, and genome. These technologies include scM&T-seq (Angermueller, C. et al. Parallel single-cell sequencing links transcriptional and epigenetic heterogeneity. Nat. Methods 13, 229-232 (2016)), scMT-seq (Hu, Y. et al. Simultaneous profiling of transcriptome and DNA methylome from a single cell. Genome Biol. 17 (2016)), scTrio-seq2 (Bian, S. et al. Single-cell multiomics sequencing and analyses of human colorectal cancer. Science 362, 1060-1063 (2018)), Smart-RRBS (Gu, H. et al. Smart-RRBS for single-cell methylome and transcriptome analysis. Nat. Protocols 16, 4004-4030 (2021)), scNMT-seq (Clark, S. J. et al. scNMT-seq enables joint profiling of chromatin accessibility DNA methylation and transcription in single cells. Nat. Commun. 9 (2018)), and scNOMeRe-seq (Wang, Y. et al. Single-cell multiomics sequencing reveals the functional regulatory landscape of early embryos. Nat. Commun. 12 (2021)). These multi-omics methods vary in terms of CpG site coverage, genomic feature characterization, sequencing depth, and the number of genes jointly covered by the methylome and transcriptome at the gene level, highlighting their diverse applications in different research scenarios. However, so far, there has been no single-cell multi-omics technology capable of karyotyping analysis, thus preventing the epigenetic characterization of chromosomal instability in tumor cells. Summary of the Invention

[0005] To solve the above technical problems, that is, the existing single-cell multi-omics technology cannot perform karyotype analysis, and thus cannot achieve the epigenetic feature analysis of chromosomal instability in tumor cells, the present invention provides a method for analyzing the epigenetic features of CIN (chromosomal instability) in tumor cells. This method can simultaneously measure DNA methylation, gene expression, and karyotype of single cells, and on this basis, analyze the correlation between DNA methylation, gene expression, and karyotype heterogeneity in tumor cells, which helps to more deeply understand the mechanism of tumorigenesis and evolution.

[0006] The specific technical solution of the present invention is as follows:

[0007] A method for analyzing the epigenetic features of CIN in tumor cells, comprising:

[0008] 1) Extract RNA and genomic DNA from a single tumor cell; prepare the RNA into an RNA library and then sequence it; perform MspI digestion reaction, bisulfite conversion, random primer amplification, and library amplification on the genomic DNA in sequence, and sequence it after obtaining the DNA library;

[0009] 2) After preprocessing the DNA library sequencing data, align it to the genomic reference sequence, identify differentially methylated regions, and perform karyotype analysis. Divide the tumor cells into different subgroups according to karyotype heterogeneity among single cells;

[0010] 3) After preprocessing the RNA library sequencing data, align it to the genomic reference sequence and quantify the gene expression level;

[0011] 4) Perform correlation analysis according to the results obtained in 2) and 3).

[0012] By using the method of the present invention, it is possible to simultaneously measure DNA methylation, gene expression, and karyotype of single cells, filling the gap in karyotype analysis of the existing single-cell multi-omics sequencing technology, thereby more intuitively and reliably distinguishing the differences between the same type of tumor cells (for example, HeLa cells were used as an example in the embodiment), identifying chromosomal karyotype heterogeneity, and clarifying the relationship between this phenomenon and epigenetic (DNA methylation) and transcriptome (gene expression) factors. Therefore, the method of the present invention can be used to explore the epigenetic features of chromosomal instability (CIN) in tumor cells, helps to more deeply understand the mechanism of tumorigenesis and evolution, and is expected to develop new treatment strategies for CIN in tumors.

[0013] The combination of different single-cell omics methods can analyze the biological characteristics of cells at different omics levels, study the correlation and synergy among various omics, and produce an effect of 1 + 1 > 2 in research such as cell subset classification and cell regulation. The difficulty lies in the fact that the scope jointly covered by multiple omics is too small, affecting the comprehensiveness of the analysis. The genomic scope and the number of genes that each single-cell single omics can cover are extremely limited, and the scope jointly covered by multiple omics technologies is even scarcer.

[0014] To overcome the above problems, a special DNA methylation library construction method is designed in the present invention. It has a high CpG coverage rate, genomic coverage rate, and read depth, can effectively capture the epigenomic landscape within the entire genome, and shows higher coverage in CGI, CGI shore, exon, intron, 5'UTR, 3'UTR, SINE, LINE, and LTR regions, as well as proximal promoters, the first intron, other introns, and distal non-coding regions. Moreover, it shows a more uniform coverage in regulatory regions and non-coding regions and on all chromosomes. The uniform coverage and higher read depth help to more accurately estimate the methylation level, improve the reliability and comprehensiveness of chromosome karyotyping, are beneficial to analyzing the CIN epigenetic characteristics of tumor cells at the multi-omics level, and more accurately reveal cell heterogeneity. In addition, in terms of gene expression analysis, the method of the present invention also has a high gene coverage rate, which is beneficial to improving the accuracy of the analysis results.

[0015] Preferably, in step 1), the specific process of preparing the DNA library includes:

[0016] 1.1) Adding an MspI digestion reaction solution to the obtained genomic DNA for digestion reaction;

[0017] 1.2) Adding a CT conversion reagent to the product obtained in 1.1) for bisulfite conversion reaction;

[0018] 1.3) Adding random primers and dNTP to the product obtained in 1.2) for random primer amplification reaction;

[0019] 1.4) Adding exonuclease I to the product obtained in 1.3) for digestion reaction to digest the primers, and then purifying the product;

[0020] 1.5) Adding universal PCR primers and index primers to the product obtained in 1.4) for library amplification reaction, and purifying the product.

[0021] Further, step 1.3) specifically includes: adding Random Hexamer Primer and dNTP to the product obtained in 1.2), reacting at 60-65°C for 30-40 min, then pausing the reaction, adding Klenow fragment, and heating to 35-37°C at a rate of 0.5-1.5°C / 1 s, and continuing the reaction for 30-40 min.

[0022] Further, the sequence of the Random Hexamer Primer is as shown in SEQ ID NO.1, specifically as follows: 5’-TACACGACGCTCTTCCGATCTNNNNNN-3’.

[0023] Further, in step 1.5), the specific process of the amplification reaction includes: reacting at 90-95°C for 3-5 min, and then performing 12-15 cycles of the following reactions: reacting at 96-98°C for 20-30 s, reacting at 60-65°C for 30-35 s, reacting at 70-72°C for 1-2 min; after completing the above cycles, continuing to react at 70-72°C for 3-5 min.

[0024] Preferably, in step 1), after completing the MspI digestion reaction, terminal repair and A-tailing reaction are carried out, and then bisulfite conversion is carried out.

[0025] Further, the specific process of the terminal repair and A-tailing reaction includes: adding Klenow fragment and dNTP to the product after completing the MspI digestion reaction to carry out terminal repair and A-tailing reaction.

[0026] Preferably, in step 1), the method for extracting RNA and genomic DNA from a single tumor cell is to lyse the cells, adsorb genomic DNA using magnetic beads, and separate the magnetic beads adsorbed with genomic DNA from the supernatant containing RNA by centrifugation; during the preparation of the DNA library, before performing the MspI digestion reaction, genomic DNA resuspension and protease digestion are carried out first.

[0027] Preferably, in step 2), the specific process of the karyotype analysis includes:

[0028] 2.1) Using AneuFinder software, inputting the sorted BAM file and the reference genome index; obtaining the base number of the reads aligned to the corresponding reference genome position, and estimating the copy status of each gene fragment accordingly; using AneuFinder software to draw a chromosome overview map;

[0029] 2.2) Calculating the chromosome ploidy.

[0030] Further, in step 2.2), the formula for calculating chromosome ploidy is as follows:

[0031]

[0032] Where: P is the chromosome ploidy, s is the ploidy characteristic of the chromosome segment (somy - state) of each segment, s > 1, and L is the length of the segment corresponding to the state of each chromosome segment (somy - state).

[0033] Further, in step 2.1), the pileups method of PySAM is used to obtain the number of bases of the reads aligned to the corresponding reference genome position.

[0034] Preferably, step 4) is specifically: according to the results obtained in 2) and 3), analyze the correlation among DNA methylation differences, gene expression differences, and karyotype heterogeneity among different subgroups.

[0035] Further, the specific process of step 4) includes:

[0036] 4.1) According to the results obtained in step 2) and 3), analyze the correlation between DNA methylation differences and gene expression differences among different subgroups;

[0037] 4.2) According to the results obtained in step 2) and 3), identify the genes that have differential methylation in the promoter region and may be related to karyotype differences among different subgroups, and then perform progression - free survival (PFS) analysis on these genes.

[0038] Further, in step 4.1), the MethGET program is used to analyze the correlation between DNA methylation differences and gene expression differences among different subgroups.

[0039] Preferably, in step 2), the specific process of the pretreatment includes: removing the sequencing adapters and their reverse complementary sequences, trimming the low - quality bases at the 3' end until the average base quality in the sliding window reaches a Phred score > 20, and discarding the paired reads with a length less than 30 bp after trimming.

[0040] Preferably, in step 3), the specific process of the pretreatment includes: removing the sequencing adapters, amplification primers, and low - quality bases at the read ends.

[0041] Preferably, in step 3), the gene expression level includes read count and reads per kilobase per million mapped reads (RPKM).

[0042] Preferably, in steps 2) and 3), the genomic reference sequence is the human hg38 genomic reference sequence.

[0043] Preferably, in step 1), according to the MATQ protocol, the obtained RNA is prepared into an RNA library.

[0044] Compared with the prior art, the present invention has the following advantages:

[0045] (1) The method of the present invention can simultaneously measure DNA methylation, gene expression and karyotype of single cells, identify karyotype heterogeneity among different subgroups of the same type of tumor cells, and clarify the relationship between them and epigenetic and transcriptomic factors, which will help to more deeply understand the mechanism of tumorigenesis and evolution, and is expected to develop new treatment strategies for CIN in tumors.

[0046] (2) In the present invention, a special DNA methylation library construction method is designed, which has higher coverage and read depth in DNA methylome analysis, and has a relatively uniform coverage at different positions on chromosomes and different chromosomes, which will help to improve the reliability and comprehensiveness of the analysis results. Description of the Drawings

[0047] Figure 1 is the CpG methylation rate of scMulOmi at different read depths.

[0048] Figure 2 is the methylation rate near the TSS and TES sites measured by scMulOmi in multiple HeLa cells.

[0049] Figure 3 is the chromosome overview map of HeLa, GM12878 and K562 cell lines.

[0050] Figure 4 is the karyotype status of two HeLa cell subgroups. Among them, "scMulOmi_HeLa1" is the HeLa subgroup in the high ploidy state, and "scMulOmi_HeLa2" is the HeLa subgroup in the low ploidy state.

[0051] Figure 5 is the epigenetic and transcriptomic landscape of the HeLa cell subgroup. Among them, Figure 5 A is the heat map visualization of DNA methylation and gene expression data between g1 and g2 of the HeLa cell subgroup; Figure 5 B is the gene-level correlation between DNA methylation changes and gene expression changes between g1 and g2 of the HeLa cell subgroup. The red dots indicate the differential genes of DNA methylation and gene expression (bivariate Gaussian mixture model; p value < 10 -4 ). Figure 5 C is the scatter plot and fitting curve of DNA methylation and related gene expression. Figure 5D represents the correlation between promoter methylation level (y-axis) and gene expression value (x-axis), and the correlation coefficient R between promoter methylation and expression level is -0.111.

[0052] Figure 6 This is a Kaplan-Meier plot for the progression-free survival (PFS) analysis of patients with cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC). Among them, the upper curve is the "Low 57 Signatures Group", which is a subgroup of HeLa cells with lower methylation values of 57 differentially methylated genes; the upper curve is the "High 57 Signatures Group", which is a subgroup of HeLa cells with higher methylation values of 57 differentially methylated genes. Detailed implementation manners

[0053] The present invention will be further described below in conjunction with embodiments.

[0054] General embodiment

[0055] A method for analyzing the epigenetic characteristics of chromosomal instability (CIN) in tumor cells, comprising:

[0056] 1) Extract RNA and genomic DNA from a single tumor cell; prepare the RNA into an RNA library and then sequence it; perform MspI digestion reaction, bisulfite conversion, random primer amplification, and library amplification on the genomic DNA in sequence, and sequence it after obtaining a DNA library;

[0057] 2) After preprocessing the DNA library sequencing data, align it to the genomic reference sequence, identify differentially methylated regions, and perform karyotype analysis. Divide the tumor cells into different subgroups according to karyotype heterogeneity between single cells;

[0058] 3) After preprocessing the RNA library sequencing data, align it to the genomic reference sequence and quantify the gene expression level;

[0059] 4) Perform correlation analysis according to the results obtained in 2) and 3).

[0060] As a specific implementation manner, in step 1), according to the MATQ protocol, the obtained RNA is prepared into an RNA library.

[0061] As a specific implementation manner, in step 1), the specific process of preparing the DNA library includes:

[0062] 1.1) Add an MspI digestion reaction solution to the obtained genomic DNA for digestion reaction;

[0063] 1.2) Add a CT conversion reagent to the product obtained in 1.1) for bisulfite conversion reaction;

[0064] 1.3) Add random primers and dNTPs to the product obtained in 1.2), and perform a random primer amplification reaction;

[0065] 1.4) Add exonuclease I to the product obtained in 1.3), perform an enzymatic digestion reaction to digest the primers, and then purify the product;

[0066] 1.5) Add universal PCR primers and index primers to the product obtained in 1.4), perform a library amplification reaction, and purify the product.

[0067] As a specific embodiment, in step 1), after completing the MspI enzymatic digestion reaction, perform end repair and A-tailing reactions, and then perform bisulfite conversion; the specific process of the end repair and A-tailing reactions includes: adding Klenow fragment and dNTPs to the product after completing the MspI enzymatic digestion reaction, and performing end repair and A-tailing reactions.

[0068] As a specific embodiment, in step 1), the method for extracting RNA and genomic DNA from a single tumor cell is to lyse the cells, adsorb genomic DNA using magnetic beads, and separate the magnetic beads adsorbed with genomic DNA from the supernatant containing RNA by centrifugation; during the preparation of the DNA library, before performing the MspI enzymatic digestion reaction, first resuspend the genomic DNA and perform protease digestion.

[0069] As a specific embodiment, in step 2), the specific process of the pretreatment includes: removing the sequencing adapters and their reverse complementary sequences, trimming the low-quality bases at the 3' end until the average base quality within the sliding window reaches a Phred score > 20, and discarding the paired reads with a length less than 30 bp after trimming.

[0070] As a specific embodiment, in step 2), the specific process of the karyotype analysis includes:

[0071] 2.1) Use the AneuFinder software, input the sorted BAM file and the reference genome index; use the pileups method of PySAM to obtain the number of bases of the reads aligned to the corresponding reference genome positions, and estimate the copy status of each gene fragment accordingly; use the AneuFinder software to draw a chromosome overview map;

[0072] 2.2) Calculate the chromosome ploidy according to the following formula:

[0073]

[0074] Where: P is the chromosome ploidy, s is the ploidy characteristic of the chromosomal segment (somy - state) of each segment, s > 1, and L is the length of the segment corresponding to the state of each chromosomal segment (somy - state).

[0075] As a specific implementation manner, in step 3), the specific process of the pretreatment includes: removing sequencing adapters, amplification primers, and low - quality bases at the read ends.

[0076] As a specific implementation manner, in step 3), the gene expression level includes read count and reads per kilobase per million mapped reads (RPKM).

[0077] As a specific implementation manner, in steps 2) and 3), the human genome reference sequence is the human hg38 genome reference sequence.

[0078] As a specific implementation manner, the specific process of step 4) includes:

[0079] 4.1) According to the results obtained in steps 2) and 3), use the MethGET program to analyze the correlation between DNA methylation differences and gene expression differences among different subgroups;

[0080] 4.2) According to the results obtained in steps 2) and 3), identify genes that have differential methylation in the promoter region and may be related to karyotype differences among different subgroups, and then perform progression - free survival (PFS) analysis on these genes. Specific embodiments

[0082] The present invention will be described below through specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and not to limit the scope of the present invention. Without departing from the spirit and scope of the inventive concept, changes and advantages that can be conceived by those skilled in the art are included in the present invention, and the scope of protection of the present invention is defined by the appended claims and any equivalents thereof.

[0083] Unless otherwise defined, all technical terms and scientific terms used in the present invention have the same meaning as commonly understood by those of ordinary skill in the art to which the present disclosure belongs. The raw materials and equipment used in the present invention are conventional raw materials and equipment in the art, which can be obtained from conventional commercial channels without special instructions; the methods used in the present invention are conventional methods in the art without special instructions.

[0084] In the following embodiments, the single - cell multi - omics sequencing technology used in the present invention is referred to as "scMulOmi".

[0085] Example 1: Cell and data sources

[0086] In the following examples, the cell and data sources used are as follows:

[0087] The HeLa-CCL2 cell line was obtained from the American Type Culture Collection (ATCC). The GM12878 cells were purchased from CoBioer Company in China. The identity of the cell line was verified by comparing the short tandem repeat (STR) typing of the HeLa-CCL2 cell line with 16 STR typings in the German Collection of Microorganisms and Cell Cultures (DSMZ) database, which contains 7 standard reference STR typing markers specified by ATCC for HeLa-CCL2 identification. Public bulk whole-genome bisulfite sequencing (WGBS) data (accession number GSE131098) of HeLa cells, single-cell reduced representation bisulfite sequencing (scXRBS-seq) data (accession number GSE149954) of GM12878 and K562 cells, and DIRECT-seq data (accession number GSE240579) of K562 cells were obtained from the NCBI Gene Expression Omnibus (GEO; https: / / www.ncbi.nlm.nih.gov / geo / ).

[0088] Example 2: Steps of the scMulOmi method

[0089] After extracting DNA and RNA from single cells, preparing DNA and RNA libraries using specific methods and sequencing, scMulOmi analyzes DNA methylation, gene expression, karyotype, and chromosome ploidy of single cells by sequence alignment. The specific steps are as follows:

[0090] 2.1 Library preparation and sequencing

[0091] Single cells were sorted using FACS and individually transferred to 200 μL PCR tubes containing 3 μL of cell lysis buffer. The lysis buffer consisted of 1× RT buffer (Invitrogen), 0.5% NP40, 5 mM DTT, 2 U RNase OUT (Invitrogen), and 0.2 μL of magnetic beads (Invitrogen, cat#65001). The cells were lysed on ice for 10 minutes and then vortexed vigorously for 30 seconds to ensure thorough mixing. The tubes were then placed on a magnetic stand for 5 minutes to separate the magnetic beads containing intact cell nuclei from the supernatant containing RNA transcripts.

[0092] Carefully transfer the supernatant to a new PCR tube, and perform RNA library preparation according to the MATQ-seq protocol (Sheng, K., Cao, W., Niu, Y., Deng, Q. & Zong, C. Effective detection of variation in single-cell transcriptomes using MATQ-seq. Nat. Methods 14, 267-270 (2017)). To evaluate library preparation efficiency and enable quantification of absolute transcription levels, an exogenous RNA control substance (ERCC spike-ins, Ambion) was added at a dilution ratio of 1:10,000 to 1:100,000.

[0093] The genomic DNA retained on the magnetic beads was subjected to DNA library preparation according to the following steps:

[0094] (1) Resuspension of genomic DNA and protease treatment:

[0095] Resuspend the DNA attached to the magnetic beads with 5 μL of water (let stand for 5 min), and incubate at 65 °C for 5 min. Take out 5 μL of the supernatant, add 0.5 μL of 20 mg / mL protease solution (Qiagen), react at 50 °C for 3 h, and then treat at 75 °C for 30 min (inactivation) to obtain product I.

[0096] (2) MspI digestion:

[0097] Prepare the MspI digestion reaction solution (total volume 10 μL) according to the following formula:

[0098] Product I, 5 μL;

[0099] DNase / RNase-free water, 2.5 μL;

[0100] 20 U / μL MspI solution, 1 μL;

[0101] 10×CutSmart buffer, 1 μL;

[0102] 1 pg / μL Lambda DNA, 0.5 μL.

[0103] Incubate this reaction solution at 37 °C for 4 h, then at 70 °C for 20 min, and then store at 4 °C. Obtain product II.

[0104] (3) End repair / A-tailing reaction

[0105] Prepare the end repair / A-tailing mixture according to the following recipe. After vortexing, centrifuge at 9000 g for 1 minute at 4°C: 5 U / μL Klenow fragment, exo-, 1 μL;

[0106] 10× CutSmart buffer, 0.4 μL;

[0107] End repair dNTP mixture (5 mM each), 2.6 μL.

[0108] Add 4 μL of the end repair / A-tailing mixture to each Product II sample to bring the total volume of each sample to 14 μL. After vortexing, centrifuge at 9000 g for 1 min at 4°C.

[0109] Incubate the mixture at 30°C for 30 min in a thermal cycler, then react at 37°C for 40 min, and then react at 75°C for 15 min. Finally, hold at 4°C to obtain Product III. After completion, immediately place the tube on ice.

[0110] (4) Bisulfite conversion:

[0111] Using the Zymo-EZ DNA Methylation-Gold Kits D5005 kit, prepare the bisulfite conversion reaction solution (total volume 85 μL) according to the following recipe:

[0112] CT Conversion Reagent, 65 μL;

[0113] Product III, 14 μL;

[0114] Water, 6 μL.

[0115] React the bisulfite conversion reaction solution at 98°C for 10 min, then react at 64°C for 4 h, and store at 4°C. Add 20 μL of elution buffer to the obtained product to elute the genomic DNA from the magnetic beads to obtain Product IV.

[0116] (5) Random primer amplification:

[0117] Prepare the amplification reaction solution (total volume 24 μL) according to the following recipe:

[0118] Product IV, 19.5 μL;

[0119] 10× NEB buffer 2, 2.5 μL;

[0120] 100 μM Random Hexamer Primer, 1.0 μL;

[0121] dNTP mixture (10 mM each), 1.0 μL.

[0122] The sequence of the above Random Hexamer Primer is shown in SEQ ID NO.1, specifically as follows: 5’-TACACGACGCTCTTCCGATCTNNNNNN-3’.

[0123] React the amplification reaction solution at 65 °C for 30 min, then pause the reaction at 4 °C. After pausing, add 1.0 μL of Klenow exo- (50 U / μL, NEB), maintain at 4 °C for 5 min, then increase the temperature to 37 °C at a rate of 1 °C every 15 s, react at 37 °C for 30 min, and store at 4 °C to obtain Product V.

[0124] (6) Digest the primer:

[0125] Add 2.0 μL of Exonuclease I (20 U / μL, NEB) to Product V and incubate at 37 °C for 1 h. Then purify using 0.8×AMPure XP beads and elute with 12 μL of water to obtain Product VI.

[0126] (7) Library amplification:

[0127] Prepare the library amplification reaction solution (total volume 14 μL) according to the following formula:

[0128] 15 μM universal PCR primer (NEB), 1.0 μL;

[0129] 15 μM PCR index primer (Index Primer), 1.0 μL;

[0130] 2×KAPA HiFi HotStart ReadyMix, 12.0 μL.

[0131] Add 12 μL of Product VII to the library amplification reaction solution, perform library amplification according to the program in Table 1, and then purify twice using AMPure XP beads and elute to obtain a solution containing the DNA library (volume 30 μL).

[0132] Table 1 Library amplification program

[0133]

[0134] The above DNA library preparation process combines MspI digestion and random primer amplification, effectively enriching CpG islands (CGIs) while taking into account the measurement of CpG sites in non-CGI regions. As a PBAT (post-bisulfite adapter tagging) methylation DNA library construction method, this method can alleviate the impact of DNA damage caused by bisulfite conversion on sequencing results.

[0135] The above-prepared RNA and DNA libraries were both sequenced in paired 150 bp mode on the Illumina HiSeq platform.

[0136] 2.2 Sequencing data processing

[0137] For single-cell DNA library sequencing data, the raw sequence read data is first trimmed of adapters and their reverse complementary sequences (4 bp). After adapter trimming, the 3'-ends of the read data are trimmed of poor-quality bases using a 4-bp sliding window until the average base quality within the window reaches a Phred score of 20 or higher. Paired reads with a length less than 30 bp after trimming are discarded. The trimmed read data is aligned to the human hg38 genome reference sequence in paired-end mode using MethylCtools (Hovestadt, V. et al. Decoding the regulatory landscape of medulloblastoma using DNA methylation sequencing. Nature 510, 537-541 (2014)). Public raw data is processed following the same pipeline. We used the mcomp method in the MOABS suite (Sun, D. et al. MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 15 (2014)) to identify differentially methylated regions (DMRs), which requires at least 3 consecutive CpG sites and an average difference in DNA methylation levels greater than 0.3 between different groups. These DMRs were visually verified manually using the Integrative Genomics Viewer IGV (Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24-26 (2011)). If there is a gene (DMG) in these DMRs with hypomethylation in its promoter region, then we will conduct a more in-depth analysis of such genes. For single-cell RNA library sequencing data, the raw sequence reads undergo quality control and preprocessing steps: sequencing adapters, amplification primers, and low-quality bases at the read ends are removed using Fastp (Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, 884-890 (2018)).The quality-controlled read data was aligned to the human reference genome hg38 using HiSAT2 (Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907-915 (2019)). Gene expression levels, including read counts and reads per kilobase per million mapped reads (RPKM), were quantified by StringTie (version 2.2.0) (Kovaka, S. et al. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 20 (2019)) with default parameters. Quality control analysis was performed using FastQC (Andrews, S. R. FastQC: A Quality Control Tool for High Throughput Sequence Data. http: / / www.bioinformatics.babraham.ac.uk / projects / fastqc (2010)).

[0138] 2.3 Karyotype analysis was performed using the processed sequencing data

[0139] Karyotype analysis was performed using the AneuFinder (Liao, W.-W. et al. MethGo: a comprehensive tool for analyzing whole-genome bisulfite sequencing data. BMC Genomics 16, S11 (2015)) software. The input data included the sorted BAM file (the DNA library sequencing data was aligned with the hg38 genome, the fragments were located in the corresponding positions of the human chromosomes, position tags were added, and then the fragments were sorted according to the order of positions, thus obtaining the sorted BAM file) and the reference genome (hg38 genome) index. In AneuFinder, the pileups method of PySAM (version 0.8.0) (Li, H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetic parameter estimation from sequencing data. Bioinformatics 27, 2987-2993 (2011)) was used to obtain the number of bases read that were aligned to the corresponding reference genome positions. Using a window size of 1 Mb, the number of bases read at all positions within each window was accumulated. Then, the copy state of each gene fragment was estimated. The "heatmapGenomewide" function in the AneuFinder software package was used to draw a chromosome overview map similar to a karyotype map, which showed the copy state of each chromosomal segment, and the copy number state of each segment was plotted in the corresponding chromosomal interval in the form of colors. Based on the estimation of the fragment duplication state, the chromosome ploidy was calculated, and the formula is as follows:

[0140]

[0141] Where:

[0142] P is the chromosome ploidy;

[0143] s is the somy-state ploidy characteristic of each fragment (s > 1);

[0144] L is the length of the fragment corresponding to each somy-state.

[0145] 2.4 Epigenetic analysis at the gene level

[0146] The MethGET program (Teng, C.-S., Wu, B.-H., Yen, M.-R. & Chen, P.-Y. MethGET: web-based bioinformatics software for correlating genome-wide DNA methylation and gene expression. BMC Genomics 21(2020)) was used to calculate the correlation of genes with both DNA methylation and transcriptome data, so as to explore the relationship between DNA methylation and gene expression at the single-cell level. The online database Gepia2 (Tang, Z., Kang, B., Li, C., Chen, T. & Zhang, Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 47, 556-560(2019)) was used for progression-free survival (PFS) analysis. Other correlation analyses were performed using the cor.test function in R version 4.3.0.

[0147] Example 3: Analysis of DNA methylation and gene expression in HeLa cells The scMulOmi method in Example 2 was used to analyze the DNA methylation and gene expression (transcriptome) of HeLa-CCL2 cells, and the performance differences between the scMulOmi method and existing methods in analyzing DNA methylation and gene expression were compared.

[0148] 3.1 Method for evaluating DNA methylation sequencing performance

[0149] The wig sum of selected CpGs was calculated by the following method:

[0150] (1) Calculate the total coverage depth St of all CpG sites (C), and the formula is as follows:

[0151]

[0152] (2) For each given threshold n, calculate the sum of depths S(n) of CpG sites (S) with a depth greater than n, and the formula is as follows:

[0153]

[0154] (3) Calculate the ratio t(n) of the sum of depths corresponding to the threshold n to the total coverage depth, and the formula is as follows:

[0155]

[0156] (4) For samples i that use the same sequencing method, calculate the average value of t(n) for these samples and the corresponding standard deviation. The formula is as follows:

[0157]

[0158] Where: N is the number of samples using the same specific sequencing method.

[0159] To reduce the impact of a few extremely high coverage values on the results, coverage values exceeding 200 are reset to this threshold.

[0160] 3.2 Performance evaluation results show that scMulOmi exhibits high performance in both DNA methylome and transcriptome analyses. In DNA methylome analysis, scMulOmi is superior to most existing single-cell bisulfite sequencing (BS-seq) methods in terms of CpG coverage, genomic coverage, and read depth (Tables 2 - 4).

[0161] The scMulOmi method generates 14.9 million to 25.9 million reads per cell, covering 1.9 million to 5.6 million CpG sites. The average bisulfite conversion rate for all samples is 0.98, ranging from 0.95 to 0.99. In addition, the gene coverage of scMulOmi in transcriptome sequencing also exceeds most existing methods and is comparable to the best-performing method (Tables 2 and 3). The number of genes detected by methylome and transcriptome sequencing ranges from 7.3k to 15.4k, with an average of 11.1k (Table 5).

[0162] Table 2 Performance comparison of different methods for single-cell multi-omics (methylome, transcriptome)

[0163]

[0164] Table 3 Performance comparison of different methods for single-cell multi-omics (methylome, transcriptome)

[0165]

[0166] Table 4 Coverage comparison of different single-cell bisulfite sequencing methods for methylated characteristic regions

[0167]

[0168]

[0169] Table 5 Effect of scMulOmi on sequencing the transcriptome and methylome of single cells simultaneously

[0170]

[0171] To evaluate the performance of scMulOmi and other single-cell bisulfite sequencing methods in terms of genome-wide CpG coverage, their coverage in different chromosomal feature regions was compared (Table 4). The results showed that scMulOmi exhibited higher coverage in CGI (75.7%), CGI shore (54.7%), exon (21.0%), intron (43.9%), 5'UTR (30.2%), 3'UTR (26.5%), SINE (14.0%), LINE (6.0%) and LTR (7.7%) regions compared to other single-cell BS-seq methods tested (including scRRBS, scXRBS, scBS and DIRECT).

[0172] In addition, the CpG methylation rates of scMulOmi at different read depths are as Figure 1 shown. It can be seen that the CpG methylation rate of this method remains highly stable at different read depths, providing a basis for accurate karyotype calculation. The methylation rates near the TSS and TES sites in multiple HeLa cells measured by scMulOmi are as Figure 2 shown. It can be seen that there are obvious changes in the methylation rate of HeLa cells near the TSS and TES sites, which is highly consistent with the existing knowledge of organism methylation, indicating that scMulOmi has high accuracy in testing methylation.

[0173] Example 4: Analysis of karyotype heterogeneity of HeLa cells by scMulOmi

[0174] Using the scMulOmi method in Example 2, karyotype analysis was performed on HeLa, GM12878 and K562 cell lines, and the performance differences between the scMulOmi method and existing single-cell BS-seq methods in karyotype analysis were compared.

[0175] The results showed significant karyotype heterogeneity in HeLa cells, in sharp contrast to the GM12878 and K562 cell lines( Figure 3 ), specifically:[[]] Figure 3Shows the chromosomal status of each chromosomal fragment (column) of individual cells (rows) from different cell lines and sequencing methods. The GM12878 and K562 cell lines exhibit minimal karyotypic differences. The chromosomal fragments of GM12878 cells mainly show a 2-copy state, while those of K562 cells are mainly in a 3-copy state. In contrast, HeLa cells show a large variation from 2 to 6 in the copy state of chromosomal fragments.

[0176] There are also significant differences among individual HeLa cells. Based on these differences, the clustering algorithm clearly divides HeLa cells into two distinct subgroups, g1 and g2. The g1 subgroup has a relatively low copy state, mainly concentrated around 3 copies, with relatively small internal differences, similar to the karyotype analysis results of the K562 cell line. In contrast, the g2 subgroup shows a higher copy state, ranging from 3 to 7, and there are also greater internal differences among individual cells ( Figure 3 and Figure 4 ). These findings highlight the large karyotypic heterogeneity in HeLa cells, which can be divided into two distinct subgroups (g1 and g2) based on their copy state differences.

[0177] The results of the analysis of the impact of data obtained by different sequencing methods on karyotype analysis and ploidy estimation are shown in Figure 3 . The karyotype analysis results of HeLa cells using scMulOmi (n = 11) and scXRBSm (n = 9) data are highly consistent, further confirming the heterogeneous karyotype characteristics of this cell line. Similarly, scXRBSm and scXRBS data also produced consistent karyotype results for GM12878 cells. Although the karyotype of HeLa cells generated by scRRBS (n = 17) data is relatively similar to that of scMulOmi and scXRBSm, the overall copy state is lower. In contrast, the copy state of HeLa cells obtained from scBS (n = 11) data is significantly lower than the previous experimental results. Since scRRBS has a low genomic coverage and the results of scBS are also biased, these two methods were excluded from subsequent epigenetic analyses. The ploidy calculation results further support this conclusion. Notably, although not directly comparable, DIRECT produced more dispersed ploidy estimates in K562 cells than scXRBS, indicating that scXRBS may provide more consistent ploidy estimates across different cell lines.

[0178] Example 5: scMulOmi analysis of the epigenetic and transcriptomic landscapes of HeLa cell subgroups To further investigate the epigenetic and transcriptomic differences between the two HeLa cell subgroups (g1 and g2) identified in the karyotype analysis of Example 4, a comprehensive analysis of DNA methylation and gene expression data was performed as described in "2.4 Epigenetic analysis at the gene level" in Example 2. The results are shown in Figure 5 。

[0179] First, we examined the association between DNA methylation changes and gene expression changes at the gene level between the g1 and g2 subgroups ( Figure 5 A). Using a bivariate Gaussian mixture model (p-value < 10 -3 ), we identified genes with significantly different methylation and expression, which are represented by red dots in the figure ( Figure 5 B). This analysis showed that by appropriately clustering HeLa cells, a large number of genes with significant differences in both DNA methylation and gene expression levels could be found.

[0180] Next, we focused on the relationship between differential methylation of CpGs in the promoter region and corresponding gene expression changes ( Figure 5 B). This quadrant plot identified a total of 2,237 differentially expressed genes based on promoter methylation levels and gene expression levels. In these four quadrants, the second quadrant representing low methylation and high expression levels identified the largest number of differentially methylated and expressed genes (715). To further explore the relationship between methylation and expression, we examined the methylation level trends around the promoter region and their association with gene expression levels ( Figure 5 C). The results clearly demonstrated the association between methylation patterns and gene expression levels, with lower methylation levels generally corresponding to higher gene expression and vice versa. Finally, we quantified the correlation between promoter methylation levels and expression levels ( Figure 5 D). The correlation coefficient (R) between promoter methylation and expression levels was -0.111, indicating a negative correlation between the two. This finding supports the view that promoter methylation is associated with gene silencing, with higher methylation levels generally corresponding to lower gene expression.

[0181] Example 6: scMulOmi analysis of epigenetic differences between HeLa cell subgroups

[0182] In this example, the epigenetic differences at the gene level between the two HeLa cell subgroups (g1 and g2) identified in the karyotype analysis of Example 4 were further investigated.

[0183] To exclude genes related to the cell cycle, we intersected the differentially methylated genes between HeLa g1 and g2 with the differentially expressed genes between GM12878 and HeLa, and identified 57 genes that were hypomethylated in the promoter region (the methylation values of each gene in the second group were lower than those in the first group) and might be related to karyotype differences.

[0184] To further characterize these 57 genes, we performed a progression-free survival (PFS) analysis. The PFS analysis showed that these 57 genes had a significant impact on the PFS of patients with cervical squamous cell carcinoma and cervical adenocarcinoma (log-rank test, p = 0.067, Figure 6 ).

[0185] As described above, the above are only the preferred embodiments of the present invention and do not impose any limitations on the present invention. Any simple modifications, changes, and equivalent transformations made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for analyzing epigenetic characteristics of chromosomal instability in tumor cells, characterized in that: include: 1) Extract RNA and genomic DNA from single tumor cells; prepare RNA into RNA library and then sequence it; The genomic DNA was subjected to MspI digestion, end repair and A-tailing reaction, sulfite conversion, amplification with random primers shown in SEQ ID NO.1, library amplification with universal PCR primers and index primers, and sequencing after obtaining the DNA library; 2) After preprocessing, the DNA library sequencing data is aligned to the genome reference sequence to identify differentially methylated regions; using AneuFinder software, input the sorted BAM file and reference genome index; obtain the number of bases read that are aligned to the corresponding reference genome position and estimate the copy status of each gene fragment; use AneuFinder software to draw a chromosome overview; multiply each s by the corresponding L and sum them, and then divide by the sum of all L to obtain the copy number of the chromosomal segment, s is the chromosome segment ploidy characteristic of each fragment, s>1, L is the fragment length corresponding to the status of each chromosome segment; tumor cells are divided into different subgroups according to the karyotype heterogeneity between single cells; 3) RNA library sequencing data is preprocessed and aligned to the genomic reference sequence to quantify gene expression levels; 4) Based on the results obtained in steps 2) and 3), analyze the correlation between DNA methylation differences and gene expression differences between different subgroups, identify genes that are differentially methylated in the promoter region between different subgroups and may be associated with karyotype differences, and then perform progression-free survival analysis on these genes.

2. The method according to claim 1, characterized in that: In step 1), the specific process of preparing the DNA library includes: 1.1) Add MspI digestion reaction solution to the obtained genomic DNA for digestion reaction, and then perform end repair and A-tailing reaction; 1.2) adding a CT conversion reagent to the product obtained in 1.1) to carry out a sulfite conversion reaction; 1.3) adding random primers and dNTPs shown in SEQ ID NO.1 to the product obtained in 1.2) to perform random primer amplification reaction; 1.4) adding exonuclease I to the product obtained in 1.3) to perform an enzyme digestion reaction to digest the primer, and then purifying the product; 1.5) Add universal PCR primers and index primers to the product obtained in 1.4) to perform library amplification reaction and purify the product.

3. The method according to claim 1 or 2, characterized in that: In step 1), the method for extracting RNA and genomic DNA from single tumor cells is to lyse the cells, adsorb the genomic DNA using magnetic beads, and separate the magnetic beads adsorbed with genomic DNA from the supernatant containing RNA by centrifugation; in the process of preparing the DNA library, the genomic DNA is resuspended and digested with protease before the MspI enzyme digestion reaction.

4. The method according to claim 1, characterized in that: In step 2.1), the pileups method of PySAM is used to obtain the number of bases of the reads that are aligned to the corresponding reference genome position.

5. The method according to claim 1, characterized in that In step 4), the MethGET program was used to analyze the correlation between DNA methylation differences and gene expression differences among different subgroups.

6. The method according to claim 1, characterized in that In step 2), the specific process of the preprocessing includes: removing the sequencing adapter and its reverse complementary sequence, trimming the low-quality bases at the 3' end until the average base quality in the sliding window reaches a level of Phred score>20, and discarding paired reads with a length of less than 30bp after trimming.

7. The method according to claim 1, characterized in that In step 3), the specific process of the pretreatment includes: removing sequencing adapters, amplification primers and low-quality bases at the end of the read.

8. The method according to claim 1, characterized in that: In step 3), the gene expression level includes read counts and copy numbers per 10 million aligned reads.

9. The method according to claim 1, characterized in that: In step 1), the obtained RNA is prepared into an RNA library according to the MATQ protocol.

Images