A method for deconvoluting DNA methylation sequencing data based on a topic model

Abstract

The present invention discloses a deconvolution method for DNA methylation sequencing data based on a topic model, belonging to the technical field of cell type deconvolution methods. This method constructs two Dirichlet distributions through the LDA algorithm, and then simulates the distribution relationships between samples and cell types, as well as between cell types and marker regions, to achieve reliable prediction of cell type composition. The purpose of constructing METRIC is to solve problems such as reduced deconvolution accuracy that may be brought about by highly sparse DNA methylation sequencing data, and to achieve efficient and reliable cell type deconvolution.

Description

Technical Field

[0001] The present invention belongs to the method of cell type deconvolution. Specifically, it relates to a deconvolution method for DNA methylation sequencing data based on a topic model. Background Art

[0002] Nucleosome-sized cell-free DNA (cfDNA) fragments released into the blood due to cell apoptosis and other reasons in the human body carry rich epigenetic information, including fragmentation, histone markers, and methylation patterns. Population-specific DNA methylation patterns can be observed in different cell types, tissues, and even cancers. Using cell-specific or tissue-specific DNA methylation patterns as biomarkers, reference maps applicable to deconvolution methods can be constructed. Deconvolution methods with reference matrices developed for DNA methylation data include MethAtlas, UXM, CelFiE, CelFEER, cfSort, etc. Due to the fact that a large number of cytosines in the genome are either not covered by sequencing or have a coverage lower than 3×, DNA methylation data usually has sparsity. The reason is that the high cost of whole genome bisulfite sequencing (WGBS) makes it very difficult to achieve sufficient depth, so it often contains a large number of CpG-deficient regions. Many published data have a sequencing depth of less than 30× and only two replicates. And reduced representation bisulfite sequencing (RRBS) lacks coverage of non-CpG dense regions. Similarly, text data in natural language processing often exhibits a high level of sparsity, and latent Dirichlet allocation (LDA) is one of the most popular topic modeling methods in text mining. It has good analytical performance for sparse data and is usually used for unsupervised topic discovery. LDA-based methods have been applied to the analysis of non-DNA methylation data. For example, cisTopic developed based on scATAC-seq, STRIDE and STdeconvolve based on spatial transcriptomics, GLDADec and GTM-decon based on RNA-seq. Although deconvolution methods for cell types applicable to DNA methylation data have been proposed, a deconvolution method friendly to sparse data for DNA methylation data has not been proposed yet. Summary of the Invention

[0003] In view of the content and related problems described in the background art, the present invention provides a deconvolution method for DNA methylation sequencing data based on a topic model. This method constructs two Dirichlet distributions through the LDA algorithm, and then simulates the distribution relationships between samples and cell types, as well as between cell types and marker regions, to achieve reliable prediction of cell type composition. The purpose of constructing METRIC is to solve problems such as reduced deconvolution accuracy that may be brought about by highly sparse DNA methylation sequencing data, and to achieve efficient and reliable cell type deconvolution.

[0004] To achieve the above object, the present invention provides the following technical solutions:

[0005] A deconvolution method for DNA methylation sequencing data based on a topic model, comprising the following steps:

[0006] Step 1: Preprocess the samples in the training set to obtain the methylation levels of the samples at each CpG site, and combine the cell type grouping information of the samples to identify the differentially methylated regions between different groups and the differential methylation scores of each region.

[0007] Step 2: Use the differentially methylated regions obtained in Step 1, combine their specificity scores to screen out cell type-specific marker regions, calculate the differential unmethylation index (DUI) of the marker regions, and construct a differential unmethylation index training matrix that can characterize cell type specificity and has the optimal number of markers.

[0008] Step 3: Regard the samples as documents, the cell types as topics, and the markers as words, and use the latent Dirichlet allocation to train the matrix obtained in Step 2 to optimize two probability distributions: the sample~cell type distribution and the cell type~marker distribution, and set the cell type labels for inferring the samples and markers.

[0009] Step 4: For the model constructed in Step 3, tune the hyperparameters α and β of the two Dirichlet distributions, and design a method for automatic assignment of cell type labels to train the cell components of the samples in the training set. After training, save the trained METRIC model and parameters.

[0010] Step 5: Use the trained METRIC model in Step 4 for prediction. Use the sequencing samples from tissues that are not included in the training set samples, mix them with leukocyte samples according to a known ratio to construct simulated test samples, extract the data of the cell type-specific marker regions identified in Step 2 and input them into the model, and the cell type component ratios consistent with the number of cell types included in the training set samples can be obtained.

[0011] Preferably, the specific implementation of step 1 adopted in the present invention is as follows: The total number of samples in the training set is N = {N1, N2, …, N j , …, N n}, N j (1 ≤ j ≤ n), and the number of cell type categories is K = {K1, K2, …, K i , …, K k}, K i (1 ≤ i ≤ k). Differentially methylated regions are identified and recognized through the wgbstools tool, and these regions have specific hypermethylation or specific hypomethylation relative to other groups.

[0012] Preferably, the specific implementation of step 2 is as follows: Considering comprehensively information such as the number of hypermethylated and hypomethylated regions with differential methylation and the differential scores of all groups, sort them within each group in descending order of the differential scores, and extract the top 25 differentially hypomethylated regions in each group as marker regions to obtain the marker region set T = {T1, T2, …, T t ,..T s}, T t (1 ≤ t ≤ s). Use wgbstools to analyze and obtain the proportion score U-score of hypomethylated reads of the training samples in each marker region. Then, for the specific marker region of each cell type, the U-score of the samples of the same cell type corresponding to it is multiplied by the differential score of the marker to obtain DUI, and then the differential non-methylation index training matrix is obtained.

[0013] Preferably, the specific implementation of step 3 is as follows: The cell types in each sample follow a multinomial distribution: , represents the cell type matrix. For the marker regions in each sample, it follows a multinomial distribution: , is the marker matrix. The distribution of cell types in sample j is obtained from the Dirichlet distribution with hyperparameter α: ~ Dirichlet(α), where the hyperparameter α controls the distribution of cell types in the sample. The distribution of marker regions in cell type i is obtained from the Dirichlet distribution with hyperparameter β: ~ Dirichlet(β), where the hyperparameter β controls the distribution of marker regions in each cell type. Then the LDA probability formula for a sample is:

[0014] .

[0015] Among them, M represents the total number of generation samples, C represents the total number of cell types, Represents the probability of the cell type distribution in sample j based on the hyperparameter α. Represents the probability of the marker distribution in cell type i based on the hyperparameter β. Represents the cell type distribution of sample j under which the cell type is assigned to sample j with a probability. Represents that in the cell type and the corresponding marker distribution under which the marker appears in sample j with a probability.

[0016] Preferably, the specific implementation of step four is as follows: The hyperparameter α in the model constructed in step three is set to auto, so that it adaptively selects the optimal value according to the data and other parameters. The hyperparameter β is a matrix of T × N, whose rows correspond to the set of cell type-specific marker regions obtained by screening, and whose columns correspond to the training set samples. The initial values in the matrix are all 0.0001. Then, for the marker region of each cell type, the values of the samples of the same cell type are set to 100000. The training outputs two result matrices. One is the distribution of the marker regions in the cell type, which describes the distribution and weights of the marker regions in the K cell types. By counting the specific cell type sources of the marker regions aggregated in each training-derived cluster, the specific cell type with the largest proportion is the label of this cluster. In the training, each cluster contains markers from the same specific cell type. The other is the distribution of the K cell types of the N samples, which describes the relative proportions of different cell types in each sample. Evaluate the two outputs in the training results, and use the entropy value to measure the degree of cell type singularity of each group. Comprehensively select the model with the smallest entropy value, and save the trained METRIC model and parameters to a file after training.

[0017] Preferably, the specific implementation of step five is as follows: Use the cell type sequencing samples from human tissues that are not included in the training set, and mix them with the white blood cell sequencing samples of normal people according to a set ratio. The mixing ratio is the true ratio, denoted as Actual = (0%, 0.01%, 0.03%, 0.1%, 0.3%, 1%, 3%, 10%, 15%, 20%, 25%, 30%). Each ratio is repeated three times to obtain a batch of simulated test samples. Then, for these simulated samples, extract the U-score of the marker region set regions identified and screened in step two, input the data into the model trained in step four, and the cell components in each simulated sample can be predicted. Extract the predicted scores corresponding to the true cell types of the simulated sample mixture from the prediction results, denoted as Predict. Pearson Correlation Coefficient (PCC), Coefficient of Determination (R²), and Root Mean Square Error (RMSE) are common indicators for evaluating the performance of statistical models. This method uses PCC to measure the strength of the linear relationship between two variables, R² to quantify the degree of explanation of the model for data variability, and RMSE to evaluate the difference between the predicted values and the actual values of the model. At the same time, in order to comprehensively evaluate the performance of the model on all indicators, this method defines a new indicator AccuracyScore(AS)=(rank(PCC)+rank(R 2 )+rank(RMSE)) / 3 to evaluate the accuracy of the model's deconvolution prediction. These indicators together provide a comprehensive evaluation of the model performance. By comparing the comprehensive performance of METRIC with several other methods on these four evaluation indicators, it shows that METRIC has good deconvolution performance.

[0018] A computer-readable storage medium, on which executable instructions are stored, and when the instructions are executed by a processor, the processor implements the above method.

[0019] An electronic device, comprising: one or more processors; a memory for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the above method.

[0020] Beneficial effects:

[0021] Compared with the existing deconvolution methods based on DNA methylation data, the advantage of the present invention is that it is a new deconvolution method based on the topic model. METRIC can simultaneously analyze the cell composition ratio and the weights of cell-specific DNA methylation markers, and has good interpretability and reliability. For data with high sparsity, it has good deconvolution performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a schematic diagram of the METRIC model structure.

[0023] Figure 2 It is a result diagram of the METRIC deconvolution simulation test sample - mammary epithelial cells.

[0024] Figure 3 It is a result diagram of the METRIC deconvolution simulation test sample - hepatocytes.

[0025] Figure 4 It is a result diagram of the METRIC deconvolution simulation test sample - villous trophoblast cells. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0026] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. However, the following embodiments are only for explaining the present invention, and the protection scope of the present invention should include all the contents of the claims. Moreover, through the description of the following embodiments, those skilled in the art can fully implement all the contents of the claims of the present invention.

[0027] Embodiment

[0028] The present invention will be described below in conjunction with the attached Figure 1-2 , Table 1, and examples. The examples herein are only for explaining the present invention and do not limit the present invention.

[0029] Figure 1The model structure diagram of the DNA methylation sequencing data deconvolution method METRIC based on the topic model is shown. A represents the methylation map, B represents feature extraction, and C represents the deconvolution model. First, DNA methylation sequencing data samples of various cell types from different human tissues are collected from public data sets, and then differential methylation analysis is performed to obtain cell type-specific marker regions. These markers are used to construct a marker pool to eliminate the effects of the order of markers on the genome. The top 25 markers in each group are screened out, and a differential non-methylation score matrix is constructed and input into the model for training. After the model is trained, the distribution of each sample will be output, including the composition ratio of different cell types and the clustering results of the marker region. Each cluster represents a cell type, in which each marker has a different weight. Then, simulated samples are generated according to known proportions through deconvolution to verify and compare the model performance. Finally, METRIC can be applied to new samples for deconvolution analysis.

[0030] Figure 2 , Figure 3 and Figure 4 The deconvolution results of METRIC applied to three groups of simulated test samples are shown. 575 RRBS sequencing samples were collected from 21 different public data sets, including 13 different cell types. After downloading the raw data of these samples, data quality control, alignment to the hg19 reference genome, methylation extraction and other analyses were performed, and then these samples were used for model training and testing. The following are the specific steps of the analysis:

[0031] Step 1: Simulate test sample generation: Among the collected RRBS test set samples, 6 samples of three cell types were randomly selected: mammary epithelial cells 2. Hepatocytes 3. Villous trophoblast 1. Mix them with a leukocyte RRBS sequencing sample respectively. The mixing ratios were set to 0%, 0.01%, 0.03%, 0.1%, 0.3%, 1%, 3%, 10%, 15%, 20%, 25%, and 30%, respectively, and recorded as Actual. Each mixing ratio was repeated three times, and a total of 6 12 3 = 216 simulated test samples.

[0032] Step 2. Training set data processing: Remove the 6 samples selected in Step 1. Use the remaining samples in the RRBS samples to construct a training set. Perform differential methylation analysis on these samples to identify 13 cell type-specific DNA methylation marker regions. After sorting according to the methylation difference level, select the top 25 markers for each cell type, obtaining a total of 325 markers to form a marker region set. Construct a differential non-methylation index matrix as the input for model training. In addition, extract the non-methylated read ratios of 216 test samples in the marker region set to form a test matrix for model prediction.

[0033] Step 3. Train METRIC: Input the training data matrix constructed in Step 2 into METRIC for training. At this time, the total number of training samples N = 569, the number of cell types K = 13, and the number of markers is T = 325. Set the hyperparameter α to auto during training, and the hyperparameter β is a matrix of T N. Set the value of the marker sample in the marker region to 100000, and other values to 0.00001. passes = 20, iterations = 200. After training, save the trained model and parameters.

[0034] Step 4. Use the trained model for prediction: Use the trained model to deconvolute the cell type composition of the simulated test samples to obtain the deconvolution results of each sample, including the composition ratio of each cell type. For the simulated samples composed of a mixture of mammary epithelial cells and white blood cells, extract the predicted values of the cell components of the mammary epithelial cells to form Predict. The same processing is applied to hepatocytes and syncytiotrophoblasts.

[0035] Step 5. METRIC model performance evaluation: Evaluate the deconvolution results of the simulated samples of 3 cell types. Take the average of the results of 2 mammary epithelial cells, and similarly take the average of the results of 3 simulated samples of hepatocytes to obtain the average Predict value. Calculate the three evaluation scores between Actual - Predict: Pearson correlation coefficient, coefficient of determination R 2 , root mean square error RMSE. Plot the line charts of Actual and Predict, and mark the values of the three coefficients on the graph. Figure 2 It shows that this method has a good deconvolution effect on the simulated test samples. The mean value of PCC is 0.99, and the mean value of R 2 is 0.72, and the mean value of RMSE is 0.05.

[0036] Table 1 Performance comparison table of METRIC and deconvolution methods of cfSort, UXM, and MethAtlas

[0037]

[0038] Table 1 shows the deconvolution comparison results of three other deconvolution methods on simulated test samples. Using three methods, cfSort, MethAtlas, and UXM, the model performance was evaluated under two different comparison scenarios. The following are the specific steps for deconvolving simulated test samples with the three methods:

[0039] (1) Deconvolving simulated test samples using cfSort: Convert the test samples into the.tfrecords input format required by cfSort, and use two models trained by cfSort, DNN1 (Deep Neural Network 1) and DNN2 (Deep Neural Network 2), to deconvolve the test samples, and calculate the correlation index scores between the actual values and the predicted values: PCC (Pearson correlation coefficient), R 2 (coefficient of determination), RMSE (root mean square error), and AS value (accuracy score).

[0040] (2) Deconvolving simulated test samples using MethAtlas: First, extract the DNA methylation levels of the training samples in the marked regions, group them by cell type and take the average to obtain the deconvolution reference matrix, and then use MethAtlas for deconvolution, and calculate the correlation index scores between the actual values and the predicted values: PCC, R 2 , RMSE, AS.

[0041] (3) Deconvolving simulated test samples using UXM: First, use the default settings of UXM to deconvolve the simulated test samples. Second, use all the training set samples to construct a new reference matrix that meets the requirements of UXM, and then deconvolve the simulated test samples again, and calculate the correlation index scores between the actual values and the predicted values: PCC, R 2 , RMSE, AS.

[0042] (4)Compare METRIC with cfSort, UXM, and MethAtlas according to two preset scenarios. Scenario 1: Compare with the two tools cfSort and UXM that perform deconvolution analysis using default settings. Since these two methods do not include trophoblasts in the reference matrix, there is no component ratio for this cell type. From the results in Table 1, it can be seen that METRIC has better deconvolution performance compared to the two methods cfSort and UXM. Scenario 2: Rebuild the reference matrix that conforms to UXM and MethAtlas using the training data and then perform deconvolution analysis on the simulated test samples. From Table 1, it can be seen that METRIC has intermediate deconvolution performance compared to the other two methods. The comprehensive results show that the deconvolution model trained using the training queue samples collected by the present invention to construct the training matrix can improve the deconvolution accuracy.

[0043] The above are only specific embodiments of the present application, enabling those skilled in the art to understand or implement the present application. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application will not be limited to these embodiments shown herein, but rather will be accorded the widest scope consistent with the principles and novel features claimed herein.

Claims

1. A deconvolution method for DNA methylation sequencing data based on a topic model, characterized in that Including the following steps: Step 1: Preprocess the samples in the training set, obtain the methylation levels of the samples at each CpG site, and combine the cell type grouping of the samples to identify differentially methylated regions and differential methylation scores; Step 2: Screen cell type-specific marker regions and calculate the differential non-methylation index of the marker regions; comprehensively consider the number and differential scores of differentially hypermethylated and differentially hypomethylated regions in all groups, sort them within each group from largest to smallest according to the differential scores, and extract the top 25 differentially hypomethylated regions within each group as marker regions to obtain the marker region set T={T1,T2,…,T t ,..T s}, Tt(1≤t≤s); use wgbstools to analyze and obtain the proportion score U-score of hypomethylated reads of the training samples in each marker region. Then, for the specific marker region of each cell type, the U-score of the samples of the same cell type corresponding to it is multiplied by the differential score of the marker to obtain the differential non-methylation index training matrix; Step 3: Treat the samples as documents, the cell types as topics, and the markers as words, use the matrix obtained in Step 2 to train with Latent Dirichlet Allocation, optimize two probability distributions: sample~cell type distribution and cell type~marker distribution, and set the cell type labels for inferring samples and markers; Step 4: Optimize the hyperparameters α and β, train the cell components of the training set samples, and after training is completed, save the METRIC model and parameters; the hyperparameter α in the model constructed in Step 3 is set to auto to adaptively select the optimal value according to the data and other parameters, and the hyperparameter β is a matrix of T × N, where the rows correspond to the set of cell type-specific marker regions obtained by screening, the columns correspond to the training set samples, and the initial values in the matrix are all 0.0001. Then, for the marker regions of each cell type, the values of the samples of the same cell type are set to 100000; two result matrices are output from the training. One is the distribution of marker regions in the cell type, which describes the distribution and weights of the marker regions in the K cell types. By statistically analyzing the specific cell type sources of the marker regions aggregated in each training-derived cluster, the specific cell type with the largest proportion is the label of this cluster. In the training, each cluster contains markers from the same specific cell type. The other is the distribution of the K cell types in the N samples, which describes the relative proportions of different cell types in each sample. Evaluate the two outputs in the training results, measure the cell type homogeneity of each group using the entropy value, and comprehensively select the model with the smallest entropy value. After training, save the trained METRIC model and parameters to a file; Step 5: Use the METRIC model trained in Step 4 for prediction. Use the sequencing samples from tissues that are not included in the training set samples, mix them with leukocyte samples in a known ratio to construct simulated test samples, extract the data of the cell type-specific marker regions identified in Step 2 and input them into the model, and the cell type component ratios consistent with the number of cell types included in the training set samples can be obtained.

2. The deconvolution method for DNA methylation sequencing data based on the topic model according to claim 1, wherein: In step one, the total number of samples in the training set is N = {N1, N2, …, N j , …, N n}, N j (1 ≤ j ≤ n), and the number of cell type categories is K = {K1, K2, …, K i , …, K k}, K i (1 ≤ i ≤ k). Differentially methylated regions are identified through the wgbstools tool, and these regions have specific hypermethylation or specific hypomethylation relative to other groups.

3. The deconvolution method for DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step three, the cell types in each sample follow a multinomial distribution: , denotes the cell type matrix, which follows a multinomial distribution for the labeled regions in each sample: , is the label matrix; the distribution of cell types in sample j is derived from a Dirichlet distribution with hyperparameter α: ∽Dirichlet(α), where the hyperparameter α controls the distribution of cell types in the sample; the distribution of labeled regions in cell type i is derived from a Dirichlet distribution with hyperparameter β: ∽Dirichlet(β), where the hyperparameter β controls the distribution of labeled regions in each cell type; then the LDA probability formula for a sample is: ； Where M represents the total number of surrogate samples, and C represents the total number of cell types. represents the probability of the cell type distribution in sample j based on the hyperparameter α; represents the probability of the marker distribution in cell type i based on the hyperparameter β; represents the cell type distribution of sample j under which the cell type is assigned to sample j; represents that in the cell type and the corresponding marker distribution under which the marker appears in sample j.

4. A deconvolution method for DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In Step 5, use the cell type sequencing samples from human tissues that are not included in the training set, mix them with the leukocyte sequencing samples of normal people in a set ratio, and the mixing ratio is the true ratio, denoted as Actual=(0%, 0.01%, 0.03%, 0.1%, 0.3%, 1%, 3%, 10%, 15%, 20%, 25%, 30%). Repeat each ratio three times to obtain a batch of simulated test samples. Then, for these simulated samples, extract the U-score of the marker region set identified and screened in Step 2, input the data into the model trained in Step 4, and the cell components in each simulated sample can be predicted. Extract the predicted scores corresponding to the true cell types of the simulated sample mixture from the prediction results, denoted as Predict; Use the Pearson correlation coefficient, coefficient of determination, root mean square error, and accuracy score as the indicators to evaluate the performance of the model.

5. A deconvolution method for DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step five, PCC is used to measure the strength of the linear relationship between two variables, R 2 is used to quantify the degree to which the model explains the data variability, and RMSE is used to evaluate the difference between the predicted values and the actual values of the model. Among them, R 2 is the coefficient of determination, and RMSE is the root mean square error.

6. The deconvolution method for DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step five, a new index, the accuracy score (AS), is defined as AccuracyScore(AS) = (rank(PCC) + rank(R 2 )) + rank(RMSE)) / 3 to evaluate the accuracy of the model's deconvolution prediction.

7. A computer-readable storage medium, characterized in that, It stores executable instructions, and when the instructions are executed by a processor, the processor implements the method according to any one of claims 1 to 6.

8. An electronic device, characterized in that, Including: One or more processors; A memory for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the method according to any one of claims 1 to 6.

Images