Method for analyzing CIN epigenetic characteristics of tumor cells

Abstract

The invention relates to the technical field of single-cell multi-omics sequencing, and discloses a method for analyzing CIN epigenetic characteristics of tumor cells, which comprises the following steps: preparing a single-cell RNA library and a DNA library, and sequencing; the method comprises the following steps: respectively preprocessing DNA and RNA sequencing data, comparing to obtain a human genome reference sequence, identifying a differential methylation region, carrying out cell karyotype analysis, dividing tumor cells into different subgroups, and quantifying a gene expression level; and analyzing the correlation among DNA methylation difference, gene expression difference and karyotype heterogeneity among different subgroups. According to the method, DNA methylation, gene expression and karyotype of a single cell can be measured at the same time, the correlation among DNA methylation, gene expression and karyotype heterogeneity of tumor cells is analyzed on the basis, and the mechanism of tumor generation and evolution can be understood more deeply.

Description

Technical Field

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

[0002] Chromosomal instability (CIN) is characterized by an increased frequency of chromosome number gains and losses between cells, manifested as heterogeneity in the karyotype. CIN is known to play a key role in the development and evolution of cancer and has been extensively documented in various cancers. Despite the strong association between CIN and cancer, their precise relationship remains a topic of intense debate.

[0003] The limitations of existing single-cell analysis tools have restricted the study 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 that there is significant heterogeneity in karyotype and chromosome stability between different HeLa subtypes, different generations, and the same clone. Although many researchers have documented the heterogeneity of HeLa chromosomes, their underlying mechanisms have not been explored in depth. Notably, in the HeLa-CCL2 subtype, the differences in chromosome number between 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 underlying 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 make up for the key gaps in the current understanding of karyotypic heterogeneity.

[0004] Recent advances in multi-omics sequencing technologies have made it possible to analyze the DNA methylome, transcriptome, and genome simultaneously. 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 DNAmethylome from a single cell. Genome Biol.17(2016)), scTrio-seq2(Bian,S.etal.Single-cell multiomics sequencing and analyzes of human colorectalcancer.Science 362,1060-1063(2018)), Smart-RRBS(Gu,H.etal.Smart-RRBS for single-cell methylome and transcriptome analysis.Nat.Protocols 16,4004-4030(2021)), scNMT-seq(Clark,SJet al.scNMT-seq enables joint profiling ofchromatin accessibility DNA methylation and transcription in single cells.Nat.Commun.9(2018)) and scNOMeRe-seq(Wang,Y.et al.Single-cell multiomicssequencing 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 covered by the methylation group and transcriptome at the gene level, highlighting their diverse applications in different research scenarios. However, to date, there is no single-cell multi-omics technology that can perform karyotype analysis, and therefore it is impossible to achieve epigenetic feature analysis of chromosomal instability in tumor cells. Summary of the invention

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

[0006] The specific technical scheme of the present invention is: A method for analyzing the epigenetic characteristics of CIN in tumor cells, comprising: 1) Extract RNA and genomic DNA from a single tumor cell; prepare RNA into an RNA library and then sequence it; perform MspI restriction digestion, sulfite conversion, random primer amplification and library amplification on the genomic DNA in sequence to obtain a DNA library and then sequence it; 2) After preprocessing, the DNA library sequencing data is aligned to the genomic reference sequence to identify the differentially methylated regions and perform cell karyotype analysis to divide tumor cells into different subgroups based on 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) Perform correlation analysis based on the results obtained in 2) and 3).

[0007] The method of the present invention can simultaneously measure DNA methylation, gene expression and karyotype of a single cell, which fills the gap in karyotype analysis of existing single-cell multi-omics sequencing technology, thereby more intuitively and reliably distinguishing the differences between tumor cells of the same type (for example, HeLa cells are 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 characteristics of chromosomal instability (CIN) of tumor cells, which helps to more deeply understand the mechanism of tumor occurrence and evolution, and is expected to develop a new treatment strategy for CIN in tumors.

[0008] The combination of different single-cell omics methods can analyze the biological characteristics of cells at different omics levels, study the correlation and synergy between omics, and produce a 1+1>2 effect in the study of cell subpopulation classification and cell regulation. The difficulty lies in the fact that the scope that multiple omics can cover together is too small, which affects the comprehensiveness of the analysis. The genome range and number of genes that can be covered by each single-cell single-omics are extremely limited, and the coverage of multiple omics technologies together is even scarcer.

[0009] In order to overcome the above problems, a special DNA methylation library construction method is designed in the present invention, which has high CpG coverage, genome coverage and reading depth, and can effectively capture the epigenetic landscape within the whole genome, and shows higher coverage in CGI, CGI shore, exon, intron, 5'UTR, 3'UTR, SINE, LINE and LTR regions, as well as proximal promoter, first intron, other introns and distal non-coding regions, and shows more uniform coverage in regulatory regions and non-coding regions and on all chromosomes. Uniform coverage and higher reading depth help to more accurately estimate methylation levels, improve the reliability and comprehensiveness of chromosome karyotyping, and are conducive to analyzing the epigenetic characteristics of tumor cell CIN at the multi-omics level, and more accurately revealing cell heterogeneity. In addition, in terms of gene expression analysis, the method of the present invention also has a higher gene coverage, which is conducive to improving the accuracy of the analysis results.

[0010] Preferably, in step 1), the specific process of preparing the DNA library includes: 1.1) Adding MspI digestion reaction solution to the obtained genomic DNA for digestion 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 to the product obtained in 1.2) to perform a 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.

[0011] 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, heating to 35-37° C. at a rate of 0.5-1.5° C. / 1 s, and continuing the reaction for 30-40 min.

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

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

[0014] Preferably, in step 1), after the MspI digestion reaction is completed, end repair and A-tailing reaction are performed, and then sulfite conversion is performed.

[0015] Furthermore, the specific process of the end repair and A-tailing reaction includes: adding Klenow fragment and dNTP to the product after the MspI enzyme digestion reaction, and performing end repair and A-tailing reaction.

[0016] Preferably, 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 the genomic DNA from the supernatant containing the 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.

[0017] Preferably, in step 2), the specific process of the cell karyotype analysis includes: 2.1) Using AneuFinder software, input the sorted BAM file and reference genome index; obtain the number of bases of the reads that are aligned to the corresponding reference genome position, and estimate the copy status of each gene fragment based on this; use AneuFinder software to draw a chromosome overview map; 2.2) Calculate chromosome ploidy.

[0018] Further, in step 2.2), the formula for calculating chromosome ploidy is as follows: Where: P is the chromosome ploidy, s is the chromosome segment (somy-state) ploidy characteristics of each fragment, s>1, L is the fragment length corresponding to each chromosome segment (somy-state) state.

[0019] Furthermore, in step 2.1), the pileups method of PySAM is used to obtain the number of bases read that are aligned to the corresponding reference genome position.

[0020] Preferably, step 4) specifically comprises: analyzing the correlation among DNA methylation differences, gene expression differences and karyotype heterogeneity among different subgroups based on the results obtained in 2) and 3).

[0021] Furthermore, the specific process of step 4) includes: 4.1) Based on the results obtained in steps 2) and 3), analyze the correlation between DNA methylation differences and gene expression differences among different subgroups; 4.2) Based on the results obtained in steps 2) and 3), genes that are differentially methylated in the promoter region between different subgroups and may be associated with karyotype differences are identified, and then these genes are analyzed for progression-free survival (PFS).

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

[0023] Preferably, 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.

[0024] Preferably, 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.

[0025] Preferably, in step 3), the gene expression level includes read counts and copies per 10 million paired reads (RPKM).

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

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

[0028] Compared with the prior art, the present invention has the following advantages: (1) The method of the present invention can simultaneously measure DNA methylation, gene expression and karyotype of single cells, identify karyotype heterogeneity between different subgroups of the same tumor cells, and clarify its relationship with epigenetic and transcriptomic factors, which will help to more deeply understand the mechanism of tumor occurrence and evolution, and is expected to develop new treatment strategies for CIN in tumors.

[0029] (2) The present invention designs a special DNA methylation library construction method, which has higher coverage and reading depth in DNA methylation group analysis, and has a more uniform coverage at different chromosome positions and on different chromosomes, which will help improve the reliability and comprehensiveness of the analysis results. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0031] Figure 2 The methylation rates near TSS and TES sites in multiple HeLa cells measured by scMulOmi.

[0032] Figure 3 Chromosome overview of HeLa, GM12878 and K562 cell lines.

[0033] Figure 4 The karyotype status of two HeLa cell subgroups, among which "scMulOmi_HeLa1" is a high ploidy HeLa subgroup, and "scMulOmi_HeLa2" is a low ploidy HeLa subgroup.

[0034] Figure 5 The epigenetic and transcriptomic landscape of HeLa cell subgroups. Figure 5 A is a heat map visualization of DNA methylation and gene expression data between HeLa cell subgroups g1 and g2; Figure 5 B shows the gene-level correlation between DNA methylation changes and gene expression changes between HeLa cell subgroups g1 and g2. Red dots indicate differentially expressed genes in terms 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 5 D is 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.

[0035] Figure 6 Kaplan-Meier plot for progression-free survival (PFS) analysis of patients with cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC). The upper curve is the "Low 57Signatures Group", which is a HeLa cell subgroup with lower methylation values ​​of 57 differentially methylated genes; the upper curve is the "High 57Signatures Group", which is a HeLa cell subgroup with higher methylation values ​​of 57 differentially methylated genes. DETAILED DESCRIPTION

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

[0037] Overall embodiment A method for analyzing the epigenetic characteristics of CIN in tumor cells, comprising: 1) Extract RNA and genomic DNA from a single tumor cell; prepare RNA into an RNA library and then sequence it; perform MspI restriction digestion, sulfite conversion, random primer amplification and library amplification on the genomic DNA in sequence to obtain a DNA library and then sequence it; 2) After preprocessing, the DNA library sequencing data is aligned to the genomic reference sequence to identify the differentially methylated regions and perform cell karyotype analysis to divide tumor cells into different subgroups based on 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) Perform correlation analysis based on the results obtained in 2) and 3).

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

[0039] As a specific implementation, in step 1), the specific process of preparing the DNA library includes: 1.1) Adding MspI digestion reaction solution to the obtained genomic DNA for digestion 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 to the product obtained in 1.2) to perform a 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.

[0040] As a specific implementation method, in step 1), after completing the MspI enzyme digestion reaction, end repair and A-tailing reaction are performed, and then sulfite conversion is performed; the specific process of the end repair and A-tailing reaction includes: adding Klenow fragment and dNTP to the product after completing the MspI enzyme digestion reaction, performing end repair and A-tailing reaction.

[0041] As a specific embodiment, in step 1), the method for extracting RNA and genomic DNA from single tumor cells is to lyse the cells, and use magnetic beads to adsorb the genomic DNA, and separate the magnetic beads adsorbed with the genomic DNA from the supernatant containing the RNA by centrifugation; in the process of preparing the DNA library, before performing the MspI enzyme digestion reaction, the genomic DNA is resuspended and digested with protease.

[0042] As a specific implementation, 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.

[0043] As a specific implementation, in step 2), the specific process of the cell karyotype analysis includes: 2.1) Using AneuFinder software, input the sorted BAM file and reference genome index; using PySAM's pileups method to obtain the number of bases of the reads that are aligned to the corresponding reference genome position, and estimate the copy status of each gene fragment; using AneuFinder software to draw a chromosome overview map; 2.2) Calculate chromosome ploidy according to the following formula: Where: P is the chromosome ploidy, s is the chromosome segment (somy-state) ploidy characteristics of each fragment, s>1, L is the fragment length corresponding to each chromosome segment (somy-state) state.

[0044] As a specific implementation, 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.

[0045] As a specific embodiment, in step 3), the gene expression level includes read counts and copies per 10 million paired reads (RPKM).

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

[0047] As a specific implementation, the specific process of step 4) includes: 4.1) Based on the results obtained in steps 2) and 3), the MethGET program was used to analyze the correlation between DNA methylation differences and gene expression differences among different subgroups; 4.2) Based on the results obtained in steps 2) and 3), genes that are differentially methylated in the promoter region between different subgroups and may be associated with karyotype differences are identified, and then these genes are analyzed for progression-free survival (PFS). Specific embodiments The present invention is described below by specific examples. It should be understood that these examples are only used to illustrate the present invention and are not used 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 thought of by those skilled in the art are included in the present invention, and the attached claims and any equivalents thereof are the protection scope of the present invention.

[0049] Unless otherwise defined, all technical terms and scientific terms used in the present invention have the same meanings as those commonly understood by ordinary technicians in the field to which the present disclosure belongs. Unless otherwise specified, the raw materials and equipment used in the present invention are conventional raw materials and equipment in the field and can be obtained from conventional commercial channels; unless otherwise specified, the methods used in the present invention are conventional methods in the field.

[0050] In the following examples, the single-cell multi-omics sequencing technology used in the present invention is referred to as “scMulOmi”.

[0051] Example 1: Cell and data sources In the following examples, the cells and data sources used are as follows: The HeLa-CCL2 cell line was obtained from the American Cell Culture Collection (ATCC). GM12878 cells were purchased from CoBioer, China. The identity of the HeLa-CCL2 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 the 7 standard reference STR typing markers specified by ATCC for HeLa-CCL2 identification. We obtained the public batch whole-genome bisulfite sequencing (WGBS) data of HeLa cells (accession number GSE131098), single-cell reduced representation bisulfite sequencing (scXRBS-seq) data of GM12878 and K562 cells (accession number GSE149954), and DIRECT-seq data of K562 cells (accession number GSE240579) from the NCBI Gene Expression Omnibus Database (GEO; https: / / www.ncbi.nlm.nih.gov / geo / ).

[0052] Example 2: Steps of the scMulOmi method scMulOmi extracts DNA and RNA from single cells, uses a specific method to prepare DNA and RNA libraries and sequence them, and then analyzes DNA methylation, gene expression, karyotype and chromosome ploidy of single cells through sequence alignment. The specific steps are as follows: 2.1 Library preparation and sequencing Single cells were sorted using FACS and transferred individually into 200 μL PCR tubes containing 3 μL cell lysis buffer. The lysis buffer consisted of 1×RT buffer (Invitrogen), 0.5% NP40, 5 mM DTT, 2U RNase OUT (Invitrogen), and 0.2 μL magnetic beads (Invitrogen, cat#65001). The cells were lysed on ice for 10 minutes and then shaken 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 nuclei from the supernatant containing RNA transcripts.

[0053] The supernatant was carefully transferred to a new PCR tube and RNA library preparation was performed 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)). In order to evaluate the efficiency of library preparation and enable quantitative absolute transcription levels, exogenous RNA control substances (ERCC spike-ins, Ambion) were added at a dilution ratio of 1:10000 to 1:100000.

[0054] The genomic DNA retained on the magnetic beads is then used to prepare the DNA library according to the following steps: (1) Resuspend genomic DNA and treat with protease: The DNA attached to the magnetic beads was resuspended in 5 μL of water (standing for 5 min) and incubated at 65°C for 5 min. 5 μL of the supernatant was taken out, 0.5 μL of 20 mg / mL protease solution (Qiagen) was added, and after reacting at 50°C for 3 h, it was treated at 75°C for 30 min (inactivation) to obtain product I.

[0055] (2) MspI digestion: Prepare the MspI digestion reaction solution (total volume 10 μL) according to the following formula: Product I, 5 μL; DNase / RNase-free water, 2.5μL; 20U / μL MspI solution, 1μL; 10× CutSmart buffer, 1 μL; 1pg / μL Lambda DNA, 0.5μL.

[0056] The reaction solution was reacted at 37° C. for 4 h, then at 70° C. for 20 min, and then stored at 4° C. to obtain product II.

[0057] (3) End repair / A-tailing reaction Prepare the end repair / A-tailing mixture according to the following recipe, vortex to mix, and centrifuge at 9000 g for 1 min at 4°C: 5 U / μL Klenow fragment, exo-, 1 μL; 10× CutSmart buffer, 0.4 μL; End repair dNTP mix (5 mM each), 2.6 μL.

[0058] Add 4 μL of end repair / A-tailing mix to each Product II sample to bring the total volume of each sample to 14 μL. Vortex to mix and centrifuge at 9000 g for 1 min at 4°C.

[0059] After the mixture was incubated at 30°C for 30 min in a thermal cycler, it was reacted at 37°C for 40 min, and then at 75°C for 15 min. Finally, it was kept at 4°C to obtain product III. After completion, the tube was immediately placed on ice.

[0060] (4) Sulfite conversion: Use Zymo-EZ DNA Methylation-Gold Kits D5005 to prepare the sulfite conversion reaction solution (total volume 85 μL) according to the following formula: CT Conversion Reagent, 65 μL; Product III, 14 μL; Water, 6 μL.

[0061] The sulfite conversion reaction solution was reacted at 98° C. for 10 min, then at 64° C. for 4 h, and stored at 4° C. 20 μL of elution buffer was added to the obtained product to elute the genomic DNA from the magnetic beads to obtain product IV.

[0062] (5) Random primer amplification: Prepare the amplification reaction solution (total volume 24 μL) according to the following formula: Product IV, 19.5 μL; 10× NEB buffer 2, 2.5 μL; 100 μM Random Hexamer Primer, 1.0 μL; dNTP mixture (10 mM each), 1.0 μL.

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

[0064] The amplification reaction solution was reacted at 65°C for 30 minutes, and then the reaction was paused at 4°C. After the pause, 1.0 μL Klenow exo- (50U / μL, NEB) was added, and after maintaining at 4°C for 5 minutes, the temperature was raised to 37°C at a rate of 1°C / 15s, and the reaction was carried out at 37°C for 30 minutes. The product V was obtained.

[0065] (6) Digestion primer: 2.0 μL of exonuclease I (20 U / μL, NEB) was added to product V and incubated at 37°C for 1 h. Then, the product was purified using 0.8×AMPure XP beads and eluted with 12 μL of water to obtain product VI.

[0066] (7) Library amplification: Prepare the library amplification reaction solution (total volume 14 μL) according to the following formula: 15 μM universal PCR primer (NEB), 1.0 μL; 15 μM PCR Index Primer, 1.0 μL; 2×KAPA HiFi HotStart ReadyMix, 12.0μL.

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

[0068] Table 1 Library amplification procedure 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. This method is a PBAT (post-bisulfite adapter tagging) methylated DNA library construction method that can alleviate the impact of DNA damage caused by sulfite conversion on sequencing results.

[0069] The RNA and DNA libraries prepared above were sequenced in 150 bp pairs on the Illumina HiSeq platform.

[0070] 2.2 Sequencing data processing For single cell DNA library sequencing data, the original sequence read (read) data first removes the joint and its reverse complementary sequence (4bp). After the joint is removed, by using a 4bp sliding window, the bases with poor quality at the 3' end of the read data are trimmed until the average base quality in the window reaches a Phred score of more than 20. After trimming, paired reads with a length less than 30bp will be discarded. The trimmed read data are 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)). After downloading the public raw data, the same process is followed. 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 between different groups greater than 0.3. These DMRs were manually visually verified by the integrated genome browser IGV (Robinson, J Tetal. Integrative genomics viewer. Nat. Biotechnol. 29, 24-26 (2011)). If there is a gene (DMG) in these DMRs that has hypomethylation in its promoter region, we will perform a more in-depth analysis of this gene. For single-cell RNA library sequencing data, the raw sequence reads underwent quality control and preprocessing steps: Fastp (Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, 884-890 (2018)) was used to remove low-quality bases in sequencing adapters, amplification primers, and read ends.The quality-controlled read data were aligned to the human reference genome hg38 using HiSAT2 (Kim, D., Paggi, JM, Park, C., Bennett, C. & Salzberg, SL Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907-915 (2019)). Gene expression levels, including read counts and copies per 10 million aligned 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 used FastQC (Andrews, SR FastQC: A Quality Control Tool for High Throughput Sequence Data. http: / / www.bioinformatics.babraham.ac.uk / projects / fastqc (2010)).

[0071] 2.3 Using processed sequencing data for karyotype analysis Karyotype analysis was performed using AneuFinder (Liao, W.-W. et al. MethGo: a comprehensive tool for analyzing whole-genome bisulfite sequencing data. BMC Genomics 16, S11 (2015)). The input data included the sorted BAM file (DNA library sequencing data was aligned with the hg38 genome, each fragment was located in the corresponding position of the human chromosome, the position label was added, and then the fragments were sorted in the order of the position, that is, the sorted BAM file was obtained) and the reference genome (hg38 genome) index. In AneuFinder, the pileups method of PySAM (version 0.8.0) (Li, HA 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 position. Using a window size of 1Mb, the number of bases read at each position in the window is accumulated. Then the copy status of each gene segment is estimated. The "heatmapGenomewide" function in the AneuFinder software package is used to draw a chromosome overview similar to a karyotype diagram, which shows the copy status of each segment of the chromosome, and the copy number status of each segment is plotted in color on the corresponding chromosome interval. Based on the estimation of the replication status of the segment, the chromosome ploidy is calculated, and the formula is as follows: in: P is the chromosome ploidy; s is the somy-state ploidy characteristic of each fragment (s>1); L is the length of the fragment corresponding to each somy-state state.

[0072] 2.4 Epigenetic analysis at the gene level 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, thereby exploring the relationship between DNA methylation and gene expression at the single-cell level. Progression-free survival (PFS) analysis was performed using 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)). Other correlation analyses were performed using the cor.test function in R version 4.3.0.

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

[0074] 3.1 DNA methylation sequencing performance evaluation method The wig sum of selected CpGs was calculated by the following method: (1) Calculate the total coverage depth St of all CpG sites (C) using the following formula: (2) For each given threshold n, calculate the sum of the depths S(n) of CpG sites (S) with a depth greater than n, using the following formula: (3) Calculate the ratio t(n) of the sum of the depths corresponding to the threshold n to the total coverage depth. The formula is as follows: (4) For sample i using the same sequencing method, calculate the mean value and corresponding standard deviation of t(n) of these samples The formula is as follows: Where: N is the number of samples using the same specific sequencing method.

[0075] To mitigate the impact of a few extremely high coverage values ​​on the results, coverage values ​​above 200 were reset to this threshold.

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

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

[0078] Table 2 Performance comparison of different methods of single-cell multi-omics (methylome, transcriptome) Table 3 Performance comparison of different methods of single-cell multi-omics (methylome, transcriptome) Table 4 Comparison of coverage of methylation feature regions by different single-cell bisulfite sequencing methods Table 5 Effect of scMulOmi on simultaneous sequencing of transcriptome and methylome of a single cell 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 characteristic regions was compared (Table 4). The results showed that compared with other single-cell BS-seq methods tested (including scRRBS, scXRBS, scBS and DIRECT), scMulOmi showed higher coverage in CGI (75.7%), CGIshore (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.

[0079] In addition, the CpG methylation rates of scMulOmi at different read depths are as follows Figure 1 As shown in the figure, it can be seen that this method maintains a highly stable CpG methylation rate at different read depths, providing a basis for accurate calculation of karyotypes. The methylation rates near TSS and TES sites in multiple HeLa cells measured by scMulOmi are shown in the figure. Figure 2 As shown, it can be seen that the methylation rate of HeLa cells changes significantly near the TSS and TES sites, which is highly consistent with the existing knowledge about methylation in organisms, indicating that scMulOmi has high accuracy in testing methylation.

[0080] Example 4: Analysis of HeLa cell karyotype heterogeneity using scMulOmi The scMulOmi method in Example 2 was used to perform karyotype analysis on HeLa, GM12878 and K562 cell lines, and the performance difference between the scMulOmi method and the existing single-cell BS-seq method in karyotype analysis was compared.

[0081] The results showed that HeLa cells have significant heterogeneity in karyotype, in sharp contrast to GM12878 and K562 cell lines ( Figure 3 ), specifically: Figure 3 The chromosome status of each chromosome segment (column) of individual cells (rows) from different cell lines and sequencing methods is shown. GM12878 and K562 cell lines show minimal karyotype differences, with GM12878 cells showing mainly 2-copy states for chromosome segments and K562 cells showing mainly 3-copy states for segments. In contrast, HeLa cells show a large variation in the copy states of chromosome segments, from 2 to 6.

[0082] HeLa cells also have significant differences between individual cells. Based on these differences, the clustering algorithm clearly divides HeLa cells into two different subgroups, g1 and g2. The copy state of the g1 subgroup is relatively low, mainly concentrated around 3 copies, with small internal differences, which is similar to the karyotype analysis results of the K562 cell line. The g2 subgroup shows a higher copy state, ranging from 3 to 7, and the internal differences between individual cells are also greater ( Figure 3 and Figure 4 These findings highlight the enormous karyotypic heterogeneity in HeLa cells, which can be distinguished into two clear subgroups (G1 and G2) based on differences in their copy states.

[0083] 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 results of karyotyping of HeLa cells using scMulOmi (n = 11) and scXRBSm (n = 9) data were highly consistent, further confirming the heterogeneous karyotypic characteristics of this cell line. Similarly, scXRBSm and scXRBS data also produced consistent karyotypic results for GM12878 cells. Although the karyotypes produced by scRRBS (n = 17) data for HeLa cells were relatively similar to those of scMulOmi and scXRBSm, the overall copy state was lower. In contrast, the copy state of HeLa cells derived from scBS (n = 11) data was significantly lower than previous experimental results. Due to the low genomic coverage of scRRBS and the biased results of scBS, these two methods were excluded from subsequent epigenetic analysis. The results of ploidy calculations further supported this conclusion. It is worth noting that although they cannot be directly compared, the ploidy estimates produced by DIRECT in K562 cells are more dispersed than those of scXRBS, suggesting that scXRBS may provide more consistent ploidy estimates across different cell lines.

[0084] Example 5: scMulOmi analysis of epigenetic and transcriptomic landscapes of HeLa cell subgroups To further study the differences in epigenetics and transcriptomes between the two HeLa cell subgroups (g1 and g2) found in the karyotype analysis of Example 4, a comprehensive analysis of DNA methylation and gene expression data was performed. The method is shown in "2.4 Epigenetic analysis at the gene level" in Example 2, and the results are shown in Figure 5 .

[0085] First, we examined the association between changes in DNA methylation at the gene level and changes in gene expression between the G1 and G2 subgroups ( Figure 5 A). Using a bivariate Gaussian mixture model (p value < 10 -3 ), we identified significantly differentially methylated and expressed genes, represented by red dots in the figure ( Figure 5B) This analysis shows that by properly clustering the HeLa cells, a large number of genes can be found that are significantly different at both the DNA methylation and gene expression levels.

[0086] Next, we focused on the relationship between differential methylation of CpG in the promoter region and changes in corresponding gene expression ( Figure 5 B). This quadrant plot identified a total of 2237 differentially expressed genes based on promoter methylation levels and gene expression levels. Among the four quadrants, the second quadrant, representing low methylation and high expression levels, identified the most differentially methylated and expressed genes (715). To further explore the relationship between methylation and expression, we examined the trend of methylation levels around the promoter region and its association with gene expression levels ( Figure 5 C). The results clearly show 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 level is -0.111, indicating a negative correlation between the two. This finding supports the view that promoter methylation is associated with gene silencing, and higher methylation levels often correspond to lower gene expression.

[0087] Example 6: scMulOmi analysis of epigenetic differences between HeLa cell subgroups In this example, the epigenetic differences at the gene level between the two HeLa cell subgroups (g1 and g2) found in the karyotype analysis of Example 4 were further studied.

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

[0089] 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 effect on PFS in patients with cervical squamous cell carcinoma and cervical adenocarcinoma (log-rank test, p = 0.067, Figure 6 ).

[0090] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any way. Any simple modification, change and equivalent transformation made to the above embodiment based on the technical essence of the present invention still falls within the protection scope of the technical solution of the present invention.

Claims

1. A method for analyzing the epigenetic characteristics of CIN 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 restriction digestion, sulfite conversion, random primer amplification and library amplification in sequence, and the DNA library was obtained and sequenced; 2) After preprocessing, the DNA library sequencing data is aligned to the genomic reference sequence to identify the differentially methylated regions and perform cell karyotype analysis to divide tumor cells into different subgroups based on 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) Perform correlation analysis based on the results obtained in 2) and 3).

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; 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 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), the specific process of the cell karyotype analysis includes: 2.1) Use AneuFinder software to 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 based on this; use AneuFinder software to draw a chromosome overview map; 2.2) Calculate chromosome ploidy.

5. The method according to claim 4, 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.

6. The method according to claim 1, characterized in that The specific process of step 4) includes: 4.1) Based on the results obtained in steps 2) and 3), analyze the correlation between DNA methylation differences and gene expression differences among different subgroups; 4.2) Based on the results obtained in steps 2) and 3), 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.

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

8. 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 30 bp after trimming; 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.

9. 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.

10. 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