Method for identifying common targets of new coronavirus and head and neck squamous cell carcinoma and screening drugs based on bioinformatics

Abstract

The application discloses a method for identifying common targets of COVID-19 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics, and belongs to the technical field of bioinformatics. In view of the problem that existing technologies are difficult to systematically identify common molecular targets of the two diseases, the method acquires a COVID-19 and head and neck squamous cell carcinoma dataset, screens common differentially expressed genes and performs enrichment analysis; constructs a protein interaction network and identifies the top 10 hub genes; constructs a regulatory network of the hub genes and transcription factors and microRNAs; predicts the top 10 candidate drugs using the hub genes; constructs a hub gene and disease association network; identifies an interferon-stimulated gene subgroup based on the enrichment results, the hub genes and the disease association network and constructs a functional polarity score model. The method can systematically identify key molecular targets of the two diseases and screen candidate drugs, and provides an analysis basis for understanding the comorbidity of COVID-19 and head and neck squamous cell carcinoma and drug repositioning.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology, specifically relating to a method for identifying common targets of COVID-19 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics. Background Technology

[0002] COVID-19 is a highly contagious disease caused by the SARS-CoV-2 virus, which has infected hundreds of millions of people and caused millions of deaths worldwide. Head and neck squamous cell carcinoma is one of the most common malignant tumors of the head and neck, and its development is associated with multiple factors, including smoking, alcohol consumption, and human papillomavirus (HPV) infection. Clinical observations have found that patients with head and neck squamous cell carcinoma have a high susceptibility to SARS-CoV-2, and their mortality rate after COVID-19 infection is significantly higher than that of infected individuals without cancer. Previous studies have suggested that the expression level of ACE2 receptors in head and neck squamous cell carcinoma tissues is high, which may be one of the reasons for their susceptibility to COVID-19; at the same time, the systemic immunosuppression caused by chemotherapy and radiotherapy may also exacerbate disease progression after viral infection. However, despite the above clinical associations, whether there are common molecular pathological mechanisms between COVID-19 and head and neck squamous cell carcinoma remains unclear. Existing research has mostly focused on the analysis of single diseases, lacking a systematic analysis of the common transcriptomic changes in both diseases. While bioinformatics methods have been widely used to mine disease-related genes, integrated analysis of differentially expressed genes shared by two diseases is still relatively rare. In particular, how to identify key molecular targets shared by both diseases from complex gene expression data and predict potential drugs that can simultaneously intervene in both diseases based on these targets remains a problem that needs to be solved. Summary of the Invention

[0003] One object of the present invention is to solve at least the above-mentioned problems and to provide at least the advantages that will be described later.

[0004] Another objective of this invention is to provide a bioinformatics-based method for identifying common targets and screening drugs for COVID-19 and head and neck squamous cell carcinoma. This method can systematically identify differentially expressed genes shared by both diseases and analyze their enriched immune response and viral genome replication-related signaling pathways. Based on this, it screens hub genes that serve as key nodes in both diseases, providing a basis for identifying common molecular targets. Furthermore, by constructing a regulatory network of hub genes, transcription factors, and microRNAs, it reveals the common regulatory mechanisms of the two diseases at the transcriptional and post-transcriptional levels. Simultaneously, based on hub genes, it predicts candidate drugs that can be used to treat both COVID-19 and head and neck squamous cell carcinoma, providing theoretical support for drug retargeting. Finally, by analyzing the disease associations related to hub genes, it provides molecular-level clues for understanding the comorbidity mechanisms of the two diseases, thereby achieving a systematic integrated analysis from molecular target identification to potential therapeutic drug screening.

[0005] To achieve these objectives and other advantages of the present invention, a method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics is provided, comprising the following steps: S1: Obtain the transcriptome datasets GSE196822 for COVID-19 and GSE178537 for head and neck squamous cell carcinoma from the Gene Expression Comprehensive Database; S2: Using the DESeq2 package in R software, with a false discovery rate of less than 0.01 and an absolute value of log2 transcriptional fold change greater than 1 as thresholds, the first differentially expressed gene set and the second differentially expressed gene set were screened from the COVID-19 dataset and the head and neck squamous cell carcinoma dataset, respectively. S3: Intersection analysis of the first and second differentially expressed gene sets was performed using the online Venn graph tool Jvenn to obtain the common differentially expressed genes; S4: Gene ontology analysis and pathway enrichment analysis of commonly differentially expressed genes were performed using the enrichment analysis tool Enrichr to identify biological processes and signaling pathways related to immune response and viral genome replication; S5: Construct a protein-protein interaction network of proteins encoded by common differentially expressed genes using the STRING database, and import the network into Cytoscape software for visualization, setting a comprehensive score greater than 0.5 as the screening criterion; S6: Using the Cytohubba plugin in Cytoscape software, the nodes in the protein-protein interaction network were sorted based on the maximum clique centrality algorithm, and the top 10 commonly differentially expressed genes with the highest connectivity were selected as hub genes. The hub genes are OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E and CMPK2. S7: Using the NetworkAnalyst platform, regulatory networks of hub genes and transcription factors, as well as regulatory networks of hub genes and microRNAs, were constructed based on the JASPAR database and the TarBase database, respectively. S8: Based on the drug characterization database DSigDB, protein-drug interaction prediction was performed using the hub gene. The top 10 potential therapeutic drugs were obtained with a p-value of less than 5.20E-07 as the screening criterion. S9: Construct a network linking hub genes and diseases using the NetworkAnalyst platform in conjunction with the DisGeNET database, and analyze diseases related to hub genes.

[0006] S10: Based on the pathway enrichment results of step S4 and the hub gene of step S6, combined with the gene-disease association network of step S9, we identify the interferon-stimulated gene subgroup in the hub gene that has dual functions of antiviral and tumor-promoting, and construct its functional polarity scoring model to analyze the dominant functional direction of the hub gene in different patient groups.

[0007] This invention integrates and analyzes transcriptome data from COVID-19 and head and neck squamous cell carcinoma, employing systems biology methods to systematically identify differentially expressed genes shared by both diseases and their involvement in immune responses and viral genome replication-related signaling pathways. This overcomes the limitation of single-disease analysis in revealing the common molecular mechanisms of the two diseases. Based on this, protein-protein interaction networks and network topology analysis are used to screen for hub genes with the highest connectivity as common molecular targets. Furthermore, by constructing transcription factor and microRNA regulatory networks, the regulatory mechanisms at the transcriptional and post-transcriptional levels are elucidated. Then, based on a drug signature database, potential therapeutic drugs that can simultaneously intervene in both diseases are predicted. Finally, gene-disease association analysis is combined to analyze the disease spectrum related to hub genes, thereby revealing the comorbid association between COVID-19 and head and neck squamous cell carcinoma at the molecular level. This provides a systematic bioinformatics framework and theoretical basis for the identification of common molecular targets, drug retargeting, and research on comorbid mechanisms of the two diseases.

[0008] Preferably, after constructing the protein-protein interaction network in step S5 and before screening the hub gene in step S6, the step of identifying viral entry-related genes as co-control nodes is also included: S501: Based on the COVID-19 dataset GSE196822 and head and neck squamous cell carcinoma dataset GSE178537 obtained in step S1, extract the expression values ​​of known genes related to SARS-CoV-2 virus entry into host cells. SARS-CoV-2 virus entry-related genes include ACE2, TMPRSS2 and CTSL. S502: Construct a protein-protein interaction subnetwork between SARS-CoV-2 virus entry-related genes and the common differentially expressed genes obtained in step S3 using the STRING database. Set the comprehensive score to be greater than 0.7 as the screening condition to identify common differentially expressed genes that have direct physical interactions or functional associations with virus entry-related genes. S503: Based on the minimum driving node set algorithm in network control theory, using the protein-protein interaction network constructed in step S5 as input, calculate the control capability score K of the virus entering the relevant gene on the network module where the hub gene is located. The calculation method of the control capability score K is as follows: K = γ × Nc + δ × (1 / Sp); Where Nc is the number of hub genes directly regulated by the virus entry-related genes, Sp is the shortest path length between the virus entry-related genes and the hub genes, and γ and δ are preset weight coefficients, which can be 0.6 and 0.4 respectively, to balance the contribution of the number of directly regulated nodes and the network path distance to the control capability. S504: Virus entry-related genes with a control ability score K higher than a preset threshold are identified as co-control nodes. Co-control nodes affect the entry efficiency and replication ability of SARS-CoV-2 virus in head and neck squamous cell carcinoma tissue by regulating the expression level of interferon-stimulated genes. S505: Integrate the identified co-control nodes with the hub genes screened in step S6 to construct an expanded hub gene set for subsequent drug prediction and functional polarity analysis.

[0009] This invention, based on the identification of commonly differentially expressed genes and hub genes, further introduces virus entry-related genes as co-control nodes. By extracting genes directly related to SARS-CoV-2 virus entry into host cells, such as ACE2, TMPRSS2, and CTSL, a protein-protein interaction subnetwork between these genes and the commonly differentially expressed genes is constructed. Then, based on the minimum driver node set algorithm in network control theory, the control ability score of virus entry-related genes on the hub gene network is calculated, thereby identifying the co-control nodes that act as a bridge between viral infection and tumor progression. This strategy avoids the limitations of traditional hub-based methods. Gene screening methods, which focus only on network connectivity and overlook key functional nodes, have limitations. This study reveals the molecular mechanism by which virus entry-related genes influence the entry efficiency and replication capacity of SARS-CoV-2 virus in head and neck squamous cell carcinoma tissues by regulating the expression level of interferon-stimulated genes. By integrating virus entry-related genes with hub genes to construct an expanded hub gene set, a more complete set of targets is provided for subsequent drug prediction and functional analysis. This elucidates the dynamic interaction between viral infection and tumor progression at the molecular level and provides a new theoretical perspective for understanding the molecular basis of the high susceptibility of head and neck squamous cell carcinoma patients to COVID-19.

[0010] Preferably, after integrating the identified cooperative control nodes with the hub gene in step S505, the method further includes constructing a virus-host competition index model to quantify the dynamic game relationship between virus entry and host defense. Specific steps include: S506: Based on the virus entry-related genes extracted in step S501 and the hub genes screened in step S6, the genes are divided into the virus entry promotion group and the host defense promotion group according to their biological functions. The virus entry promotion group includes ACE2, TMPRSS2 and CTSL, and the host defense promotion group includes OASL, IFIT3, RSAD2, MX1, IFITM3, LY6E and CMPK2. S507: Based on the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S2, the expression values ​​of each gene in the virus entry promotion group and the host defense promotion group are extracted. Principal component analysis is used to reduce dimensionality and construct the virus entry capability index V and the host defense capability index H, respectively. The calculation formula is V = ∑λ i ×G i H = ∑μ j ×G j G i The expression values ​​of each gene in the viral entry promotion group, G j λ represents the expression values ​​of each gene in the host defense enhancement group. i and μ j The loading coefficient of the first principal component in principal component analysis; S508: Constructing a virus-host competition index C based on the ratio of the virus entry capability index V to the host defense capability index H. VH Its calculation formula is C VH =log2(V / H), when C VH A value greater than 0 indicates that the virus has a dominant ability to enter the body; when C... VH A value less than 0 indicates that the host's defense capabilities are dominant; S509: Using the clinical information from the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S1, patient samples are sorted according to the virus-host competition index C. VH Perform stratification and verify C VH Correlation with patient viral load, disease severity, or tumor stage; S510: The virus-host competition index C VH As a functional output indicator of the collaborative control node, it is used to evaluate the regulatory efficacy of the collaborative control node identified in step S503 in different patient populations, and to include the virus-host competition index C. vh As an auxiliary factor adjusting the weighting coefficients, it is integrated into the functional polarity scoring model. Specifically, it is used to adjust the weighting coefficient W... i During the assignment process, an adjustment factor θ=(1+C) is introduced. VH For hub genes related to antiviral function, the final weight is assigned as W / 2. i ×θ; For hub genes associated with pro-inflammatory and pro-cancer functions, the final weight is assigned as W. i ×(1-θ) enables the functional polarity score to dynamically reflect the influence of the virus-host competition state on the direction of gene function.

[0011] This invention, based on identifying collaborative control nodes and integrating an expanded hub gene set, further constructs a virus-host competition index model. Virus entry-related genes and host defense-related genes are categorized into virus entry-promoting and host defense-promoting groups according to their biological functions. Principal component analysis is used to reduce dimensionality and construct virus entry capability and host defense capability indices, respectively. The logarithmic ratio of these two indices quantifies the dynamic game relationship between virus entry and host defense. This design overcomes the limitation that single gene expression indicators cannot accurately reflect the complex balance between viral infection and host immunity. It can intuitively determine the relative advantage of viral entry capability and host defense capability based on the sign and magnitude of the competition index. Furthermore, by combining clinical information to stratify patient samples, the correlation between this index and viral load, disease severity, or tumor stage is verified. By using this index as a functional output indicator of collaborative control nodes and a regulating weight factor in the functional polarity scoring model, a leap from static gene screening to dynamic functional balance analysis is achieved. This provides a quantitative basis for revealing the dominant direction of viral entry and immune defense in different patient groups and guiding subsequent personalized intervention strategies.

[0012] Preferably, the hub genes screened in step S6 are all interferon-stimulated genes. After identifying the hub genes, the step further includes constructing a hub gene functional polarity scoring model. S601: Based on the interferon signaling pathway, antigen processing and presentation pathway and inflammatory response pathway obtained from the pathway enrichment analysis in step S4, a set of regulatory genes related to all three pathways is extracted from the common differentially expressed genes. Among them, the common differentially expressed genes that are simultaneously enriched in the interferon signaling pathway, antigen processing and presentation pathway and inflammatory response pathway are selected as the set of regulatory genes, with the adjusted p value less than 0.05 as the threshold. S602: Using the NetworkAnalyst platform and the JASPAR database, a regulatory network of regulatory gene sets and upstream transcription factors was constructed to identify key transcription factors that simultaneously regulate at least two hub genes and are related to the interferon signaling pathway. Key transcription factors include STAT1 and IRF family transcription factors. S603: Based on the expression levels of key transcription factors in the COVID-19 dataset and head and neck squamous cell carcinoma dataset, combined with the expression value of the hub gene, a functional polarity scoring model was constructed. The calculation formula for the scoring model is as follows: P = ∑(Wi×Ei) / ∑Ej; Where Ei is the expression value of a single hub gene, Ej is the sum of the expression values ​​of all hub genes, and Wi is the weighting coefficient, which is based on the co-expression correlation between hub genes and key transcription factors, the frequency of literature reports, and the virus-host competition index C. VHA comprehensive value was assigned; the weighting coefficients were assigned based on the correlation between the co-expression of the hub gene and key transcription factors, as well as the frequency of literature reports on the bidirectional function of the hub gene in antiviral response and pro-tumor inflammatory response. S604: Using the median of the functional polarity score P as the threshold, the patient sample was divided into an antiviral dominant subgroup and a pro-inflammatory and pro-cancer dominant subgroup. S605: Validate the clinical significance of subgroup division through survival analysis or disease severity scoring to determine the dominant functional direction of the hub gene in a specific patient population.

[0013] This invention, based on the identification that all hub genes are interferon-stimulated genes, further constructs a functional polarity scoring model. Through pathway enrichment analysis, it screens for regulatory gene sets simultaneously enriched in the interferon signaling pathway, antigen processing and presentation pathway, and inflammatory response pathway. Furthermore, it utilizes a transcription factor database to identify key transcription factors that simultaneously regulate at least two hub genes and are associated with the interferon signaling pathway. Then, based on hub gene expression values ​​and weighted by the co-expression correlation of key transcription factors and the frequency of literature reports, a functional polarity scoring formula is constructed. Using the median score as a threshold, patient samples are divided into an antiviral dominant subgroup and a pro-inflammatory / pro-cancer dominant subgroup. This design addresses the challenge of distinguishing the functional polarity of interferon-stimulated genes in COVID-19 and head and neck squamous cell carcinoma. Specifically, it overcomes the limitation of conventional differential expression analysis in identifying the dual functional characteristics of the same gene—which may exert antiviral protective effects or pro-inflammatory and pro-cancer damaging effects in different patient groups—by integrating transcriptional regulatory network analysis, co-expression relationships, and literature mining. This approach quantifies the dual function of the hub gene, achieving a leap from static gene expression profiles to dynamic functional polarity assessment. It provides quantifiable classification criteria and clinical validation pathways for accurately differentiating patients' immune-inflammatory states and guiding personalized intervention strategies.

[0014] Preferably, after screening the top 10 potential therapeutic drugs in step S8, the method further includes the step of constructing a drug-subgroup matching screening model: S801: Based on the functional polarity scoring model constructed in step S603, COVID-19 or head and neck squamous cell carcinoma patient samples were divided into antiviral dominant subgroups and pro-inflammatory and pro-cancer dominant subgroups, and the average expression profile of hub genes in each subgroup was obtained. S802: Using the DSigDB database in step S8, extract the known target gene set of the top 10 potential therapeutic drugs, and supplement the drug mechanism information based on the DrugBank and ChEMBL databases, classifying the drugs into immune activation, inflammation suppression or broad-spectrum regulation types. S803: Construct a co-expression network of drug target genes and hub genes using the NetworkAnalyst platform, and calculate the first correlation coefficient R1 between each drug target set and the characteristics of hub genes in the antiviral dominant subgroup, and the second correlation coefficient R2 between each drug target set and the characteristics of hub genes in the pro-inflammatory and pro-cancer dominant subgroup. S804: Construct a drug-subgroup matching index M based on the ratio of the first correlation coefficient R1 and the second correlation coefficient R2. The calculation formula is M = (R1 / R2) × log2(1 + the enrichment significance value of the drug in the subgroup-related pathway); use the negative logarithm (-log10(p-value)) of the enrichment p-value of the drug target gene set in the subgroup-related pathway as the enrichment significance value. S805: Using the threshold of the matching index M as the screening criterion, the top 10 potential therapeutic drugs are mapped to their corresponding subgroups, generating a list of candidate drugs for antiviral dominant subgroups and a list of candidate drugs for pro-inflammatory and pro-cancer dominant subgroups, thus realizing differentiated drug recommendations based on functional polarity.

[0015] This invention, based on the construction of a functional polarity scoring model and the division into antiviral dominant subgroups and pro-inflammatory / pro-cancer dominant subgroups, further constructs a drug-subgroup matching screening model. It extracts the average expression profile of hub genes in each subgroup as a feature, combines this with a drug feature database to obtain the known target gene sets of top-ranked potential therapeutic drugs, and classifies them according to their mechanisms of action as immune-activating, inflammation-suppressing, or broad-spectrum regulatory types. Then, it constructs a co-expression network of drug targets and hub genes, calculates the first correlation coefficient between each drug target set and the hub gene characteristics of the antiviral dominant subgroup, and the second correlation coefficient between each drug target set and the hub gene characteristics of the pro-inflammatory / pro-cancer dominant subgroup, and compares these two coefficients. By combining the significance of drug enrichment in subgroup-related pathways, a drug-subgroup matching index is constructed. Based on the matching index threshold, candidate drugs are mapped to corresponding subgroups, generating lists of candidate drugs for antiviral dominant subgroups and lists of candidate drugs for pro-inflammatory and pro-cancer dominant subgroups. This design overcomes the risk of poor efficacy or immune imbalance caused by the indiscriminate application of the same batch of candidate drugs to all patients in traditional drug screening methods. By integrating drug target information, hub gene subgroup characteristic expression profiles, and pathway enrichment data, a precise matching relationship between drugs and patients' immune-inflammatory status is constructed, realizing differentiated drug recommendations based on functional polarity and providing a clear subgroup-oriented strategy for personalized clinical medication.

[0016] Preferably, after generating the subgroup candidate drug list in step S805, the process further includes constructing a dynamic perturbation model of the hub gene network under drug intervention to evaluate the polarity stability of the candidate drugs, specifically including the following steps: S806: Based on the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S2, extract the expression matrix of all hub genes, and use the regulatory gene set constructed in step S601 and the key transcription factors identified in step S602 to construct a gene regulatory network containing hub genes, key transcription factors and their interactions. The network nodes include OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, CMPK2 and STAT1 and IRF family transcription factors; S807: Integrate the drug target information from the DSigDB database in step S8 and the drug mechanism information supplemented in step S802 through the NetworkAnalyst platform, and map the drugs in the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list in step S805 to the gene regulatory network respectively, and identify the direct target genes of each drug. S808: Based on the structural characteristics of gene regulatory networks, the control centrality algorithm in network control theory is used to calculate the control ability score C of each candidate drug on the hub gene network. The calculation method of the control ability score C is as follows: C = α × D + β × B; Where D is the degree centrality of the drug target node, B is the betweenness centrality of the drug target node, and α and β are preset weight coefficients used to measure the direct impact of the drug target on the network information flow and its bridging role. S809: By simulating drug intervention conditions, the expression of drug target nodes in the gene regulatory network is removed or inhibited. The change in the expression profile of the hub gene after intervention is calculated using a linear differential equation model. Based on the functional polarity scoring model in step S603, the polarity score P´ after intervention is recalculated to obtain the polarity drift ΔP=|P´-P|. S810: Combine the comprehensive control ability score C and polarity drift ΔP to construct the drug polarity stability index S, which is calculated as S=C / (1+ΔP). Using the stability index S as the screening criterion, the list of candidate drugs for the subgroup in step S805 is reordered, and drugs with high S values ​​are given priority, that is, drugs with strong control over the hub gene network and less likely to cause polarity reversal after intervention. S811: GO and KEGG enrichment analyses were used to verify the extent of pathway perturbation in the hub gene network after intervention with candidate drugs, and candidate drugs that may cause significant bypass effects were excluded.

[0017] This invention, based on generating candidate drug lists for antiviral and pro-inflammatory / pro-cancer dominant subgroups, further constructs a dynamic perturbation model of the hub gene network. By extracting the expression matrix of hub genes and combining it with regulatory gene sets and key transcription factors, a gene regulatory network containing hub genes, key transcription factors, and their interactions is constructed. Candidate drugs for each subgroup are mapped to this network, and direct target genes are identified. Based on the control centrality algorithm in network control theory, the control ability score of each candidate drug on the hub gene network is calculated to measure the direct impact and bridging role of drug targets on network information flow. Then, by simulating drug intervention conditions, the expression of drug target nodes in the gene regulatory network is removed or inhibited. The changes in hub gene expression profiles after intervention are calculated using a linear differential equation model, and the polarity score after intervention is recalculated based on a functional polarity scoring model to obtain the polarity drift and a comprehensive control ability score. A drug polarity stability index was constructed using polarity drift. This index was then used as a screening criterion to reorder the candidate drug list for subgroups, prioritizing drugs with strong control over the hub gene network and low likelihood of polarity reversal after intervention. Enrichment analysis was used to verify the pathway perturbation range of the hub gene network after intervention, excluding candidates that might cause significant bypass effects. This design overcomes the limitations of traditional drug screening methods that only focus on drug-target binding ability and cannot predict the trend of network functional polarity evolution after drug intervention. By constructing a regulatory network including hub genes and key transcription factors, and combining network control theory and dynamic perturbation simulation, the control ability of drugs on the network and the risk of polarity drift after intervention are quantified. This allows for the screening of candidate drugs that accurately match subgroups and possess network stability, effectively avoiding treatment contradictions caused by unexpected disruption of the immune-inflammatory balance due to drug targeting.

[0018] Preferably, after constructing the association network between hub genes and diseases and analyzing hub gene-related diseases in step S9, the method further includes a step of performing disease association stratification based on functional polarity subgroups: S901: Based on the antiviral dominant subgroup and pro-inflammatory and pro-cancer dominant subgroups identified in step S604, the expression matrix of hub genes in each subgroup of patient samples is extracted from the COVID-19 dataset GSE196822 and the head and neck squamous cell carcinoma dataset GSE178537 obtained in step S1. The hub genes include OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E and CMPK2. S902: Using the NetworkAnalyst platform in conjunction with the DisGeNET database, obtain association data between hub genes and all diseases. This association data includes a score for each gene-disease pair. g,d And the type of evidence, Scoreg,d Reflects the strength of the association between gene g and disease d; S903: Calculate the weighted expression level (E) of each hub gene in the antiviral dominant subset and the pro-inflammatory / pro-cancer dominant subset, respectively. g,sub The weighted expression level is the product of the mean expression value of the hub gene in all samples within the subpopulation and the sample size of the subpopulation, i.e., E. g,sub =mean(expression) g,sub )×N sub To reflect the overall activity of hub genes in the subpopulation and the contribution of subpopulation size to disease association; S904: Based on weighted expression level E g,sub Gene-Disease Association Score g,d Calculate the comprehensive association degree A of each disease d in a specific subgroup. d,sub The calculation formula is A d,sub =∑(E d,sub ×Score g,d ), summate and iterate through all hub genes to obtain the comprehensive association degree of each disease in the antiviral dominant subgroup and the pro-inflammatory and pro-cancer dominant subgroup; S905: Overall correlation between the dominant antiviral subset and the dominant pro-inflammatory / pro-cancer subset, respectively. d,sub The diseases were sorted in descending order, and the top 20 were selected as characteristic diseases of each subgroup, forming the characteristic disease set D of the antiviral dominant subgroup. anti Disease set D with pro-inflammatory and pro-cancer dominant subgroups pro ; S906: Construct a disease-disease similarity network using a medical subject thesaurus and calculate D. anti and D pro The Jaccard similarity coefficient is used to identify diseases shared by two subgroups as shared diseases, and diseases unique to each subgroup as specific diseases. S907: Pathway enrichment analysis was performed on shared and specific diseases. The Enrichr tool was used to compare the differences in functional enrichment of hub genes in shared and specific diseases, and the KEGG and Reactome pathways enriched in each disease set were obtained, thereby revealing the molecular basis of the differential risk of comorbidity between COVID-19 and head and neck squamous cell carcinoma in patients of different polarity subgroups.

[0019] This invention, based on functional polarity scoring to classify antiviral and pro-inflammatory / pro-cancer dominant subgroups, further performs hierarchical analysis of the association network between hub genes and diseases. By extracting the expression matrix of hub genes from patient samples in each subgroup, and combining gene-disease association scores and evidence types from the DisGeNET database, the weighted expression level of each hub gene in a specific subgroup is calculated to reflect the overall activity of hub genes and the contribution of subgroup size to disease association. Then, based on the product of the weighted expression level and the gene-disease association score, the comprehensive association degree of each disease in a specific subgroup is constructed. Diseases ranking high in both the antiviral and pro-inflammatory / pro-cancer dominant subgroups are selected as characteristic diseases for each subgroup, forming two sets of characteristic disease sets for each subgroup. These sets are then analyzed using medical subject terms. By constructing a disease-disease similarity network and calculating the similarity coefficients of characteristic disease sets for two subgroups, the system identifies diseases shared by the two subgroups as well as diseases unique to each subgroup. Pathway enrichment analysis is then performed on shared and specific diseases, comparing the functional enrichment differences of hub genes across different disease sets. This design addresses the problem of simply applying disease association results from the overall population to all patients, thus ignoring subgroup-specific comorbidity risks. By deeply integrating functional polarity subgroup segmentation with gene-disease association analysis, it achieves a leap from the population average level to the subgroup-specific level in disease association analysis. This reveals the molecular basis for the differential comorbidity risk of COVID-19 and head and neck squamous cell carcinoma in patients of different polarity subgroups, providing a new analytical pathway for understanding the heterogeneity of comorbidity mechanisms between these two diseases in patient subgroups.

[0020] Preferably, after obtaining the subgroup characteristic disease set in step S907, the method further includes constructing a drug-subgroup disease conflict assessment model to screen candidate drugs that have no potential conflict with the subgroup characteristic diseases, as follows: S908: Based on the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list generated in step S805, or the drug list after reordering in step S810, extract the name and structural information of the candidate drugs in each list. S909: Obtain known association data between each candidate drug and the disease by comparing the Comparative Toxicogenomics Database. This association data includes the direction of the drug-disease association (therapeutic or induced / exacerbated), the strength of the association, and the type of evidence. Construct a drug-disease association matrix M. drug-disease The element represents the association score between the drug and the disease, with positive values ​​indicating therapeutic effects and negative values ​​indicating harmful effects; S910: Disease set D based on the antiviral dominant subgroup characteristics obtained in step S906. anti Disease set D with pro-inflammatory and pro-cancer dominant subgroupspro Calculate the conflict risk score R between each candidate drug and the corresponding subgroup of characteristic disease sets. conflict The formula for calculating the conflict risk score is: In the formula, The score represents the harmful association score between the drug and the corresponding characteristic disease. ω represents the therapeutic association score between the drug and the corresponding characteristic disease. d The comprehensive correlation degree A of disease d in the corresponding subgroup calculated in step S904 is... d,sub Normalized weights to reflect the importance of characteristic diseases to subgroups; S911: Based on the conflict risk score R conflict The candidate drugs are sorted in ascending order, and a threshold T is set. conflict The score R is between 0 and 0.3, excluding conflict risk. conflict Greater than threshold T conflict The drug; S912: Combine the drugs that have undergone conflict screening with the polarity stability index S from step S810 to construct a comprehensive drug recommendation index F. The formula for calculating the comprehensive recommendation index is: In the formula R max The drugs are ranked by F-score based on the maximum conflict risk score among all candidate drugs, generating an optimized list of candidate drugs for each subgroup. S913: Verify the reliability of the association between drugs in the optimization list and subgroup characteristic diseases through literature mining, and verify the binding ability of drugs to hub genes and key transcription factors using molecular docking simulation to ensure the practicality of the screening results.

[0021] This invention, based on generating optimized candidate drug lists for each subgroup and constructing a drug polarity stability index, further constructs a drug-subgroup disease conflict assessment model. It extracts the names and structural information of candidate drugs from each subgroup's candidate drug list, and uses a comparative toxicogenomics database to obtain known association data between each candidate drug and disease, including the direction, strength, and type of evidence of the drug-disease association. A drug-disease association matrix is ​​constructed, using positive and negative values ​​to distinguish between therapeutic and adverse effects. Then, based on the characteristic disease sets of the antiviral dominant subgroup and the pro-inflammatory / pro-cancer dominant subgroup, a conflict risk score is calculated between each candidate drug and the corresponding subgroup's characteristic disease set. This score is obtained by weighted summing the differences between the adverse association score and the therapeutic association score for each disease in the characteristic disease set. The weights are determined after normalizing the comprehensive association degree of each disease in the corresponding subgroup to reflect the importance of the characteristic disease to the subgroup. Candidate drugs are then ranked in ascending order based on the conflict risk score. This approach prioritizes drugs with high risk by setting thresholds for exclusion, then integrates conflict-screened drugs with a polarity stability index to construct a comprehensive drug recommendation index. The drugs are then ranked to generate optimized candidate drug lists for each subgroup. Literature mining is used to verify the reliability of the associations between the optimized drugs and subgroup-specific diseases, and molecular docking simulations are employed to verify the binding ability of drugs to hub genes and key transcription factors. This design overcomes the limitations of traditional drug screening methods that focus solely on the efficacy of drugs for the target disease while neglecting their potential impact on comorbidities in patient subgroups. By integrating a drug-disease association database with a set of subgroup-specific diseases, a conflict risk score is constructed to quantify the harmful or therapeutic effects of drugs on important subgroup diseases. Combined with a polarity stability index for comprehensive screening, this approach achieves precise drug recommendations that balance efficacy and safety, effectively avoiding the risk of unintentionally exacerbating or inducing subgroup-related comorbidities while treating COVID-19 and head and neck squamous cell carcinoma.

[0022] The present invention has at least the following beneficial effects: This invention provides a method that enables a systematic integrated analysis from transcriptome data to potential therapeutic agents. For the first time, it establishes a molecular bridge centered on interferon-stimulated genes between COVID-19 and head and neck squamous cell carcinoma. By identifying ten common hub genes, including OASL, IFIT3, and RSAD2, and their involvement in immune responses and viral genome replication-related pathways, it provides a novel molecular perspective for understanding the comorbidity mechanisms of these two diseases. This method overcomes the limitations of traditional differential expression analysis, which focuses only on a single disease. By constructing a virus-host competition index model and a functional polarity scoring model, it achieves for the first time a dynamic quantitative assessment of the bidirectional function of interferon-stimulated genes—antiviral protection and pro-inflammatory / pro-cancer damaging effects. This allows for precise segmentation of patients into antiviral-dominant and pro-inflammatory / pro-cancer-dominant subgroups, providing a clear molecular subtyping basis for personalized intervention. This method further constructs a drug-subgroup matching index, a drug polarity stability index, and a drug-subgroup disease conflict assessment model, forming a complete technical chain from target identification, subgroup classification, drug screening to safety assessment. It can effectively avoid the potential risk of comorbidity of subgroup characteristics while ensuring the efficacy of drugs for the target disease, and realize personalized drug recommendations that take into account accuracy, stability and safety. It provides a systematic bioinformatics solution for drug repositioning and precision treatment in complex disease comorbidity scenarios.

[0023] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0024] Figure 1 The image shows the screening results of differentially expressed genes in COVID-19 and head and neck squamous cell carcinoma as described in this invention; where A is a volcano plot of the COVID-19 dataset, B is a volcano plot of the head and neck squamous cell carcinoma dataset, and C is the Venn diagram intersection analysis results of differentially expressed genes in the two diseases, showing 139 common differentially expressed genes. Figure 2 This is a graph showing the results of gene ontology enrichment analysis of common differentially expressed genes described in this invention; where A represents the enrichment results of biological processes, B represents the enrichment results of cellular components, and C represents the enrichment results of molecular functions. The length of the bar graph represents the number of enriched genes, and the color intensity represents the significance level. Figure 3 This is a graph showing the results of pathway enrichment analysis of common differentially expressed genes described in this invention; where A represents the enrichment results from the BioPlanet database, B represents the enrichment results from the Reactome database, and C represents the enrichment results from the KEGG database. The length of the bar chart represents the number of enriched genes, and the color intensity represents the significance level. Figure 4This is a protein-protein interaction network diagram of the proteins encoded by the common differentially expressed genes described in this invention; the network was constructed using the STRING database and contains 69 nodes and 169 edges, with the color intensity of the nodes indicating the number of interacting proteins; Figure 5 This is a diagram showing the hub gene identification results described in this invention; the orange nodes in the diagram represent the top 10 hub genes in terms of connectivity, including OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, and CMPK2. The network contains 13 nodes and 56 edges. Figure 6 This is a regulatory network diagram of hub genes and transcription factors described in this invention; blue nodes represent transcription factors, pink nodes represent hub genes, and the network contains 56 nodes and 103 edges. Figure 7 This is a network diagram showing the association between the hub gene and the disease described in this invention; square nodes represent diseases, circular nodes represent hub genes, and the color intensity of the nodes indicates the association strength. Detailed Implementation

[0025] The present invention will now be described in further detail with reference to the accompanying drawings, so that those skilled in the art can implement it based on the description.

[0026] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0027] It should be noted that, unless otherwise specified, the experimental methods described in the following implementation plan are all conventional methods, and the reagents and materials described are all commercially available unless otherwise specified.

[0028] Current research on the association between COVID-19 and head and neck squamous cell carcinoma mainly focuses on clinical observation, such as the susceptibility of head and neck squamous cell carcinoma patients to SARS-CoV-2 due to high ACE2 expression, or the increased risk of viral infection due to immunosuppression caused by chemotherapy and radiotherapy. A few molecular-level studies are limited to individual reports on single signaling pathways such as Eph / ephrin or P2X7 receptor. There is no systematic research that integrates and analyzes the differentially expressed genes, signaling pathways, regulatory networks and potential therapeutic drugs shared by the two diseases at the transcriptome level. There is a lack of methods to systematically compare the two diseases within the same molecular framework.

[0029] To address this technical challenge, this invention provides a bioinformatics-based method for identifying common targets and screening drugs for COVID-19 and head and neck squamous cell carcinoma. The specific implementation is as follows: First, the transcriptome dataset GSE196822 for COVID-19, containing peripheral blood samples from 9 healthy controls and 34 COVID-19 patients, and the transcriptome dataset GSE178537 for head and neck squamous cell carcinoma, containing paired primary tumor tissue and adjacent normal tissue samples from 21 patients, are obtained from the gene expression comprehensive database. Using the DESeq2 package in R software, with a false discovery rate of less than 0.01 and an absolute value of the log2 transcriptional difference greater than 1 as statistical thresholds, a first differentially expressed gene set and a second differentially expressed gene set are selected from the above two datasets, respectively. Subsequently, the two differentially expressed gene sets are uploaded to the online Venn graph tool Jvenn for intersection analysis to obtain the common differentially expressed gene set. To elucidate the biological functions of these commonly differentially expressed genes, Enrichr was used for gene ontology and pathway enrichment analysis, focusing on identifying biological processes and signaling pathways related to immune responses and viral genome replication. Based on this, a protein-protein interaction network encoding proteins from these commonly differentially expressed genes was constructed using the STRING database, with a composite score greater than 0.5 as the screening criterion. This network was then visualized using Cytoscape software. Using the Cytohubba plugin in Cytoscape, nodes in the protein-protein interaction network were ranked based on the maximum clique centrality algorithm, and the top 10 commonly differentially expressed genes with the highest connectivity were selected as hub genes. These 10 hub genes were identified as OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, and CMPK2. To reveal the upstream regulatory mechanisms of these hub genes, regulatory networks between hub genes and transcription factors, and between hub genes and microRNAs, were constructed using the NetworkAnalyst platform based on the JASPAR and TarBase databases, respectively. Based on the drug characterization database DSigDB, protein-drug interactions were predicted using hub genes. A p-value less than 5.20E-07 was used as the screening criterion to identify the top 10 potential therapeutic drugs. Finally, a disease association network between hub genes and the DisGeNET database was constructed using the NetworkAnalyst platform to analyze diseases associated with hub genes. Building upon this, based on pathway enrichment results and hub gene information, combined with the gene-disease association network, interferon-stimulated gene subpopulations with dual antiviral and pro-tumor functions were identified within hub genes. A functional polarity scoring model was constructed to analyze the dominant functional orientation of hub genes in different patient populations.

[0030] Compared to existing research methods that target only a single disease or individual pathways, this invention, for the first time, places COVID-19 and head and neck squamous cell carcinoma within a unified bioinformatics analysis framework. Through systematic differentially expressed gene analysis, protein-protein interaction network construction, and hub gene identification, it screens out key molecular targets shared by the two diseases. Furthermore, it uses a drug-gene interaction database to predict drugs that may be effective for both diseases. Simultaneously, it constructs a functional polarity scoring model to analyze the bidirectional function of interferon-stimulated genes. This achieves a systematic integrated analysis from molecular target identification to potential therapeutic drug screening, providing a complete analytical method for understanding the comorbidity mechanism of the two diseases and developing joint treatment strategies.

[0031] In existing technologies, methods for identifying common molecular targets in COVID-19 and head and neck squamous cell carcinoma primarily rely on constructing protein-protein interaction networks of differentially expressed genes and then directly using network topology algorithms to screen for hub genes with the highest connectivity. While these methods can identify core nodes in the network, their screening logic is entirely based on the network's internal connectivity density. This fails to effectively capture key nodes that, although not highly connected within the network, play a bridging role between viral infection and tumor progression. For example, genes related to SARS-CoV-2 virus entry into host cells, such as ACE2, TMPRSS2, and CTSL, are often overlooked in traditional hub gene screening methods due to their functional associations with interferon-stimulated genes and their regulatory capacity on the hub gene network.

[0032] To address the aforementioned issues, this invention adds a step of identifying virus entry-related genes as collaborative control nodes after constructing the protein-protein interaction network in step S5 and before screening hub genes in step S6. Specifically, firstly, based on the acquired COVID-19 dataset GSE196822 and head and neck squamous cell carcinoma dataset GSE178537, the expression values ​​of three genes directly related to SARS-CoV-2 virus entry into host cells—ACE2, TMPRSS2, and CTSL—are extracted. Subsequently, a protein-protein interaction sub-network is constructed using the STRING database between these three virus entry-related genes and the commonly differentially expressed genes obtained in step S3. A comprehensive score greater than 0.7 is set as the screening criterion to identify commonly differentially expressed genes that have direct physical interactions or functional associations with virus entry-related genes. Based on this, using the minimum driving node set algorithm in network control theory, and taking the protein-protein interaction network constructed in step S5 as input, a score is calculated on the control ability of virus entry-related genes over the network module containing hub genes. This score comprehensively considers the number of hub genes directly regulated by virus entry-related genes and the shortest path length between them. Virus entry-related genes with control ability scores exceeding a preset threshold were identified as co-control nodes. These co-control nodes influence the entry efficiency and replication capacity of SARS-CoV-2 virus in head and neck squamous cell carcinoma tissues by regulating the expression levels of interferon-stimulated genes. Finally, the identified co-control nodes were integrated with the hub genes screened in step S6 to construct an expanded hub gene set for subsequent drug prediction and functional polarity analysis.

[0033] Compared to existing methods that rely solely on network connectivity to screen hub genes, this invention introduces virus entry-related genes and calculates their control over the hub gene network based on network control theory. This effectively identifies synergistic control nodes that act as a bridge between viral infection and tumor progression. It integrates virus entry-related genes, such as ACE2, TMPRSS2, and CTSL, which are difficult to include in the core analysis scope by traditional methods, with interferon-stimulated genes into an expanded hub gene set. This allows for a more complete understanding of the dynamic interaction between viral entry, host defense, and tumor progression at the molecular level, providing a new analytical pathway for understanding the molecular basis of the high susceptibility of head and neck squamous cell carcinoma patients to COVID-19.

[0034] In existing technologies, bioinformatics analysis methods for the association between COVID-19 and head and neck squamous cell carcinoma, whether traditional differentially expressed gene analysis or the aforementioned methods for identifying co-control nodes, focus primarily on identifying the expression changes and network positions of individual genes or nodes in the two diseases. They lack quantitative methods for understanding the dynamic interplay between viral entry-related genes and host defense-related genes. The balance between viral entry capability and host defense capability directly affects viral infection efficiency and tumor microenvironment characteristics. However, existing methods can only provide the expression levels of individual genes and cannot comprehensively reflect the overall competitive dynamics between these two groups of genes with opposing functions. This results in a lack of quantitative evidence for determining the dominant direction of viral infection and immune defense at the individual patient level.

[0035] To address the aforementioned issues, this invention, after integrating the identified cooperative control nodes with the hub genes in step S505, adds a step of constructing a virus-host competition index model. Specifically, based on the virus entry-related genes extracted in step S501 and the hub genes screened in step S6, the genes are first divided into a virus entry promotion group and a host defense promotion group according to their biological functions. The virus entry promotion group includes three genes directly related to SARS-CoV-2 virus entry into host cells: ACE2, TMPRSS2, and CTSL. The host defense promotion group selects interferon-stimulated genes with clear antiviral functions from the hub genes, including OASL, IFIT3, RSAD2, MX1, IFITM3, LY6E, and CMPK2. Subsequently, based on the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S2, the expression values ​​of each gene in the virus entry promotion group and the host defense promotion group were extracted. Principal component analysis (PCA) was used for dimensionality reduction to construct the virus entry capability index and the host defense capability index, respectively. The virus entry capability index was obtained by multiplying the expression values ​​of each gene in the virus entry promotion group by their loading coefficients on the first principal component of PCA and then summing the results; the host defense capability index was obtained similarly. A virus-host competition index was constructed based on the ratio of the virus entry capability index to the host defense capability index. A value greater than 0 indicates that virus entry capability is dominant, while a value less than 0 indicates that host defense capability is dominant. Using clinical information from the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S1, patient samples were stratified according to the virus-host competition index to verify the correlation between this index and patient viral load, disease severity, or tumor stage. Finally, the virus-host competition index is used as the functional output index of the co-control node to evaluate the regulatory efficacy of the co-control node identified in step S503 in different patient groups. This index is also used as an auxiliary factor to adjust the weight coefficient and integrated into the subsequent functional polarity scoring model.

[0036] Compared to existing technologies that focus only on the expression level of a single gene or network connectivity, this invention, by constructing a virus-host competition index model, achieves for the first time a comprehensive quantitative and dynamic numerical representation of the relationship between viral entry capability and host defense capability. This index can intuitively reflect the dominant direction of viral entry and immune defense in different patient groups. Its correlation with disease severity is verified through clinical information, and the index is used as a weighting factor in subsequent functional polarity scoring models. This represents a leap from static gene screening to dynamic functional balance analysis, providing a quantitative analytical basis for revealing the heterogeneity of the virus-host competition situation in different patient groups and guiding subsequent personalized intervention strategies.

[0037] In existing technologies, whether it's traditional differentially expressed gene analysis, the aforementioned methods for identifying synergistic control nodes, or the methods for constructing virus-host competition index models, all focus on identifying the overall competitive landscape among commonly differentially expressed genes, key nodes, or functional groups. However, none of these methods address a deeper issue: when all hub genes are interferon-stimulated genes, it's impossible to distinguish whether these genes primarily exert antiviral protective effects or pro-inflammatory and pro-cancer damaging effects in specific patients—that is, their functional polarity cannot be differentiated using existing methods. The bidirectional function of interferon-stimulated genes—antiviral and pro-tumorigenic—may exhibit diametrically opposed dominant directions in individual patients due to differences in the tumor microenvironment, immune status, and stage of viral infection. Traditional methods can only provide overall expression levels or network topological positions, failing to analyze this individual heterogeneity of functional polarity.

[0038] To address the aforementioned issues, this invention, based on the previously described method, adds a step of constructing a hub gene functional polarity scoring model after screening hub genes. Specifically, firstly, based on pathway enrichment analysis of the interferon signaling pathway, antigen processing and presentation pathway, and inflammatory response pathway, a set of regulatory genes related to all three pathways is extracted from commonly differentially expressed genes. Using an adjusted p-value less than 0.05 as a threshold, genes simultaneously enriched in all three pathways are screened as the regulatory gene set. Subsequently, a regulatory network of the regulatory gene set and upstream transcription factors is constructed using the JASPAR database on the NetworkAnalyst platform. Key transcription factors that simultaneously regulate at least two hub genes and are related to the interferon signaling pathway are identified. These key transcription factors include STAT1 and IRF family transcription factors. Based on the expression levels of these key transcription factors in the COVID-19 dataset and head and neck squamous cell carcinoma dataset, combined with the expression values ​​of hub genes, a functional polarity scoring model is constructed. This scoring model calculates the sum of the expression values ​​of individual hub genes multiplied by weighting coefficients, then divides by the sum of the expression values ​​of all hub genes. The weighting coefficients are based on a comprehensive assessment of the co-expression correlation between hub genes and key transcription factors, the frequency of literature reports, and the virus-host competition index. The weighting process emphasizes the bidirectional functional reporting of hub genes in both antiviral and pro-tumor inflammatory responses. Using the median functional polarity score as a threshold, patient samples are divided into antiviral-dominant and pro-inflammatory / pro-tumor-dominant subgroups. Finally, survival analysis or disease severity scoring is used to validate the clinical significance of the subgroup divisions, thereby determining the dominant functional direction of hub genes in specific patient populations.

[0039] Weighting coefficient W i The following rules can be used to assign values: ① Calculate the Pearson correlation coefficient r between the hub gene and key transcription factors (STAT1, IRF1). If |r|>0.6, add 0.3 to the base weight; ② Search the PubMed database and count the frequency of the gene in the literature under the themes of "antiviral" and "tumor-promoting". Use the normalized value of the frequency ratio as the correction coefficient; ③ The final weight is the sum of the above two items, normalized to the [0,1] interval.

[0040] Compared to existing technologies that only focus on the expression level or network position of interferon-stimulated genes, this invention, by integrating transcriptional regulatory network analysis, key transcription factor identification, co-expression relationships, and literature mining, has for the first time constructed a scoring model capable of quantifying the bidirectional functional polarity of interferon-stimulated genes. This model divides the same set of hub genes into antiviral dominant subgroups and pro-inflammatory / pro-cancer dominant subgroups in different patient populations, and has been clinically validated through survival analysis or disease severity scoring. This represents a leap from static gene expression profiling to dynamic functional polarity assessment, providing a quantifiable classification basis for accurately distinguishing patients' immune-inflammatory states and guiding personalized intervention strategies.

[0041] In existing technologies, methods for predicting potential therapeutic drugs based on hub genes typically recommend selected candidate drugs to all patients indiscriminately, assuming that the same group of drugs has the same applicability to all patients. This "one-size-fits-all" drug screening strategy ignores the essential differences in immune-inflammatory states among patients. Especially when the hub gene itself has a bidirectional function, the same drug may produce drastically different efficacy or even opposite adverse reactions in patients in antiviral-dominant subgroups and pro-inflammatory / pro-cancer-dominant subgroups. Current technologies lack screening models that accurately match drugs to patient subgroups.

[0042] To address the aforementioned issues, this invention, building upon the previously described method, adds a step of constructing a drug-subgroup matching screening model after identifying the top ten potential therapeutic drugs. Specifically, based on the aforementioned functional polarity scoring model, COVID-19 or head and neck squamous cell carcinoma patient samples are first divided into antiviral dominant subgroups and pro-inflammatory / pro-tumor dominant subgroups, and the average expression profile of hub genes in each subgroup is obtained as the gene characteristic of that subgroup. Subsequently, the known target gene sets of the top ten potential therapeutic drugs are extracted using the drug signature database DSigDB, and the mechanism of action information of the drugs is supplemented based on the DrugBank and ChEMBL databases, classifying these drugs into immune-activating, inflammation-suppressing, or broad-spectrum regulatory categories. A co-expression network of drug target genes and hub genes is constructed using the NetworkAnalyst platform, and the first correlation coefficient between each drug target set and the hub gene characteristics of the antiviral dominant subgroup, as well as the second correlation coefficient with the hub gene characteristics of the pro-inflammatory / pro-tumor dominant subgroup, are calculated. Based on the ratio of the first correlation coefficient to the second correlation coefficient, and combined with the enrichment significance value of the drug in the relevant pathways of the subgroup, a drug-subgroup matching index is constructed. Using the preset threshold of the matching index as the screening criterion, the top ten potential therapeutic drugs are mapped to their corresponding subgroups, generating a list of candidate drugs for antiviral dominant subgroups and a list of candidate drugs for pro-inflammatory and pro-cancer dominant subgroups, thereby achieving differentiated drug recommendations based on functional polarity.

[0043] Compared to existing technologies that indiscriminately recommend candidate drugs to all patients, this invention, by constructing a drug-subgroup matching screening model, is the first to deeply integrate patient subgroup segmentation results with drug target information, mechanisms of action, and pathway enrichment data. This model can accurately determine whether each drug is more suitable for an antiviral-dominant subgroup or a pro-inflammatory / pro-cancer-dominant subgroup based on the correlation strength between the drug target set and the hub gene characteristics of different subgroups, combined with the significant enrichment of drugs in subgroup-related pathways. This allows for the matching of the same batch of candidate drugs to corresponding patient subgroups, effectively avoiding the risks of poor efficacy or immune imbalance that may result from a "one-size-fits-all" approach to medication, and providing a clear subgroup-oriented strategy for personalized clinical medication.

[0044] In existing technologies, drug-based screening methods, whether traditional drug-target binding ability assessments or the aforementioned construction of drug-subgroup matching screening models, focus on predicting the static correlation between drugs and hub genes and the degree of matching with patient subgroups. However, they lack the ability to predict the functional polarity evolution trend of the hub gene network after drug intervention. Existing technologies can only answer whether a drug "might be effective" and "which subgroup is more suitable," but cannot answer whether the drug will cause an unexpected reversal of the immune-inflammatory balance after actual intervention. That is, although the drug can accurately match the subgroup, the complexity of network regulation during the intervention process may cause the originally antiviral dominant subgroup to shift to a pro-inflammatory and pro-cancer dominant state, thus causing treatment contradictions.

[0045] To address the aforementioned issues, this invention, based on the method described above, adds a step of constructing a dynamic perturbation model of the hub gene network under drug intervention after generating the subgroup candidate drug list. Specifically, firstly, based on the acquired COVID-19 dataset and head and neck squamous cell carcinoma dataset, the expression matrices of all hub genes are extracted. Then, using the previously constructed regulatory gene set and identified key transcription factors, a gene regulatory network containing hub genes, key transcription factors, and their interactions is constructed. The network nodes cover OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, CMPK2, and STAT1 and IRF family transcription factors. By integrating drug target information from the drug feature database and the aforementioned supplementary drug action mechanism information through the NetworkAnalyst platform, drugs from the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list are mapped to this gene regulatory network, respectively, identifying the direct target genes of each drug. Based on the structural characteristics of gene regulatory networks, a control centrality algorithm from network control theory is used to calculate the control ability score of each candidate drug on the hub gene network. This score comprehensively considers the degree centrality and betweenness centrality of drug target nodes to measure the direct impact and bridging role of drug targets on network information flow. By simulating drug intervention conditions, the expression of drug target nodes in the gene regulatory network is removed or inhibited. The changes in hub gene expression profiles after intervention are calculated using a linear differential equation model, and the polarity score after intervention is recalculated based on the aforementioned functional polarity scoring model to obtain the polarity drift. Combining the control ability score and polarity drift, a drug polarity stability index is constructed, and the candidate drug list for subgroups is reordered using this index as a screening criterion, prioritizing drugs with strong control over the hub gene network and less likely to cause polarity reversal after intervention. Alternatively, the median or a preset threshold (e.g., S>0.5) of the drug polarity stability index S value can be used as the screening criterion to sort the candidate drugs in descending order, prioritizing the top 30% of drugs. Finally, GO and KEGG enrichment analyses were used to verify the extent of pathway perturbation in the hub gene network after intervention with candidate drugs, and candidate drugs that may cause significant bypass effects were excluded.

[0046] Compared to existing technologies that only focus on the static binding ability of drugs to targets and the matching of subgroups, this invention, by constructing a regulatory network including hub genes and key transcription factors, and combining network control theory and dynamic perturbation simulation, achieves for the first time the prediction of the evolution trend of network functional polarity after intervention by candidate drugs. This model, by quantifying the drug's control over the network and the risk of polarity drift after intervention, can screen for candidate drugs that can accurately match subgroups and possess network stability. This effectively avoids treatment contradictions caused by the unexpected disruption of the immune-inflammatory balance due to drug targeting, providing a prospective assessment method for the safety and stability of clinical drug use.

[0047] In existing technologies, disease association analysis based on hub genes typically integrates the expression levels of hub genes in a whole patient sample with a gene-disease association database to screen for a spectrum of diseases associated with the hub gene as a whole, and treats these diseases as a common risk for all patients. This approach of simply applying disease association results from the overall population to all patients ignores the essential differences in immune-inflammatory states among patients. In particular, when the hub gene exhibits different functional polarities in different patient groups, such as antiviral dominance or pro-inflammatory / pro-cancer dominance, the associated disease spectrum may be drastically different, and existing technologies cannot resolve such subgroup-specific comorbidity risks.

[0048] To address the aforementioned issues, this invention, building upon the previously described method, adds a step of hierarchical analysis of disease associations based on functional polarity subgroups after constructing the association network between hub genes and diseases and analyzing hub gene-related diseases. Specifically, based on the aforementioned antiviral dominant subgroup and pro-inflammatory / pro-cancer dominant subgroup, expression matrices of hub genes from patient samples in each subgroup are extracted from the acquired COVID-19 dataset GSE196822 and head and neck squamous cell carcinoma dataset GSE178537. These hub genes include OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, and CMPK2. Using the NetworkAnalyst platform in conjunction with the DisGeNET database, association data between these hub genes and all diseases is obtained, including the score and evidence type for each gene-disease pair. The score reflects the strength of the association between the gene and the disease. For the antiviral and pro-inflammatory / pro-cancer dominant subgroups, the weighted expression levels of each hub gene in each subgroup were calculated. This weighted expression level was the product of the mean expression value of the hub gene across all samples in the subgroup and the subgroup sample size, reflecting the overall activity of the hub gene and the contribution of subgroup size to disease association. Based on the weighted expression levels and gene-disease association scores, the comprehensive association degree of each disease in a specific subgroup was calculated. This was achieved by summing the products of the weighted expression levels of all hub genes and their corresponding disease association scores, yielding the comprehensive association degree of each disease in both the antiviral and pro-inflammatory / pro-cancer dominant subgroups. The comprehensive association degrees of the two subgroups were sorted in descending order, and the top twenty diseases were selected as the characteristic-related diseases of each subgroup, forming the characteristic disease sets for the antiviral dominant subgroup and the pro-inflammatory / pro-cancer dominant subgroup. A disease-disease similarity network was constructed using a medical subject thesaurus, and the Jaccard similarity coefficients of the characteristic disease sets of the two subgroups were calculated. Diseases common to both subgroups were identified as shared diseases, and diseases unique to each subgroup were identified as specific diseases. Finally, pathway enrichment analysis was performed on shared and specific diseases. The functional enrichment differences of hub genes in shared and specific diseases were compared using the Enrichr tool to obtain the KEGG and Reactome pathways enriched in each disease set.

[0049] Compared to existing techniques that simply apply disease association results from the overall population to all patients, this invention achieves a breakthrough in disease association analysis, moving from the population-level average to the subgroup-specific level by deeply integrating functional polarity subgroup classification with gene-disease association analysis. This method can identify the unique disease sets of antiviral-dominant and pro-inflammatory / pro-cancer-dominant subgroups, and reveals the molecular basis of differences in comorbidity risk between the two subgroups through disease similarity network and pathway enrichment analysis. This provides a new analytical pathway for understanding the heterogeneity of comorbidity risk of COVID-19 and head and neck squamous cell carcinoma in patients of different polarity subgroups.

[0050] Furthermore, after obtaining the subgroup characteristic disease set in step S907, the process also includes constructing a drug-subgroup disease conflict assessment model to screen candidate drugs that have no potential conflict with the subgroup characteristic diseases, as follows: S908: Based on the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list generated in step S805, or the drug list after reordering in step S810, extract the name and structural information of the candidate drugs in each list. S909: Obtain known association data between each candidate drug and the disease through the Comparative Toxicogenomics Database. This association data includes the direction of the drug-disease association (therapeutic or induced / exacerbated), the strength of the association, and the type of evidence. Construct a drug-disease association matrix M. drug-disease The element represents the association score between the drug and the disease, with positive values ​​indicating therapeutic effects and negative values ​​indicating harmful effects; S910: Disease set D based on the antiviral dominant subgroup characteristics obtained in step S906. anti Disease set D with pro-inflammatory and pro-cancer dominant subgroups pro Calculate the conflict risk score R between each candidate drug and the corresponding subgroup of characteristic disease sets. conflict The formula for calculating the conflict risk score is: In the formula, The score represents the harmful association score between the drug and the corresponding characteristic disease. ω represents the therapeutic association score between the drug and the corresponding characteristic disease. d The comprehensive correlation degree A of disease d in the corresponding subgroup calculated in step S904 is... d,sub Normalized weights to reflect the importance of characteristic diseases to subgroups; S911: Based on the conflict risk score R conflict The candidate drugs are sorted in ascending order, and a threshold T is set. conflict The score R is between 0 and 0.3, excluding conflict risk. conflictGreater than threshold T conflict The drug; S912: Combine the drugs that have undergone conflict screening with the polarity stability index S from step S810 to construct a comprehensive drug recommendation index F. The formula for calculating the comprehensive recommendation index is: In the formula R max The drugs are ranked by F-score based on the maximum conflict risk score among all candidate drugs, generating an optimized list of candidate drugs for each subgroup. S913: Verify the reliability of the association between drugs in the optimization list and subgroup characteristic diseases through literature mining, and verify the binding ability of drugs to hub genes and key transcription factors using molecular docking simulation to ensure the practicality of the screening results.

[0051] In existing technologies, both traditional drug screening methods and the aforementioned drug-subgroup matching screening models and drug polarity stability assessments primarily focus on the potential efficacy of candidate drugs against target diseases COVID-19 and head and neck squamous cell carcinoma, as well as their suitability for patient subgroups. However, they neglect another crucial issue: when patients already have other disease risks associated with subgroup characteristics, will the candidate drug exacerbate or induce these subgroup-specific diseases? Existing technologies lack assessment mechanisms for potential conflicts between drugs and subgroup-specific diseases, potentially leading to unintended exacerbations of existing comorbidities or the induction of new diseases while treating the target disease. This is especially true when antiviral-dominant subgroups and pro-inflammatory / pro-cancer-dominant subgroups each possess different disease profiles; the adverse effects of the same drug on specific diseases may present differentiated safety risks depending on the subgroup.

[0052] To address the aforementioned issues, this invention, based on the method described above, adds a step of constructing a drug-subgroup disease conflict assessment model after obtaining the subgroup characteristic disease set. Specifically, firstly, based on the aforementioned generated lists of candidate drugs for antiviral dominant subgroups and pro-inflammatory / pro-cancer dominant subgroups, or drug lists reordered using polarity stability indices, the names and structural information of the candidate drugs in each list are extracted. Known association data between each candidate drug and the disease are obtained by comparing toxicogenomics databases. This association data includes clearly distinguishing the direction of the drug-disease association as therapeutic or aggravating, the association strength, and the type of evidence. Based on this, a drug-disease association matrix is ​​constructed, where positive values ​​represent therapeutic effects and negative values ​​represent harmful effects. Based on the aforementioned obtained sets of characteristic diseases for antiviral dominant subgroups and pro-inflammatory / pro-cancer dominant subgroups, the conflict risk score between each candidate drug and its corresponding subgroup characteristic disease set is calculated. The score is calculated by weighted summing the differences between the harmful association score and the therapeutic association score of a drug for each disease in the characteristic disease set. The weights are determined by normalizing the overall association degree of the disease within the corresponding subgroup to reflect the importance of the characteristic disease to the subgroup. Candidate drugs are sorted in ascending order based on their conflict risk scores, with a threshold of 0 to 0.3 to exclude drugs with conflict risk scores greater than the threshold. The conflict-screened drugs are then combined with the aforementioned polarity stability index to construct a comprehensive drug recommendation index. This index is used to finally rank the drugs, generating an optimized candidate drug list for each subgroup. Finally, literature mining is used to verify the reliability of the associations between the drugs in the optimized list and the characteristic diseases of the subgroups, and molecular docking simulations are used to verify the binding ability of the drugs to hub genes and key transcription factors, ensuring the practicality of the screening results.

[0053] Compared to existing technologies that only focus on the efficacy of drugs for target diseases, this invention, by integrating a drug-disease association database with a subgroup-specific disease set, is the first to construct a conflict risk scoring model capable of quantifying the harmful or therapeutic effects of drugs on important subgroup diseases. This model integrates the conflict risk score with a polarity stability index, effectively mitigating potential risks of subgroup-specific comorbidities while ensuring the drug's efficacy for the target disease. This achieves precise drug recommendations that balance efficacy and safety, providing a complete safety assessment framework for personalized medication in complex comorbidity scenarios.

[0054] Application examples This application example, based on the proposed bioinformatics analysis method, systematically studies the common molecular mechanisms and potential therapeutic drugs of COVID-19 and head and neck squamous cell carcinoma (HNSCC), and validates the key steps involved.

[0055] 1. Data Acquisition and Differentially Expressed Gene Identification: The COVID-19 dataset GSE196822 (containing peripheral blood samples from 9 healthy controls and 34 COVID-19 patients) and the HNSCC dataset GSE178537 (containing paired primary tumor tissue and adjacent normal tissue samples from 21 patients) were obtained from the GEO database.

[0056] Using the DESeq2 package in R, with FDR < 0.01 and |log2FC| > 1 as thresholds, 808 differentially expressed genes (DEGs) were identified in the COVID-19 dataset (554 upregulated and 254 downregulated). Figure 1 As shown in Figure A, 4458 DEGs were identified in the HNSCC dataset (2085 upregulated and 2373 downregulated), as follows: Figure 1 As shown in B. Intersection analysis of the two DEGs using the Jvenn tool yielded 139 common differentially expressed genes (common DEGs), as shown in Figure B. Figure 1 As shown in C. These common genes form the basis for subsequent analyses.

[0057] 2. Functional Enrichment and Pathway Analysis: Enrichment analysis of GO and KEGG pathways was performed on 139 common DEGs using Enrichr. For example... Figure 2 As shown, these genes are significantly enriched in pathways related to "immune response" and "viral genome replication," such as the "type I interferon signaling pathway," "neutrophil-mediated immunity," and "defense response against viruses." KEGG pathway analysis further indicates that, as Figure 3 As shown, these genes are associated with "cytokine-cytokine receptor interactions" and "transcriptional dysregulation in cancer," validating the crucial role of immune and virus-related pathways in both diseases.

[0058] 3. PPI Network Construction and Cooperative Control Node Identification: 139 common DEGs were imported into the STRING database to construct a PPI network (overall score > 0.5). This network contains 69 nodes and 169 edges, as shown below. Figure 4 As shown Based on this, identifying viral entry-related genes serves as a collaborative control node: Three viral genes, ACE2, TMPRSS2 and CTSL, were extracted and incorporated into the relevant genes to construct a PPI subnetwork with common DEGs.

[0059] Based on network control theory, the control ability score K of viral entry into related genes on the hub gene network was calculated. The calculation revealed that although the number of hub genes directly connected to the ACE2 gene in the network (Nc) is not large, the shortest path (Sp) with hub genes (such as IFITM3 and LY6E) is extremely short, and its control ability score K (when γ=0.6, δ=0.4) is significantly higher than the preset threshold.

[0060] ACE2 was identified as a key co-control node. This node affects the entry efficiency of SARS-CoV-2 into HNSCC tissues by regulating the expression of interferon-stimulated genes such as IFITM3. ACE2 was integrated with subsequently screened hub genes to form an expanded hub gene set.

[0061] 4. Construction of virus-host competition index model: Genes are divided into two groups according to function: virus entry promotion group (ACE2, TMPRSS2, CTSL) and host defense promotion group (OASL, IFIT3, RSAD2, MX1, IFITM3, LY6E, CMPK2).

[0062] Principal component analysis (PCA) was used to construct the virus entry capability index (V) and the host defense capability index (H). The virus-host competition index (C) was calculated based on the ratio of the two indices. VH =log2(V / H)). Analysis results show that in the GSE196822 dataset, most COVID-19 patients have C... VH A value <0 indicates dominant host defense; however, in tumor samples from GSE178537, some patients had C... VH A value >0 indicates that the virus has a dominant ability to enter the system.

[0063] HNSCC patient samples were analyzed according to C VH Stratified survival analysis showed that C VH The patient population with a >0 (i.e., those with dominant viral entry capacity) had a significantly shorter disease-free survival (P<0.05), validating the clinical relevance of this index.

[0064] 5. Hub gene screening and functional polarity scoring model construction: Using Cytohubba's MCC algorithm, the top 10 hub genes with the highest connectivity were identified from the PPI network: OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, and CMPK2. Figure 5 These core hub genes and their interaction networks are shown. These genes are all interferon-stimulated genes (ISGs).

[0065] From common DEGs, a set of regulatory genes simultaneously enriched in the "interferon signaling," "antigen processing and presentation," and "inflammatory response" pathways was screened, and key upstream transcription factors STAT1 and IRF1 were identified. Figure 6 It demonstrates the regulatory network of 46 transcription factors, including FOXC1 and STAT1, and hub genes.

[0066] Construct a functional polarity rating model: P = ∑(W i ×E i ) / ∑E j The weighting coefficient Wi is assigned based on the correlation between the hub gene and STAT1 / IRF1 co-expression and the frequency of reports on its antiviral and pro-tumor functions in the literature. Virus-host competition index C. VH As an auxiliary adjustment factor, the weighting coefficients are corrected.

[0067] HNSCC patients were divided into two subgroups based on the median P value: the antiviral dominant subgroup (high P value, dominated by host defense genes) and the pro-inflammatory and pro-cancer dominant subgroup (low P value, dominated by inflammation-related genes).

[0068] Clinical pathological feature analysis showed that patients in the pro-inflammatory and pro-cancer dominant subgroup had higher tumor stages and higher lymph node metastasis rates, further validating the clinical significance of subgroup classification.

[0069] 6. Candidate drug prediction and drug-subgroup matching screening: Based on the DSigDB database, drug prediction was performed using 10 hub genes. With p < 5.20E-07 as the threshold, the top 10 potential therapeutic drugs were obtained: sulocitidil, acetoxemamide, prenylamine, clioquinol, prochlorperazine, terfenadine, thioridazine, vanoxerine, chlorophyllin, 3'-Azido-3'-deoxythymidine.

[0070] Obtain the average expression profiles of hub genes for the dominant antiviral and pro-inflammatory / pro-tumor subgroups. Calculate the correlation coefficients R1 (with the antiviral subgroup) and R2 (with the pro-inflammatory / pro-tumor subgroup) between the target set of each drug and the hub gene characteristics of the two subgroups using network analysis.

[0071] The drug-subgroup matching index M = (R1 / R2) × log2(1 + enrichment significance) was constructed. The results showed that clioquinol and prochlorperazine had M values ​​> 1 and were matched to the dominant antiviral subgroup; while thioridazine had an M value < 1 and was matched to the dominant pro-inflammatory and pro-cancer subgroup.

[0072] 7. Drug polarity and stability assessment: A gene regulatory network containing 10 hub genes and key transcription factors (STAT1, IRF1) was constructed. Candidate drugs clioquinol and thioridazine were mapped to the network, with their direct targets being CMPK2 and OASL / IFITM3, respectively.

[0073] Calculate the drug's control ability score C on the network, and simulate drug intervention to calculate the polarity drift ΔP.

[0074] The polarity stability index S = C / (1+ΔP) was calculated. Clioquinol showed a higher S value, indicating its ability to stabilize network polarity in the antiviral dominant subpopulation and its low likelihood of inducing reversal. Thioridazine also showed a high S value, indicating its ability to stably inhibit the expression of inflammation-related genes in the pro-inflammatory and pro-tumorigenic subpopulation. Calculations showed that clioquinol had a control ability score of C = 0.72, polarity drift ΔP = 0.08, and a polarity stability index S = 0.72 / (1+0.08) = 0.667; thioridazine had C = 0.68, ΔP = 0.12, and S = 0.607. Both are above the median threshold of 0.55, therefore, they are preferred.

[0075] GO enrichment analysis showed that after clioquinol intervention, antiviral-related pathways (such as "response to the virus") remained stably enriched, with no significant bypass activation.

[0076] 8. Disease association stratification and drug conflict assessment: Combining the DisGeNET database, disease association stratification was performed on antiviral dominant subgroups and pro-inflammatory and pro-cancer dominant subgroups. Figure 7 This study demonstrated the association network between hub genes and diseases, and further identified the characteristic disease set D of the antiviral dominant subgroup. anti It mainly includes infectious diseases such as "viral pneumonia" and "influenza"; while the characteristic disease set of the pro-inflammatory and pro-cancer dominant subgroup is D. pro This mainly includes autoimmune and inflammatory diseases such as "rheumatoid arthritis" and "Crohn's disease".

[0077] Pathway enrichment analysis showed that D anti The associated diseases are enriched in the "cytokine storm" and "viral entry" pathways, while D pro The associated diseases are enriched in the IL-17 signaling pathway and the NF-kappa B signaling pathway.

[0078] A drug-subgroup disease conflict assessment model was constructed. For thioridazine, a candidate drug in the pro-inflammatory and pro-cancer subgroup, the CTD database was used to identify drug-disease conflicts.pro The term "rheumatoid arthritis" in the text has a known treatment association (positive value) and is associated with D. pro The harmful association score for diseases in China is low. Therefore, its conflict risk score R is low. conflict It is only 0.1, which is below the threshold T. conflict (0.3). Based on its polarity and stability index S, the final comprehensive drug recommendation index F was obtained, confirming that thioridazine is an optimized candidate drug for the dominant pro-inflammatory and pro-cancer subgroups.

[0079] Literature review supports the role of thioridazine in enhancing the chemosensitivity of HNSCC, and molecular docking simulations also verified its good binding ability to OASL protein.

[0080] Experimental verification results: To preliminarily verify the above analysis results, this study selected thioridazine, a candidate drug matched to the pro-inflammatory and pro-cancer dominant subgroup, and clioquinol, a candidate drug matched to the antiviral dominant subgroup, for in vitro experimental verification.

[0081] In validating the pro-inflammatory and pro-cancer dominant subsets, head and neck squamous cell carcinoma lines CAL27 and SCC-9 were treated with different concentrations of thioridazine. MTT assays showed that thioridazine significantly inhibited the proliferation of both cell lines in a dose-dependent manner. Western blotting experiments indicated that thioridazine treatment downregulated the expression of hub genes such as OASL and IFITM3, while simultaneously inhibiting the activity of the downstream PI3K / AKT signaling pathway, consistent with bioinformatics predictions of network regulation. Furthermore, the combination of thioridazine and the chemotherapeutic drug carboplatin showed a synergistic inhibitory effect.

[0082] In the validation of the dominant antiviral subpopulation, A549 cells were used, and the regulatory role of clioquinol was observed by inducing an interferon response via poly(I:C). The results showed that clioquinol synergistically enhanced the expression of interferon-stimulated genes, consistent with the predicted antiviral function of the subpopulation.

[0083] The above in vitro experimental results preliminarily validate the reliability of this analytical method in identifying candidate drugs that match the characteristics of the subgroup.

[0084] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details and illustrations shown and described herein.

Claims

1. A method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics, characterized in that, Includes the following steps: S1: Obtain the transcriptome datasets GSE196822 for COVID-19 and GSE178537 for head and neck squamous cell carcinoma from the Gene Expression Comprehensive Database; S2: Using the DESeq2 package in R software, with a false discovery rate of less than 0.01 and an absolute value of log2 transcriptional fold change greater than 1 as thresholds, the first differentially expressed gene set and the second differentially expressed gene set were screened from the COVID-19 dataset and the head and neck squamous cell carcinoma dataset, respectively. S3: Intersection analysis of the first and second differentially expressed gene sets was performed using the online Venn graph tool Jvenn to obtain the common differentially expressed genes; S4: Gene ontology analysis and pathway enrichment analysis of commonly differentially expressed genes were performed using the enrichment analysis tool Enrichr to identify biological processes and signaling pathways related to immune response and viral genome replication; S5: Construct a protein-protein interaction network of proteins encoded by common differentially expressed genes using the STRING database, and import the network into Cytoscape software for visualization, setting a comprehensive score greater than 0.5 as the screening criterion; S6: Using the Cytohubba plugin in Cytoscape software, the nodes in the protein-protein interaction network were sorted based on the maximum clique centrality algorithm, and the top 10 commonly differentially expressed genes with the highest connectivity were selected as hub genes. The hub genes are OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E and CMPK2. S7: Using the NetworkAnalyst platform, regulatory networks of hub genes and transcription factors, as well as regulatory networks of hub genes and microRNAs, were constructed based on the JASPAR database and the TarBase database, respectively. S8: Based on the drug characterization database DSigDB, protein-drug interactions were predicted using the hub gene. The top 10 potential therapeutic drugs were obtained with a p-value of less than 5.20E-07 as the screening criterion. S9: Construct a network linking hub genes and diseases using the NetworkAnalyst platform in conjunction with the DisGeNET database, and analyze diseases related to hub genes; S10: Based on the pathway enrichment results of step S4 and the hub gene of step S6, combined with the gene-disease association network of step S9, we identify the interferon-stimulated gene subgroup in the hub gene that has dual functions of antiviral and tumor-promoting, and construct its functional polarity scoring model to analyze the dominant functional direction of the hub gene in different patient groups.

2. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 1, characterized in that, After constructing the protein-protein interaction network in step S5 and before screening hub genes in step S6, the process also includes identifying viral entry-related genes as co-control nodes. S501: Based on the COVID-19 dataset GSE196822 and head and neck squamous cell carcinoma dataset GSE178537 obtained in step S1, extract the expression values ​​of known genes related to SARS-CoV-2 virus entry into host cells. SARS-CoV-2 virus entry-related genes include ACE2, TMPRSS2 and CTSL. S502: Construct a protein-protein interaction subnetwork between SARS-CoV-2 virus entry-related genes and the common differentially expressed genes obtained in step S3 using the STRING database. Set the comprehensive score to be greater than 0.7 as the screening condition to identify common differentially expressed genes that have direct physical interactions or functional associations with virus entry-related genes. S503: Based on the minimum driving node set algorithm in network control theory, using the protein-protein interaction network constructed in step S5 as input, calculate the control capability score K of the virus entering the relevant gene on the network module where the hub gene is located. The calculation method of the control capability score K is as follows: K = γ × Nc + δ × (1 / Sp); Where Nc is the number of hub genes directly regulated by viral entry-related genes, Sp is the shortest path length between viral entry-related genes and hub genes, and γ and δ are preset weighting coefficients. S504: Virus entry-related genes with a control ability score K higher than a preset threshold are identified as co-control nodes. Co-control nodes affect the entry efficiency and replication ability of SARS-CoV-2 virus in head and neck squamous cell carcinoma tissue by regulating the expression level of interferon-stimulated genes. S505: Integrate the identified co-control nodes with the hub genes screened in step S6 to construct an expanded hub gene set for subsequent drug prediction and functional polarity analysis.

3. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 2, characterized in that, Step S505, after integrating the identified cooperative control nodes with the hub gene, also includes constructing a virus-host competition index model to quantify the dynamic game relationship between virus entry and host defense. Specific steps include: S506: Based on the virus entry-related genes extracted in step S501 and the hub genes screened in step S6, the genes are divided into the virus entry promotion group and the host defense promotion group according to their biological functions. The virus entry promotion group includes ACE2, TMPRSS2 and CTSL, and the host defense promotion group includes OASL, IFIT3, RSAD2, MX1, IFITM3, LY6E and CMPK2. S507: Based on the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S2, the expression values ​​of each gene in the virus entry promotion group and the host defense promotion group are extracted. Principal component analysis is used to reduce dimensionality and construct the virus entry capability index V and the host defense capability index H, respectively. The calculation formula is V = ∑λ i ×G i H = ∑μ j ×G j G i The expression values ​​of each gene in the viral entry promotion group, G j λ represents the expression values ​​of each gene in the host defense enhancement group. i and μ j The loading coefficient of the first principal component in principal component analysis; S508: Constructing a virus-host competition index C based on the ratio of the virus entry capability index V to the host defense capability index H. VH Its calculation formula is C VH =log2(V / H), when C VH A value greater than 0 indicates that the virus has a dominant ability to enter the body; when C... VH A value less than 0 indicates that the host's defense capabilities are dominant; S509: Using the clinical information from the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S1, patient samples are sorted according to the virus-host competition index C. VH Perform stratification and verify C VH Correlation with patient viral load, disease severity, or tumor stage; S510: The virus-host competition index C VH As a functional output indicator of the collaborative control node, it is used to evaluate the regulatory efficacy of the collaborative control node identified in step S503 in different patient populations, and to include the virus-host competition index C. vh As an auxiliary factor for adjusting the weight coefficients, it is integrated into the functional polarity scoring model.

4. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 3, characterized in that, The hub genes screened in step S6 are all interferon-stimulated genes. After identifying the hub genes, the next step is to construct a hub gene functional polarity scoring model. S601: Based on the interferon signaling pathway, antigen processing and presentation pathway and inflammatory response pathway obtained from the pathway enrichment analysis in step S4, a set of regulatory genes related to all three pathways is extracted from the common differentially expressed genes. Among them, the common differentially expressed genes that are simultaneously enriched in the interferon signaling pathway, antigen processing and presentation pathway and inflammatory response pathway are selected as the set of regulatory genes, with the adjusted p value less than 0.05 as the threshold. S602: Using the NetworkAnalyst platform and the JASPAR database, a regulatory network of regulatory gene sets and upstream transcription factors was constructed to identify key transcription factors that simultaneously regulate at least two hub genes and are related to the interferon signaling pathway. Key transcription factors include STAT1 and IRF family transcription factors. S603: Based on the expression levels of key transcription factors in the COVID-19 dataset and head and neck squamous cell carcinoma dataset, combined with the expression value of the hub gene, a functional polarity scoring model was constructed. The calculation formula for the scoring model is as follows: P = ∑(Wi×Ei) / ∑Ej; Where Ei is the expression value of a single hub gene, Ej is the sum of the expression values ​​of all hub genes, and Wi is the weighting coefficient, which is based on the co-expression correlation between hub genes and key transcription factors, the frequency of literature reports, and the virus-host competition index C. VH A comprehensive value was assigned; the weighting coefficients were assigned based on the correlation between the co-expression of the hub gene and key transcription factors, as well as the frequency of literature reports on the bidirectional function of the hub gene in antiviral response and pro-tumor inflammatory response. S604: Using the median of the functional polarity score P as the threshold, the patient sample was divided into an antiviral dominant subgroup and a pro-inflammatory and pro-cancer dominant subgroup. S605: Validate the clinical significance of subgroup division through survival analysis or disease severity scoring to determine the dominant functional direction of the hub gene in a specific patient population.

5. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 4, characterized in that, After screening the top 10 potential therapeutic drugs in step S8, the next step is to construct a drug-subgroup matching screening model: S801: Based on the functional polarity scoring model constructed in step S603, COVID-19 or head and neck squamous cell carcinoma patient samples were divided into antiviral dominant subgroups and pro-inflammatory and pro-cancer dominant subgroups, and the average expression profile of hub genes in each subgroup was obtained. S802: Using the DSigDB database in step S8, extract the known target gene set of the top 10 potential therapeutic drugs, and supplement the drug mechanism information based on the DrugBank and ChEMBL databases, classifying the drugs into immune activation, inflammation suppression or broad-spectrum regulation types. S803: Construct a co-expression network of drug target genes and hub genes using the NetworkAnalyst platform, and calculate the first correlation coefficient R1 between each drug target set and the characteristics of hub genes in the antiviral dominant subgroup, and the second correlation coefficient R2 between each drug target set and the characteristics of hub genes in the pro-inflammatory and pro-cancer dominant subgroup. S804: Construct the drug-subgroup matching index M based on the ratio of the first correlation coefficient R1 and the second correlation coefficient R2. The calculation formula is M = (R1 / R2) × log2(1 + the enrichment significance value of the drug in the subgroup-related pathway). S805: Using the threshold of the matching index M as the screening criterion, the top 10 potential therapeutic drugs are mapped to their corresponding subgroups, generating a list of candidate drugs for antiviral dominant subgroups and a list of candidate drugs for pro-inflammatory and pro-cancer dominant subgroups, thus realizing differentiated drug recommendations based on functional polarity.

6. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 5, characterized in that, After generating the subgroup candidate drug list in step S805, the process also includes constructing a dynamic perturbation model of the hub gene network under drug intervention to assess the polarity stability of the candidate drugs. This specifically includes the following steps: S806: Based on the COVID-19 dataset and head and neck squamous cell carcinoma dataset obtained in step S2, extract the expression matrix of all hub genes, and use the regulatory gene set constructed in step S601 and the key transcription factors identified in step S602 to construct a gene regulatory network containing hub genes, key transcription factors and their interactions. The network nodes include OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E, CMPK2 and STAT1 and IRF family transcription factors; S807: Integrate the drug target information from the DSigDB database in step S8 and the drug mechanism information supplemented in step S802 through the NetworkAnalyst platform, and map the drugs in the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list in step S805 to the gene regulatory network respectively, and identify the direct target genes of each drug. S808: Based on the structural characteristics of gene regulatory networks, the control centrality algorithm in network control theory is used to calculate the control ability score C of each candidate drug on the hub gene network. The calculation method of the control ability score C is as follows: C = α × D + β × B; Where D is the degree centrality of the drug target node, B is the betweenness centrality of the drug target node, and α and β are preset weight coefficients used to measure the direct impact of the drug target on the network information flow and its bridging role. S809: By simulating drug intervention conditions, the expression of drug target nodes in the gene regulatory network is removed or inhibited. The change in the expression profile of the hub gene after intervention is calculated using a linear differential equation model. Based on the functional polarity scoring model in step S603, the polarity score P´ after intervention is recalculated to obtain the polarity drift ΔP=|P´-P|. S810: Combine the comprehensive control ability score C and polarity drift ΔP to construct the drug polarity stability index S, which is calculated as S=C / (1+ΔP). Using the stability index S as the screening criterion, the list of candidate drugs for the subgroup in step S805 is reordered, and drugs with high S values ​​are given priority, that is, drugs with strong control over the hub gene network and less likely to cause polarity reversal after intervention. S811: GO and KEGG enrichment analyses were used to verify the extent of pathway perturbation in the hub gene network after intervention with candidate drugs, and candidate drugs that may cause significant bypass effects were excluded.

7. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 6, characterized in that, After constructing the association network between hub genes and diseases and analyzing hub gene-related diseases in step S9, the process also includes a step of hierarchical analysis of disease associations based on functional polarity subgroups: S901: Based on the antiviral dominant subgroup and pro-inflammatory and pro-cancer dominant subgroups identified in step S604, the expression matrix of hub genes in each subgroup of patient samples is extracted from the COVID-19 dataset GSE196822 and the head and neck squamous cell carcinoma dataset GSE178537 obtained in step S1. The hub genes include OASL, IFIT3, RSAD2, OAS3, MX1, OAS1, HERC5, IFITM3, LY6E and CMPK2. S902: Using the NetworkAnalyst platform in conjunction with the DisGeNET database, obtain association data between hub genes and all diseases. This association data includes a score for each gene-disease pair. g,d And the type of evidence, Score g,d Reflects the strength of the association between gene g and disease d; S903: Calculate the weighted expression level (E) of each hub gene in the antiviral dominant subset and the pro-inflammatory / pro-cancer dominant subset, respectively. g,sub The weighted expression level is the product of the mean expression value of the hub gene in all samples within the subpopulation and the sample size of the subpopulation, i.e., E. g,sub =mean(expression) g,sub )×N sub To reflect the overall activity of hub genes in the subpopulation and the contribution of subpopulation size to disease association; S904: Based on weighted expression level E g,sub Gene-Disease Association Score g,d Calculate the comprehensive association degree A of each disease d in a specific subgroup. d,sub The calculation formula is A d,sub =∑(E d,sub ×Score g,d ), summate and iterate through all hub genes to obtain the comprehensive association degree of each disease in the antiviral dominant subgroup and the pro-inflammatory and pro-cancer dominant subgroup; S905: Overall correlation between the dominant antiviral subset and the dominant pro-inflammatory / pro-cancer subset, respectively. d,sub The diseases were sorted in descending order, and the top 20 were selected as characteristic diseases of each subgroup, forming the characteristic disease set D of the antiviral dominant subgroup. anti Disease set D with pro-inflammatory and pro-cancer dominant subgroups pro ; S906: Construct a disease-disease similarity network using a medical subject thesaurus and calculate D. anti and D pro The Jaccard similarity coefficient is used to identify diseases shared by two subgroups as shared diseases, and diseases unique to each subgroup as specific diseases. S907: Pathway enrichment analysis was performed on shared and specific diseases. The Enrichr tool was used to compare the differences in functional enrichment of hub genes in shared and specific diseases, and the KEGG and Reactome pathways enriched in each disease set were obtained, thereby revealing the molecular basis of the differential risk of comorbidity between COVID-19 and head and neck squamous cell carcinoma in patients of different polarity subgroups.

8. The method for identifying common targets of SARS-CoV-2 and head and neck squamous cell carcinoma and screening drugs based on bioinformatics according to claim 7, characterized in that, After obtaining the subgroup-characteristic disease set in step S907, the process also includes constructing a drug-subgroup disease conflict assessment model to screen candidate drugs that have no potential conflict with the subgroup-characteristic diseases, as follows: S908: Based on the antiviral dominant subgroup candidate drug list and the pro-inflammatory and pro-cancer dominant subgroup candidate drug list generated in step S805, or the drug list after reordering in step S810, extract the name and structural information of the candidate drugs in each list. S909: Obtain known association data between each candidate drug and the disease by comparing toxicogenomics databases. This association data includes the direction, strength, and type of evidence of the drug-disease association, and a drug-disease association matrix M is constructed. drug-disease The element represents the association score between the drug and the disease, with positive values ​​indicating therapeutic effects and negative values ​​indicating harmful effects; S910: Disease set D based on antiviral dominant subgroup characteristics obtained in step S906 anti Disease set D with pro-inflammatory and pro-cancer dominant subgroups pro Calculate the conflict risk score R between each candidate drug and the corresponding subgroup of characteristic disease sets. conflict The formula for calculating the conflict risk score is: In the formula, The score represents the harmful association score between the drug and the corresponding characteristic disease. ω represents the therapeutic association score between the drug and the corresponding characteristic disease. d The comprehensive correlation degree A of disease d in the corresponding subgroup calculated in step S904 is... d,sub Normalized weights to reflect the importance of characteristic diseases to subgroups; S911: Based on the conflict risk score R conflict The candidate drugs are sorted in ascending order, and a threshold T is set. conflict The score R is between 0 and 0.3, excluding conflict risk. conflict Greater than threshold T conflict The drug; S912: Combine the drugs that have undergone conflict screening with the polarity stability index S from step S810 to construct a comprehensive drug recommendation index F. The formula for calculating the comprehensive recommendation index is: ; In the formula R max The drugs are ranked by F-score based on the maximum conflict risk score among all candidate drugs, generating an optimized list of candidate drugs for each subgroup. S913: Verify the reliability of the association between drugs in the optimization list and subgroup characteristic diseases through literature mining, and verify the binding ability of drugs to hub genes and key transcription factors using molecular docking simulation to ensure the practicality of the screening results.

Images