Method and system for identifying downstream motion factors for regulating brown fat metabolism and storage medium

Abstract

The invention relates to a downstream motion factor identification method and system for regulating brown fat metabolism and a storage medium, and relates to the technical field of bioinformatics, and the method comprises the following steps: (1) processing a sample and obtaining transcriptome original data; (2) difference analysis; (3) time sequence analysis; (4) enrichment analysis; (5) carrying out weighted gene co-expression network analysis; (6) machine learning; according to the method, downstream metabolism-promotion factors of brown adipocytes are identified by using RNA next-generation sequencing results, overall characteristics of transcriptomes are obtained through difference analysis, enrichment analysis and time sequence analysis, data are further reduced by using weighted gene co-expression network analysis and machine learning, and key regulatory genes are screened out; according to the invention, key regulation molecules capable of promoting metabolism of brown adipocytes can be found, and a new biological target is provided for treating obesity.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and particularly to a method, system, and storage medium for identifying downstream exercise factors that regulate brown fat metabolism. Background Art

[0002] Brown adipocytes are emerging therapeutic targets for metabolic diseases. After genetically edited brown adipocytes are transplanted into obese mice, their obesity and metabolic phenotypes can be significantly improved. Exercise can affect the activation and function of brown fat through multiple mechanisms, including enhancing mitochondrial activity, increasing fatty acid oxidation, regulating heat consumption, etc. Nevertheless, the specific molecular mechanisms by which exercise regulates brown fat, especially the downstream factors and signaling pathways involved, have not been fully revealed. In recent years, with the development of high-throughput RNA second-generation sequencing technology, transcriptome-based analysis has become a powerful tool for studying the biological processes of brown fat. RNA sequencing can systematically identify genes whose expression changes under specific stimuli (such as exercise), providing rich information for revealing the regulation of exercise on brown fat. However, how to screen out exercise-related downstream factors from complex gene expression data and further verify their roles in brown fat function remains a challenge. Therefore, there is an urgent need for an efficient and accurate method to identify downstream factors of brown fat related to exercise, so as to promote the potential of brown fat in the treatment of obesity and metabolic diseases. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method, system, and storage medium for identifying downstream exercise factors that regulate brown fat metabolism. The purpose is to find key molecules that can promote brown adipocyte metabolism and provide new biological targets for the treatment of obesity.

[0004] The technical solution of the present invention to solve the above technical problems is as follows:

[0005] In the first aspect, a method for identifying downstream exercise factors that regulate brown fat metabolism includes the following steps:

[0006] (1) Sample processing and obtaining of raw transcriptome data:

[0007] Divide brown adipocytes into a control group and an experimental group. Administer a blank solvent to the control group (CON) to obtain samples of the control group at different time points after the administration treatment. Administer isoproterenol (ISO) or / and COX2 inhibitor NS-398 (NS) to the experimental group to obtain samples of the experimental group at different time points after the administration treatment;

[0008] Perform second-generation RNA sequencing on the samples of the control group at different time points after drug administration and the samples of the experimental group at different time points after drug administration to obtain the original transcriptome data of the control group and the original transcriptome data of the experimental group;

[0009] (2) Differential analysis:

[0010] Process the original transcriptome data of the control group and the original transcriptome data of the experimental group to obtain the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points; perform differential analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain differentially expressed genes;

[0011] (3) Temporal analysis:

[0012] Perform temporal analysis on the gene expression matrix of the experimental group at different time points, cluster according to the gene expression trend to obtain clusters, and perform functional annotation on the clusters;

[0013] (4) Enrichment analysis:

[0014] Perform GSEA (Gene Set Enrichment Analysis) enrichment analysis on the differentially expressed genes to obtain the GSEA analysis results;

[0015] Perform GSVA (Gene Set Variation Analysis) enrichment analysis on the differentially expressed genes to obtain the GSVA analysis results;

[0016] Perform GO and KEGG enrichment analysis on the differentially expressed genes to obtain the GO and KEGG analysis results;

[0017] Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways;

[0018] (5) Weighted gene co-expression network (WGCNA) analysis:

[0019] Perform weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain WGCNA key genes;

[0020] (6) Machine learning:

[0021] Perform Lasso regression feature screening on the differentially expressed genes, the genes in the clusters, and the WGCNA key genes to obtain Lasso candidate feature genes;

[0022] Perform random forest feature (RF) screening on the Lasso candidate feature genes to obtain the RF key genes;

[0023] Take the intersection of the WGCNA key genes, Lasso candidate feature genes, and RF key genes, which is the key regulatory molecule that promotes brown adipocyte metabolism.

[0024] Among them, isoproterenol (ISO) is used to activate ADRβ3; ADRβ3 (Adrenoceptor beta 3) is a G protein-coupled receptor that is mainly expressed in adipocytes and participates in the regulation of lipid metabolism, especially promoting lipolysis, that is, the decomposition of triglycerides in adipocytes; when epinephrine or norepinephrine binds to ADRβ3, a series of signal transduction pathways will be activated, ultimately leading to fat breakdown and energy release.

[0025] The COX2 inhibitor NS-398 (NS) is a selective cyclooxygenase-2 (COX-2) inhibitor. Among them, COX-2 (cyclooxygenase 2) is an enzyme necessary for the synthesis of prostaglandins (PGs) and is also a key rate-limiting enzyme in the initial step of PGs synthesis; its main function is to catalyze the conversion of arachidonic acid into prostaglandins.

[0026] Brown adipocytes are differentiated from the C3H10T1 / 2 cell line. The differentiation step is to induce the C3H10T1 / 2 cell line to differentiate into brown adipocytes using a medium containing insulin, dexamethasone, and isobutylmethylxanthine. The differentiation time is 10 days, and the differentiation efficiency is observed under the microscope.

[0027] The beneficial effects of the present invention are as follows: The present invention uses the results of RNA second-generation sequencing to identify the downstream promoting factors of brown adipocyte metabolism. Through differential analysis, enrichment analysis, and time-series analysis, the overall characteristics of the transcriptome are obtained. WGCNA and machine learning (Lasso regression, random forest) are used to further reduce the data and screen out the key genes; the present invention helps to find the key regulatory molecules that can promote brown adipocyte metabolism and provides new biological targets for the treatment of obesity.

[0028] On the basis of the above technical solutions, the present invention can also be improved as follows.

[0029] Further, in step (1), the experimental groups include the ISO group, the NS group, and the ISO-NS combined group;

[0030] The ISO group was given isoproterenol at a final concentration of 10 μM; the NS group was given the COX2 inhibitor NS-398 at a final concentration of 100 μM; the ISO-NS combined group was given isoproterenol and the COX2 inhibitor NS-398 at a final concentration of 10 μM for isoproterenol and 100 μM for the COX2 inhibitor NS-398.

[0031] In step (1), the different time points after the administration treatment include any time point between 0 and 24 h after the administration treatment.

[0032] The beneficial effects of adopting the above further scheme are as follows: The ISO group and the NS group were established to effectively compare the metabolic effects of brown adipocytes after simulating exercise after COX2 inhibition; samples of the control groups at different time points after the administration treatment (0, 3, 6, 12, 24 h) were collected to perform transcriptome sequencing of the cells, in order to conveniently observe the differential gene changes of time factors in brown adipocytes after different drug treatments.

[0033] Furthermore, in step (2), differential analysis was performed using DESeq2 or edgeR software; the FDR threshold was set to be lower than 0.05, and fold change ≥ 1.5.

[0034] Furthermore, the time series analysis was performed using the Masigpro package and the two-step regression method (StepwiseRegression) for time series analysis.

[0035] Furthermore, the time series analysis includes the following specific steps:

[0036] (3-1) Construct an input data frame containing the gene expression matrix of the experimental groups at different time points and the experimental group information (experimental group and time point);

[0037] (3-2) First, fit a global regression model to screen for differential genes related to time changes, and set the screening threshold to FDR < 0.05; then fit a differential model between experimental groups to identify differential genes related to time changes with significant differences between experimental groups, and perform k-means clustering according to the screened differential genes related to time changes with significant differences between experimental groups, and group them according to the expression trend to obtain clustering groups;

[0038] (3-3) For each of the clustering groups, perform GO and KEGG functional enrichment analysis using the ClusterProfiler package, and set the screening threshold to FDR < 0.05 for functional annotation.

[0039] Further, in step (4), the GSEA enrichment analysis includes the following specific steps: constructing an input framework containing differentially expressed genes and their significance statistical values (log2FC and Padj), evaluating the gene set enrichment using the GSEA method, setting the screening threshold as FDR < 0.05, and obtaining the GSEA analysis result;

[0040] In step (4), the GSVA enrichment analysis includes the following specific steps: using the differentially expressed genes as input, evaluating the enrichment degree of each group on a specific gene set, calculating the GSVA enrichment score (ES), performing a statistical test on the GSVA enrichment score, setting the threshold for screening significant gene sets as FDR < 0.05, and obtaining the GSVA analysis result;

[0041] In step (4), the GO and KEGG enrichment analysis includes the following steps: constructing an input framework containing differentially expressed genes and significance statistical values, performing enrichment analysis using the ClusterProfiler package, and setting the screening threshold as FDR < 0.05; classifying and statistically analyzing the GO and KEGG enrichment results, extracting significant functions and pathways, and obtaining the GO and KEGG analysis results.

[0042] Further, step (5) includes the following specific steps:

[0043] (5-1) Input data construction: Using the gene expression matrices of the control group at different time points and the gene expression matrices of the experimental group at different time points as input, constructing a data framework containing the gene expression information of all groups;

[0044] (5-2) Soft threshold screening: Using the pickSoftThreshold function in the WGCNA package to screen for suitable soft threshold parameters (power); selecting the soft threshold range as 9 according to the scale-free network characteristic standard (Scale-Free Topology Fit Index);

[0045] (5-3) Network construction and module division: Constructing a weighted gene co-expression network based on the selected soft threshold, calculating the weighted correlation coefficient between genes; using the Dynamic Tree Cut algorithm to divide genes into modules, classifying genes with similar expression patterns into the same module; setting the minimum number of genes in a module as 30, and merging modules with a module similarity > 0.75 to obtain the final co-expression module;

[0046] (5-4) Co-expression module and phenotype correlation analysis: Calculating the characteristic gene (moduleeigengene, ME) of each co-expression module, and performing a correlation analysis with the phenotype data of the experimental group and time points; screening modules significantly correlated with the phenotype, and setting the significance threshold as p < 0.05;

[0047] (5-5) Key gene screening: In the modules significantly related to phenotypes, screening is performed through gene module membership (MM) and the correlation between genes and phenotypes (gene significance, GS); the screening criteria are genes with MM > 0.8 and GS > 0.2, and the key genes in the module are extracted, which are the WGCNA key genes.

[0048] Furthermore, in step (6), the Lasso regression feature screening includes the following specific steps: First, use the Lasso regression model to perform feature screening on the differentially expressed genes, the genes of the clustering groups, and the WGCNA key genes, with the standardized gene expression matrix and the binary classification labels of the experimental groups as inputs; determine the optimal value of the regularization parameter λ through 10-fold cross-validation, and screen out the genes with non-zero regression coefficients as candidate feature genes, which are the Lasso candidate feature genes;

[0049] In step (6), the random forest feature (RF) screening includes the following specific steps: Construct a random forest classification model, calculate the importance scores (Mean Decrease Accuracy and Mean Decrease Gini) of each Lasso candidate feature gene, and screen out the key genes with an importance threshold > 0.01; the model parameters are set as the number of decision trees (n_estimators) 500, the maximum tree depth (max_depth) 10, and the minimum number of samples for splitting (min_samples_split) 2.

[0050] In the second aspect, the system adopted by the method for identifying downstream exercise factors regulating brown fat metabolism includes:

[0051] Module for sample processing and obtaining of raw transcriptome data: Divide brown fat cells into a control group and an experimental group, administer a blank solvent to the control group (CON) to obtain samples of the control group at different time points after the administration treatment, and administer isoproterenol (ISO) or / and COX2 inhibitor NS-398 (NS) to the experimental group to obtain samples of the experimental group at different time points after the administration treatment;

[0052] Perform RNA second-generation sequencing on the samples of the control group at different time points after the administration treatment and the samples of the experimental group at different time points after the administration treatment to obtain the raw transcriptome data of the control group and the raw transcriptome data of the experimental group;

[0053] Differential analysis module: Process the original transcriptome data of the control group and the original transcriptome data of the experimental group to obtain the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points; perform differential analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain differentially expressed genes at different time points;

[0054] Time series analysis module: Perform time series analysis on the gene expression matrix of the experimental group at different time points, cluster according to the gene expression trend to obtain clusters, and perform functional annotation on the clusters;

[0055] Enrichment analysis module: Perform GSEA (Gene Set Enrichment Analysis) enrichment analysis on the differentially expressed genes to obtain GSEA analysis results;

[0056] Perform GSVA (Gene Set Variation Analysis) enrichment analysis on the differentially expressed genes to obtain GSVA analysis results;

[0057] Perform GO and KEGG enrichment analysis on the differentially expressed genes to obtain GO and KEGG analysis results;

[0058] Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways;

[0059] Weighted gene co-expression network (WGCNA) analysis module:

[0060] Perform weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain WGCNA key genes;

[0061] Machine learning module: Perform Lasso regression feature screening on the differentially expressed genes, the genes in the clusters, and the WGCNA key genes to obtain Lasso candidate feature genes;

[0062] Perform random forest feature (RF) screening on the Lasso candidate feature genes to obtain RF key genes;

[0063] Take the intersection of the WGCNA key genes, the Lasso candidate feature genes, and the RF key genes, which is the key molecule promoting brown adipocyte metabolism.

[0064] In a third aspect, a computer-readable storage medium storing a computer program, where the computer program causes a computer to execute the steps (2) to (6). Brief Description of the Drawings

[0065] Figure 1 This is the experimental result graph of Oil Red O staining and glycerol quantification under different treatment conditions after culturing brown adipose cells of the present invention; wherein, A is the time-sequential graph of Oil Red O staining of lipid droplets in brown adipose cells; B is the time-sequential change graph of the Oil Red O staining area of lipid droplets in brown adipose cells; C is the comparison graph of the area under the curve (AUC) of the lipid droplet area change graph in brown adipose cells; D is the graph of the change in relative glycerol content in the culture medium of brown adipose cells; E is the comparison graph of the area under the curve of the relative glycerol content change curve in brown adipose cells;

[0066] Figure 2 This is the visualization graph of the differential expression gene analysis and time-sequential analysis results of the present invention; wherein, A is the heat map of differentially expressed genes in brown adipose cells at each time point after ISO and ISO_NS treatments; B is the difference in the time-sequential differential gene expression trends of 4 clusters; C is the bubble graph of the functional enrichment analysis of differentially expressed genes after 6h of drug treatment; D is the bubble graph of the functional enrichment analysis of differentially expressed genes after 12h of drug treatment; E is the time-sequential principal component analysis graph of differentially expressed genes in Cluster 2; F is the Wenn graph of differentially expressed genes in brown adipose cells after 6h and 12h of treatment; G is the functional enrichment graph of differentially expressed genes at 6h that are up-regulated after ISO treatment and down-regulated after ISO_NS treatment; H is the functional enrichment graph of differentially expressed genes at 12h that are up-regulated after ISO treatment and down-regulated after ISO_NS treatment;

[0067] Figure 3 This is the schematic diagram of WGCNA, machine learning and external chip sequencing verification of the present invention; wherein, A is the selection of the optimal soft threshold for the WGCNA module; B is the dendrogram of clustering based on highly correlated feature genes; C is the heat map of the correlation between the WGCNA module and traits; D is the venn diagram of genes up-regulated in ISO vs CON, down-regulated in ISO_NS vs ISO and module genes; E is the screening of genes based on the random forest model; F is the top ten genes ranked by the random forest model feature screening; G is the screening of genes based on the Lasso regression model and cross-validation method; H is the screening of genes based on Lasso regression analysis; I is the Venn diagram according to the results of random forest feature screening and Lasso regression feature screening, with 3 interactive genes; J is the expression trend of the Ywhah gene among different treatment groups of brown adipose cells; K is the external chip sequencing (Fu et al) verifying that the expression of the Ywhah gene increases in brown adipose tissue after exercise;

[0068] Figure 4This is the enrichment analysis result graph of the proteins that the molecules screened by the affinity purification-mass spectrometry (AP-MS) verification of the present invention can bind; among them, A is the KEGG signaling pathway graph; B is the molecular component GO graph; C is the biological process GO graph; D is the molecular function GO graph; E is the correlation between the expression level of Ywhah and the GSVA score; F is the GSVA score trend graph of different treatment groups of brown adipocytes. Detailed implementation mode

[0069] The principles and features of the present invention are described below. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention. For those without specific technical or conditions noted in the examples, the techniques or conditions described in the literature in this field are followed, or the product specifications are followed. For reagents or instruments without the manufacturer noted, they are all conventional products that can be purchased through regular channels.

[0070] Example

[0071] 1. Identification method of downstream exercise factors regulating brown fat metabolism

[0072] This example relates to an identification method of downstream exercise factors regulating brown fat metabolism, including the following steps:

[0073] (1) Sample treatment and acquisition of raw transcriptome data:

[0074] The brown adipocytes are divided into a control group and an experimental group. The control group (CON) is treated with a blank solvent to obtain samples of the control group at different time points after the administration treatment. The experimental group is treated with isoproterenol (ISO) or / and COX2 inhibitor NS-398 (NS) to obtain samples of the experimental group at different time points after the administration treatment;

[0075] The samples of the control group at different time points after the administration treatment and the samples of the experimental group at different time points after the administration treatment are subjected to RNA second-generation sequencing to obtain the raw transcriptome data of the control group and the raw transcriptome data of the experimental group;

[0076] It includes the following specific steps:

[0077] (1-1) Cultivation, differentiation and treatment of brown adipocytes:

[0078] ① The cultivation steps of C3H10T1 / 2 cells include:

[0079] The C3H10T1 / 2 cells (provided by the Chinese Academy of Sciences Cell Bank, Shanghai, China) are placed in DMEM / F12 medium containing 10% fetal bovine serum (10091148, Gibco) at 37 °C and 5% CO 2Cultured under culture conditions; when the cell confluence reached 70%-80%, the fresh medium was replaced until the confluence reached 90%-100%.

[0080] ② The steps of brown adipocyte differentiation include:

[0081] Transfer the C3H10T1 / 2 cells with a confluence of 90%-100% to the differentiation medium and culture for 4 days. The differentiation medium is DMEM / F12 medium supplemented with 10% fetal bovine serum, 20 nM insulin (I2643, Sigma), and 1 nM 3,3′,5-triiodo-L-thyronine (T2877, Sigma), and the medium is changed daily; on the 4th day, transfer the cells to the brown adipocyte induction medium and culture for 2 days. The cell induction medium is the differentiation medium supplemented with 2 μg / mL dexamethasone (D1756, Sigma), 0.5 mM isobutylmethylxanthine (I5879, Sigma), 0.125 mM indomethacin (I7378, Sigma), and 1 μM rosiglitazone (R2408, Sigma); after the induction culture, transfer the cells back to the differentiation medium and change the medium every other day until the 10th day of culture.

[0082] ③ Grouping after differentiation

[0083] The fully differentiated brown adipocytes were divided into the following four groups:

[0084] Control group (CON): Treated with an equal amount of DMSO and used as an experimental control; Isoproterenol group (ISO) (I5627, Sigma): Added with isoproterenol (ISO) at a final concentration of 10 μM to activate the β3-adrenergic receptor (ADRβ3); NS-398 group (NS): Added with NS-398 (COX2 inhibitor) (HY-13913, MCE) at a final concentration of 100 μM to inhibit the arachidonic acid metabolic pathway; Combined intervention group (ISO-NS): Treated with a final concentration of 10 μM ISO and a final concentration of 100 μM NS-398 simultaneously.

[0085] ④ Sample collection

[0086] Cell samples and culture supernatants were collected at different time points (0, 3, 6, 12, and 24 hours) after treatment for RNA extraction and detection of lipolysis-related indicators respectively; the cell samples were treated with TRIzol reagent (15596026, Thermo Fisher) and stored at -80 °C; the culture supernatants were centrifuged at 1000 rpm for 5 minutes to remove impurities and then stored in a -80 °C refrigerator as well.

[0087] The above steps are aimed at preparing sequencing samples and experimental verification to explore downstream motility factors.

[0088] (1-2) Second-generation RNA sequencing:

[0089] ① The steps for obtaining RNA-seq data include:

[0090] Through high-throughput sequencing technology, obtain the RNA sequence data of the target sample, perform paired-end sequencing using the Illumina NovaSeq6000 platform, with a read length of 150 bp (PE150); the generated sequencing data should reach no less than 30 M reads per sample to ensure sufficient coverage depth; the preliminary processing steps include removing low-quality sequences, adapter sequences, and possible sequencing contamination sequences, and using the Trim Galore software for cleaning.

[0091] ② The steps for data quality control include:

[0092] Use the FastQC software to perform quality detection on the sequencing data, mainly evaluating the Q30 ratio (should be ≥ 85%), GC content, and base quality distribution; the cleaned data should meet the requirements of downstream analysis, and the qualified data is used for subsequent alignment.

[0093] ③ The steps for calculating gene expression levels include:

[0094] Use the HISAT2 alignment tool to align the cleaned RNA-seq reads to the reference genome Mus musculus GRCm38.p5; perform the alignment using default parameters, and the alignment results include the total alignment rate, unique alignment rate, and non-alignment rate; ensure that the total alignment rate reaches ≥ 85% to guarantee the reliability of data analysis.

[0095] ④ The steps for statistical analysis of gene expression levels include:

[0096] Use the featureCounts software to statistically analyze the expression levels of each gene based on the alignment results; the generated gene expression level data is represented in the form of counts (count), which is used for subsequent differential analysis and functional enrichment analysis.

[0097] ⑤ The steps for data normalization include:

[0098] To remove technical biases between samples and differences in sequencing depth, use the DESeq2 package to perform normalization processing on the count data; the normalized gene expression level data is used for subsequent differential expression analysis, time series analysis, and network construction. Obtain the gene expression matrix.

[0099] The above steps aim to ensure the data quality and the reliability of analysis results throughout the process from RNA-seq raw data to the normalized gene expression matrix, providing an accurate expression data basis for subsequent research.

[0100] (2) Differential analysis:

[0101] Process the raw transcriptome data of the control group and the raw transcriptome data of the experimental group to obtain the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points; perform differential analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain differentially expressed genes;

[0102] The differential analysis includes the following specific steps:

[0103] (2-1) The differential analysis processing steps include:

[0104] Use the count value expression data as input, and use the DESeq2 package to perform differential analysis on the gene expression between different experimental groups to obtain differentially expressed genes; the criteria for screening differentially expressed genes are: log2 Fold Change (log2FC) > 1.5 indicates significantly upregulated genes; log2 Fold Change (log2FC) < -1.5 indicates significantly downregulated genes; adjusted p-value (Padj) < 0.05, which is used to control the False Discovery Rate (FDR); use the PCA method to observe the sample distribution.

[0105] The finally screened differentially expressed genes will be used for subsequent functional enrichment analysis and network construction.

[0106] (2-2) The normalization steps include:

[0107] To reduce the influence of technical noise and experimental batches, perform normalization processing on the data of differentially expressed genes, and use the vst (variance-stabilizing transformation) method of the DESeq2 package to perform variance stability transformation on the expression count data to generate a normalized expression matrix.

[0108] (3) Time series analysis:

[0109] Perform time series analysis on the gene expression matrix of the experimental group at different time points, cluster according to the gene expression trend to obtain a clustering group, and perform functional annotation on the clustering group;

[0110] The time series analysis includes the following specific steps:

[0111] Using the standardized expression data as input, construct an input data framework that includes a gene expression matrix and experimental design information (experimental groups and time points), and perform time series analysis using a two-step regression method (Stepwise Regression): In the first step, fit a global regression model to screen genes related to time changes, and set the screening threshold to FDR < 0.05; in the second step, fit a differential model between experimental groups to identify genes with significant differences between experimental groups. Perform k-means clustering on the selected genes, group the genes according to their expression trends, and obtain clustering groups. Perform GO and KEGG functional enrichment analysis on the genes in each clustering group to reveal the potential biological functions and involved signaling pathways of genes with similar expression trends. Use the ClusterProfiler package for enrichment analysis, and set the screening threshold to FDR < 0.05. Use the built-in function of the maSigPro package to generate dynamic change curves of gene expression trends to visually display the gene expression patterns between different treatment groups.

[0112] (4) Enrichment analysis:

[0113] Perform GSEA (Gene Set Enrichment Analysis) enrichment analysis on the differentially expressed genes to obtain GSEA analysis results;

[0114] Perform GSVA (Gene Set Variation Analysis) enrichment analysis on the differentially expressed genes to obtain GSVA analysis results;

[0115] Perform GO and KEGG enrichment analysis on the differentially expressed genes to obtain GO and KEGG analysis results;

[0116] Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways;

[0117] The enrichment analysis includes the following specific steps:

[0118] The enrichment analysis takes the standardized genes and expression levels as input for GSEA, GSVA, GO, and KEGG analysis.

[0119] (4-1) GSEA (Gene Set Enrichment Analysis) analysis: Based on the standardized gene expression matrix, construct an input framework that includes a gene list and its significance statistical values (log2FC and Padj), and use the GSEA method to evaluate gene set enrichment. Set the screening threshold to FDR < 0.05, and the significant gene sets include genes involved in key biological processes and signaling pathways.

[0120] (4-2) GSVA (Gene Set Variation Analysis) analysis: Using the standardized gene expression data as input, evaluate the enrichment degree of each sample in a specific gene set, and calculate the GSVA enrichment score (ES). Conduct a statistical test on the GSVA scores, and the threshold for screening significant gene sets is FDR < 0.05 to reveal the functional differences between samples in different groups.

[0121] (4-3) GO and KEGG enrichment analysis:

[0122] Perform functional enrichment analysis on the selected genes, and analyze Biological Process (BP), Cellular Component (CC), Molecular Function (MF), and KEGG signaling pathways respectively; construct an input framework containing the gene list and significance statistical values, and use the ClusterProfiler package for enrichment analysis. The screening threshold is FDR < 0.05; classify and statistically analyze the GO and KEGG enrichment results, and extract significant functions and pathways. Use the ClusterProfiler package to generate visualization charts of enrichment analysis (such as bubble charts and bar charts), and combine the gene expression patterns to reveal the dynamic changes of enriched pathways and functions.

[0123] (4-4) Integration:

[0124] Finally, integrate the results of GSEA, GSVA, GO, and KEGG analyses, extract significantly enriched gene sets and pathways, provide a basis for studying the exercise regulation mechanism, and obtain the GO enrichment analysis results and / or KEGG pathway enrichment analysis results.

[0125] (5) Weighted gene co-expression network (WGCNA) analysis:

[0126] Perform weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain WGCNA key genes;

[0127] The weighted gene co-expression network (WGCNA) analysis includes the following specific steps:

[0128] (5-1) Input data construction: Using the standardized gene expression matrix as input, construct a data framework containing the gene expression information of all samples for co-expression network analysis

[0129] (5-2) Soft threshold screening:

[0130] Use the pickSoftThreshold function in the WGCNA package to screen for suitable soft threshold parameters (power); select a soft threshold range of 9 according to the scale-free network characteristic standard (Scale-Free Topology Fit Index) to ensure that the network conforms to the scale-free distribution characteristic;

[0131] (5-3) Network construction and module division:

[0132] Construct a weighted gene co-expression network based on the selected soft threshold, and calculate the weighted correlation coefficient between genes; use the Dynamic Tree Cut algorithm to divide genes into modules, and group genes with similar expression patterns into the same module; set the minimum number of genes in a module to 30, and merge modules with a module similarity > 0.75 to obtain the final co-expression module;

[0133] (5-4) Module and phenotype correlation analysis:

[0134] Calculate the module eigengene (ME) of each module, and perform a correlation analysis with phenotype data such as experimental groups and time points; screen modules significantly correlated with phenotypes, and set the significance threshold to p < 0.05;

[0135] (5-5) Key gene screening:

[0136] In modules significantly correlated with phenotypes, screen through module membership (MM) and gene-phenotype correlation (gene significance, GS); the screening criteria are genes with MM > 0.8 and GS > 0.2, and extract the key genes in the module for subsequent functional verification;

[0137] (6) Machine learning:

[0138] Perform Lasso regression feature screening on the differentially expressed genes, genes in the clustering groups, and WGCNA key genes to obtain Lasso candidate feature genes;

[0139] Perform random forest feature (RF) screening on the Lasso candidate feature genes to obtain RF key genes;

[0140] Take the intersection of WGCNA key genes, Lasso candidate feature genes, and RF key genes, which are the key regulatory molecules promoting brown adipocyte metabolism.

[0141] The machine learning steps include:

[0142] (6-1) Input data construction:

[0143] Based on the gene lists selected by differential analysis, time series analysis, and weighted gene co-expression network analysis, a data framework containing gene expression matrices and experimental group information is constructed for the training and feature screening of machine learning models.

[0144] (6-2) Lasso regression feature screening:

[0145] Lasso regression feature screening includes: First, use the Lasso regression model to screen genes for features. Using the standardized gene expression matrix and the binary classification labels of experimental groups as inputs, determine the optimal value of the regularization parameter λ through 10-fold cross-validation, and select genes with non-zero regression coefficients as candidate feature genes. The goal of the Lasso regression model is to eliminate redundant variables and only retain genes that make important contributions to group discrimination. The screening result is a preliminary set of feature genes that are potentially important for experimental group discrimination.

[0146] (6-3) Random forest feature screening:

[0147] Use the random forest (RF) algorithm to further optimize the candidate genes screened by Lasso. Construct a random forest classification model, calculate the importance scores (Mean Decrease Accuracy and Mean Decrease Gini) of each candidate feature gene, and screen for key genes with an importance threshold > 0.01. The model parameters are set as the number of decision trees (n_estimators) 500, the maximum tree depth (max_depth) 10, and the minimum number of samples for splitting (min_samples_split) 2.

[0148] (6-4) Finally, integrate the screening results of Lasso and random forest to generate a list of key genes, and generate a bar chart of gene importance through the ggplot2 visualization tool.

[0149] Take the intersection of the results of WGCNA, Lasso regression, and RF for further verification.

[0150] 2. External data verification

[0151] (1) External animal experiment data verification

[0152] The animal experiments with the same experimental design as this study using external data were carried out for verification, and it was found that the expression trends of the screened molecules were the same after exercise. The external data is the single-cell sequencing dataset of brown adipose tissue (GSE207706) cited in this article: BURL R B, RONDINI E A, WEI H, et al. Deconstructing cold-induced brown adipocyte neogenesis in mice[J]. eLife, 2022, 11.

[0153] (2) Verification of AP-MS data

[0154] Since the characteristics of the screened molecules are binding proteins, the AP-MS technique was continued to verify that the proteins bound by this molecule are indeed related to the process of exercise improving diseases. The AP-MS data comes from SEGAL D, MAIER S, MASTROMARCO G J, et al. A central chaperone-like role for 14-3-3 proteins in human cells[J]. Molecular cell, 2023, 83(6): 974-93.e15.

[0155] (2-1) Obtain the proteins bound by the screened molecules from the supplementary materials of Segal et al., 2023. After removing low-quality entries, use the GO and KEGG enrichment analyses in step (4-3) to annotate the functions of these proteins, and generate the enrichment results of these proteins through the ggplot2 visualization tool.

[0156] 3. Procedures and verification results

[0157] (1) Figure 1 The results showed that cell culture was successful. The administration verification found that ISO promoted the lipid metabolism of brown adipocytes, but the NS treatment inhibited the promoting effect of ISO. Oil red O staining ( Figure 1 A-C) and glycerol quantification ( Figure 1 D and E) both supported this result. ISO enhanced the function of brown adipocytes by promoting lipid metabolism and reducing lipid droplets in brown adipose tissue, but this process could be blocked by COX2 inhibitors.

[0158] (2) Differential analysis at each time point ( Figure 2 A) and time-series differential gene analysis ( Figure 2 B), KEGG enrichment analysis of time-series differential genes ( Figure 2 C and D) and PCA distribution ( Figure 2E); The number of differentially expressed genes that were upregulated after ISO treatment but downregulated after ISO_NS treatment at the 6h and 12h time points ( Figure 2 F) Perform KEGG enrichment analysis separately ( Figure 2 G and H). The time-series differentially expressed genes were concentrated in stress-related pathways, cytoskeleton pathways, and fatty acid metabolism pathways. The genes downstream of ISO-COX2 were concentrated in cytoskeleton pathways and pathways such as PI3K-AKT.

[0159] (3) Through WGCNA ( Figure 3 A-C), random forest ( Figure 3 D-F), Lasso regression ( Figure 3 G and H), regression to locate the top 10 downstream genes ( Figure 3 F), Venn diagram ( Figure 3 I) identified and screened the most important gene as Ywhah, and the expression trend of the Ywhah gene in each group of brown adipocytes was statistically analyzed ( Figure 3 J). The intersection of the results of multiple algorithms confirmed the robustness of the selected genes, and the Ywhah gene is a key link in exercise promoting brown fat metabolism.

[0160] (4) External chip data verified that exercise increased the expression level of Ywhah in brown adipose tissue ( Figure 3 K), AP-MS data verification ( Figure 4 A-F), the molecule Ywhah selected by this process can also be verified in animal experiments, and the proteins bound by this molecule are closely related to the phenotypes of interest.

[0161] 4. The system used for the method of identifying downstream exercise factors regulating brown fat metabolism

[0162] The system used for the method of identifying downstream exercise factors regulating brown fat metabolism includes:

[0163] Sample processing and acquisition module of transcriptome raw data: Divide brown adipocytes into a control group and an experimental group, administer a blank solvent to the control group (CON) to obtain samples of the control group at different time points after the administration treatment, and administer isoproterenol (ISO) or / and COX2 inhibitor NS-398 (NS) to the experimental group to obtain samples of the experimental group at different time points after the administration treatment;

[0164] Perform RNA second-generation sequencing on the samples of the control group at different time points after the administration treatment and the samples of the experimental group at different time points after the administration treatment to obtain the transcriptome raw data of the control group and the transcriptome raw data of the experimental group;

[0165] Differential analysis module: Process the original transcriptome data of the control group and the original transcriptome data of the experimental group to obtain the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points; perform differential analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain differentially expressed genes at different time points;

[0166] Time series analysis module: Perform time series analysis on the gene expression matrix of the experimental group at different time points, cluster according to the gene expression trend to obtain clusters, and perform functional annotation on the clusters;

[0167] Enrichment analysis module: Perform GSEA (Gene Set Enrichment Analysis) enrichment analysis on the differentially expressed genes to obtain GSEA analysis results;

[0168] Perform GSVA (Gene Set Variation Analysis) enrichment analysis on the differentially expressed genes to obtain GSVA analysis results;

[0169] Perform GO and KEGG enrichment analysis on the differentially expressed genes to obtain GO and KEGG analysis results;

[0170] Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways;

[0171] Weighted gene co-expression network (WGCNA) analysis module:

[0172] Perform weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain WGCNA key genes;

[0173] Machine learning module: Perform Lasso regression feature screening on the differentially expressed genes, the genes in the clusters, and the WGCNA key genes to obtain Lasso candidate feature genes;

[0174] Perform random forest feature (RF) screening on the Lasso candidate feature genes to obtain RF key genes;

[0175] Take the intersection of the WGCNA key genes, the Lasso candidate feature genes, and the RF key genes, which is the key molecule promoting brown adipocyte metabolism.

[0176] 5. A computer-readable storage medium storing a computer program

[0177] A computer-readable storage medium storing a computer program, the computer program causing a computer to execute the steps (2) to (6).

[0178] In summary, the present invention uses the results of RNA second-generation sequencing to identify the downstream metabolic-promoting factors of brown adipocytes. Through differential analysis, enrichment analysis, and time-series analysis, the overall characteristics of the transcriptome are obtained. WGCNA and machine learning (Lasso regression, random forest) are used to further reduce the data and screen out key genes. The present invention helps to find key molecules that can promote the metabolism of brown adipocytes and provides new biological targets for the treatment of obesity.

[0179] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A method for identifying downstream motor factors that regulate brown fat metabolism, characterized in that: The steps include: (1) Sample processing and acquisition of raw transcriptome data: The brown fat cells are divided into a control group and an experimental group, the control group is treated with a blank solvent to obtain samples of the control group at different time points after the treatment, and the experimental group is treated with isoproterenol and / or COX2 inhibitor NS-398 to obtain samples of the experimental group at different time points after the treatment; Performing RNA next generation sequencing on samples of the control group at different time points after the drug administration and samples of the experimental group at different time points after the drug administration to obtain original transcriptome data of the control group and original transcriptome data of the experimental group; (2) Difference analysis: The original transcriptome data of the control group and the original transcriptome data of the experimental group are processed to obtain gene expression matrices of the control group at different time points and gene expression matrices of the experimental group at different time points; the gene expression matrices of the control group and the gene expression matrices of the experimental group are differentially analyzed to obtain differentially expressed genes; (3) Timing analysis: Performing time series analysis on the gene expression matrix of the experimental groups at different time points, clustering according to the gene expression trend to obtain cluster groups, and performing functional annotation on the cluster groups; (4) Enrichment analysis: Performing GSEA enrichment analysis on the differentially expressed genes to obtain GSEA analysis results; Performing GSVA enrichment analysis on the differentially expressed genes to obtain GSVA analysis results; Performing GO and KEGG enrichment analysis on the differentially expressed genes to obtain GO and KEGG analysis results; Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways; (5) Weighted gene co-expression network analysis: Performing weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain WGCNA key genes; (6) Machine Learning: Performing Lasso regression feature screening on the differentially expressed genes, the genes of the cluster group and the WGCNA key genes to obtain Lasso candidate feature genes; The Lasso candidate feature genes are screened by random forest features to obtain RF key genes; The intersection of WGCNA key genes, Lasso candidate characteristic genes and RF key genes is the key regulatory molecule that promotes brown fat cell metabolism.

2. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: In step (1), the experimental groups include ISO group, NS group and ISO-NS combined group; The ISO group was treated with isoproterenol, and the final concentration of the isoproterenol was 10 μM; the NS group was treated with the COX2 inhibitor NS-398, and the final concentration of the COX2 inhibitor NS-398 was 100 μM; the ISO-NS combination group was treated with isoproterenol and the COX2 inhibitor NS-398, and the final concentration of the isoproterenol was 10 μM, and the final concentration of the COX2 inhibitor NS-398 was 100 μM; In step (1), different time points after the drug administration include any time point between 0-24 hours after the drug administration.

3. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: In step (2), the difference analysis was performed using DESeq2 or edgeR software; the FDR threshold was set to be less than 0.05 and fold change ≥ 1.

5.

4. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: The time series analysis was performed using the Masigpro package and a two-step regression method.

5. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 4, characterized in that: The timing analysis includes the following specific steps: (3-1) Constructing an input data framework containing the gene expression matrix of experimental groups at different time points and experimental group information; (3-2) First, fit the global regression model to screen differentially expressed genes related to temporal changes, and set the screening threshold to FDR < 0.05; then fit the differential model between experimental groups to identify differentially expressed genes with significant temporal changes between experimental groups, perform k-means clustering based on the screened differentially expressed genes with significant temporal changes between experimental groups, and group them according to expression trends to obtain cluster groups; (3-3) The ClusterProfiler package was used to perform GO and KEGG functional enrichment analysis on each cluster group, with a screening threshold of FDR < 0.05 for functional annotation.

6. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: In step (4), the GSEA enrichment analysis includes the following specific steps: constructing an input framework including differentially expressed genes and their significant statistical values, using the GSEA method to evaluate the gene set enrichment, setting the screening threshold to FDR<0.05, and obtaining the GSEA analysis results; In step (4), the GSVA enrichment analysis includes the following specific steps: taking differentially expressed genes as input, evaluating the enrichment degree of each group on a specific gene set, calculating the GSVA enrichment score, performing a statistical test on the GSVA enrichment score, screening a significant gene set with a threshold of FDR<0.05, and obtaining a GSVA analysis result; In step (4), the GO and KEGG enrichment analysis includes the following steps: constructing an input framework including differentially expressed genes and significant statistical values, using the ClusterProfiler package for enrichment analysis, and the screening threshold is FDR<0.05; classifying and counting the GO and KEGG enrichment results, extracting significant functions and pathways, and obtaining GO and KEGG analysis results.

7. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: Step (5) includes the following specific steps: (5-1) Input data construction: Using the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points as input, a data framework containing the gene expression information of all groups is constructed; (5-2) Soft threshold screening: Use the pickSoftThreshold function in the WGCNA package to screen suitable soft threshold parameters; select the soft threshold range as 9 according to the scale-free network characteristic standard; (5-3) Network construction and module division: A weighted gene co-expression network was constructed based on the selected soft threshold, and the weighted correlation coefficients between genes were calculated; the genes were divided into modules using the dynamic cutting tree algorithm, and genes with similar expression patterns were divided into the same module; the minimum number of module genes was set to 30, and modules with module similarity > 0.75 were merged to obtain the final co-expression module; (5-4) Analysis of association between co-expression modules and phenotypes: Calculate the characteristic genes of each co-expression module and perform correlation analysis with the phenotypic data of the experimental groups and time points; screen modules that are significantly correlated with the phenotype, with the significance threshold set at p < 0.05; (5-5) Key gene screening: In the modules with significant correlation with phenotypes, screening was performed based on gene module members and the correlation between genes and phenotypes; the screening criteria were genes with MM>0.8 and GS>0.2, and the key genes in the modules were extracted, namely the WGCNA key genes.

8. The method for identifying downstream motor factors that regulate brown fat metabolism according to claim 1, characterized in that: In step (6), the Lasso regression feature screening includes the following specific steps: first, using the Lasso regression model to perform feature screening on the differentially expressed genes, the genes of the cluster group and the WGCNA key genes, with the standardized gene expression matrix and the binary classification labels of the experimental group as input; determining the optimal value of the regularization parameter λ\lambdaλ through 10-fold cross validation, and screening out genes with non-zero regression coefficients as candidate feature genes, namely, Lasso candidate feature genes; In step (6), the random forest feature screening includes the following specific steps: constructing a random forest classification model, calculating the importance score of each Lasso candidate feature gene, and screening key genes with an importance threshold > 0.01; the model parameters are set to 500 decision trees, a maximum tree depth of 10, and a minimum sample split number of 2.

9. A system used in the method for identifying downstream motor factors that regulate brown fat metabolism according to any one of claims 1 to 8, characterized in that: include: Sample processing and transcriptome raw data acquisition module: brown fat cells are divided into a control group and an experimental group, the control group is treated with a blank solvent to obtain samples of the control group at different time points after the treatment, and the experimental group is treated with isoproterenol and / or COX2 inhibitor NS-398 to obtain samples of the experimental group at different time points after the treatment; Performing RNA next generation sequencing on samples of the control group at different time points after the drug administration and samples of the experimental group at different time points after the drug administration to obtain original transcriptome data of the control group and original transcriptome data of the experimental group; Differential analysis module: processing the original transcriptome data of the control group and the original transcriptome data of the experimental group to obtain the gene expression matrix of the control group at different time points and the gene expression matrix of the experimental group at different time points; performing differential analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain differentially expressed genes at different time points; Time series analysis module: perform time series analysis on the gene expression matrix of the experimental group at different time points, cluster according to the gene expression trend to obtain cluster groups, and perform functional annotation on the cluster groups; Enrichment analysis module: performing GSEA enrichment analysis on the differentially expressed genes to obtain GSEA analysis results; Performing GSVA enrichment analysis on the differentially expressed genes to obtain GSVA analysis results; Performing GO and KEGG enrichment analysis on the differentially expressed genes to obtain GO and KEGG analysis results; Integrate the GSEA analysis results, the GSVA analysis results, and the GO and KEGG analysis results to obtain significantly enriched gene sets and pathways; Weighted gene co-expression network analysis module: Performing weighted gene co-expression network analysis on the gene expression matrix of the control group and the gene expression matrix of the experimental group to obtain key genes; Machine learning module: performing Lasso regression feature screening on the differentially expressed genes, the genes of the cluster group and the key genes to obtain Lasso candidate feature genes; The Lasso candidate feature genes are screened by random forest features to obtain RF key genes; The intersection of WGCNA key genes, Lasso candidate characteristic genes and RF key genes is the key regulatory molecule that promotes brown fat cell metabolism.

10. A computer-readable storage medium storing a computer program, characterized in that: The computer program enables a computer to execute steps (2) to (6) according to any one of claims 1 to 8.

Images