Method for screening breast cancer biomarkers based on weighted gene co-expression network analysis and application

Abstract

The invention discloses a method for screening breast cancer biomarkers based on weighted gene co-expression network analysis and application. High-expression keratin KRT5 and low-expression ubiquitin coupling enzyme UBE2C in breast cancer patients are screened out and are closely related to good prognosis of the high-expression keratin KRT5 and the low-expression ubiquitin coupling enzyme UBE2C, and the KRT5 and the UBE2C are effective molecular markers for breast cancer early-stage diagnosis and prognosis analysis. In-vitro experiments show that the development process of the breast cancer can be effectively inhibited by promoting the expression of the KRT5 gene and inhibiting the expression of the UBE2C gene. According to the invention, early-stage diagnosis and prognosis analysis of breast cancer can be rapidly carried out, a new direction is provided for diagnosis, treatment and prognosis of breast cancer, and development of individualized treatment is promoted.

Description

Technical Field

[0001] The present invention relates to the field of biomedical technology, and more specifically to a screening method for breast cancer biomarkers based on weighted gene co-expression network analysis and its application. Background Art

[0002] Breast cancer (BRCA) is a type of epithelial tumor and a very common cancer. Female cancer patients often suffer from this disease, with more than 42,000 cases reported in 2020. Primary BRCA itself is not fatal; but BRCA cells can easily metastasize to other organs, including the brain, lungs, liver and bones. From the perspective of age of onset, the incidence of BRCA in my country has increased year by year since around the age of 20, reaching a peak around the age of 50. Of course, the global breast cancer mortality rate has gradually declined with the popularization of new technologies, new treatment options and new treatment strategies. However, the trend of declining mortality in BRCA is not obvious, especially in China's vast rural and underdeveloped areas. Although there have been some progress in diagnosis and treatment strategies, data show that nearly 60% of BRCA cases are diagnosed in the late stage, and the mortality rate is very high. The stage of breast cancer development is inseparable from its accurate prognosis. In short, the earlier the tumor is discovered in the patient's body, the higher the probability of survival within five years. Statistics from the International Cancer Agency show that the five-year relative survival rate of breast cancer patients is about 90%, and the five-year survival rate of patients with early invasive cancer is about 85%. However, the survival rate after invasive cancer metastasis is very different. The five-year survival rate of patients with distant metastasis has dropped to less than 28%. 25% to 50% of BRCA patients will eventually develop fatal metastatic tumors, usually decades after the primary tumor was diagnosed and removed.

[0003] Although significant progress has been made in understanding the progression and prognosis of BRCA, the current research shows that the occurrence and development mechanism of BRCA is not very clear, and there are many influencing factors. Some clinical pathological factors, such as tumor size, axillary lymph node status, and pathology, are used to determine the development direction of breast cancer patients, but these factors alone are often not used as a basis for selecting treatment strategies and predicting BRCA prognosis. Some unhealthy lifestyles are related to the onset of BRCA and will increase the incidence rate to a certain extent, such as overnutrition, obesity, high-fat diet, excessive drinking, etc. BRCA is a multifactorial disease caused by gene mutations. There are many gene changes that can affect the relevant processes of BRCA and can also indicate the prognosis of BRCA. Finding biological markers for breast cancer is extremely important for evaluating its progression and prognosis, which will help patients with early treatment and improve survival rates.

[0004] Therefore, how to screen and analyze breast cancer biomarkers is an urgent problem that technicians in this field need to solve. Summary of the invention

[0005] In view of the problems existing in the prior art, the purpose of the present invention is to provide a screening method and application of breast cancer biomarkers based on weighted gene co-expression network analysis. The potential biomarkers keratin (KRT5) and ubiquitin conjugating enzyme (UBE2C) screened by the method of the present invention are described. The close correlation between high-expression keratin (KRT5) and low-expression ubiquitin conjugating enzyme (UBE2C) can be used to quickly perform early diagnosis and prognosis analysis of breast cancer. The present invention provides a new direction for the diagnosis, treatment and prognosis of breast cancer, and promotes the development of "individualized treatment".

[0006] In order to achieve the above object, the present invention provides a method for screening breast cancer biomarkers based on weighted gene co-expression network analysis; specifically comprising the following steps:

[0007] 1) Data download and preprocessing:

[0008] Download genes and expression levels related to breast cancer from the TCGA database and the GEO database to obtain the TCGA-BRCA database; filter the obtained database to obtain the GSE83591 database;

[0009] Specifically: TCGA cancer genome database (https: / / portal.gdc.cancer.gov / ) and GEO high-throughput sequencing database (https: / / www.ncbi.nlm.nih.gov / gds) were used to download genes and expression levels related to breast cancer; there were 1222 BRCA samples in the database, including 1109 breast cancer cases and 113 normal tissues; the symbol script and the corresponding ensemble ID were converted to gene names, and genes with cpm (per million) ≥ 1 were retained; after using the function filtering in the edgeR software package, a total of 14902 genes with RPKM values ​​were calculated by dividing the gene count by the gene length for the subsequent analysis; GSE83591 consists of 79 tumor samples and 30 paired normal tissues from BRCA patients, and was analyzed by GPL14550 Agilent-028004SurePrint G3 Human GE8x60K Microarray (ProbeName Version). In the end, a list of 22070 genes was selected for subsequent analysis;

[0010] 2) Data analysis:

[0011] The gene co-expression network of TCGA-BRCA database and GSE83591 database was constructed by using WGCNA software package. Specifically, the genes screened were clustered by traits, and the correlation between the positive correlation modules of the clusters and the clinical traits was analyzed; the correlation between the negative correlation modules of the clusters and the clinical traits was analyzed; finally, the genes with correlation coefficients greater than 99% were clustered.

[0012] 3) Screening for overlapping genes:

[0013] The differentially expressed genes DEGs in the TCGA-BRCA database and the GSE83591 database were screened and gradually converted into visual volcano maps and heat maps. The software package was then used to screen the overlapping genes in the important expression module genes in DEGs and WGCNA and present them;

[0014] 4) Determine overlapping gene functions:

[0015] The clusterProfiler software package was used to determine the functions of overlapping genes, with the standard adjusted to p < 0.05, and bar graphs and bubble graphs of GO and KEGG gene functional pathways were drawn;

[0016] 5) Construction of protein interaction PPI and screening of Hub genes:

[0017] Predict the interaction between overlapping genes and proteins, construct the network PPI of proteins, screen genes with a score greater than 0.9 for preservation, build the model using Cytoscape, use the maximum cluster centrality algorithm to obtain the central nodes, and select the top ten as hub genes;

[0018] 6) Verify the value of Hub genes:

[0019] The R software package ggpubr was used to verify the expression levels of Hub genes in different clinical periods to ensure authenticity and reliability;

[0020] The R software packages survminer and survival were then used to analyze and screen genes with a correlation coefficient of ≥99% between overall survival and Hub genes, and the relationship between Hub genes and disease-free survival in BRCA patients was analyzed, among which genes KRT5 and UBE2C with P < 0.01 were retained;

[0021] 7) Correlation analysis of breast cancer samples:

[0022] Clinical analysis of genes KRT5 and UBE2C was performed to determine the intensity of gene expression in each stage and the relationship between the four stages according to the four stages of breast cancer;

[0023] 8) Verify the protein expression of Hub gene:

[0024] The KRT5 gene was obtained from the human protein atlas database to obtain relevant sample cells, and the differences between normal samples and tumor samples were compared and analyzed to observe the expression changes of the KRT5 gene in tumor samples;

[0025] 9) Obtain metabolic pathways of key genes:

[0026] Metabolic pathway search was performed for the UBE2C gene.

[0027] Preferably, in step 1), the TCGA database includes 1222 BRCA samples, and the BRCA samples include: 1109 breast cancers and 113 normal tissues;

[0028] The GSE83591 database consists of 9 tumor samples and 30 paired normal tissues in the BRCA database, and a list of 22,070 genes was determined by screening with the GPL14550 Agilent-028004 SurePrint G3 Human GE 8x60K Microarray analysis.

[0029] Preferably, step 3) screening overlapping genes specifically includes:

[0030] 31) The logFC filtering condition was set to ≥2.5% and the corrected p value was 0.01 to obtain 3377 DEGs in TCGA and 6784 DEGs in GEO;

[0031] 32) The DEGs were then converted into a visual volcano map using the R package ggplot2 and then into a heat map using the pheatmap package;

[0032] 33) R software package VennDiagram was then used to extract overlapping genes from DEGs and important expression module genes after WGCNA and present them as a Venn diagram.

[0033] Furthermore, the R package limma (linear model of microarray data) was used to find differentially expressed genes (DEGs) between BRCA and normal tissues; limma was used to screen DEGs in the TCGA-BRCA and GSE83591 datasets respectively; 3377 DEGs in TCGA and 6784 DEGs in GEO were found by setting logFC filtering conditions and modified p-value filtering conditions; the R package ggplot2 was used to convert the DEGs of the TCGA-BRCA and GSE83591 datasets into a visual volcano map, and then the pheatmap package was used to convert them into a heat map; the R package VennDiagram was then used to extract overlapping genes from the DEGs and the important expression module genes after WGCNA, and presented as a Venn diagram.

[0034] Preferably, the GO annotation in step 4) includes biological processes, cellular components and molecular functions.

[0035] Preferably, the specific process of step 5) constructing protein interaction PPI and screening Hub genes includes:

[0036] 51) Use the online analysis tool STRING to effectively predict the interactions between overlapping gene proteins to construct the protein network PPI;

[0037] 52) Genes with a score greater than 0.9 were selected for preservation, and the model was constructed using Cytoscape. The maximum cluster centrality algorithm was used to find the central node; the top ten MCC values ​​were taken as Hub genes.

[0038] Furthermore, a graph theory algorithm (Bron-Kerbosch algorithm) is used to find all the largest gene clusters in the graph: the first step is to find the gene points that have been added to the current maximum cluster; the second step is to find the gene points that may be added; the third step is to re-determine the gene points that have completed the maximum cluster counting; the fourth step is to take the intersection of the gene points screened out in the second and third steps, and put the intersection gene cluster into the first step gene cluster, so as to ensure that the gene points obtained in the first step are the largest gene cluster.

[0039] Then the centrality of the nodes in each largest cluster is calculated, and finally the sum of the centralities of these largest clusters is compared to find the cluster with the highest centrality.

[0040] Preferably, the verification in step 6) is specifically as follows: using the R software package ggpubr to classify the screened overlapping genes into high, medium and low categories, and corresponding them with the clinical stages of breast cancer T1, T2, T3 and T4 according to the correlation coefficient ≥99%, and then presenting them in the form of a box plot.

[0041] Preferably, the clinical analysis in step 7) is performed using RStudio's limma package, ggplot2, and ggpubr, combined with RPKM gene expression data.

[0042] Preferably, step 9) is to use the Mapper tool in KEGG to search the metabolic pathway of the UBE2C gene.

[0043] Another object of the present invention is to provide an application, specifically the application of the screening method in clinical research on breast cancer, and more specifically the application in early diagnosis and prognosis analysis of breast cancer based on gene expression levels.

[0044] A further object of the present invention is to use the screening method in the preparation of breast cancer drugs.

[0045] In summary, the technical effects produced by the technical solution provided by the present invention include at least:

[0046] The present invention not only confirms the relationship between genes KRT5 and UBE2C and breast cancer to a certain extent. At the same time, it confirms the correlation between genes KRT5 and UBE2C and breast cancer progression, prognosis, etc., filling the gap in current research. Applying genes KRT5 and UBE2C to breast cancer diagnosis and prognosis can save patient treatment time and cost and improve the accuracy of diagnosis. Further converting them into applications such as kits, test reagents, gene chips, etc. can better play their practicality. It can be seen that the present invention can achieve "individualized treatment" to a certain extent, and at the same time promote the development of "precision medicine". BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0048] Figure 1 :GO enrichment analysis diagram;

[0049] Figure 2.1 : OS survival analysis chart;

[0050] Figure 2.2 : Disease-free survival analysis chart;

[0051] Figure 3.1 For the KRT5 gene, enter the HPA human protein image database to search for its related expression status images, and select the sample number CAB000129;

[0052] Figure 3.2 For the KRT5 gene, enter the HPA human protein image database to search for its related expression status images and compare the images. The sample number is CAB011464;

[0053] Note: The same lowercase letters are not significant (p>0.05), and the different lowercase letters are significant (p<0.05). DETAILED DESCRIPTION

[0054] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0055] Example 1

[0056] This embodiment provides a method for screening breast cancer biomarkers based on weighted gene co-expression network analysis, specifically:

[0057] 1WGCNA Analysis

[0058] 1.1 Data Source

[0059] The TCGA cancer genome database (https: / / portal.gdc.cancer.gov / ) and the GEO high-throughput sequencing database (https: / / www.ncbi.nlm.nih.gov / gds) were used to download the clinical information gene expression profile of BRCA. There are 1222 BRCA samples in the database, including 1109 breast cancers and 113 normal tissues. The symbol script and the corresponding ensemble ID were used to convert them into gene names.

[0060] Genes with cpm (per million) ≥ 1 were retained to ensure the effectiveness of screening genes associated with breast cancer. After filtering using the function of correlation>99% in the edgeR software package, a total of 14,902 genes with RPKM values ​​were calculated by dividing the gene count by the gene length for subsequent analysis. GSE83591 consists of 79 tumor samples and 30 paired normal tissues from BRCA patients and is analyzed by GPL14550 Agilent-028004SurePrint G3Human GE 8x60KMicroarray (Probe Name Version). In the end, a list of 22,070 genes was selected for subsequent analysis.

[0061] 1.2 Analysis process

[0062] 1.2.1 Identification of key co-expression using WGCNA modules

[0063] The gene expression data maps of TCGA-BRCA and GSE83591 were constructed into a gene co-expression network using the WGCNA package in R. Cluster analysis of gene expression and correlation analysis of modules and clinical traits were performed using the WGCNA package.

[0064] 1.2.2 Differential expression analysis and intersection with module genes

[0065] The R package limma (linear model of microarray data) was used to find differentially expressed genes (DEGs) between BRCA and normal tissues. Limma was used to screen DEGs in the TCGABRCA and GSE83591 datasets. By setting the logFC filter condition to ≥2.5% and the corrected p value to 0.01, 3377 DEGs in TCGA and 6784 DEGs in GEO were found. The R package ggplot2 was used to convert the DEGs of the TCGA-BRCA and GSE83591 datasets into a visual volcano map, and the pheatmap package was used to convert them into a heat map. The R package VennDiagram was then used to extract overlapping genes from the DEGs and the important expression module genes after WGCNA, and presented as a Venn diagram.

[0066] 1.2.3 Functional annotation of overlapping genes (GO enrichment and KEGG pathway)

[0067] In order to explore the functions of the overlapping genes screened out, the R package clusterProfiler was used to explore the functions between genes, and the cutoff standard was adjusted p < 0.05.

[0068] GO annotation includes biological process (BP), cellular component (CC) and molecular function (MF). Draw bar charts and bubble charts of GO and KEGG gene function pathways.

[0069] 1.2.4 Construction of protein interaction PPI and screening of Hub genes

[0070] Among the 65 overlapping genes obtained, the online analysis tool STRING (https: / / string-db.org / ) can be used to make effective predictions on protein interactions and to construct a protein network PPI. Genes with a score greater than 0.9 were selected for preservation, and then the model was constructed using Cytoscape, using the maximum cluster centrality algorithm (MCC) to find the central node (a total of 10 nodes). The top ten MCC values ​​were taken as hub genes, and the MCC value was calculated by CytoHubba, a plug-in in Cytoscape.

[0071] 1.2.5 Validation of key gene expression patterns and prognostic value

[0072] To verify the reliability of the Hub gene, the R software package ggpubR was used to divide the genes screened in the previous step into three categories of high, medium, and low, and the correlations were compared with the four clinical stages T1, T2, T3, and T4 (correlation coefficient ≥ 99%);

[0073] Among them, T1: the maximum diameter of the tumor is less than 2 cm;

[0074] T2; the maximum diameter of the tumor is 2 to 5 cm;

[0075] T3: The maximum diameter of the tumor exceeds 5 cm;

[0076] T4: Regardless of tumor size, the tumor invades the chest wall or skin (including inflammatory breast cancer);

[0077] The data were then presented in a box plot. The R software packages survminer and survival were used to analyze the correlation between overall survival (OS) and hub genes (correlation coefficient ≥ 99%), and the online tool GEPIA2 was used to perform a correlation analysis on the relationship between the above hub genes and disease-free survival (DFS) in BRCA patients (correlation coefficient ≥ 99%). The data with P < 0.01 were statistically significant and could be retained as reliable records.

[0078] 1.2.6 Clinical analysis of two hub genes

[0079] Clinical analysis was performed on the two most relevant genes KRT5 and UBE2C identified by survival analysis. RStudio's limma package, ggplot2, and ggpubr were used in combination with RPKM gene expression data. Breast cancer was divided into four stages, and the intensity of gene expression in each stage and the relationship between the four stages were explored.

[0080] 1.2.7 Validation of survival-related Hub gene protein expression in the HPA database

[0081] HPA (Human Protein Atlas) is a valuable database from which a large amount of human-specific tissue and cell transcriptomics and proteomics data can be obtained. The KRT5 gene was searched for relevant sample cells in the Human Protein Atlas Database (https: / / www.proteinatlas.org), and the differences between normal samples and tumor samples were compared and analyzed to check the expression status; such as comparing the expansion of mammary ducts, the number of epithelial cells in the glandular ducts, the number of myoepithelial cells, whether the cells were arranged in an orderly manner, whether the nuclei were atypical, and whether the cells had mitotic figures.

[0082] 1.2.8 Use KEGG Mapper to query the metabolic pathways of key genes

[0083] The Mapper tool in KEGG was used to search the metabolic pathways of the UBE2C gene screened previously. KEGG (www.kegg.jp) was version 102.0 updated on April 1, 2022.

[0084] Example 2

[0085] Establishment of induced model

[0086] 2.1 Experimental animals and materials

[0087] Sixty 4-week-old female Sprague-Dawley (SD) rats weighing (60±5) g were purchased from the Experimental Animal Center of Xiangya School of Medicine, Central South University, and were specific pathogen-free (SPF) grade. Rats were housed in a standard animal laboratory, 8 rats per cage, and had free access to food and water. The room temperature was controlled at 18-25°C, and the relative humidity was 30%-50%. The laboratory floor and the bottom cover of the mouse cage were treated with 2% sodium thiosulfate phosphate buffer.

[0088] N-Methyl-N-nitrosourea (1-methyl-1-nitrosourer, MNU) was purchased from Sigma.

[0089] 2.2 Experimental Grouping

[0090] SD rats were injected intraperitoneally with MNU 25 mg / kg to establish the model, and then randomly divided into a control group and an experimental group according to the random number table method, with 30 rats in each group.

[0091] 2.3 Morphology and light microscopy observation

[0092] After 4 weeks (early stage), 8 weeks (middle stage) and 16 weeks (late stage) of modeling, the rats were killed and the mammary gland and tumor tissues were taken, and the tumor weight and tumor number of rats in each group were recorded. The mammary gland and tumor tissues were fixed with neutral formalin solution and paraffin sections were stained with HE for histopathological diagnosis, and malignant tumors were confirmed when there were more than 5 nuclear division phases in every 10 high-power fields (400 times).

[0093] 2.4 Detection of KRT5 and UBE2C expression

[0094] The mammary tissues of the model rats were collected at different stages of early, mid and late periods, crushed and centrifuged, and the contents of KRT5 and UBE2C were determined using ELISA kits (purchased from Shanghai Hengyuan Biotechnology Co., Ltd.).

[0095] 3WACNA analysis results

[0096] 3.1 Construction of weighted gene co-expression module

[0097] The R package WGCNA was used to construct gene co-expression networks with the data of TCGABRCA and GSE83591, respectively. Different colors were assigned to different modules, and 9 modules in TCGA-BRCA and 15 modules in GSE83591 were determined. A heat map was drawn to show the differential association between normal and tumor in the module. It showed that the turquoise module in TCGA-BRCA and the yellow module in GSE83591 had the highest correlation with normal tissue (turquoise: r = 0.77, p = 8.999999999999999e-236; yellow: r = 0.83, p = 3e-29). The volcano map of gene expression differences was drawn using the set parameters. The genes with greater differences were distributed at both ends, so the data points in the upper left and upper right corners were often more meaningful for biological research.

[0098] 3.2 Identification of genes between differentially expressed and co-expressed modules

[0099] P<0.05. Using the limma software package, 3377 DEGs and 6784 DEGs were found to be dysregulated in tumor tissues in the TCGA dataset and GSE83591 dataset, respectively. 7321 and 3445 co-expressed genes were found in the turquoise module of the TCGA dataset and the yellow module of GSE83591, respectively. A total of 65 overlapping genes were extracted by taking the intersection method to verify the co-expression network.

[0100] 3.3 Functional enrichment analysis of 65 genes

[0101] To further understand the potential functions of the 65 genes that overlapped with the differentially expressed gene list and the two co-expression modules, clusterProfiler, enrichplot was used after screening GO enrichment analysis, as shown in Figure 1Several enriched gene sets are displayed. From the bar chart, the longer the bar, the more genes with this function, and vice versa. The redder the color, the more significant the difference in the results, and the more statistically significant it is. From the bubble chart, the size of the circle indicates the number of genes, and the depth of the red represents the size of the P value, that is, the significance of the difference. The biological processes (BP) of the 65 genes are mainly concentrated in mitotic nuclear division, negative regulation of hydrolase activity, and enzyme inhibitor activity. In the cell component analysis, these genes are mainly involved in myosin filaments, blood microparticles, and spindle microtubules. In the molecular function (MF) analysis, peptidase regulation activity and enzyme inhibitor activity are related to the 65 genes. At the same time, combined with KEGG analysis, it can be concluded that these genes mainly control cell mitosis, and the control method is to regulate enzyme activity and regulate microtubule binding proteins.

[0102] 3.4 Construction of protein interaction PPI network and identification of hub genes

[0103] The hub genes selected from the PPI network were selected using the MCC algorithm of the CytoHubba plugin in Cytoscape. According to MCCsores, the ten genes with the highest scores, including FOS, CONB2, SPAG5, UBE2C, NUSAP1, BIRC5, CDCA8, PTTG1, ALB, and ACTG1, were selected as hub genes. By comparing the number of gene connections, it was found that the two genes ALB and ACTGI had little interaction with other genes and individuals, so KRT5 and KRT6B, which are also two nodes, were used to replace the position of their core genes.

[0104] 3.5 Key gene expression patterns, prognostic value, and protein expression verification

[0105] After screening out 10 hub genes (FOS, CONB2, SPAG5, UBE2C, NUSAP1, BIRC5, CDCA8, PTTG1, KRT5, KRT6B) through the CytoHubba plug-in, OS and DFS analysis was performed on the 10 hub genes to study the prognostic value of hub genes in BRCA patients. Figure 2.1The horizontal and vertical axes represent the survival time and survival rate, respectively. As time goes by, the patient's survival rate will gradually decrease. To make a conclusion, it is necessary to make a significance judgment based on the size of the P value. Among the 10 hub genes, the higher expression of KRT5 in BRCA patients is associated with increased overall survival (OS) and is the most significant gene. UBE2C is the most significant gene in the disease-free survival analysis and is negatively correlated with survival rate ( Figure 2.2 ). For other genes, the P values ​​were large and there was no significant difference, so no analysis was performed for the time being.

[0106] 3.6 Clinical analysis of two hub genes

[0107] The redundant parts of the clinical data were deleted, and the results were obtained after staging. The horizontal axis represents the stage of breast cancer, and the vertical axis represents the relative expression. The horizontal line is the P value. It can be seen that KRT5 has a large expression difference between stage one and stage four, but it has always been at a higher state, which corresponds to the OS survival analysis. The high expression of the whole process has a certain positive effect on the patient's survival status. On the other hand, UBE2C is only significant compared with the later stage, and it is less and less significant from stage two to later periods, which means that the expression of UBE2C decreases in later patients. This is also consistent with the results of the disease-free survival analysis, where patients in the later stage have low expression.

[0108] 3.7 Comparison of HPA samples

[0109] For the KRT5 gene, we searched the HPA human protein image database for relevant expression status images and made comparisons. We selected sample number CAB000129 ( Figure 3.1 ) and CAB011464( Figure 3.2 ). The difference can be clearly observed. Figure 3.1 This is a detection diagram of KRT5 protein expression in breast cells by immunohistochemistry. The tumor image was taken from a 37-year-old female patient with patient number 2174 in the database, and the normal image was taken from a 39-year-old female patient with patient number 1801 in the database. UBE2C selected a sample with number CAB011464, and the tumor image was taken from a 41-year-old female patient with number 2898 ( Figure 3.2 A in the figure), the normal image is a 45-year-old woman with the ID number 3544 ( Figure 3.2 B).

[0110] 3.8 Gene UBE2C affects ubiquitin-mediated protein hydrolysis

[0111] By affecting ubiquitin-mediated protein hydrolysis, the result was a pathway numbered hsa04120, which showed that UBE2C, as an effective gene affecting the production of ubiquitin-conjugating enzymes, together with genes such as UBE2A, determines the production of ubiquitin-conjugating enzymes. Thereafter, together with E1 (ubiquitin activating enzyme) and E3 (ubiquitin ligase), ubiquitin is added to the protein to be degraded.

[0112] 4 Animal testing results

[0113] 4.1 Pathological observation

[0114] The animal model of breast cancer precancerous lesions is characterized by atypical mammary hyperplasia. Its pathological features include the expansion of mammary ducts and the formation of adjacent masses. In these glandular ducts, the number of epithelial cell layers increases significantly, while the number of myoepithelial cell layers decreases relatively. The shapes and sizes of cells vary and are arranged in disorder. Sometimes, cells can be seen squeezing or overlapping each other. The nuclei are atypical, and mitotic figures are occasionally seen. These pathological changes provide important experimental models and theoretical basis for the occurrence and development of breast cancer.

[0115] 4.2 KRT5 and UBE2C expression detection results

[0116] Evaluation of KRT5 and UBE2C expression in breast tumors. Compared with the negative control group, the expression of KRT5 in the model group breast cancer mice at different stages was significantly reduced, while the expression of UBE2C was significantly increased. In addition, the expression of KRT5 in the early stage of breast cancer mice was 5 times that of the middle stage and 5 times that of the late stage; the expression of UBE2C in the early stage of breast cancer mice was 1 / 4 times that of the middle stage and 1 / 9 times that of the late stage. The results showed that the high expression of KRT5 and the low expression of UBE2C in the early stage of breast cancer mice were closely related to different stages of the breast.

[0117] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0118] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the concept or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for screening breast cancer biomarkers based on weighted gene co-expression network analysis, comprising the following steps: 1) Data download and preprocessing: Download genes and expression levels related to breast cancer from the TCGA database and the GEO database to obtain the TCGA-BRCA database; filter the obtained database to obtain the GSE83591 database; 2) Data analysis: The gene co-expression network of TCGA-BRCA database and GSE83591 database was constructed by WGCNA software package. Specifically, the genes screened were clustered into traits, and the correlation analysis between the positive correlation modules of the clusters and the clinical traits was performed; the correlation analysis between the negative correlation modules of the clusters and the clinical traits was performed; finally, the genes with correlation coefficients greater than 99% were clustered. 3) Screening for overlapping genes: The differentially expressed genes DEGs in the TCGA-BRCA database and the GSE83591 database were screened and gradually converted into visual volcano maps and heat maps. The software package was then used to screen the overlapping genes in the important expression module genes in DEGs and WGCNA and present them; 4) Determine overlapping gene functions: The clusterProfiler software package was used to determine the functions of overlapping genes, with the standard adjusted to p < 0.05, and bar graphs and bubble graphs of GO and KEGG gene functional pathways were drawn; 5) Construction of protein interaction PPI and screening of Hub genes: Predict the interaction between overlapping genes and proteins, construct the network PPI of proteins, screen genes with a score greater than 0.9 for preservation, build the model using Cytoscape, use the maximum cluster centrality algorithm to obtain the central nodes, and select the top ten as hub genes; 6) Verify the value of Hub genes: The R software package ggpubr was used to verify the expression levels of Hub genes in different clinical periods to ensure authenticity and reliability; The R software packages survminer and survival were then used to analyze and screen genes with a correlation coefficient of ≥99% between overall survival and Hub genes, and the relationship between Hub genes and disease-free survival in BRCA patients was analyzed, among which genes KRT5 and UBE2C with P < 0.01 were retained; 7) Correlation analysis of breast cancer samples: Clinical analysis of genes KRT5 and UBE2C was performed to determine the intensity of gene expression in each stage and the relationship between the four stages according to the four stages of breast cancer; 8) Verify the protein expression of Hub gene: The KRT5 gene was obtained from the human protein atlas database to obtain relevant sample cells, and the differences between normal samples and tumor samples were compared and analyzed to observe the expression changes of the KRT5 gene in tumor samples; 9) Obtain metabolic pathways of key genes: Metabolic pathway search was performed for the UBE2C gene.

2. A method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 1) The TCGA-BRCA database includes 1222 BRCA samples, including: 1109 breast cancers and 113 normal tissues; The GSE83591 database consists of 9 tumor samples and 30 paired normal tissues in the BRCA database, and a list of 22,070 genes was determined by screening with the GPL14550 Agilent-028004 SurePrint G3 Human GE 8x60K Microarray analysis.

3. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 3) Screening for overlapping genes specifically includes: 31) The logFC filter condition was set to ≥2.5 and the corrected p value was 0.01 to obtain 3377 DEGs in TCGA and 6784 DEGs in GEO; 32) The DEGs were then converted into a visual volcano map using the R package ggplot2 and then into a heat map using the pheatmap package; 33) R software package VennDiagram was then used to extract overlapping genes from DEGs and important expression module genes after WGCNA and present them as a Venn diagram.

4. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 4) The GO annotation includes biological process, cellular component and molecular function.

5. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 5) Constructing protein interaction PPI and screening Hub genes. The specific process includes: 51) Use the online analysis tool STRING to effectively predict the interactions between overlapping gene proteins to construct the protein network PPI; 52) Genes with a score greater than 0.9 were selected for preservation, and the model was constructed using Cytoscape. The maximum cluster centrality algorithm was used to find the central node; the top ten MCC values ​​were taken as hub genes.

6. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: The verification in step 6) is specifically as follows: using the R software package ggpubr to classify the screened overlapping genes into high, medium and low categories, and corresponding them with the clinical stages of breast cancer T1, T2, T3 and T4 according to the correlation coefficient ≥99%, and then presenting them in the form of a box plot.

7. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 7) The clinical analysis was performed using RStudio's limma package, ggplot2, and ggpubr, combined with RPKM gene expression data.

8. The method for screening breast cancer biomarkers based on weighted gene co-expression network analysis according to claim 1, characterized in that: Step 9) is to use the Mapper tool in KEGG to search the metabolic pathway of the UBE2C gene.

9. The screening method according to any one of claims 1 to 8, characterized in that: The screening method is placed on application in clinical research of breast cancer.

10. The screening method according to any one of claims 1 to 8, characterized in that: Application of the screening method in the preparation of breast cancer drugs.

Images