Identifying Potential Gene Defect Patterns Related to COVID-19 Based on Pharmacological and Bioinformatics Analysis for Lung Adenocarcinoma

Introduction

Studies reported that cancer patients are more susceptible to coronavirus disease 2019 (COVID-19) and may develop severe illness due to suppressed immunity, a higher frequency to health-care centers, and older ages.¹ In addition, cancer patients infected with COVID-19 are more likely to deteriorate.² A normal health care and management of cancer patients could be disrupted in the pandemic as a result of the limitations and re-arrangement of medical resources. Clinical investigations also supported the conclusion that cancer patients with COVID-19 are associated with worse outcomes than with cancer only.^3,4 The mortality of cancer patients with COVID-19 is approximately 25% to 30%.^5,6 Moreover, patients who have recently received therapeutic management show a higher risk to develop worse outcomes, such as complications or even death.⁷

Lung cancer is a malignant tumor posing great threat to human health. In 2020, the number of lung cancer cases was about 2.21 million in the world, ranking the second among all new cancer cases of that year and the highest in all cancer-related mortality. It is estimated that in 2040, the number of lung cancer cases worldwide will reach about 3.6 million.^8,9 Study showed that lung cancer is the most common cancer type among patients infected with COVID-19.¹⁰ Lung cancer patients face a higher death risk because of the possibilities such as smoking-induced pulmonary injury and insufficiency of respiratory healthcare during the pandemic.^10,11 To alleviate the adverse effects resulted from COVID-19, professional organizations have proposed a series of recommendations to guide decision-making for managing patients with COVID-19. In China, the traditional Chinese medicine (TCM) plays a positive role in inhibiting COVID-19 progression in patients either in mild or severe condition.^12–14 TCM protocol in diagnosis and treatment protocol for COVID-19 (trial version 8) developed by the National Health Commission in China proposes nine different TCM prescriptions for different stages of COVID-19 infection.¹⁵

As TCM is effective in treating COVID-19 and the unfavorable situations of cancer patients with COVID-19, we, therefore, identified TCM drug targets related to both COVID-19 and lung adenocarcinoma (LUAD).^16–18 By using data from public databases, pharmacologic analysis and bioinformatics analysis, we want to identify the main COVID-19 target genes both associated with TCM and LUAD. These target genes may be new molecular targets for treating lung cancer patients with COVID-19. Furthermore, based on target genes, we also would like to investigate novel molecular subtypes capable of guiding targeted therapies such as immunotherapy for LUAD patients.

Materials and Methods

Data Source and Data Preprocessing

The Cancer Genome Atlas (TCGA)-LUAD dataset including RNA-seq data and clinical information of LUAD samples was downloaded from TCGA database on June 30, 2021. For RNA-seq data, samples without survival information or follow-up information were excluded. Ensembl ID was transferred into gene symbols. Median expression value was taken when one gene had multiple gene symbols. After preprocessing, 472 LUAD samples in the TCGA-LUAD dataset remained. Data of gene mutation and genomic variation including single nucleotide variation (SNV) and copy number variation (CNV) were independently downloaded from TCGA database. The GSE31210 dataset was downloaded from Gene Expression Omnibus (GEO) database on June 30, 2021. In the GSE31210 dataset, samples without survival information or follow-up information were excluded. Probes were transferred into gene symbols. One probe corresponded to multiple genes was excluded. Median value was taken when multiple probes corresponded to one gene. After preprocessing, 226 LUAD samples remained.

Identification of Differentially Expressed Genes Associated with COVID-19

COVID-19-related genes with a total of 1566 genes were downloaded from GeneCards database (https://www.genecards.org) in June 30, 2021. According to “gencode.v32.annotation.gff3” file, protein-coding genes (PCGs) were annotated. Limma R package was applied to identify differentially expressed genes (DEGs) between tumor samples and normal samples in the TCGA-LUAD dataset.¹⁹ |log2(fold change)| >1 and false discovery rate (FDR) <0.05 were determined as the threshold to screen DEGs. Venn analysis was used to show DEGs associated with COVID-19.

Identification of COVID-19 Target Genes

According to the TCM protocol in diagnosis and treatment protocol for COVID-19 (trial version 8), nine TCM prescriptions were obtained. For each prescription, main active ingredients were identified using the Traditional Chinese Medicine Database and Analysis Platform (TCMSP, https://tcmsp-e.com/). Effective molecules of active ingredients and drug targets were screened under the criteria of oral bioavailability (OB) >30% and drug-likeness (DL) >0.18. UniProt database (https://www.uniprot.org/) was used to annotate drug targets. CytoScape (V3.8.0)²⁰ was introduced to develop a network among TCM, molecules and drug targets related to COVID-19. NetworkAnalyzer tool in CytoScape was employed to analyze topology for each interaction. The top 10 molecules and drug targets of FPQXZ were listed based on the degree and BetweennessCentrality. Finally, drug targets interacting with molecules and COVID-19 were confirmed as COVID-19 target genes.

Construction of Protein–Protein Interaction (PPI) Networks for COVID-19 Target Genes

Interactions among the 26 COVID-19 target genes were predicted by String online tool (https://cn.string-db.org/) and visualizing by CytoScape. NetworkAnalyzer tool in CytoScape was used to calculate the degree and determine the importance of one gene in the network. Under conditions of degree cutoff = 2, node score cutoff = 0.2, k-core = 2 and max depth = 100, MCODE and Haircut algorithms were applied to identify gene modules.

Functional Analysis

WebGestalt, which is a popular tool for interpreting functional categories based on the enrichment of gene sets, was employed to annotate gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.²¹ The top 10 enriched terms and pathways were visualized.

Assessment of Immune Features

ESTIMATE was employed to calculate immune score of tumor microenvironment (TME) in the TCGA-LUAD dataset based on immune-related gene sets.²² CIBERSORT was used to calculate the enrichment score for the 22 immune cells.²³ Pearson correlation rank analysis was performed to assess the relation between the immune score and COVID-19 target genes, between the enrichment of immune cells and COVID-19 target genes.

Mutation Analysis of 26 COVID-19 Target Genes

Mutect2 algorithm was applied to detect SNVs, including missense, nonsense, splice site mutations, frame-shift insertions and deletions.²⁴ CNVs were identified by Segment_Mean, where Segment_Mean >0.2 indicates amplification, −0.2 < Segment_Mean <0.2 indicates diploid, and Segment_Mean < −0.2 indicates deletion.

Correlation Analysis Between COVID-19 Target Genes and Cell Death

Gene sets of cell death including pyroptosis, ferroptosis and necroptosis were downloaded from Molecular Signatures Database v7.4 (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/). Single sample gene set enrichment analysis (ssGSEA) was conducted to calculate the enrichment score of the three types of cell death for each sample.²⁵ Pearson correlation analysis was used to evaluate the association between COVID-19 target genes and cell death.

Unsupervised Consensus Clustering

Unsupervised consensus clustering in ConsensusClusterPlus R package is a popular method for an efficient and automatic clustering of samples based on gene expression patterns into different groups with different molecular features.²⁶ Expression profiles were centralized and normalized using scale. PAM algorithm and Spearman was the basis to conduct 500 times of bootstraps, with 90% samples in the TCGA-LUAD dataset included in each bootstrap. Cluster number k was set from 2 to 10. Cumulative distribution function (CDF) and consensus matrix were used to confirm the optimal cluster number.

Prediction of Sensitivity to Immunotherapy and Chemotherapy

To predict the response of different molecular subtypes to immunotherapy, a treatment dataset IMvigor210 was introduced as a reference and SubMap analysis was used for comparing the similarity of expression profiles.²⁷ P values between IMvigor210 and the TCGA-LUAD or the GSE31210 represent the similarity, where lower P values indicate a higher similarity. pRRophetic R package was applied to assess the biochemical half-maximal inhibitory concentration (IC50) of chemotherapeutic drugs including cisplatin, erlotinib, rapamycin, sunitinib and MG-132.²⁸

Statistical Analysis

Student's t-test was conducted to detect statistical significance between the two groups. Log-rank test was used in Kaplan–Meier survival analysis. P < 0.05 was defined as significant. ns, no significance. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Statistical analysis was performed in R software (4.1.1). Parameters of tools without introduction were set as default.

Results

Identifying Dysregulated Genes Related to Both COVID-19 and LUAD

Firstly, we compared gene expression profiles between normal and tumor samples in the TCGA-LUAD dataset. Two thousand six hundred and seventeen DEGs (1012 upregulated and 1605 downregulated) were screened under the threshold of |log2 (fold change)| >1 and FDR < 0.05 (Figure 1A). Two hundred and ninety-one genes (167 downregulated and 124 upregulated) were found in the intersection of the 2617 DGEs with COVID-19-related genes (Figure 1B), indicating that these 291 genes may be potential drug targets for treating LUAD patients with COVID-19. To further understand the function of these genes, we analyzed enriched KEGG pathways with respect to upregulated and downregulated genes, respectively. Upregulated genes were significantly enriched to 5 KEGG pathways (FDR < 0.05), including drug metabolism, biosynthesis of amino acids, Fanconi anemia pathway, Interleukin (IL)-17 signaling pathway, and pyrimidine metabolism (Figure 1C). Downregulated genes were significantly enriched to 50 KEGG pathways, and the top 5 enriched pathways, namely, Malaria, Staphylococcus aureus infection, Rheumatoid arthritis, AGE-RAGE signaling pathway in diabetic complications, and hematopoietic cell lineage, were visualized (Figure 1D).

Figure 1 Identification of COVID-19-related target genes. (A) Volcano plot of upregulated (blue) and downregulated (red) genes in the TCGA-LUAD dataset. Vertical dashed lines indicate |log2(fold change)| = 1. Horizontal dashed line indicates FDR < 0.05. (B) Venn plot among downregulated genes (purple), upregulated genes (red) and COVID-19-related genes (green). (C and D) Enriched KEGG pathways on upregulated (C) and downregulated (D) COVID-19-related genes. Different colors indicate different pathways, and each radiation lines indicate one gene in one pathway.

Constructing Networks of “TCM-Chemical Compounds-Target Genes” in Combination COVID-19 with LUAD

According to the TCM protocol in diagnosis and treatment protocol for COVID-19 (trial version 8), we obtained nine TCM prescriptions. For each TCM, the main active ingredients were analyzed using TCMSP. The effective molecules and drug targets were screened with criteria of OB ≥ 30% and DL ≥ 0.18. Drug targets were annotated by UniProt database. Next, CytoScape was applied to develop networks among drug targets, TCMs and ingredients grouped by each TCM prescription. For YDBFZ, 12 TCMs were included, and a number of effective ingredients of 12 TCMs were presented in the left of the network (Figure 2). A large number of targets of these molecules were visualized, and 18 targets were related to COVID-19, as shown in the right of the network. Networks of other eight TCM prescriptions were also presented, and COVID-19-related genes were identified (Supplementary Figures S1–S8). According to the degree of these molecules and COVID-19-related genes in the networks, the top 10 of them such as FPQXZ were screened (Tables 1 and 2). The one-to-one correspondence among molecules, COVID-19 targets, and FPQXZ were listed (Supplementary Table S1).

Table 1 The Top 10 Effective Molecules of FPQXZ Associated with COVID-19-Related Target Genes

Table 2 The Top 10 COVID-19-Related Target Genes Associated with FPQXZ

Figure 2 A network among YDBFZ, molecules of YDBFZ, and COVID-19-related target genes. Purple indicates TCMs of YDBFZ. Green indicates effective molecules in TCMs. Blue indicates target genes of molecules. Light Orange indicates COVID-19-related target genes.

Genomic Features and Functional Analysis of COVID-19 Target Genes

In total, we identified 26 genes associated with both COVID-19 and TCMs. In tumor samples, these genes were noticeably mutated, particularly, the missense mutation type was the majority (Figure 3A). CD163 gene, which was reported as a marker of M2 macrophages, presented the highest mutations among others (Figure 3B). Apart from SNVs, a large number of CNVs, such as SELP (52.06% CNVs), IL6 (43.32 CNVs), and IGFBP3 (40.00% CNVs), were observed in tumor samples (Supplementary Table S2). The expression of these genes significantly changed with amplifications and/or deletions (Supplementary Figure S9).

Figure 3 Mutation features of 26 COVID-19 target genes in the TCGA-LUAD dataset. (A) Mutation numbers of different types of mutations including splice site, nonsense, missense, frame-shift insertions and deletions. (B) Mutation counts of 26 COVID-19 target genes.

To understand the function of these identified COVID-19 target genes, WebGestaltR was used to annotate GO terms and KEGG pathways for them. A total of 1610, 61 and 128 terms of biological process, cellular component, and molecular function were annotated, respectively, with FDR < 0.05, and the top 10 of each group were visualized (Figure 4A–C). In addition, we identified 156 KEGG pathways (FDR < 0.05). Tumor-related pathways such as prostate cancer, tumor necrosis factor (TNF) signaling pathway and pathways in cancer were annotated, and IL-17 signaling pathway related to the immune system was also identified (Figure 4D).

Figure 4 Functional analysis of 26 COVID-19 target genes in the TCGA-LUAD dataset. (A–C) The top 10 enriched terms of biological process, cellular component and molecular function. (D) The top 10 enriched KEGG pathways. Size indicates gene counts. FDR, false discovery rate.

The Relation Between COVID-19 Target Genes and Immune Feature

As COVID-19 target genes were reported to be closely related to immunity, we further investigated their correlations in the TCGA-LUAD dataset through ESTIMATE and CIBERSORT tools. ESTIMATE score was calculated for each gene. Spearman correlation analysis revealed that ESTIMATE score was correlated with the expression of most COVID-19 target genes (Figure 5A), especially with ADRB2 (R = 0.40, P = 2.31e-21), ALOX5 (R = 0.65, P = 2.61e-63), Interleukin-1beta (IL1B) (R = 0.48, P = 1.82e-30), CD163 (R = 0.64, P = 2.05e-60), and SELP (R = 0.42, P = 6.56e-23). Eleven of 26 genes showed significant positive correlations with ESTIMATE score, indicating that COVID-19 target genes may be closely associated with the regulation of tumor immune microenvironment. Furthermore, we introduced CIBERSORT tool to calculate enrichment score of 22 immune cells for each sample. Apart from CD4 naive T cells, other 21 immune cells all obtained corresponding scores. Using Pearson correlation analysis between enrichment score and the expression of 16 genes, a heatmap was showing strong correlations between the two was plotted (Figure 5B). In general, plasma cells, follicular helper T cells, regulatory T cells and activated NK cells were negatively correlated with COVID-19 target genes, while neutrophils were positively correlated with COVID-19 target genes. Different immune cells manifested different correlation patterns with the 26 genes, suggesting that COVID-19 target genes played different roles on these immune cells.

Figure 5 Correlation analysis between immune score and COVID-19 target genes in the TCGA-LUAD dataset. (A) Spearman correlation analysis between immune score and the expression of 26 COVID-19 target genes. Vertical axis represents immune score and horizontal axis represents gene expression. (B) Pearson correlation analysis between the enrichment of 21 immune cells and the expression of 26 COVID-19 target genes. Colors from red to blue indicate negative to positive correlations. Colors from green to purple indicate P values from high to low. *P < 0.05, **P < 0.01, ***P < 0.001.

Identifying Key COVID-19 Target Genes and Pathways Through PPI Analysis

To better understand the interactions among these 26 COVID-19 target genes, String online tool was used to analyze their correlations and CytoScape was used to visualize the interactions (Figure 6A). Based on the degree of each gene and the analysis of MCODE module, two clusters were identified to be the key components by Haircut algorithm (Figure 6B). Specifically, cluster 1 included IL6, IL1B, ICAM1, PLAU, IGFBP3, LDLR and ALOX5, and cluster 2 included TIMP1, CAV1, MMP3, SAA1, JUN, SELP, CXCL2, and IL1A. The remained 8 genes were dissociative from the two clusters. For the 15 genes within the two clusters, we enriched 31 KEGG pathways including tumor-related pathways (such as NF-kappa B signaling, NOD-like receptor signaling, TNF signaling pathways, and transcriptional misregulation in cancer) and immune-related pathways (such as cytokine–cytokine receptor interaction, IL-17 signaling and Th17 cell differentiation) (Figure 6C).

Figure 6 Analysis of interactions among 26 COVID-19 target genes. (A) A PPI network among the 26 genes. Higher dot size indicates more interactions with other genes. (B) Three groups (clust1, clust2 and genes outside clusters) identified by MOCODE. (C) Enrichment analysis of genes in clust1 and clust2. Red indicates genes in clust1 and clust2. Blue dots indicate enriched KEGG pathways.

COVID-19 Target Genes are Associated with Cell Death

Subsequently, we tried to know if there was a relation between COVID-19 target genes and cell death. Three classes of genes related to cell death were obtained from MsigDB database, including pyroptosis, necroptosis and ferroptosis. The enrichment score of three types of cell death was calculated by ssGSEA. Spearman correlation analysis showed that 25 of 26 genes (except for ALDH2) were markedly correlated with at least one type of cell death (Figure 7). Especially, 13 of 26 genes, including ARDB2, PLAU, ICAM1, IL6, MMP3, CAV1, IL1B, THBD, IL1A, TIMP1, CD163, EPHP2 and SAA1, were significantly related to either one type of cell death (Figure 7 and Supplementary Figures S10–S12).

Figure 7 Pearson correlation analysis between 26 COVID-19 target genes and cell death including ferroptosis, necroptosis and pyroptosis. Blue indicates negative correlation and red indicates positive correlation. *P < 0.05, **P < 0.01, ***P < 0.001.

Constructing Two Molecular Subtypes Based on 26 COVID-19 Target Genes

As 26 COVID-19 target genes were highly involved in tumor-related pathways, immune response and cell death in LUAD, we then attempted to construct molecular subtypes based on their expression. ConsensusClusterPlus was used to conduct unsupervised consensus clustering. Cluster number k was set from 2 to 10, and CDF and consensus matrix were used to determine the optimal cluster number (Figure 8A and B). In the TCGA-LUAD dataset, samples were clearly clustered into the two groups when k = 2 (Figure 8C). Survival analysis showed that the two groups (Clust1 and Clust2) had differential overall survival with Clust1 had a better prognosis in the TCGA-LUAD dataset (P = 0.01, Figure 8D). In another independent dataset (the GSE31210), the same results were presented (P = 0.00026, Figure 8E).

Figure 8 Unsupervised consensus clustering based on 26 COVID-19 target genes. (A and B) Consensus CDF curve and relative change in area under CDF curve when cluster number k = 2–10. (C) Consensus matrix when k = 2. (D and E) Kaplan–Meier survival analysis of Clust1 and Clust2 in the TCGA-LUAD (D) and the GSE31210 (E) datasets. Log-rank test was conducted.

TME of Two Molecular Subtypes

To understand whether 26 COVID-19 target genes had an effect on TME, we analyzed the proportion of 22 immune cells in the two groups. As a result, 12 of 21 immune cells (naïve CD4 T cells with no expression were excluded) had differential enrichment scores in the TCGA-LUAD dataset (P < 0.05, Figure 9A and B). No significant difference in CD8 T cells was shown between the two groups. However, M0 and M1 macrophages were more enriched in Clust2, and M2 macrophages were more enriched in Clust1. Within 10 oncogenic pathways, 7 pathways were differentially enriched between Clust1 and Clust2 (P < 0.05, Figure 9C). Interestingly, except for the cell cycle and PI3K signaling pathway, other 5 pathways were all lower enriched in Clust2, which seemed contradictory with its unfavorable prognosis.

Figure 9 TME features of two molecular subtypes in the TCGA-LUAD dataset. (A) A heatmap of the proportion of 22 immune cells in two subtypes. (B) The enrichment score of the 22 immune cells in two subtypes. (C) The enrichment score of 10 oncogenic pathways. (D) Expression of 47 immune checkpoints in two subtypes. (E and F) Expression of chemokines (E) and chemokine receptors (F) in the two subtypes. Student t-test was conducted between the two groups. ns, no significance. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

To find out the reason for the above observations, we further assessed 47 immune checkpoints identified to play an important role in cell–cell communication in TME. The result showed that nearly half of the immune checkpoints were differentially expressed between two Clust1 and Clust2 (P < 0.05, Figure 9D). Notably, PDCD1 (PD-1), CD274 (PD-L1), IDO1 and LAG3, which were considered as immunosuppressive factors, were significantly upregulated in Clust2, and this may contribute to a worse prognosis of Clust2. In addition, we evaluated the expression of chemokines and chemokine receptors. It was observed that most of chemokines were more expressed in Clust2 compared with Clust1 (Figure 9E). However, chemokine receptors apart from CCR10 were more enriched in Clust1 in 8 differentially expressed receptors (Figure 9F). We also observed similar results on the GSE31210 dataset (Supplementary Figure S13), demonstrating that COVID-19 target genes could affect TME.

Sensitivity of Two Subtypes to Immunotherapy and Chemotherapy

As the TME was different in Clust1 and Clust2, we next explored if there was a difference in immunotherapy as well. SubMap analysis was applied to compare the similarity of expression patterns between the non-treated dataset and the treated dataset (IMvigor210). The result showed that the expression pattern of Clust1 was more similar to the anti-PD-1-treated samples with a status of complete response (CR) and partial response (PR) (P = 0.091), while Clust1 was more similar to stable disease (SD) and progressive disease (PD) (P = 0.001, Figure 10A). It could be speculated that Clust1 may benefit more from immunotherapy than Clust2. In addition, we assessed the sensitivity of the two groups to five chemotherapeutic drugs, including cisplatin, erlotinib, rapamycin, sunitinib and MG-132. Compared to Clust1, Clust2 showed a lower estimated IC50 of the five drugs, indicating that Clust2 was more sensitive to chemotherapy (P < 0.01, Figure 10B–F). The results were further verified in the GSE31210 dataset (Figure 10G–L), indicating that two subtypes had differential responses to targeted therapies.

Figure 10 Differential sensitivity of the two subtypes to immunotherapy and chemotherapy. (A) SubMap analysis for comparing expression patterns between the TCGA-LUAD and IMvigor210 datasets grouped by CR/PR and SD/PD. (B–F) Comparison of estimated IC50 of five chemotherapeutic drugs between Clust1 and Clust2 in the TCGA-LUAD dataset. (G) SubMap analysis for comparing expression patterns between the GSE31210 and IMvigor210 datasets grouped by CR/PR and SD/PD. (H–L) Comparison of estimated IC50 of five chemotherapeutic drugs between Clust1 and Clust2 in the GSE31210 dataset. The gray dots in the figure represent tumor samples. Student t-test was conducted between the two groups. **P < 0.01, ****P < 0.0001.

Discussion

The present identified 26 COVID-19 target genes related to TCM, COVID-19 and LUAD. Functional analysis showed that these genes were significantly enriched in tumor-related pathways including IL-17 signaling, TNF signaling, prostate cancer and pathways in cancer. Numerous studies showed that IL-17 plays an important role in tumor development, especially in TME regulation.^29,30 IL-17 secreted from type 17 cells mediates the recruitment of immunosuppressive cells dependent on the collaboration with chemokines.³⁰ Moreover, IL-17 signaling in TME can lead to epithelial-mesenchymal transition (EMT) via STAT3 and NF-κB signaling pathways in LUAD.^31,32 Those data indicated that the 26 genes may be involved in the development of LUAD.

As for the relation between the 26 genes and TME, Pearson correlation analysis revealed that many genes were strongly correlated with immune infiltration and the proportion of certain immune cells, especially ADRB2, ALOX5, IL1B, CD163, and SELP. IL1B and CD163 have been widely reported to be associated with tumor development and TME. A high expression of IL1B could increase the risk of non-small cell lung cancer (NSCLC).³³ CD163 is considered as one of the markers of M2 macrophages related to worse outcomes.^34,35 Our result also showed a positive correlation between CD163 expression and M2 macrophages. PPI analysis presented that IL6 and IL1B showed the most interactions with other genes, indicating that they played central roles in regulation network among these genes.

Combining the above two results, we found that IL1B may be an important gene, the expression level of which could largely affect the interactions in the network and TME. Th17 cells are one of the major cells secreting IL-17. IL1B, a pro-inflammatory cytokine, drives Th17 differentiation together with transforming growth factor-beta (TGF-β), IL-6, and IL-21.³⁶ Chen et al have found that lung cancer microparticles (L-MPs) induce macrophages to release IL1B, which has a tumor-promoting effect.³⁷ Mature IL1B is released from immune cells through pyroptosis-mediated cell death.³⁸ Active IL1B promotes tumor cell growth and metastasis through vascular endothelial growth factor (VEGF), prostaglandin E2 (PGE2) and TGF-β.³⁹ Furthermore, through binding with IL-1, IL1B can activate MAPK and NF-κB pathways and induce downstream effects such as releasing granulocyte-macrophage colony-stimulating factor related to the activation of myeloid-derived suppressor cells (MDSCs) and M2 macrophages.⁴⁰ Therefore, aberrant IL1B expression and signaling can alter the formation of TME and then promote tumor cell invasion and migration.

In addition, the close relation between cell death and 26 target genes were illustrated. Among the three types of cell death, necroptosis and pyroptosis were obviously shown to be positively correlated with the 26 target genes, especially PLAU, ALOX5, IL1B, IL1A and CD163. As indicated before, IL1B is one of the main products of pyroptosis,⁴¹ and IL-1β secretion can also be triggered by necroptosis,⁴² which is in accordant with our results. To some extent, cell death can affect the modulation of TME through inducing IL-1β. Kaplanov et al have found that IL1B inhibition can restore an effective immune response in breast cancer, and has a synergetic effect when combined with anti-PD-1 therapy in mice model.⁴³ In IL1B-deficient mice, a high anti-tumor activity was shown, such as increased activated CD8 T cells and dendritic cells.⁴³ We speculated that IL1B may be also a potential target for treating LUAD, but additional experiments are needed in future studies.

Based on the expression of 26 target genes, we constructed two novel molecular subtypes showing distinct prognosis. Notably, the two subtypes manifested differential TME, including the proportion of immune cells, expression of immune checkpoints, chemokines and chemokine receptors. The observations indirectly proved the role of 26 target genes in TME modulation and tumor development. Importantly, the two subtypes showed differential sensitivity to the immune checkpoint blockade and chemotherapy, which may guide clinical treatments for LUAD patients with COVID-19. Those data may prove worthwhile in developing new strategies for combating COVID-19 for LUAD.

Conclusions

Although potential COVID-19 targets predicted for LUAD therapy were identified from a large number of samples applying bioinformatics approaches, some limitations still exist in this study. Firstly, the sample lacked clinical follow-up information; thus, factors such as the presence of other health conditions were not considered during the identification of the biomarkers. Secondly, the results obtained by bioinformatics analysis alone were not convincing enough, which requires further experimental verifications. Genetic and experimental studies with larger sample sizes and experimental validation are also needed. In conclusion, this study proposed 26 potential COVID-19 targets, particularly IL1B, which may serve as a therapeutic drug target for developing effective molecular drugs for LUAD patients with COVID-19. The effective molecules in TCM interacting with IL1B may be considered as potential drugs, and can also guide the exploration of new therapeutic drugs for LUAD patients. In addition, the two new developed molecular subtypes can provide guidance for target therapies in clinical practice.