2.1 Screening of active compounds of EH-ASA herb pairs

The active compounds of EH-ASA were obtained from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (Ru et al., 2014) (https://tcmsp-e.com), containing all herbs registered in Chinese Pharmacopoeia, which came into being to promote the integration of modern medicine and traditional medicine in drug research, and further development of TCM. It provides a total of twelve related properties on absorption, distribution, metabolism, and excretion to evaluate herbal medicines with potential biological efficacy, among which oral bioavailability (OB) and drug-likeness (DL) are leading indicators to appraise pharmacokinetics of a particular ingredient and whether it takes effect. We know that the higher the OB and DL scores, the better the component's pharmacological activity is. Thus, we set the screening conditions as OB > 30 % (Xu et al., 2012) and DL > 0.18 (Tao et al., 2013) in our study.

2.2 Predicting of targets of EH-ASA herb pairs

After getting the active ingredients of EH-ASA, we downloaded its corresponding target protein, namely “Target name” in TCMSP, and then the “Target name” was standardized as “Gene name” in the Universal Protein (UniProt) database, and the duplications were removed (https://www.uniprot.org) (UniProt, 2021). Finally, the gene symbol related to EH-ASA was obtained.

2.3 The construction of the Herb-Compound -target network

We used the Cytoscape software (Version.3.6.2) (Shannon et al., 2003) to construct the network of active compounds corresponding target genes from herbs in TCMSP that satisfied the threshold for OB and DL, so that we could observe the core relationship of the compound-target-herb network. The CytoNCA plug-in was used to calculate each node parameter in the network. Deeper color and larger node areas represented more robust connectivity of the node and its neighbors, suggesting the effective critical compound and key candidate target.

2.4 Identification of COVID-19 differentially expressed genes

We browsed the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013) with “COVID-19” as the search term and finally selected the dataset of GSE180226 as following subject. The dataset contains 23 cases of lung tissue gene expression, pathology, and plasma protein as clinical sample information, of which 20 samples were from the disease group and others were from the control. The screening of differentially expressed genes (DEGs) between the two groups and the standardization of DEGs were performed using the “limma” package of R language software (Version 3.6.3, https://www.r-project.org/). To identify genes with significant differences in expression, the screening criteria were set as fold change (FC) value (|log2FC|) > 2 and P value < 0.05 (Li et al., 2022).

The standardization results were shown in box plots (Supplementary Fig. 1A–B), while the significantly expressed genes of the GSE180226 dataset were visualized with a heatmap and volcano plot by the “ggplot2” package in R software (Ito and Murphy, 2013). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

2.5 Construction of protein–protein interaction (PPI) network

The Venn online mapping tool (https://bioinformatics.psb.ugent.be/cgi-bin/liste/Venn/) was used to intersect the DEGs and EH-ASA herb pairs targets, namely common genes, which were regarded as potential targets against COVID-19, and were prepared for constructing PPI network. The list of genes obtained above was input into the STRING database (https://string-db.org/), with the species set as “homo sapiens”, the minimum required interaction score set as more than or equal to 0.4 (Szklarczyk et al., 2019), free targets hidden, and the rest of the parameters remained default. The related parameters in the PPI network were saved in “tsv” file format and then imported into Cytoscape software to visualize the relationship between targets of EH-ASA and COVID-19 genes. Subsequently, the nodes revealed critical genes according to the degree value calculated by the “NetworkAnalyzer” plug-in.

2.6 Functional enrichment analysis of common genes

To further illuminate the biological mechanism of active compounds of EH-ASA herb pairs on COVID-19, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted. First, we imported common genes into the DAVID platform (https://david.ncifcrf.gov/) for GO and KEGG analyses. Then, the functional annotation tables were downloaded and arranged in ascending order according to P value. Notably, the GO analysis examined the following three aspects: biological process (BP), cellular component (CC), and molecular function (MF). Then, we visualized the top 10 common genes of topological analysis and significant terms of GO annotation (top 10 of each category) and KEGG analysis (top 20). Finally, the relationship between hub genes and possible mechanisms was visualized using the function of “SankeyNetwork” of the “Goplot” package of R language software with a Sankey diagram. The P value < 0.05 was considered a threshold for significant terms.

2.7 Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) is a widely used method to discover the Gene Characteristic Modules (GCMs) and key genes by connecting gene networks to clinical traits realized by the scale-free topology criterion. The characteristic module is the first essential element of a GCM, which could be used to indicate the whole expression of the GCM. In this study, we applied the “WGCNA” package of R software (Langfelder and Horvath, 2008) to construct the WGCNA co-expression network of the DEGs from the GSE180226 dataset, with healthy control and COVID-19 as clinical traits. Firstly, the outlier cases were filtered by hierarchical clustering analysis. Secondly, the pickSoft Threshold function was conducted to measure the soft threshold, also called soft power and adjacencies. So far, we have obtained an adjacency matrix, which was subsequently performed into a topological overlap matrix (TOM). The dynamic tree-cutting means was used to determine the modules in the co-expression network. We then appraised the correlation between gene modules and COVID-19 by gene significance (GS) and module membership (MM). The higher the values, the closer the connection between genes and clinical characteristics. At last, we merged the genes of the key module with the most excellent Pearson correlation coefficient and DEGs, namely hub genes, for further investigation.

2.8 Functional enrichment analysis of hub genes

We picked up the hub genes by crossing the genes in the key module obtained by the WGCNA method with DEGs, which were input into the DAVID online analysis platform. The results were visualized as circle plots through the “GOplot” package in R software. In addition, we exerted the gene set enrichment analysis (GSEA) to assess the differences in the stimulation between the COVID-19 samples and phenotypes. GSEA on DEGs was selected with P < 0.05 as a filtration standard to identify the association between hub genes and significant enrichment terms (Subramanian et al., 2005).

2.9 Identification of signature genes

After a series of filtration, we initially acquired common genes via the intersection of EH-ASA herb pairs targets and DEGs of the GSE180226 dataset; next, we also got hub genes of COVID-19 by merging DEGs and genes of the characteristic module with the highest correlation coefficient. Now, we were going to identify the candidate signature genes by combining common genes with hub genes, regarded as signature genes of EH-ASA that could treat COVID-19 with high probability. The diagnostic accuracy of screened signature genes was evaluated by employing the area under the curve (AUC) of the receiver operating characteristic curve (ROC). An AUC was more significant than 0.7, which manifested preferable diagnostic performance.

2.10 Molecular docking

The binding situation and interaction force between the protein target, usually receptor, and small molecule compound, usually ligand, could be predicted and obtained with molecular docking analysis (Pinzi and Rastelli, 2019). In this study, the active compounds from EH-ASA obtained through topological analysis and signature genes acquired by WGCNA were finally picked for molecular docking analysis. Furthermore, we selected another two protein targets, including apin-like protease (PLpro) and main protease (Mpro), from the COVID-19 Docking Server website (https://ncov.schanglab.org.cn/index.php) (Kong et al., 2020), which were considered as promising targets for discovering novel anti-viral drugs. Then, the 2D structure of drug ligands and protein receptors were separately downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (Wang et al., 2017) and Protein Data Bank (PDB) database (https://www.rcsb.org/) (Wang et al., 2012). Small molecule energy minimization was performed to optimize conformation using PyMol 2.4.0 software. Autodock Vina software was conducted to delete water, add perhydro, calculate charges, detect and choose a torsion tree, an active docking centre via the grid box function, and then processed the molecular docking with default parameters (Trott and Olson, 2010). Finally, the binding energy of docking results was shown in a heatmap. We also chose the best affinity as the final docking conformation and visualized it in Pymol 2.4.0 software (Seeliger and de Groot, 2010) and LigPlot 2.2.8 software.

3.1 The screening of active compounds and targets of EH-ASA herb pairs

A total of 42 qualified chemical compounds of EH-ASA were collected in TCMSP, among which Stigmasterol and (+)-catechin were the mutual ingredients of EH and ASA, and the details of these screened constituents are shown in Table 1. Subsequently, we identified the targets corresponding to the above compounds in the UniProt KB database. After discarding the repetition, we finally obtained 200 potentially therapeutic targets of EH-ASA. Whereas compounds MOL001506 from EH and MOL010922, MOL002372 and MOL003410 from ASA could not relate to their targets, so we removed these four chemicals and then built a network about the interaction of active compounds (C) and targets (T) of EH and ASA (drug, D).

As shown in Fig. 2, it was found that there was a one-to-one correspondence between 42 active constituents and 200 targets of EH-ASA by interconnection in this complex C-T-D network, which contained a total of 238 nodes and 580 edges, demonstrating that EH and ASA could directly or indirectly exert the therapeutic influence on a mass of targets.

Close

Fig. 2

C-T-D network of EH-ASA herb pairs.

Export to PPT

3.2 Identification of different expression genes related to COVID-19

Compared with the healthy samples, there were 548 genes with significant differences in expression of the subjects of COVID-19 total, among which 137 genes were highly up-regulated, while 411 genes were very down-regulated in the GSE180226 dataset. A heat map and a volcano diagram were drawn with the screened DEGs to demonstrate the distribution of these genes in this dataset between the control and COVID-19 group using “limma”, “heatmap”, and “ggplot2” packages in R software (Fig. 3A, B).

Close

Fig. 3

Results of DEGs in GSE180226 related to COVID-19. (A) Heatmap. (B) Volcano plot. The red dots represent upregulated genes, the blue dots represent downregulated genes, and the grey dots represent stable genes. FC, fold change. (C) Venn diagram of the common genes of DEGs with targets of EH-ASA herb pairs. EH, Ephedra Herba; ASA, Armeniacae Semen Amarum.

Export to PPT

In addition, the therapeutic effects of EH-ASA on COVID-19 needed to be validated. Thus, As shown in Fig. 3C, after combining the 200 targets corresponding to active constituents of 548 DEGs of the GSE180226 dataset via limma difference analysis, there were a total of 74 overlap genes identified, namely common genes, which were regarded as the potential targets in the treatment of COVID-19 and prepared for the subsequent study.

3.3 Protein-protein interaction (PPI) network construction of common genes

As shown in Fig. 4A, the PPI network was constructed with common genes, including 74 nodes and 1320 edges. To determine the interactive relationships of DEGs with active targets, we applied these 74 common genes to the STRING online platform and formed a PPI network through Cytoscape software. The result showed that multiple targets had strong interactions with each other, such as tumour protein P53 (TP53), Serine/threonine-protein kinase (AKT1), tumour necrosis factor (TNF), and recombinant caspase 3 (CASP3).

Close

Fig. 4

The results of PPI. (A) PPI network of the common genes. The targets were ranked in order based on their degree values, larger dots, and darker colors indicating greater significance. (B) Heatmap of 23 samples of GSE180226 and its corresponding DEGs of common genes.

Export to PPT

Moreover, we elucidated the interaction between targets of EH-ASA and DEGs of the GSE180226 dataset related to COVID-19. Therefore, the expression of 74 genes that were necessary from these 23 COVID-19 samples was selected, and a heatmap was then created using “circlize” and “ComplexHeatmap” packages in R software (Fig. 4B).

3.4 Functional enrichment analysis of common genes

We also applied the common genes to GO and KEGG enrichment analysis to clarify the active mechanism. As shown in Fig. 5A-B, the BP category extensively focused on cellular response to reactive oxygen species, ageing and process of apoptosis, and positive regulation of transcription from the RNA polymerase II promoter. CC was mainly involved in the cytoplasm, mitochondria, and RNA polymerase II transcription factor complex. Moreover, the targets were widely aggregated in transcription factor binding and protein binding in MF. When referred to the KEGG analysis, the most common targets were enriched in the TNF signaling pathway, Influenza A, and Toll-like receptor signaling pathway. The complex overlapping relationship between genes and pathways was also shown in the Sankey diagram, in which MAPK3 is related to most KEGG pathways, suggesting a crucial status in associating the various cross-talk with different pathways.

Close

Fig. 5

Sankey diagram of functional enrichment analysis of common genes. (A) GO; (B) KEGG. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Export to PPT

3.5 Construction of co-expression network and selection of key module

The WGCNA of COVID-19 subjects and healthy controls in the GSE180226 dataset was conducted. We finally established a scale-free co-expression network through hierarchical clustering and dynamic tree cut. A scale-free index of 0.85 and relatively good mean connectivity identified the soft threshold power as β = 19 (Fig. 6A, B). The cluster gene dendrogram is shown in Fig. 6C-E. Finally, a total of 17 differently coloured modules were obtained, among which the MEsinna3 module, including 1226 genes, was considered to be the most significantly associated with COVID-19 (r = 0.9, P < 1e-200) (Fig. 6F, G). As shown in Fig. 5H, a total of 189 mapping genes between DEGs and genes in the Mesinna3 module were shown in a Venn diagram.

Close

Fig. 6

The results of WGCNA of GSE180226. (A) The soft thresholding power of scale independence. (B) The soft thresholding power of mean connectivity values. (C) Heatmap of the sample cluster. (D) Clustering of module eigengenes of correlation diagram of modules. (E) Merged dynamic tree cut diagram of all DEGs according to the cluster. The different color stripes show the results of the monolithic analysis in automation. (F) The heatmap of the correlation between module characteristic genes and phenotypes. The vertical coordinate is the correlation coefficient of feature modules, and the darker color indicates a more significant correlation. (G) The scatter diagram of module characteristic genes and phenotypes. (H) Venn diagram of DEGs with genes in the key module.

Export to PPT

3.6 Functional enrichment analysis of hub genes

To further discover the pathogenesis of the GSE180226 dataset related to COVID-19, we performed GO and KEGG enrichment analysis with 189 hub genes (Fig. 7A-D). The BP annotation implied that positive regulation of regulatory T cell differentiation, response to epidermal growth factor, and UMP biosynthetic process were mainly enriched. Also, lateral element, brush border, and terminal bouton were primarily collected in CC. In MF, transcription coactivator binding, amino acid binding, and transcription cofactor binding were mainly involved. As for KEGG analysis, the above hub genes were primarily focused on the IL-17 signaling pathway, Human T-cell leukaemia virus one infection, Type I diabetes mellitus, and so on. In addition, the GSEA results suggested that terms of fat digestion and absorption, glycosaminoglycan degradation, linoleic acid metabolism, maturity onset diabetes of the young and mucin type O-glycan biosynthesis were significantly up-regulated in COVID-19 samples; meanwhile, autophagy, cysteine and methionine metabolism, legionellosis, proteasome and ribosome were activated in control samples (Fig. 7E, F).

Close

Fig. 7

Circle plot of functional enrichment analysis of hub genes. (A) BP; (B) CC; (C) MF; (D) KEGG. (E) GSEA in COVID-19 group; (F) GSEA in control group. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. UMP, uridine monophosphate.

Export to PPT

3.7 Diagnostic accuracy of signature genes in COVID-19

We further performed difference analysis to determine the diagnostic efficacy of signature genes screened through the PPI and co-expression networks of WGCNA. The result showed that these genes significantly distinguished the subjects infected with COVID-19 from those who were healthy, indicating that these genes from EH-ASA might perform a good remedy for COVID-19 (Fig. 8A-C). Among them, ESR1 and IL2 were increased in COVID-19 samples (P < 0.01), while NFKBIA increased in control samples significantly (P < 0.01). Furthermore, the AUC of the signature gene was 1.00 of ESR1, 0.967 of IL2, and 0.967 of NFKBIA, respectively (Fig. 8D-F), suggesting that these three signature genes might have good diagnostic accuracy in detecting COVID-19 and might be against the virus.

Close

Fig. 8

The performance of signature genes. The boxplot of gene expression of ESR1(A), IL2(B) and NFKBIA(C). The ROC of diagnostic accuracy of ESR1(D), IL2(E) and NFKBIA(F). Con, Control; Treat, COVID-19.

Export to PPT

3.8 Molecular docking

The scores of molecular docking between the main active compound of EH-ASA and signature genes were shown in Fig. 8. It is known that the lower the binding energy, the stronger the docking affinity is, indicating that the compounds were closer to the targets and were more effective (Pantsar and Poso, 2018). The corresponding PDB IDs were ESR1 (1L2I), IL2 (4NEJ), NFKBIA (1NFI), Mpro (6LU7), and PLpro (6WUU), respectively, and the details of main active compounds for molecular docking were shown in Table 2.

Close

Fig. 9

Heatmap of Molecular docking.

Export to PPT

Table 2 Details of main active compounds for molecular docking.

MOL ID	Compounds	Pubchem ID	Molecular formula	Source	2D Structure
MOL002311	Glycyrol	5320083	C21H18O6	Armeniacae Semen Amarum
MOL004841	Licochalcone B	5318999	C16H14O5	Armeniacae Semen Amarum
MOL004908	Glabridin	124052	C20H20O4	Armeniacae Semen Amarum
MOL005017	Phaseol	44257530	C20H16O5	Armeniacae Semen Amarum
MOL000006	Luteolin	5280445	C15H10O6	Ephedra Herba
MOL000098	Quercetin	5280343	C15H10O7	Ephedra Herba
MOL004328	Naringenin	932	C15H12O5	Ephedra Herba
MOL010788	Leucopelargonidin	3286789	C15H14O6	Ephedra Herba

The combined scores of the top six conformations were in the rankings, Glycyrol-ESR1, Glabridin-ESR1, Glabridin-PLpro, Phaseol-ESR1, Luteolin-ESR1, Glabridin-Mpro, suggesting that the above-mentioned active compounds of EH-ASA might possess an excellent affinity with COVID-19 related genes.

In addition, hydrogen bonds and hydrophobic contacts were recorded to meticulously observe the interactions between protein and ligand, which were considered the most common type of interactions (Ferreira de Freitas and Schapira, 2017). As shown in Fig. 10A-L, Glycyrol and Glabridin are bound to ESR1 with the amino acid residues GLU-353 and ARG-394, forming three hydrogen bonds and two hydrophobic contacts. Glabridin bound to PLpro with the amino acid residues LYS-306, forming three hydrogen bonds and two hydrophobic contacts. Phaseol is tied to ESR1 with the amino acid residues ARG-394 and GLU-353, forming three hydrogen bonds and two hydrophobic contacts. Luteolin bound to ESR1 with the amino acid residues ARG-394, GLU-353, LEU-387, LEU-346, and GLY-521, forming three hydrogen bonds and nine hydrophobic contacts. Glabridin bound to Mpro with the amino acid residues GLU-71 and MET-17, forming three hydrogen bonds and four hydrophobic contacts.

Close

Fig. 10

Results of Molecular docking. Molecular docking models(3D and 2D) with the lower binding energy of sites of Glycyrol with ESR1(A, B); Site of Glabridin with ESR1(C, D); Site of Glabridin with PLpro(E, F); Site of Phaseol with ESR1(G, H); Site of Luteolin with ESR1(I, J); Glabridin with Mpro(K, L). GLU, Glutamate; ARG, Arginine; LYS, Lysine; LEU, Leucine; GLY, Glycine; MET, Methionine; PLpro, Papain-Like protease; Mpro, Main protease; Green dotted lines, hydrogen bonds; Red combs, hydrophobic contacts.

Export to PPT