DNA methylation sequencing data deconvolution method based on topic model

Abstract

The invention discloses a DNA methylation sequencing data deconvolution method based on a topic model, and belongs to the technical field of cell type deconvolution methods. According to the method, two Dirichlet distributions are constructed through the LDA algorithm, so that the distribution relation between the sample and the cell type and the distribution relation between the cell type and the marked area are simulated, and reliable cell type composition prediction is achieved. The purpose of constructing the METRIC is to solve the problem that the deconvolution precision is reduced possibly caused by DNA methylation sequencing data with relatively high sparsity, and efficient and reliable cell type deconvolution is realized.

Description

Technical Field

[0001] The present invention belongs to a cell type deconvolution method, and specifically relates to a DNA methylation sequencing data deconvolution method based on a topic model. Background Art

[0002] Nucleosome-sized cfDNA (cell-free DNA) fragments released into the blood due to apoptosis and other reasons in the human body carry rich epigenetic information, including fragmentation, histone marks 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 suitable for deconvolution methods can be constructed. Deconvolution methods with reference matrices developed for DNA methylation data include: MethAtlas, UXM, CelFiE, CelFEER, cfSort, etc. DNA methylation data are usually sparse because a large number of cytosines in the genome are either not covered by sequencing or have a coverage of less than 3×. 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 loss regions. Many published data have a sequencing depth of less than 30× and only two repeats. Reduced Representation Bisulfite Sequencing (RRBS) lacks coverage of non-CpG dense regions. Similarly, text data in natural language processing often exhibit 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 often 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, etc. Although deconvolution methods suitable for cell types of DNA methylation data have been proposed, deconvolution methods that are friendly to sparse data for DNA methylation data have not yet been proposed. Summary of the invention

[0003] In view of the content and related problems described in the background technology, the present invention provides a deconvolution method for DNA methylation sequencing data based on a topic model. The method constructs two Dirichlet distributions through the LDA algorithm, thereby simulating the distribution relationship 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 the problem of reduced deconvolution accuracy caused by highly sparse DNA methylation sequencing data, and to achieve efficient and reliable cell type deconvolution.

[0004] In order to achieve the above-mentioned purpose, the present invention provides the following technical solutions:

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

[0006] Step 1: Preprocess the samples in the training set to obtain the methylation level of the samples at each CpG site, combine the cell type grouping information of the samples, 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 and 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: Consider samples as documents, cell types as topics, and tags as words. Use latent Dirichlet allocation to train the matrix obtained in step 2 to optimize two probability distributions: sample-cell type distribution and cell type-tag distribution. Set the cell type labels of the inferred samples and tags.

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

[0010] Step 5. Use the METRIC model trained in step 4 to make predictions. Use sequencing samples from tissues that are not included in the training set samples and mix them with white blood cell samples in a known ratio to construct simulated test samples. Extract the cell type-specific marker region data identified in step 2 and input them into the model to obtain the cell type component ratio that is consistent with the number of cell types contained in the training set samples.

[0011] As a preferred embodiment, the specific implementation method of step 1 adopted by the present invention is: the total number of samples in the training set is N={N1,N2,…,N j ,…,N n},N j (1≤j≤n), the number of cell types is K={K1,K2,…,K i ,…,K k},K i (1≤i≤k), the differentially methylated regions were identified by the wgbstools tool. These regions were specifically hypermethylated or hypomethylated relative to other groups.

[0012] As a preferred embodiment, the specific implementation method of step 2 is: comprehensively consider the number of differentially methylated high-methylated and low-methylated regions, difference scores and other information of all groups, sort the groups from large to small according to the difference scores, extract the first 25 differentially low-methylated regions in each group as marker regions, and obtain the marker region set T={T1,T2,…,T t ,..T s},T t (1≤t≤s). wgbstools was used to analyze the proportion of low methylation reads in each marker region of the training sample, U-score. Then, for each cell type-specific marker region, the corresponding U-score of the sample of the same cell type was multiplied by the difference score of the marker to obtain DUI, and then the differential unmethylation index training matrix was obtained.

[0013] Preferably, the specific implementation of step 3 is: the cell types in each sample follow a multinomial distribution: , Represents a cell type matrix. For each labeled region in a sample, it follows a multinomial distribution: , is the labeling matrix. The distribution of cell types in sample j is given by 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 given by 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:

[0014] .

[0015] Among them, M represents the total number of 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 Next, cell type The probability assigned to sample j. Expressed in cell types and the corresponding marker distribution Next, mark The probability of occurrence in sample j.

[0016] As a preferred embodiment, the specific implementation of step 4 is: the hyperparameter α in the model constructed in step 3 is set to auto, so that it can adapt to the data and other parameters to select the optimal value. The hyperparameter β is a T The N matrix has rows corresponding to the cell type-specific marker region sets obtained by screening, and columns corresponding to the training set samples. The initial values ​​in the matrix are all 0.0001. Then, for each cell type marker region, the value of the same cell type sample is set to 100000. The training outputs two result matrices, one is the distribution of marker regions in cell types, which describes the distribution and weight of marker regions in K cell types. By counting the specific cell type sources of the marker regions clustered in each training cluster, the specific cell type with the largest proportion is the label of this cluster. In the training, each cluster is a marker from the same specific cell type. The other is the distribution of K cell types in N samples, which describes the relative proportion of different cell types in each sample. Evaluate the two outputs in the training results, use entropy to measure the degree of cell type homogeneity of each group, and comprehensively select the model with the smallest entropy value. After training, save the METRIC model and parameters of the training number to the file.

[0017] Preferably, the specific implementation of step 5 is: use a cell type sequencing sample from human tissue that is not included in the training set, mix it with a normal white blood cell sequencing sample according to a set ratio, the mixing ratio is the real ratio, recorded as Actual = (0%, 0.01%, 0.03%, 0.1%, 0.3%, 1%, 3%, 10%, 15%, 20%, 25%, 30%), and 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 obtained by identification and screening in step 2, input the data into the model trained in step 4, and predict the cell components in each simulated sample. Extract the prediction score corresponding to the real cell type mixed with the simulated sample from the prediction results, recorded as Predict. Pearson Correlation Coefficient (PCC), Coefficient of Determination (R²), Root Mean Square Error (RMSE) are commonly used 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 to which the model explains the variability of the data, and RMSE to evaluate the difference between the model's predicted value and the actual value. 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 is used to evaluate the accuracy of the model's deconvolution prediction. These indicators together provide a comprehensive evaluation of the model's performance. By comparing the comprehensive performance of METRIC with several other methods on these four evaluation indicators, it is shown that METRIC has good deconvolution performance.

[0018] A computer-readable storage medium stores executable instructions, which, when executed by a processor, enable the processor to implement the above method.

[0019] An electronic device comprises: one or more processors; and 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 present invention has the advantages that it is a new deconvolution method based on topic model, METRIC can simultaneously analyze the cell composition ratio and the weight 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 Schematic diagram of the METRIC model structure.

[0023] Figure 2 This is the result of the METRIC deconvolution simulation test sample - breast epithelial cells.

[0024] Figure 3 This is the result of the METRIC deconvolution simulation test sample - liver cells.

[0025] Figure 4 This is the result of the METRIC deconvolution simulation test sample - villous trophoblast. DETAILED DESCRIPTION

[0026] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. However, the following embodiments are limited to explaining the present invention, and the protection scope of the present invention should include the entire contents of the claims, and through the description of the following embodiments, those skilled in the art can fully implement the entire contents of the claims of the present invention.

[0027] Example

[0028] The following is combined with Figure 1-2 , Appendix 1, and examples are used to illustrate the present invention. The examples herein are only used to explain 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 sample to construct the training set, perform differential methylation analysis on these samples, identify 13 cell type-specific DNA methylation marker regions, sort by methylation difference level, screen the top 25 markers for each cell type, and obtain a total of 325 markers to form a marker region set, and construct a differential non-methylation index matrix as model training input. In addition, extract the non-methylated reads ratio of 216 test samples in the marker region set to form a test matrix for model prediction.

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

[0034] Step 4: Use the trained model to make predictions: Use the trained model to deconvolve the cell type composition of the simulated test sample to obtain the deconvolution results of each sample, including the composition ratio of each cell type. For the simulated sample composed of a mixture of mammary epithelial cells and white blood cells, extract the predicted value of the cell components of mammary epithelial cells to form Predict, and the same treatment is used for hepatocytes and villous trophoblasts.

[0035] Step 5. METRIC model performance evaluation: Evaluate the deconvolution results of the simulated samples of the three cell types, average the results of the two breast epithelial cells, and similarly average the results of the simulated samples of the three liver cells to obtain the average Predict value, and calculate the Pearson correlation coefficient and determination coefficient R between the three evaluation scores of Actual-Predict. 2 , root mean square error RMSE, draw Actual and Predict line graphs, and mark the values ​​of the three coefficients on the graph. Figure 2 This shows that the proposed method has a good deconvolution effect on the simulated test samples, with a PCC mean of 0.99 and R 2 The mean of is 0.72, and the mean of RMSE is 0.05.

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

[0037] Table 1 shows the deconvolution comparison results of the other three deconvolution methods on the simulated test samples. The model performance is evaluated in two different comparison scenarios using the three methods of cfSort, MethAtlas, and UXM. The following are the specific steps of the deconvolution simulation test samples of the three methods:

[0038] (1) Use cfSort deconvolution to simulate test samples: Convert the test samples into the .tfrecords input format required by cfSort, use the two models DNN1 (deep neural network 1) and DNN2 (deep neural network 2) trained by cfSort to deconvolve the test samples, and calculate the correlation index scores between the actual value and the predicted value: PCC (Pearson correlation coefficient), R 2 (coefficient of determination), RMSE (root mean square error), AS value (accuracy score).

[0039] (2) Use MethAtlas deconvolution to simulate test samples: First, extract the DNA methylation level of the training samples in the marked region, group them by cell type and take the average to obtain the deconvolution reference matrix, then use MethAtlas to perform deconvolution and calculate the correlation index scores between the actual value and the predicted value: PCC, R 2 , RMSE,AS.

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

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

[0042] The above description is only a specific implementation of the present application, so that those skilled in the art can understand or implement the present application. Various modifications to these embodiments will be apparent 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 the embodiments shown herein, but will conform to the widest range consistent with the principles and novel features applied for herein.

Claims

1. A DNA methylation sequencing data deconvolution method based on a topic model, characterized in that: The following steps are involved: Step 1: Preprocess the samples in the training set to obtain the methylation level of the samples at each CpG site, and identify the differentially methylated regions and differentially methylated scores based on the cell type grouping of the samples; Step 2: Screen cell type-specific marker regions, calculate differential unmethylation indexes of the marker regions, and construct a differential unmethylation index training matrix; Step 3: Consider the sample as a document, the cell type as a topic, and the tag as a word. Use latent Dirichlet allocation to train the matrix obtained in step 2 to optimize two probability distributions: sample-cell type distribution and cell type-tag distribution. Set the cell type labels of the inferred sample and tag. Step 4: Tune the hyperparameters α and β, train and predict the cell components of the training set samples, and save the METRIC model and parameters after training. Step 5. Use the METRIC model trained in step 4 to make predictions. Use sequencing samples from tissues that are not included in the training set samples and mix them with white blood cell samples in a known ratio to construct simulated test samples. Extract the cell type-specific marker region data identified in step 2 and input them into the model to obtain the cell type component ratio that is consistent with the number of cell types contained in the training set samples.

2. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 1, the total number of samples in the training set is N={N1,N2,…,N j ,…,N n },N j (1≤j≤n), the number of cell types is K={K1,K2,…,K i ,…,K k },K i (1≤i≤k), the differentially methylated regions were identified by the wgbstools tool. These regions were specifically hypermethylated or hypomethylated relative to other groups.

3. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 2, the number and difference scores of differentially hypermethylated and differentially hypomethylated regions in all groups are comprehensively considered, and the groups are sorted from large to small according to the difference scores. The first 25 differentially hypomethylated regions in each group are extracted as marker regions, and the marker region set T = {T1, T2, …, T t ,..T s }, Tt (1≤t≤s); wgbstools was used to analyze the proportion of low-methylation reads in each marker region of the training sample, U-score. Then, for each cell type-specific marker region, the corresponding U-score of the sample of the same cell type was multiplied by the difference score of the marker to obtain the differential unmethylated index training matrix.

4. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 3, the cell types in each sample follow a multinomial distribution: , represents the cell type matrix, which follows a multinomial distribution for the labeled regions in each sample: , is the labeling matrix; The distribution of cell types in sample j is given by 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 given by 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: ； Among them, M represents the total number of 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 distribution of markers in cell type i based on the hyperparameter β; Represents the cell type distribution of sample j Next, cell type The probability of being assigned to sample j; Expressed in cell types and the corresponding marker distribution Next, mark The probability of occurrence in sample j.

5. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 4, the hyperparameter α in the model constructed in step 3 is set to auto, so that it can adapt to the data and other parameters to select the optimal value, and the hyperparameter β is a T The matrix is ​​N, whose rows correspond to the cell type-specific marker region sets obtained by screening, and the columns correspond to the training set samples. The initial values ​​in the matrix are all 0.0001. Then, for each cell type marker region, the value of the same cell type sample is set to 100000; the training outputs two result matrices, one is the distribution of marker regions in cell types, which describes the distribution and weight of marker regions in K cell types. By counting the specific cell type sources of the marker regions clustered in each training 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 K cell types of N samples, which describes the relative proportion of different cell types in each sample. The two outputs in the training results are evaluated, and the entropy value is used to measure the degree of cell type uniformity of each group. The model with the smallest entropy value is comprehensively selected. After training, the trained METRIC model and parameters are saved to the file.

6. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step five, cell type sequencing samples from human tissues that are not included in the training set are used and mixed with normal white blood cell sequencing samples according to a set ratio. The mixing ratio is the actual 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, the U-score of the marker region set obtained by identification and screening in step two is extracted, and the data is input into the trained model in step four to predict the cell components in each simulated sample. The prediction score corresponding to the actual cell type mixed with the simulated sample is extracted from the prediction result, denoted as Predict. The Pearson correlation coefficient, determination coefficient, root mean square error and accuracy score are used as indicators to evaluate model performance.

7. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 5, PCC is used to measure the strength of the linear relationship between two variables, R 2 to quantify the extent to which the model explains the variability of the data, and RMSE to assess the difference between the model's predicted values ​​and the actual values, where R 2 is the coefficient of determination, and RMSE is the coefficient of determination.

8. The method for deconvolution of DNA methylation sequencing data based on a topic model according to claim 1, characterized in that: In step 5, a new indicator AccuracyScore(AS)=(rank(PCC)+rank(R 2 )+rank(RMSE)) / 3 is used to evaluate the accuracy of the model’s deconvolution prediction.

9. A computer-readable storage medium, characterized in that: Executable instructions are stored thereon, and when the instructions are executed by a processor, the processor implements the method according to any one of claims 1 to 8.

10. An electronic device, characterized in that: include: 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 8.

Images