Construction method and detection device of malignant lung nodule screening gene marker and screening model

Abstract

The application provides a malignant pulmonary nodule screening gene marker, a screening model construction method and a detection device. Through high-throughput sequencing results, 16bp tumor new short sequences of DNA differences between benign and malignant lung nodules with high risk in image performance are obtained, a model is constructed by using automatic machine learning, benign and malignant lung nodules with high risk in image performance are predicted, noninvasive accurate diagnosis of malignant pulmonary nodules is realized, and unnecessary benign pulmonary nodule resection surgery is reduced.

Description

Technical Field

[0001] This invention relates to the early screening of benign and malignant pulmonary nodules that appear to be of high risk on imaging, and belongs to the field of molecular biomedical technology. Background Technology

[0002] Lung cancer is one of the most frequently diagnosed cancers worldwide. High-risk groups include individuals over 65 years of age who possess one or more of the following risk factors: heavy smoking, a history of smoking, family history of lung cancer, previous chest radiation therapy, and exposure to carcinogens. Because early-stage lung cancer often lacks obvious symptoms, patients are generally diagnosed in the middle to late stages (stage III and IV). However, numerous studies have shown that lung cancer patients diagnosed earlier have higher survival rates. Patients diagnosed in stage I lung cancer have a 13 times higher five-year survival rate than those diagnosed in stage IV lung cancer. Therefore, early detection and diagnosis of lung cancer are crucial for improving the survival rate of lung cancer patients.

[0003] Low-dose computed tomography (LDCT) is currently the most common diagnostic method for lung tumors. Surgical resection of imaging-associated lung nodules can effectively reduce lung cancer mortality by 20%-39%. However, approximately 15%-35% of lung nodules initially identified as high-risk on LDCT are later confirmed to be pathologically harmless after surgical resection. Therefore, imaging examinations have certain limitations. Diagnosing malignant lung tumors solely based on imaging results increases the risk of unnecessary surgeries, complications, and increased medical expenses for patients. Therefore, the determination of benign or malignant lung nodules in individuals identified as high-risk for lung cancer solely based on imaging is particularly important. Summary of the Invention

[0004] The technical problem to be solved by this invention is that the existing technology lacks non-invasive detection methods in the process of diagnosing benign and malignant lung cancer nodules, which leads to unnecessary surgery and increases the burden on patients.

[0005] The technical solution of this patent provides a method for performing WGS sequencing on plasma sample cfDNA. By analyzing the high-throughput sequencing results, fragment information is obtained, and differential DNA fragment length ratios (fragment size ratio), breakpoint motifs (5' end breakpoint sequence proportion), 1 MB-bin copy number variants, 16 bp neomers, and nucleosome coverage patterns are analyzed. Four algorithms—elastic net logistic regression (GLM), extreme gradient boosting (XGBoost), random forest, and neural network—are used to construct a multi-feature, multi-algorithm integrated model through automated machine learning, achieving non-invasive and accurate diagnosis of malignant lung nodules.

[0006] The specific technical solution is as follows: Application of genetic markers in the preparation of screening reagents for malignant pulmonary nodules; The genetic markers mentioned include: The first biomarker: the proportion of short reads and the proportion of long reads in different windows of the cfDNA fragment aligned to the reference genome; The second biomarker: the proportion of m bases from the 5' end of the reference genome that different types of cfDNA fragments are aligned to; The third biomarker: the copy number in different windows on the chromosome in WGS data; Fourth biomarker: Percentage of new short sequences in tumors; Fifth marker: nucleosome coverage pattern.

[0007] The first biomarker is obtained through the following steps: aligning the read data results to the reference genome, dividing the reference genome into multiple windows, and obtaining the proportion of short reads and ultra-long reads within each window.

[0008] The second marker is obtained through the following steps: taking the m bases at the 5' end of the read data as a set of base fragments, and obtaining the proportion of each base fragment in the total fragments.

[0009] The third biomarker is obtained through the following steps: the reference genome is divided into multiple windows, and copy number data of WGS data in different windows on chromosomes 1-22 are obtained respectively.

[0010] The fourth marker is obtained through the following steps: A set A of short sequences of length 16bp is generated by exhaustive search; all short sequences of length 16bp are exhaustively searched in the human reference gene sequence, and the data of set B is removed from set A and defined as invalid. WGS sequencing results of samples from different cancer types were obtained from a cancer database, and multiple base substitution mutations were extracted. Based on the location of the base substitutions, a set C of short invalid sequences containing these base substitutions was found from invalid sequences. Obtain base substitution mutations with a frequency greater than 0.01 in East Asian populations; based on the position of the base substitution, find a set D of invalid short sequences containing these base substitutions from invalid sequences; exclude invalid sequences from set D from set C and define them as new short sequences; Count the number of samples in the sample that can read any new short sequence. Then, for each new short sequence, search for the number of samples containing these new short sequences. Calculate the ratio of the number of samples with each new short sequence to the total number of samples that can read any new short sequence.

[0011] The fifth marker is obtained through the following steps: Transcription factors were obtained from the GTRD database, and those transcription factors that did not have known transcription sites in the CIS-BP database were excluded. The range of -5kb to +5kb near the transcription site of the obtained transcription factor is used as a window to obtain fragments of 100-220bp in length that can be aligned to these windows. The read data in the window are then subjected to GC correction and sequencing depth smoothing to obtain the coverage pattern curve of each transcription factor. For each transcription factor, the following three features are obtained, which together serve as the nucleosome coverage pattern: 1) For all transcription sites of a transcription factor, calculate the average depth from the upper 1kb to the lower 1kb of the transcription site; 2) For the obtained coverage pattern curve, obtain the amplitude value of the curve trough as the center depth of the transcription factor; 3) Perform a fast Fourier transform on the obtained coverage pattern curve to obtain the amplitude value of the highest point of the nucleosome amplitude signal.

[0012] The method for constructing a screening model for malignant pulmonary nodules includes the following steps: Step 1: Extract and sequence cfDNA from samples in the positive and control groups to obtain read data; Step 2: Align the read data results to the reference genome, divide the reference genome into multiple windows, and obtain the proportion of short reads and long reads in each window as the first feature set; Step 3: Take the m bases at the 5' end of the read segment data as a base fragment set, and obtain the proportion of each base fragment in the whole fragment as the second feature set; Step 4: Divide the reference genome into multiple windows and obtain the copy number data within each window as the third feature set; Step 5: The ratio of the number of samples that read the new 16bp short sequence to the total number of samples that can read any new short sequence is used as the fourth feature set. Step 6: The nucleosome coverage pattern features of the selected transcription factors are used as the fifth feature set; Step 7: Use the first, second, third, fourth and fifth feature sets together as initial feature values, input them into the classification model as model feature vectors, and use the benign or malignant nature of lung nodules as the output value to train the model and obtain the early screening model.

[0013] Step 3 includes: Step 3-1: Divide the reference genome into multiple windows and obtain the number of long reads and short reads within each window. Step 3-2: Standardize the number of short segments and the number of long segments in all windows of Step 3-1, and use the ratio of the standardized number of short segments to the number of long segments as the first feature value.

[0014] In step 3-1, the window size is 5Mb, and a total of 541 windows are divided.

[0015] The short read segment refers to a length of 100-150bp, and the long read segment is 151-220bp.

[0016] In step 3, m is 4.

[0017] In step 4, the window size is 1Mb, and a total of 2475 windows are created.

[0018] In step 5, the steps for obtaining the fourth feature set are as follows: Step 5-1: Generate a set A of short sequences with a length of 16bp using an exhaustive search method; exhaustively search for all short sequences with a length of 16bp in the human reference gene sequence, and remove the data from set B from set A, defining them as invalid substrings; Step 5-2: Obtain WGS sequencing results of samples from different cancer types from the cancer database, and extract the base substitution mutations that occur multiple times; Based on the position of the base substitution, find the set C of short invalid sequences containing these base substitutions from the invalid sequences; Step 5-3: Obtain base substitution mutations with a frequency greater than 0.01 in East Asian populations; based on the position of the base substitution, find the set D of invalid short sequences containing these base substitutions from the invalid sequences; exclude the invalid sequences in set D from set C and define them as new short sequences; Step 5-4: Count the number of samples in the sample that can read any new short sequence. Then, for each new short sequence, search for the number of samples containing these new short sequences. Use the ratio of the number of samples of each new short sequence to the total number of samples that can read any new short sequence as the fourth feature set of the model.

[0019] The cancer database mentioned is the PCAWG database.

[0020] Different types of cancer refer to colorectal cancer, lung cancer, breast cancer, stomach cancer, prostate cancer, and liver cancer.

[0021] Base substitution mutations in East Asian populations were obtained using the gnomAD database.

[0022] Step 6 includes: Step 6-1: Obtain transcription factors from the GTRD database and exclude transcription factors that do not have known transcription sites in the CIS-BP database. Step 6-2: Using the range of -5kb to +5kb near the transcription site of the transcription factor obtained in Step 6-1 as a window, obtain fragments of length 100-220bp that can be aligned to these windows. Perform GC correction and sequencing depth smoothing on the read data in the window in sequence to obtain the coverage pattern curve of each transcription factor. Step 6-3: For each transcription factor, the following three features are obtained, which together serve as nucleosome coverage pattern features: 1) For all transcription sites of a transcription factor, calculate the average depth from the upper 1kb to the lower 1kb of these transcription sites; 2) For the obtained coverage pattern curve, obtain the amplitude value of the curve trough as the center depth of the transcription factor; 3) Perform a fast Fourier transform on the obtained coverage pattern curve to obtain the amplitude value of the highest point of the nucleosome amplitude signal.

[0023] In step 7, the classification model includes the following steps: Step 7-1: Input the first, second, third, fourth and fifth feature sets into different classifier models respectively, train the models, and obtain one or more optimal classifier models for the first, second, third, fourth and fifth feature sets respectively; Step 7-2: Perform secondary ensemble training on the optimal classifier models of the first, second, third, fourth, and fifth feature sets obtained in Step 7-1 to construct an ensemble classifier model.

[0024] The different classifier models mentioned are selected from elastic net logistic regression (GLM), extreme gradient boosting (XGBoost), random forest, and deep learning neural network (NN).

[0025] In the second set training, a generalized linear model (GLM), an extreme gradient boosting Xgboost model, or a deep learning regression model are used.

[0026] A device for detecting malignant pulmonary nodules, comprising: The sequencing module is used to extract and sequence cfDNA from samples in the positive and control groups to obtain read data and WGS sequencing data. The first feature acquisition module is used to align the read data results to the reference genome, divide the reference genome into multiple windows, and obtain the proportion of short reads and ultra-long reads within each window range, which are used as the first feature set. The second feature acquisition module is used to take the m base data at the 5' end of the read segment data as a base fragment set, and obtain the proportion of each base fragment in the whole fragment as the second feature set. The third feature acquisition module is used to divide the reference genome into multiple windows and obtain the copy number data of WGS data in different windows on the chromosome, as the third feature set; The fourth feature acquisition module uses the ratio of the number of samples that read a new 16bp short sequence to the total number of samples that can read any new short sequence as the fourth feature set. The fifth feature acquisition module analyzes the nucleosome coverage pattern features of the selected transcription factors as the fifth feature set; The prediction module is used to train the classification model by taking the first, second, third, fourth and fifth feature sets as initial feature values ​​and inputting them as model feature vectors. The benign or malignant nature of the lung nodules is used as the output value. Attached Figure Description

[0027] Figure 1 This is a schematic diagram of the model building process; Figure 2 This is a difference distribution chart of the feature with the largest contribution value among all features; Figure 3 The AUC curves are for models using individual features on the training set and models using all features. Figure 4 It is the AUC curve of the model using all features on the validation set; Figure 5 This is a distribution of the predicted scores of the classifiers after all models are set on the training set; Figure 6 It is a distribution of the predicted scores of the classifiers after the entire model set is applied to the validation set. Detailed Implementation

[0028] The calculation method in this invention is described in detail below: This invention first requires the extraction, library construction, and sequencing of cfDNA from blood samples. The extraction and library construction methods are not particularly limited and can be adapted from existing extraction methods. The sequencing process can utilize existing sequencing technologies to obtain the base information of the cfDNA.

[0029] The purpose of the model in this patent is to distinguish between benign lung nodules and malignant lung nodules. The samples are categorized: patients with lung nodules identified as high-risk on LDCT imaging are classified as a control group, while those identified as having benign lung nodules based on subsequent postoperative pathology are classified as a positive group.

[0030] The dataset used in the model building process of this invention is as follows: Methods for extraction and sequencing of plasma cfDNA samples: Before LDCT imaging diagnosis, a liquid biopsy was performed on the patient. A 10ml whole blood sample was collected from the patient using a purple blood collection tube (EDTA anticoagulant tube), and plasma was immediately centrifuged (within 2 hours). The plasma was then frozen at -80°C and transferred to the laboratory for analysis. Upon arrival at the laboratory, ctDNA was extracted from the plasma sample using the QIAGEN plasma DNA extraction kit according to the manufacturer's instructions. After library construction of the collected cfDNA samples, WGS 5x sequencing was performed. After obtaining the sequencing data, the data was aligned to the human reference genome to obtain the base information of the corresponding reads.

[0031] The model establishment process of this patent is mainly as follows: Step 2: Extract and sequence cfDNA from samples in the positive and control groups to obtain reads; Step 3: Align the read data results to the reference genome to obtain the number of reads in different length intervals within different window ranges on the reference genome. The ratio of the number of reads of different lengths is used as the first feature value. Step 4: Align the read data results to the reference genome to obtain the position of the 5' end of the read on the reference genome; obtain the sequence data of m bp bases upstream and downstream of the position as a set of base fragments; use the proportion of each base fragment in the total fragments as the second feature value. Step 5: Divide the reference genome into multiple windows and obtain the copy number data within each window range as the third feature value; Step 6: The ratio of the readings of the 16bp new short sequence to the total readings is used as the fourth feature; Step 7: Analyze the nucleosome coverage pattern of the selected transcription factors as the fifth feature; Step 8: Input the model feature vectors of the positive group and the control group into the first layer model, and select the model with the best AUC for each feature. A total of 25 feature models are selected in the first layer model. Step 9: Input the 25 models selected in Step 9 into the second layer of the integrated model, and output the top three integrated models in terms of AUC. Take the average of the predicted probability results output by these three integrated models as the final judgment result. There are five types of feature values ​​in this patent, which are described in detail below: 1. DNA Fragment Size Ratios (FSR) The ratio of DNA fragment sizes reflects the distribution characteristics of cfDNA read lengths in benign and malignant tumors. Machine learning is used to build predictive models based on the ratio of short to long DNA fragments, specifically for benign lung nodules and malignant lung nodules. The cfDNA read length data were obtained as follows: In the aligned BAM (Bio-Aligned Unit-Analog Transformer), the quality, length, and alignment position information of each read were recorded. The human reference genome was the hg19 sequence provided by the University of California, Santa Cruz (UCSC). The human reference genome was divided into 541 windows of 5 Mb each. The number of short reads (100 bp-151 bp) and long reads (151 bp-220 bp) in each window were counted. Based on the statistical results of the number of reads in all windows, the number of each type of read was standardized, i.e., standardized value = (original value – mean) / standard deviation. This resulted in 541 sets of read length data.

[0032] 2. Percentage of sequence reads at the 5' end breakpoint (Breakpoint Motifs, BPM) The human reference genome is a DNA double helix structure linked by complementary base pairing hydrogen bonds. During normal aging and cancer progression, the pH of the cellular environment changes, disrupting these complementary hydrogen bonds and causing breaks. Because the base sequences at the break points differ, the proportion of sequences containing information from different break points also varies. Data collection method: The alignment BAM records the basic information and alignment location of each read. The 5' end of each read is located 4 bp to the left and right of the break point in the human reference genome sequence coordinates. The number of 8 bp sequences at each break point (a total of 4**8=65536 types) is counted, thus calculating the proportion of reads at each of the 65536 break points. For example, the proportion of the AAAAAAAA read = the number of AAAAAAAA reads / the total number of reads at all break points.

[0033] 3. 1 Mb window copy number variation (CNV) Copy number changes are highly correlated with individual cancers. Although it is possible to distinguish between cancers by detecting copy number changes in some cancer-related genes or specific genomic regions, there are still other rare or unknown genes or regions that can provide information on potential copy number changes. Data Collection Method: First, WGS data from 30 healthy individuals were collected. Reference genes on chromosomes 1-22 were divided into non-overlapping windows of 1 Mb length. The read depth within each window was calculated using bedtools coverage for each sample. Corrections were made based on the GC content and average alignment capability records (UCSC BigWig file) for each window. The median depth of 30 healthy individuals within each window was used as a representative, resulting in a population control baseline of 2475 window read depths. For each test sample, 2475 individual read depths were also obtained. Using a Hidden Markov Model (HMM) and the population control baseline depth for each window, the logarithm of copy number change for each window was constructed, i.e., log2 (corrected and normalized depth of the test sample / corrected and normalized depth of the population baseline), thus obtaining the copy number change information for each test sample.

[0034] 4.16bp novel short tumor sequences (16bp Neomers, NEO) Nullomers are short DNA sequences that do not exist in the human genome. 16bp novel tumor sequences are a subset of nullomers, specifically referring to nullomers of 16bp length that are repeatedly found in the genome of tumor tissues but do not appear in the human genome.

[0035] In this patent, the feature value is obtained in the following way: First, an exhaustive search is used to generate a set A of all possible short sequences of length 16 bp. Then, in the human reference gene sequence (hg19 version), using a 1 bp sliding window, an exhaustive search algorithm is used to exhaustively search for all sets B of short sequences of length 16 bp and their occurrence counts. Short sequences of length 16 bp that appear only in set A are defined as invalid sequences.

[0036] This patent focuses on analyzing the PCAWG database (https: / / dcc.icgc.org / pcawg) to obtain the WGS mutation results of 2577 patients with six different types of cancer (colorectal cancer, lung cancer, breast cancer, stomach cancer, prostate cancer, and liver cancer). From these, 977 frequently occurring base substitutions (appearing at least twice) were extracted. Based on the location of the base substitution, an exhaustive method was used to extract the set C of all possible invalid subsequences containing that base substitution.

[0037] Using the gnomAD (https: / / gnomad.broadinstitute.org / ) database, mutation sites with a frequency greater than 0.01 in the East Asian population were collected. Based on the location of these mutation sites, a set D of short sequences containing these sites was identified from the null subsets, serving as a common null subset for the East Asian population. Null sequences from set D that appeared in set C were then excluded, resulting in 4616 novel 16bp tumor-related short sequences. For these 4616 novel short sequences, the number of samples in the dataset that could read any of the 4616 novel short sequences was first counted. Then, for each novel short sequence, the number of samples containing these sequences was searched. The ratio of the number of samples containing each novel short sequence to the total number of samples that could read any of the 4616 novel short sequences (a total of 4616 ratio values) was used as the fourth feature of the model.

[0038] 5. Nucleosome coverage patterns (NCP) Transcription factors were selected from the GTRD database (https: / / gtrd.biouml.org / #!) (v21.12), excluding those with known transcription sites not found in the CIS-BP database (http: / / cisbp.ccbr.utoronto.ca / ) (v2.00), resulting in 334 transcription factors with more than 10,000 highly matched sites.

[0039] For the transcription factors obtained above, a window of -5kb to +5kb near the transcription site of these target transcription factors was used to obtain fragments of length 100-220bp that could be aligned to these windows. For these fragments, GC correction was performed, and a Savitzky-Golay filter with a polynomial power of 3 was used to smooth the curve for sequencing depth, resulting in the final coverage pattern curve for each transcription factor.

[0040] After obtaining the above coverage pattern curves, three features were extracted for each transcription factor: 1) For all transcription sites of a transcription factor, calculate the average depth from the upper 1kb to the lower 1kb of these transcription sites; 2) The depth of the center of transcription factors; 3) Perform a fast Fourier transform on the coverage pattern curve obtained above to obtain the amplitude value of the highest point of the nucleosome amplitude signal.

[0041] These three features are used together as the feature values ​​of the nucleosome coverage pattern.

[0042] Through the above data acquisition, initial data vectors for these five data categories can be obtained. Next, corresponding calculation methods are designed. In this patent, a conventional classifier algorithm can be used to classify the data using the aforementioned feature values. These feature values ​​are input into the classifier, with the probability of malignant pulmonary nodules as the output. Alternatively, known classifiers can be optimized and adjusted. The classifier models used for optimization in this patent include the following four types. Furthermore, during classifier optimization, sub-models with different model parameters under the same model are simultaneously generated for selection. The four main models include: 1. Elastic net logistic regression (GLM) model Elastic network regression models are a common class of algorithms in machine learning. They are models that fit generalized linear simulations by penalizing maximum likelihood, combining L2 regularization from ridge regression and L1 regularization from LASSO regression. The regularization path is calculated using the lasso or elastic network penalty value at the grid of regularization parameter λ, thus solving the overfitting problem in regression. The hyperparameter alpha controls the distribution of regularization L1 and L2.

[0043] 2. Extreme gradient boosting (XGBoost) It is an optimization algorithm for additive models based on the ensemble concept of Gradient Boosting Decision Tree (GBDT). It uses second-order Taylor expansion to optimize the loss function and improve computational accuracy, uses regularization terms to simplify the model and avoid overfitting, and employs a Blocks storage structure for parallel computation.

[0044] The XGboost learner used in this patent is a tree model. Increasing the tree depth increases the tree's complexity, which can lead to better model training but may also cause overfitting. The hyperparameter `max_depth` controls the maximum depth of the tree, and the hyperparameter `min_rows` controls the minimum number of samples per leaf node.

[0045] 3. Random Forest Random forest is a powerful classification and regression tool for high-dimensionality and multicollinearity scenarios. Given a dataset, random forest can randomly extract information to generate a tree of decision trees to aid in classification or regression. It performs node splitting attribute sampling, repeatedly extracting attributes until no further splitting is possible. Finally, it combines all the splitting attribute results to obtain the final prediction. The complexity of the tree is controlled by hyperparameters such as `max_depth` and `min_rows`.

[0046] 4. Deep Learning Neural Network Model (NN) A neural network consists of inputs, weights, biases or thresholds, and outputs. If the output of any single node exceeds a specified threshold, that node is activated, and the data is sent to the next layer of the network. Each node in the input layer performs point-to-point computation with each node in the hidden layer, using weighted summation and activation methods. Each value calculated by the hidden layers is then used to perform computation with the output layer in the same way. It has advantages such as high classification accuracy, strong parallel distributed processing capabilities, and strong distributed storage and learning capabilities.

[0047] The deep learning neural network model used in this patent is a multi-layer feedforward neural network. The feedforward network structure consists of an input layer at the beginning, hidden layers in the middle, and an output layer at the end. Multiple, complex hidden layers can be present. Network signals are transmitted unidirectionally from the input layer to the output layer, and neurons in each layer are only connected to neurons in the layer preceding them to receive information from the previous layer. The deep learning neural network uses stochastic gradient descent for optimization and training. It uses data from the training dataset to estimate the error gradient of the model's current state and then uses the error to update the model's weights. The amount of weight update during training is the learning rate (epsilon), which controls the speed at which the model adapts to the problem. The learning rate decay exponent (rho) also controls the model's learning rate; both are important configurable hyperparameters in neural network training.

[0048] In addition, this patent also employs a Random Grid Search Parameters algorithm to optimize the model.

[0049] Random search is a common method for hyperparameter optimization in machine learning. Random search involves randomly sampling parameter values ​​within a specific range of model parameters and selecting the optimal combination of parameters from multiple samples. Instead of trying all possible combinations, this method selects a specific number of random combinations of a random value for each hyperparameter. Compared to manual tuning and web search for model hyperparameter tuning, random search can achieve good results with fewer searches, providing a more efficient solution (especially when dealing with a large number of parameters).

[0050] In the implementation of the four algorithms in this patent, the following specific algorithm types were used (or please specify which algorithm toolkit was called). The hyperparameters of the four algorithms used in this patent are shown in the table below: After obtaining the initial data from 247 patients with malignant pulmonary nodules and 60 patients with benign pulmonary nodules, the fragment size ratio (FSR) statistics were used as input values ​​(each sample's input vector included 541 feature values ​​composed of fragment size ratios). Four classification models were used to classify malignant and benign pulmonary nodule samples. During the selection process, the parameters and structure of these four models were varied through random hyperparameter search, serving as sub-models for data training and model suggestions. Finally, the feature was selected. The five optimal sub-models were selected using the AUC curve of the model's training set as the classification performance indicator. Similarly, after collecting the proportion of DNA fragment 5' end break sequence reads from patients with benign and malignant pulmonary nodules, the proportion of DNA fragment 5' end break sequence reads (65,536 types) was used as input. Four classification models were used to classify malignant and benign pulmonary nodule samples based on five features, and the five optimal sub-models for each feature were selected (the specific model optimization and hyperparameter adjustment process is the same as above). Similarly, copy number variation (2,475 types), new short sequences (4,616 types), and nucleosome coverage patterns (1,002 types) were also used as input, and four separate models were used for classification, with the five optimal sub-models selected for each feature (the specific model optimization process is the same as above). Through the above calculation process, a total of 5 × 5 = 25 model calculation results were obtained. In each calculation, the contribution value of each feature vector to the classification result can be obtained.

[0051] The five optimal models selected for each feature (out of a total of 25 models) are shown in the table below: The top-ranked feature variables and their contributions for each feature in the optimal model are as follows: 1. DNA Fragment Size Ratio (FSR) Deep Learning Neural Network Model (NN): 2. 5' end breakpoint sequence (BPM) read elastic network regression model (GLM): 3. Copy Number Variation (CNV) in Deep Learning Neural Network Models (NN) 4.16bp Novel Neo-Original (NEO) XgBoost Model for Tumors: 5. Nucleosome Coverage Mode (NCP) XgBoost Model: To further improve the predictive performance of the classifier, a secondary stacking training was performed on the results of 25 trained models. Stacking is an ensemble learning technique that stacks the results of 25 underlying weak classifiers (1... st Meta-learning is performed again on the -level basemodel (2) nd Stacking (a type of meta-learning algorithm) collects the features of each underlying classifier and finds the optimal integration method to improve the model's predictive performance. Ultimately, the training algorithms used in this patent's Stacking are the Generalized Linear Model (GLM), the Extreme Gradient Boosting (XGBoost) model, and the Deep Learning Regression model.

[0052] The importance of the three best ensemble models and the variable models is shown in the table below: The optimal stacking model with the highest AUC is the Generalized Linear Model (GLM). It establishes the relationship between the expected value of the response variable and the linear combination of predictor variables through a connection function, transforming 25 training models into the final linear equation: ALL Stacked = Intercept + A*CNV model 1 + B*CNV model 2 + C*CNV model 3 + D*CNV model 4 + E*CNV model 5 + F* BPM model 1 + G* BPM model 2 + H*BPM model 3 + I*BPM model 4 + J*BPM model 5 + K*FSR model 1 + L* FSR model 2 + M*FSR model 3 + N*FSR model 4 + O* FSR model 5 + P* NEO model 1 + Q* NEO model 2 + R* NEO model 3 + S* NEO model 4 + T* NEO model 5 + U* NCP model 1 + V* NCP model 2 + W* NCP model 3 + X* NCP model 4 +Y*NCP model 5 Where Intercept and AI are parameters of the linear equation.

[0053] The specific coefficients are as follows: Each feature showed some predictive performance under different training algorithms, and training a single feature with a secondary set improved the predictive performance of that feature. The final prediction AUC for the training set reached 0.9474, and the prediction AUC for the validation set reached 0.931, with a sensitivity of 85% and a specificity of 98.7%.

Claims

1. A method for constructing a screening model for malignant pulmonary nodules, characterized in that, Includes the following steps: Step 1: Extract and sequence cfDNA from samples in the positive and control groups to obtain read data; Step 2: Analyze the nucleosome coverage patterns of the selected transcription factors as a feature set; Step 3: Use the feature set as the model feature vector input into the classification model, and use the benign or malignant nature of lung nodules as the output value to train the model and obtain the early screening model; Step 2 includes: Step 2-1: Obtain transcription factors from the GTRD database and exclude transcription factors that do not have known transcription sites in the CIS-BP database; Step 2-2: Use the range of -5kb to +5kb near the transcription site of the obtained transcription factor as a window to obtain fragments of 100-220bp in length that can be aligned to the window. Perform GC correction and sequencing depth smoothing on the read data in the window in sequence to obtain the coverage pattern curve of each transcription factor. Steps 2-3: For each transcription factor, obtain the following three features, which together serve as the nucleosome coverage pattern: 1) For all transcription sites of the transcription factor, calculate the average depth from the upper 1kb to the lower 1kb of the transcription site; 2) For the obtained coverage pattern curve, obtain the amplitude value of the curve trough, which serves as the center depth of the transcription factor; 3) Perform a fast Fourier transform on the obtained coverage pattern curve to obtain the amplitude value of the highest point of the nucleosome amplitude signal. The classification model is selected from one or more of the following: elastic network regression model, extreme gradient boosting model, random forest model, and deep learning neural network model.

2. The construction method according to claim 1, characterized in that, In step 2-1, transcription factors were screened from the GTRD database v21.12, and transcription factors that did not have known transcription sites in the CIS-BP database v2.00 were excluded, resulting in 334 transcription factors with more than 10,000 highly matched sites.

3. The construction method according to claim 1, characterized in that, The sequencing depth smoothing process in step 2-2 uses a Savitzky-Golay filter with a polynomial power of 3, and the nucleosome coverage pattern features contain 1002 feature values.

4. The construction method according to claim 1, characterized in that, The classification model is optimized using random search hyperparameters.

5. The construction method according to any one of claims 1-3, characterized in that, The following feature set is further added, and this feature set, together with the nucleosome coverage pattern, is used as the initial feature value input into the classification model for training: First feature set: the proportion of short reads and the proportion of long reads in different windows of the cfDNA fragment aligned to the reference genome; Second feature set: the proportion of m bases in the cfDNA fragment aligned to the 5′ end of the reference genome in all base fragments; The third feature set: the copy number of WGS data in different windows of the chromosome; The fourth feature set: the proportion of new short tumor sequences of 16bp.

6. The construction method according to claim 5, characterized in that, The first feature set, the second feature set, and the third feature set are obtained respectively in the following ways: (1) The reference genome was divided into multiple 5Mb windows, totaling 541 windows. The number of short reads with a length of 100-150bp and the number of long reads with a length of 151-220bp were obtained within each window. The ratio of the standardized number of short reads to the number of long reads was used as the first feature set. (2) Take the base sequences of 4 bp upstream and downstream of the 5′ end of the read segment as the base fragment set, and count the proportion of 65,536 8 bp breakpoint sequences in all fragments as the second feature set; (3) Divide chromosomes 1-22 into 2475 windows with 1Mb non-overlapping, and obtain the copy number data in each window as the third feature set.

7. The construction method according to claim 6, characterized in that, The fourth feature set is obtained in the following manner: Exhaustive search is used to generate the set A of all short sequences with a length of 16bp; The set B of all short sequences with a length of 16bp was obtained by exhaustively searching through the human reference gene sequence, and the invalidation component was obtained by removing set B from set A. 977 base substitution mutations that occurred at least twice were extracted from the WGS sequencing results of 2577 patients with colorectal cancer, lung cancer, breast cancer, gastric cancer, prostate cancer and liver cancer in the PCAWG database, and a set of null subsequences C was obtained based on the base substitution mutations. The base substitution mutations with a frequency greater than 0.01 in East Asian populations were obtained from the gnomAD database, and the set of invalid short sequences D was obtained accordingly. After excluding set D from set C, 4616 new 16bp tumor-related short sequences were obtained; The ratio of the number of samples of each new short sequence to the total number of samples that can read any new short sequence is used as the fourth feature set.

8. The construction method according to claim 7, characterized in that, The training of the model includes: The first, second, third, fourth and fifth feature sets are respectively input into different classifier models for training, and one or more optimal classifier models are obtained for each feature set. The optimal classifier model is then subjected to a second ensemble training to construct an ensemble classifier model.

9. The construction method according to claim 8, characterized in that, The secondary ensemble training employs a generalized linear model, an extreme gradient boosting model, or a deep learning regression model, and takes the average of the predicted probability results output by the top three ensemble models in terms of AUC ranking as the final judgment result.

Images