A method for establishing a cervical disease progression prediction model based on low-depth WGS

Abstract

The present invention discloses a method for establishing a cervical disease progression prediction model based on low-depth whole-genome sequencing technology, which uses cervical scraping or vaginal swabs to collect exfoliated cells, extract DNA, and perform whole-genome sequencing to obtain the original offline data of low-depth whole-genome sequencing of the DNA sample; the above-mentioned original data is subjected to standard quality control, and then compared with the human reference genome sequence and the repetitive sequences are marked; the genome instability index is calculated by the following formula:; In addition, the present invention also relates to a method for constructing a cervical disease prediction model; the present invention is conducive to improving the accuracy and sensitivity of detection, reducing the cost of detection, and improving the accessibility of cervical disease screening and prediction, while contributing to the early diagnosis of cervical disease and the formulation of treatment decisions, and improving the management and treatment effects of patients with cervical disease.

Description

Technical Field

[0001] The present invention belongs to the field of medical-assisted screening, and specifically involves using low-depth whole-genome sequencing technology to analyze the sequencing results, calculate the genomic instability index, establish a prediction model for cervical disease progression, and further predict cervical disease progression. Background Art

[0002] Cervical cancer is one of the most common cancers among women worldwide, and its incidence continues to rise in developing countries. The 2018 GLOBOCAN report showed that 85% of cervical cancer deaths occur in less developed countries. These low- and middle-income countries have 18 times the rate of cervical cancer deaths compared to developed countries. China and India account for over one-third of global cervical cancer cases. The five-year survival rate for patients with early-stage cervical cancer is 92%. However, cervical cancer often presents insidiously, often asymptomatic in its early stages, and is often discovered in its advanced stages, with limited treatment options and a poor prognosis. Data from the 2023, first edition of the NCCN Clinical Practice Guidelines for Cervical Cancer indicate that, due to the lack of access to screening for the majority of people in developing countries, approximately 275,000 people die from cervical cancer each year. As the only malignant tumor with a clearly defined cause that is both preventable and controllable, early detection of cervical cancer is crucial for its intervention and treatment.

[0003] Genome instability (GI) refers to the tendency or tendency for abnormal genetic changes to occur in the genome of a cell or individual. Genomic instability mainly includes chromosomal instability, point mutations and nucleotide alterations, copy number variation, etc. Increased genomic instability may lead to the development and progression of cancer cells. Many types of cancer are associated with genomic instability, including cervical cancer, breast cancer, and colon cancer. Copy number variation (CNV) generally refers to genomic structural variations in DNA fragments longer than 1Kb in the genome compared to the genomic reference sequence, including DNA deletions, insertions, inversions, translocations, and / or duplications at the microscopic and submicroscopic levels. Copy number variation is one of the important pathogenic factors of many human diseases (such as tumors, genetic diseases, nervous system and autoimmune diseases, etc.).

[0004] Research has shown that HPV is a key pathogen that promotes the development of cervical cancer. However, approximately 90% of HPV infections are cleared by the immune system within 1-2 years, and less than 2% of women will develop cervical cancer. Influencing factors include HPV infection and host genomic alterations, which can contribute to the progression of precancerous lesions into invasive cancer. Furthermore, multiple studies have shown that copy number gains or losses often inactivate tumor suppressor genes in the cancer genome or affect oncogene expression, playing a key role in tumor development. For example, Mittal et al. found that patients with persistent HPV infection and a high genomic copy number were at increased risk of progressing from low-grade to high-grade cervical lesions. Loharamtaweethong et al. found that increased PD-L1 gene copy number was associated with locally advanced cervical cancer. Compared to traditional cytology, HPV testing offers higher sensitivity and reproducibility, but HPV screening results are associated with a higher rate of false positives. Therefore, the use of a highly accurate and cost-effective method to detect copy number variations is crucial for cervical cancer screening and the prediction of cervical lesion grade.

[0005] At present, the main methods used to detect cell genome instability are: (1) banding karyotype analysis technology, also known as karyotype analysis; (2) molecular hybridization technology, including fluorescence in situ hybridization (FISH), chromosome microarray analysis technology (CMA), etc.; (3) polymerase chain reaction related technologies, including polymerase chain reaction (PCR), real-time quantitative PCR (qPCR), restriction fragment length polymorphism, microsatellite DNA polymorphism and denaturing high performance liquid chromatography (DHPLC), etc.; (4) gene sequencing technology, including first-generation sequencing (Sanger sequencing), second-generation sequencing (NGS), also known as high-throughput sequencing, among which second-generation sequencing includes whole genome sequencing (WGS), genome copy number variation sequencing (CNV-seq), whole exome sequencing (WES) and single cell sequencing.

[0006] Among them, karyotype analysis technology has always been considered the "gold standard" for diagnosing chromosomal aberrations, but its detection cycle is long and the resolution is low, and it cannot detect CNVs below 5Mb; CMA, a technology used for whole-genome CNVs detection, has high technical costs and low throughput. In addition, due to the limited coverage of the chip probes used in CMA, some pathogenic genomic copy number variations (pCNVs) may not be detected; single-cell sequencing relies on gene amplification, so there are still major challenges in sequencing bias and sufficient genome coverage.

[0007] The genomic copy number variation sequencing method involved in this invention is a high-throughput sequencing method that performs low-depth whole-genome sequencing (WGS) on DNA, aligns the sequencing results with the human reference genome sequence, and uses bioinformatics analysis to identify CNVs in the sample. It can detect chromosomal aneuploidy mosaicism as low as 5%. Compared with existing detection technologies, low-depth whole-genome copy number variation sequencing is cost-effective, simple to operate, has a wide detection range, and requires a low DNA sample size for testing.

[0008] Based on this, the present invention provides a highly accurate and cost-effective method for detecting genomic instability by equipping low-depth whole-genome sequencing technology with a suitable model, thereby predicting the progression of cervical disease. Summary of the Invention

[0009] To address the low sensitivity, poor reproducibility, and low accuracy of existing detection technologies, a method for detecting genomic instability for cervical disease assessment is provided using exfoliated cells collected via cervical scraping or vaginal swabs. Another object of the present invention is to provide a method for constructing a cervical disease progression prediction model.

[0010] To achieve the above object, the technical solution of the present invention is as follows:

[0011] The present invention provides a method for determining the genome instability index based on low-depth whole-genome sequencing technology, comprising the following steps:

[0012] Step 1: Collect cervical exfoliated cells through cervical scraping or vaginal swab;

[0013] Step 2: DNA samples were extracted from the collected cervical exfoliated cells and whole-genome sequencing was performed using the PE150 library construction method. The average sequencing depth was 0.5-1X; where X is the ratio of the total number of sequenced bases (bp) to the genome size (Genome);

[0014] Step 3: Obtain the raw data of low-depth whole-genome sequencing of DNA samples;

[0015] Step 4: The raw data were subjected to standard quality control (Q30>0.90), and then compared with the standard human reference genome sequence and repeated sequences were marked;

[0016] Step 5: Calculate the genomic instability index using the following formula:

[0017]

[0018] Wherein, n is 5547, i is a natural number, and GII represents the genomic instability index of the sample in step 2 above.

[0019] In the embodiments of the present invention, the sequencing method of the above-mentioned samples is low-depth whole-genome sequencing, which can reduce the cost of whole-genome sequencing on the one hand; on the other hand, compared with CMA gene chip detection, low-depth whole-genome sequencing has a larger detection range and higher accuracy.

[0020] In the embodiments of the present invention, the above-mentioned "standard human reference genome" refers to the standard human reference genome sequence of GRCh38 in the NCBI database. The sequence of the standard human reference genome can also be obtained from gene databases such as UCSC, Ensemble, GENCODE and RefSeq.

[0021] In an embodiment of the present invention, the method for aligning the sequencing results with the human reference genome sequence is BWA-MEM. The alignment can also be performed using methods known to professionals in the field, such methods include but are not limited to: Needleman-Wunsch global alignment algorithm, Smith-Waterman local alignment algorithm, MSA multiple sequence alignment algorithm, BLAST fast alignment algorithm, Pairwise alignment algorithm and FASTA alignment algorithm, etc.

[0022] In step 5, the genomic instability index is calculated, which specifically includes the following steps:

[0023] S1-1: Compare the result sequence of the original offline data with the standard human reference genome sequence;

[0024] S1-2: Divide the human reference genome into 5547 non-overlapping sliding windows, each with a size of 500 kb;

[0025] S1-3: Remove blacklist regions in non-overlapping continuous windows, such as chromosome centromeres and repeat regions, and remove X chromosomes and Y chromosomes at the same time;

[0026] S1-4: Use CNVkit software, a Python and command-line software toolkit for studying CNVs, to calculate the copy number of each region and then calculate the genomic instability index.

[0027] In an embodiment of the present invention, the multiple DNA samples are from multiple individuals, including cancer samples and non-cancerous samples, the DNA samples are from collected cervical exfoliated cells, and the DNA samples are derived from specially processed tissue samples.

[0028] In the embodiments of the present invention, blacklist regions refer to regions in the genome that are abnormal or have high signals in any next-generation sequencing experiment. Excluding these regions can better analyze functional genomic data.

[0029] In this embodiment of the present invention, low-depth whole-genome sequencing is used to detect copy number variations in chromosomes. The main principle is to perform whole-genome sequencing on sample DNA using next-generation sequencing technology, compare the sequencing results with the human reference genome, and use bioinformatics analysis to identify possible chromosomal abnormalities in the sample. The sequencing technology of this invention can be performed on the following platforms: MGISEQ-2000 (MGI), NextSeq 500 (Illumina), and Ion Torrent (Thermo Life).

[0030] The present invention also provides a method for constructing a cervical disease progression prediction model, comprising the following steps:

[0031] Step A: The raw data obtained using low-depth whole-genome sequencing was reviewed. The model used CNVkit software to analyze the Log2 ratio value of each window, removed features with a large number of missing values ​​in the window, and filled the missing values ​​with the mean. The reviewed data was then divided into two parts proportionally: a training set (80%) for model training and a test set (20%) for testing model performance.

[0032] Step B: Perform supervised machine learning on the training set using the "sklearn" library, use random forest to rank feature importance, use Lasso regression to further analyze and screen the features, adjust the parameters and structure, optimize the random forest algorithm, build a random forest model, and obtain prediction values;

[0033] Step C: Using the "matplotlib" Python library to draw an ROC curve based on the predicted values ​​determined above, the random forest prediction model is evaluated based on the test set data;

[0034] Step D: Using K-fold cross validation to improve the generalization ability of the training set and prevent the random forest model from overfitting due to excessive complexity;

[0035] Step E: The predictive ability of the random forest model was evaluated by sensitivity, specificity, AUC, and accuracy indicators to construct an accurate and robust prediction model for cervical disease progression assessment.

[0036] In the embodiment of the present invention, the supervised machine learning algorithms used include: random forest model and Lasso regression. In step B, further analyzing and screening the features using Lasso regression specifically includes the following steps:

[0037] S2-1: First, a window with statistical differences among three groups of samples was screened by rank sum test or variance analysis, where the three groups of samples were: CC (cervical cancer), HSIL (high-grade squamous intraepithelial lesion), and YZ_LSIL (inflammation / low-grade squamous intraepithelial lesion);

[0038] S2-2: Select the feature importance through the random forest and select the top b features;

[0039] S2-3: Take the intersection of the above features a and b to obtain c features;

[0040] S2-4: Based on the above c features, further analyze and screen the features through Lasso regression to obtain d, and add the GII calculated in claim 1 to finally obtain e features, e=d+GII.

[0041] In an embodiment of the present invention, based on the e features finally obtained above, the sensitivity is used as the vertical axis to represent the true positive rate, and the specificity is used as the horizontal axis to represent the false positive rate. An ROC curve is drawn to screen the optimal threshold value to obtain the corresponding sensitivity, specificity and accuracy in distinguishing the level of cervical disease.

[0042] In the embodiments of the present invention, Sklearn (also known as Scikit-learn) is a free software machine learning library for the Python programming language, which includes many commonly used machine learning algorithms, preprocessing techniques, model selection and evaluation tools, etc.; Lasso regression is a compression estimation method based on the idea of ​​reducing the set of variables.

[0043] In an embodiment of the present invention, the indicators for evaluating the predictive ability of the prediction model include sensitivity and specificity. Indicators that can be used to evaluate the predictive ability of the model include, but are not limited to, sensitivity, specificity, precision, and F1 value.

[0044] In the embodiment of the present invention, given the current clinical tendency to perform surgical treatment on patients with HSIL or above, the prediction results of CC and HSIL are combined into CC_HSIL to determine CC_HSIL and non-CC_HSIL.

[0045] In an embodiment of the present invention, the three-category model for predicting cervical disease is: CC, HSIL, and LSIL_YZ.

[0046] In the embodiments of the present invention, cervical precancerous lesions are clinically divided into three grades: CIN1 (very mild / mild dysplasia), CIN2 (moderate dysplasia), and CIN3 (severe dysplasia and carcinoma in situ). These grades are based on the proportion of cells with abnormal nuclear mitotic figures within the cervical squamous epithelium. HSIL, as described in the present invention, encompasses both CIN2 and CIN3 lesions. LSIL encompasses CIN1. YZ refers to inflammation of the cervix, including inflammation of the ectocervicovaginal portion (i.e., the lower end of the cervix leading to the vagina) and inflammation of the endocervical mucosa.

[0047] In an embodiment of the present invention, based on low-depth whole-genome sequencing technology, the genomic instability index is calculated and a prediction model for cervical disease progression is established, which reduces costs and is conducive to large-scale application.

[0048] The present invention also provides an application of a prediction model for cervical disease progression assessment established by the above method in a tool for predicting cervical disease progression.

[0049] The beneficial effects of the present invention are:

[0050] 1. To address the problems of low throughput and limited probe coverage in karyotyping and CMA, the present invention uses a prediction model constructed by combining the genomic instability index and low-depth whole-genome sequencing technology to improve the accuracy and sensitivity of detection overall.

[0051] 2. It is simple and convenient to operate by testing a small amount of exfoliated cell samples from cervical scraping or vaginal swab.

[0052] 3. The present invention uses low-cost sequencing technology to assess genomic instability, which reduces detection costs and improves the accessibility of cervical disease screening and prediction.

[0053] 4. Based on low-depth whole-genome sequencing technology, the present invention provides a method for predicting the progression of cervical disease. This method helps in the early diagnosis of cervical disease and the formulation of treatment decisions, improves the management and treatment effects of patients with cervical disease, and is of great significance. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 : This is a ROC curve diagram of the threshold value based on the training set feature screening in Example 1 of the present invention;

[0055] Figure 2 This is the ROC curve diagram of LSIL_YZ in Example 1 of the present invention. DETAILED DESCRIPTION

[0056] The present invention will be further described below with reference to the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0057] The definitions of the main abbreviations and key terms involved in this invention are as follows:

[0058] WGS: Whole Genome Sequencing;

[0059] GLOBOCAN: Global Cancer Statistics, a global cancer epidemiology database established by the World Health Organization;

[0060] CNV: Copy Number Variation;

[0061] GI: Genome Instability;

[0062] CC: Cervical Cancer;

[0063] YZ: YanZheng, inflammation;

[0064] HSIL: High Grade Squamous Intraepithelial Lesion, high grade squamous intraepithelial lesion;

[0065] LSIL: Low Grade Squamous Intraepithelial Lesion, low grade squamous intraepithelial lesion;

[0066] CIN: Cervical Intraepithelial Neoplasia;

[0067] AUC: Area Under the Curve, the area under the ROC curve and the coordinate axis;

[0068] FISH: Fluorescence in Situ Hybridization;

[0069] CMA: Chromosome Microarray Analysis, chromosome microarray analysis;

[0070] PCR: Polymerase Chain Reaction;

[0071] NGS: Next Generation Sequencing;

[0072] NCBI: National Center for Biotechnology Information, a public database of standard human reference genome sequences;

[0073] ROC curve: Receiver Operating Characteristic Curve, receiver operating characteristic curve;

[0074] True Positive Rate: True positive rate;

[0075] False Positive Rate: False positive rate.

[0076] In Example 1 of the present invention, the 474 samples were derived from the First Affiliated Hospital of University of Science and Technology of China (Anhui Provincial Hospital). Example 1

[0077] A method for constructing a cervical disease progression prediction model based on low-depth whole-genome sequencing.

[0078] 1. Calculation of genomic instability index:

[0079] 474 cervical cell samples were collected by cervical scraping, and DNA of the samples was extracted for whole genome library construction. The library construction mode was PE150, and the samples were sequenced at a low depth of 0.5-1X. The formula

[0080] Obtain GII.

[0081] 2. Prediction model establishment process:

[0082] Step 1: Low-depth whole-genome sequencing technology was used to sequence 474 samples to obtain the original offline data and review them. The features were 5547 windows. The Log2 ratio value of each window was obtained by analysis using CNVkit software. Features with a large number of missing values ​​were removed (for example, if a feature was missing in 450 samples, it was removed), and the missing values ​​were filled with the mean. Finally, 5069 windows remained, which were divided into training and test sets in an 8:2 ratio.

[0083] Step 2: Use the "sklearn" library to perform supervised machine learning on the training set, use random forest to rank the feature importance, and use Lasso regression to further analyze and screen the features. Finally, 25 features are obtained, and then the prediction values ​​are obtained based on the model. The ROC curve is drawn based on the above prediction values. Figure 1 .

[0084] Step 3: To prevent the prediction model from being too complex and overfitting, K-fold cross validation with CV=10 is used to improve the generalization ability of the above training set.

[0085] Step 4: Based on the 25 features, the predictive ability of the above model was evaluated through indicators such as sensitivity and specificity, and finally a three-classification model was established.

[0086] 3. Specific steps for feature screening in step 2 above, "further analysis and screening of the features using Lasso regression":

[0087] S2-1: First, 1657 windows with statistical differences were screened in CC, HSIL, and YZ_LSIL groups by rank sum test or variance analysis.

[0088] S2-2: Random forest is used to sort the features and select the top 1000 features.

[0089] S2-3: Take the intersection of the above 1657 features and 1000 features to obtain 544 features.

[0090] S2-4: The above 544 features were further analyzed and screened by Lasso regression to obtain 24 features, and the GII was added to finally obtain 25 features, see Table 1.

[0091] Table 1. 25 feature importance rankings

[0092] .

[0093] 4. Based on the above training set characteristics, the optimal threshold is selected. The sensitivity is used as the vertical axis to represent the true positive rate, and the specificity is used as the horizontal axis to represent the false positive rate. The ROC curve is plotted. Figure 1 .

[0094] See also Figure 1When the CC prediction probability threshold was set at 0.4, the model achieved the highest accuracy in distinguishing CC from non-CC, with a sensitivity and specificity of 79% and 84%, respectively, and an AUC of 0.87. When the HSIL prediction probability threshold was set at 0.34, the model achieved the highest accuracy in distinguishing HSIL from non-HSIL, with a sensitivity and specificity of 81% and 76%, respectively. When the LSIL_YZ prediction probability threshold was set at 0.38, the model achieved the highest accuracy in distinguishing LSIL_YZ from non-LSIL_YZ, with a sensitivity and specificity of 89% and 78%, respectively. Therefore, this model demonstrated the best predictive performance for LSIL_YZ and non-LSIL_YZ, with an AUC of 0.88.

[0095] Currently, the clinic tends to obtain the prediction results of CIN2+. In view of this, in this embodiment 1, lesions of CIN2 and CIN2+ are judged as non-LSIL_YZ, and those below CIN2 are judged as LSIL_YZ. First, the ROC curve is drawn based on the extracted sensitivity and specificity, see Figure 2 On the test set, the AUC of LSIL_YZ and non-LSIL_YZ is about 0.88. Then based on the threshold of 0.38, the accuracy of LSIL_YZ and non-LSIL_YZ judgment on the test set is 85.3%.

[0096] In the external test set, based on a threshold of 0.38, the accuracy of distinguishing LSIL_YZ and non-LSIL_YZ in the test set is 75%.

[0097] It should be noted that the above content merely illustrates the technical idea of ​​the present invention and cannot be used to limit the scope of protection of the present invention. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications all fall within the scope of protection of the claims of the present invention.

Claims

1. A method for constructing a prediction model for evaluating cervical disease progression, characterized in that: The following steps are involved: Step A: Divide the original offline data into two parts: 80% training set and 20% test set; Step B: Perform supervised machine learning on the training set using the "sklearn" library, use random forest to rank feature importance, use Lasso regression to further analyze and screen the features, adjust the parameters and structure, optimize the random forest algorithm, build a random forest model, and obtain prediction values; Step C: Use the "matplotlib" Python library to draw the ROC curve based on the predicted values ​​determined above, and evaluate the random forest prediction model based on the data of the test set; Step D: Using K-fold cross validation to improve the generalization ability of the training set and prevent the random forest model from overfitting due to excessive complexity; Step E: Evaluate the predictive ability of the model using sensitivity, specificity, AUC, and accuracy metrics to construct an accurate and robust prediction model for cervical disease progression assessment. The original offline data in step A is the original offline data obtained by low-depth whole-genome sequencing. The low-depth whole-genome sequencing method is a method for determining the genome instability index based on low-depth whole-genome sequencing technology, comprising the following steps: Step 1: Collect cervical exfoliated cells through cervical scraping or vaginal swab; Step 2: DNA samples were extracted from the collected cervical exfoliated cells and whole genome sequencing was performed using the PE150 library construction method, with an average sequencing depth of 0.5-1X; Step 3: Obtain the raw data of low-depth whole-genome sequencing of DNA samples; Step 4: The raw data is subjected to standard quality control, then aligned with the standard human reference genome sequence and repetitive sequences are marked; Step 5: Calculate the genomic instability index using the following formula: ； Where n is 5547, i is a natural number, and GII represents the genomic instability index of the sample in step 2 above; In step B, Lasso regression is used to further analyze and screen the features, specifically including: S2-1: First, a window with statistical differences among three groups of samples was screened by rank sum test or variance analysis, wherein the three groups of samples were: cervical cancer, high-grade squamous intraepithelial lesion, and inflammation / low-grade squamous intraepithelial lesion; S2-2: Select the feature importance through the random forest and select the top b features; S2-3: Take the intersection of the above features a and b to obtain c features; S2-4: Based on the above c features, further analyze and screen the features through Lasso regression to obtain d, and add the GII calculated in step 5 to finally obtain e features, e=d+GII.

2. The method for constructing a prediction model for cervical disease progression assessment according to claim 1, characterized in that: Step 5 specifically includes the following steps: S1-1: Compare the result sequence of the original offline data with the standard human reference genome sequence; S1-2: Divide the standard human reference genome into n non-overlapping continuous windows of 500 kb in order of arrangement; S1-3: Remove the blacklist regions in non-overlapping continuous windows and remove both X and Y chromosomes; S1-4: Use CNVkit to calculate the copy number of each region and then calculate the genomic instability index.

3. Use of a prediction model for cervical disease progression assessment established by the method according to any one of claims 1 to 2 in the preparation of a tool for predicting cervical disease progression assessment.

Images