Transcriptome profiles of STAD patients were downloaded from The Cancer Genome Atlas (TCGA) database in November 2021. The differential genes in STAD were screened and obtained using the ‘limma’ package of R-language Bioconductor. Furthermore, genes related to COVID-19 were selected from the OMIM database, Genecard database, KEGG and NCBI databases. Finally, these genes were compared to obtain the overlapping genes in STAD and COVID-19 (Li et al., 2020; Li et al., 2021a).

VD3-related genes were collected from accessible online tools, such as Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP), Swiss Target Prediction, TargetNet and Drugbank. Protein-gene name conversion was then performed using the protein database UniProt. Finally, a list of related genes was obtained by removing duplicates of these genes.

The correlation of STAD/COVID-19-related genes with survival in STAD/COVID-19 patients was analyzed in the ‘survival package’ of the R package. Prognostic analysis of STAD/COVID-19 patients was performed using univariate and multivariate Cox proportional hazards regression. In addition, we analyzed the different characteristics of patients with STAD/COVID-19 (Fisher and Lin, 1999).

Using R-language packages, including ‘ClusterProfiler', ‘org.Hs.eg.Db’, ‘pathview’, and ‘GOplot’, etc. for enrichment analysis and visualization of intersecting genes by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO enrichment analysis annotated the mechanism of action of target genes in three parts: biological process (BP), cellular composition (CC) and molecular function (MF), and analyzed the top ten entries of each part according to their p values. KEGG enrichment analysis was performed to annotate the signaling pathways in which the target genes were involved.

Protein-protein interaction network (PPI) and TSV files were obtained by importing the intersecting genes into the STRING database. Next, the top five proteins of the network in terms of attribute ranking were obtained using the MCC calculation method of the cytoHubba plugin in Cytoscape software (Chin et al., 2014).

The protein structures of COVID-19 associated proteins and target proteins were obtained from the PDB database (https://www.rcsb.org/) (Zardecki et al., 2022). The molecular structures of VD3 and Calcitriol were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (Li et al., 2021a). The corresponding proteins and small molecules were processed and docked using Discovery studio-2019 (DS). To ensure the accuracy of the docking prediction results, we used DeepSite to predict the possible docking pockets of the target proteins. DeepSite is freely available at www.playmolecule.org (Jiménez et al., 2017). DeepSite considers various molecular descriptors related to the protein through the 3D-deep convolutional neural networks (DCNNs) validated using an extensive test set based on over 7,000 proteins from the scPDB database (Kellenberger et al., 2006; Jiménez et al., 2017). The cartesian coordinates of the center of the binding pocket found by DeepSite were used to define site in the DS with a radius of 12. And for other parameters, defaults were used. CDOCKER was used for small molecular-protein docking, Top Hit was set to 10, Pose Cluster Radius was set to 0.5. Finally, the interaction and binding modes of the active compounds were visualized and analyzed using PyMoL-2.1.0 and DS.

The flow of the study was exhibited in Figure 1. From The Cancer Genome Atlas (TCGA) database, we obtained 32 normal samples and 375 tumor samples. When duplicated genes were excluded from all samples, the obtained genes were processed with the egdeR package for normalization and differential expression analysis. An adjusted p-value of < 0.05 and a fold change value higher than 1.5 was used to indicate genes that were significantly differentially expressed. We obtained 1,439 differentially expressed genes (DEGs) after analysis and mapped the volcano plot (Figure 2A). After that, we collected 2,537 genes associated with COVID-19 from the Genecard, OMIM, KEGG and NCBI databases by the method of network pharmacology. Comparing these two gene clusters, we identified 137 intersection genes for STAD and COVID-19 (Figure 2B). Of them, 71 were upregulated, and the others were downregulated (Figure 2C).

To understand the correlation between STAD/COVID-19-related genes and prognosis, univariate and multifactorial Cox analyses were performed on 137 DEGs. First, univariate Cox analysis identified 15 genes that were significantly (p < 0.05) associated (Table 1). However, multivariate Cox analysis identified 4 of these genes, including KIT, HBA1, SERPINE1 and FKBP10 (Table 2). Results in table 2 showed that the 4 genes were bad prognostic factors. In addition, hazard ratios (HR) and 95% confidence intervals (CI) for risk scores from univariate Cox (uni-Cox) regression was 1.21 and 1.136–1.29, respectively, (p < 0.001). In comparison, the values of multivariate Cox (multi-Cox) regression were 1.189 and 1.115–1.268, respectively, (p < 0.001) (Supplementary Figures S1A,B). In addition, we found two other independent prognostic parameters: age (1.033 and 1.015–1.052; p < 0.001) and stage (1.586 and 1.115–1.268; p < 0.001) (Supplementary Figure S1B).

We investigated the correlation between the expression of the four genes obtained in Table 2 and the prognosis of STAD patients through the GEPIA online platform (http://gepia2.cancer-pku.cn/#index). We found that STAD patients with lower SERPINE1 levels exhibited significantly longer overall survival (OS) and disease-free survival (DFS) than those with higher SERPINE1 levels (p < 0.01). Likewise, STAD patients with lower levels of FKBP10 or KIT showed a better prognosis (p < 0.05, Figure 3). The above results indicated that SERPINE1, FKBP10 and KIT were closely related to the prognosis of STAD patients.

We downloaded the 3D molecular structure and SMILES [CC(C)CCC​C(C)C1CCC2C1(CCC​C2 = CC = C3CC(CCC3 = C)O)C] of VD3 from PubChem. The PubChem CID of VD3 is 5280795. The chemical formula and molecular weight of VD3 are C₂₇H₄₄O and 384.6 g/mol, respectively. Then, a portion of the major target genes of VD3 were searched from the Swiss Target Prediction databases and TargetNet web server. Meanwhile, we also collected a portion of VD3-related target genes from GeneCards databases and TCMSP databases. After name conversion of the targets and removing duplicates, 2274 VD3-related target genes were obtained. Identification of 60 VD3 intersecting genes against STAD/COVID-19 after overlap of STAD/COVID-19-related genes with VD3-related target genes (Figure 4A, Supplementary Table S1). By the same method, we identified 52 calcitriol anti- STAD/COVID-19 intersecting genes Supplementary Figure S2A.

For a further, more comprehensive and macroscopic understanding of the biological functions of the intersection genes, we performed GO functional enrichment analysis and KEGG pathway enrichment analysis on these 60 genes. The results showed that VD3 affected 477 biological processes (BP), 32 cellular components (CC), 56 molecular functions (MF) (p < 0.05, Supplementary Table S2), involving response to lipopolysaccharide, response to molecule of bacterial origin, cell chemotaxis, leukocyte migration, and other biological processes; vesicle lumen, secretory granule lumen, cytoplasmic vesicle lumen, and other cellular components; receptor ligand activity, signaling receptor activator activity, and other molecular functions (Figure 4B). The 10 biological process terms with the lowest adjusted p values in each category are shown in the dot plot. And the above results implied the possible main functions of the 60 intersecting genes.18 signaling pathways associated with all intersection genes were obtained by the KEGG pathway enrichment analysis (p < 0.05, Supplementary Table S3), and shown in dot plot according to the gene ratios (Figure 4C). The results showed 60 intersecting genes focused on transcriptional misregulation in cancer, coronavirus disease-COVID-19, IL-17 signaling pathway and cytokine-cytokine receptor interaction pathway. The above results implied the main pathways involved in the 60 intersecting genes. Notably, SERPINE1 and FKBP10 genes were found in the hematopoietic cell lineage (hsa04640), acute myeloid leukemia (hsa05221) and complement and coagulation cascades (hsa04610) pathways. Taken together, these results demonstrated that VD3 and calcitriol (Supplementary Figures S2B,C) act on COVID-19 and STAD through multiple pathways, multiple targets, and overall cooperation.

Determination of the 60 intersecting target gene-mediated PPI-networks of VD3 against COVID19/STAD using STRING analysis (Figure 5A). Ranked network centrality nodes using the Maximal Clique Centrality (MCC) algorithm of the cytoHubba plugin with 60 intersection genes input to Cytoscape software. Figure 5B showed the connectivity of the top five hub nodes in the network, the darker the node color means the higher the score. The results showed that the top five core target genes were IL1A, CXCL8, ALB, CSF2 and SERPINE1. Furthermore, the hub nodes appeared in 16 pathways in Figure 4C except for neutrophil extracellular trap formation and renin-angiotensin system pathways (Supplementary Table S3).

Molecular docking was performed to visualize the patterns of interactions between the VD3 and the identified STAD/COVID-19 target proteins (Figure 5C). We selected the SARS-COV-2 Spike RBD/ACE2 complex (PDB:6LZG) (Wierbowski et al., 2021) and main protease complex crystal structure of COVID-19 (PDB:5R84) (Li et al., 2021a) for molecular docking analysis with VD3. These crystal structures of the COVID-19 were collected from the PDB database (http://www.rcsb.org). In order to predict the VD3-protein binding site, we performed binding pocket simulations of the target structures in this study by referring to the methods applied in the literature for natural product binding mode studies (Singh et al., 2021). In this paper, pockets were predicted using DeepSite and pocket locations were visualized using PyMoL (Supplementary Figure S4). Supplementary Table S4 shows the scores for each prediction pocket.

Two cavity coordinates were detected by DeepSite in SARS-CoV-2 RBD/ACE2 (Supplementary Figure S4A). We believed that interference in protein-protein interactions occurs mainly on the protein surface, and the key amino acid residue (Y41, K353, N501, Y505) reported in the literature that affects the binding of Spike (S) glycoprotein to ACE2 is located near cavity 2 (Supplementary Figure S5) (Lan et al., 2020), so we selected cavity 2 for molecular docking. Our results showed that the hydroxyl group of VD3 forms hydrogen bond with D405 of SARS-CoV-2 RBD/ACE2, and hydrophobic interactions were also found. (Figure 6A). Notably, the binding pocket is located in the protein-protein interaction region of the SARS-COV-2 Spike RBD/ACE2 complex. Based on the docking results, we hypothesized that VD3 may play an inhibitory role in the binding of the Spike glycoprotein of SARS-CoV-2 to the host angiotensin-converting enzyme 2 (ACE2) receptor. We also did the molecular docking of calcitriol with SARS-CoV-2 RBD/ACE2. The results showed that the hydroxyl groups of calcitriol form hydrogen bonds with G504 and Q325 of SARS-CoV-2 RBD/ACE2 (Supplementary Figure S6A). Four cavity coordinates were detected by DeepSite in SARS-CoV-2 Mpro. We selected the highest scoring cavity 2 for molecular docking (Supplementary Figure S4B). The results showed that the hydroxyl group of VD3 forms hydrogen bond with E166 of Mpro, and the hydrophobic interactions were also found. (Figure 6B). Similarly, we did the molecular docking of calcitriol with Mpro. The results showed that the hydroxyl groups of calcitriol could also form hydrogen bonds with Q189, T190, T26 of Mpro (Supplementary Figure S6B). Our results also showed that the amino acid residue VD3 bound to Mpro was consistent with the amino acid residue bound the proto-ligand 2-cyclohexyl-N-(pyridin-3-yl) acetamide (GWS) to Mpro, which was E166 (Supplementary Figure S6D).

In addition, we analyzed the possible binding between VD3 and the top 5 STAD/COVID-19 core target proteins (IL1A, IL8, ALBU, CSF2, PAI1) identified in the STRING analysis. Similarly, we also considered the binding of calcitriol. Our result showed that VD3 only bind to PAI1, as does calcitriol. Four cavity coordinates were detected by DeepSite in PAI1. We selected the highest scoring cavity 1 for molecular docking (Supplementary Figure S4C). The hydroxyl group of VD3 was bound to Q322 via 1 hydrogen bonding to PAI1 (PDB: 1A7C). And hydrophobic interactions were also found (Figure 6C). The hydroxyl groups of calcitriol were bound to N167, Q322, R163 and L163 of PAI1 by hydrogen bonding (Supplementary Figure S6C). The scores of each docking posture were listed in Table 3.