Effects of Bacille Calmette Guerin (BCG) vaccination during COVID-19 infection.

The coronavirus disease 2019 (COVID-19) is caused by the infection of highly contagious severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), also known as the novel coronavirus. In most countries, the containment of this virus spread is not controlled, which is driving the pandemic towards a more difficult phase. In this study, we investigated the impact of the Bacille Calmette Guerin (BCG) vaccination on the severity and mortality of COVID-19 by performing transcriptomic analyses of SARS-CoV-2 infected and BCG vaccinated samples in peripheral blood mononuclear cells (PBMC). A set of common differentially expressed genes (DEGs) were identified and seeded into their functional enrichment analyses via Gene Ontology (GO)-based functional terms and pre-annotated molecular pathways databases, and their Protein-Protein Interaction (PPI) network analysis. We further analysed the regulatory elements, possible comorbidities and putative drug candidates for COVID-19 patients who have not been BCG-vaccinated. Differential expression analyses of both BCG-vaccinated and COVID-19 infected samples identified 62 shared DEGs indicating their discordant expression pattern in their respected conditions compared to control. Next, PPI analysis of those DEGs revealed 10 hub genes, namely ITGB2, CXCL8, CXCL1, CCR2, IFNG, CCL4, PTGS2, ADORA3, TLR5 and CD33. Functional enrichment analyses found significantly enriched pathways/GO terms including cytokine activities, lysosome, IL-17 signalling pathway, TNF-signalling pathways. Moreover, a set of identified TFs, miRNAs and potential drug molecules were further investigated to assess their biological involvements in COVID-19 and their therapeutic possibilities. Findings showed significant genetic interactions between BCG vaccination and SARS-CoV-2 infection, suggesting an interesting prospect of the BCG vaccine in relation to the COVID-19 pandemic. We hope it may potentially trigger further research on this critical phenomenon to combat COVID-19 spread.

Introduction

The World Health Organization (WHO) declared a global pandemic on March 11, 2020 for the Coronavirus disease 2019 (COVID-19) [1], caused by the highly contagious novel coronavirus (SARS-CoV-2) infection [2]. It was first detected in Wuhan, China in December 2019, although its epidemiological origin yet remains debatable. Eventually, it has spread in many countries throughout the world very quickly [3]. Since the disease is a novel one, the knowledge about its underlying mechanism is still sparse. Numerous studies have already manifested strong and consistent evidences regarding the impact of various disease conditions during COVID-19 [[4], [5], [6]], including cardiovascular diseases [7], malignancies [8], chronic kidney diseases [9], Chronic obstructive pulmonary disease (COPD) [10], type II diabetes [11] and many more.

The world population is at the highest risk due to the COVID-19. As on September 16, 2021, over 225 million confirmed SARS-CoV-2 infected cases have been reported in more than 217 countries and regions with ≈4.64 million deaths according to WHO (https://covid19.who.int/), and the grievous thing is that its mortality rate is increasing day-by-day. However, the infection pattern as well as its severity and mortality are shown some salient disparity in different regions indicating the influence of social norms and healthcare strategies. Till now, 80 vaccines against SARS-CoV-2 are on clinical trials on humans and 23 among them are at the final phase. Fortunately, 7 vaccines have been already approved for full use [12]. But still, the vaccines are not available and affordable for everyone all over the world, only 5.53 billion vaccine doses have been administered till September 16, 2021 (https://covid19.who.int/). Moreover, the efficacy of these vaccines are subjects for verification through long-term follow-up. In these circumstances, until an effective and affordable vaccine has been available for all, adaptation of existing and safe vaccines that reinforces the immunity system may be beneficial. This strategy suggesting the protective impact of the Bacillus Calmette-Guérin (BCG) vaccination on the intensity of the COVID-19 has gained considerable research focus. Therefore, it is of significant importance to conduct a rigorous system-level study if the BCG vaccination can boost immune response during COVID-19 infection and reduce its mortality risk.

Many countries all over the world have been using the BCG-vaccination to fight against tuberculosis (TB), organised through their national TB programs. It is obtained from Mycobacterium bovis isolation and currently it is the most widely used but amongst the most controversial vaccines. The BCG is an attenuated variant of a Mycobacterium bovis, which is firmly identified with Mycobacterium tuberculosis, the operator liable for TB. As one of the most widely used vaccines throughout the world, BCG has also been reported to reduce infant mortality due to infections other than TB [13,14]. BCG vaccine bolsters the inherent immunity system and thus protects from a wide range of other infections. For example, it is routinely used in the treatment of bladder cancer [15] and also reduced the respiratory syncytial virus infections [16]. Wardhana et al. has demonstrated its preventive impact on respiratory tract infections in elderly people [17], whereas a clinical trial evidenced protective effect against pneumonia in tuberculin-negative senior individuals [18]. Inspired by this evidence, it has been hypothesised that BCG vaccination might alleviate the severity and fatality of SARS-CoV-2 infection and thus provoke quick rescue [19,20]. Various studies are being under clinical trial to evaluate the effect of BCG vaccination on COVID-19 pandemic (for example, NCT04379336, NCT04537663, NCT04475302, NCT04327206 etc. on clinicaltrials.gov). All these evidences raise research need to investigate the influence of BCG vaccination on COVID-19 at the genetic level that has not been done yet.

Availability of high throughput technologies to analyse large-scale transcriptomic data have excelled these methodologies as promising tools in the biomedical research field [[21], [22], [23], [24], [25]]. Genetic inspection into the transcriptomic data yields better insight into the molecular pathogenesis of the SARS-CoV-2 infection and its related complications that includes idiopathic pulmonary fibrosis (IPF) [26], pulmonary arterial hypertension [27,28], common cancers [29], cardiovascular, hypertensive disorders [30] and psychiatric disorders [31]. This study aims to explore the genetic interaction of BCG vaccination on the COVID-19 through investigating the coexisting differentially expressed genes (DEGs), shared molecular pathways induced by those DEGs and their protein-protein interactome. The underlying analytical approach for this study is depicted in Fig. 1 . We used the shared DEGs to identify hub-genes, regulatory factors, potential drug targets and putative chemical agents. The findings could help to fight against the COVID-19 pandemic [32].

Fig. 1

Schematic diagram outlining the workflow of our proposed approach. (A) To conduct differential expression analysis, we have designed three individual experiments for each of the datasets. In those experiments, the case conditions were SARS-CoV-2 infection and BCG non-vaccination (2 datasets), and the control conditions were healthy status and BCG vaccination (2 datasets) respectively. (B) Common DEGs were then identified for both health conditions. (C) Biological functions of these DEGs were assessed and therapeutic targets were found by PPI analysis. (D) Functional enrichment analysis was performed with GO and cell signalling pathway databases. (E) Regulatory elements and possible comorbidities were determined. (F) Putative drug candidates and chemical agents were identified using curated databases. (G) All the gained results were validated through an extensive literature review.

Methods

Data

To identify the relationship between BCG vaccination and SARS-CoV-2 infection, we have analysed gene expression microarray and RNA-Seq transcriptomic data. In this study, we have collected two datasets from National Center for Biotechnology Information Gene Expression Omnibus (NCBI-GEO) with accession numbers GSE108363 and GSE90748, and one dataset from EBI array express with accession number E-MTAB-8871. GSE90748 is an RNA-Seq data and the other two are microarray transcriptomic data. The GSE90748 study has measured the expression profiles of 15 samples with 5 replicates of BCG injected and 10 replicates of non-injected lessions generated by high throughput sequencing technique, namely Illumina HiSeq 2000 by comparing [33]. GSE108363 is an expression profiling data by array using Illumina HumanHT-12 platform for 2 cohorts of blood samples infected with Mycobacterium tuberculosis and the same with Mycobacterium bovis BCG [34]. The E-MTAB-8871 study was designed to study the comparative gene expression profiles derived from human peripheral blood mononuclear cells (PBMCs) of 10 healthy control and 23 SARS-CoV-2 infected patients using the NanoString Human Immunology Panel [35].

Deferential expression analysis

RNA-Seq, the next-generation sequencing technology measures the gene expression with a high level of accuracy and mitigates many limitations of microarrays. Using this high-throughput sequencing technology and global trasncriptiomic analyses, we have designed three individual experiments for each of the datasets. In those experiments, the case conditions were SARS-CoV-2 infection and BCG non-vaccination (2 datasets), and the control conditions were healthy status and BCG vaccination (2 datasets) respectively [Fig. 1 A]. To identify DEGs associated with the respective conditions, we have used an R Bioconductor package DESeq2 [36], which identifies DEGs based on the Negative Binomial (also known as Gamma-Poisson distribution. Moreover, we have applied the R Bioconductor package, namely limma [[37], [38], [39], [40]] for the microarray dataset analysis to obtain the dys-regulated genes. To negate the errors introduced in the preparation and analysis of microarray data due to diverse operational set-ups and experimental system, the transcriptomic data were transformed using z-score normalisation defined aswhere is the expression data for gene and sample (for both case and control), and are respectively the mean and standard deviation of the expression levels for gene considering all the samples.

After obtaining the DEGs for each disease condition, we have selected the significant genes by setting the threshold level of the absolute value of log fold change 1 and FDR-adjusted (false discovery rate) -value .

Protein-protein interaction analysis

Proteins exhibit physical contact with each other in a cell or in a living organism indicating some biochemical events, typically function as some molecular processes within a cell, and thereby form a protein-protein interaction (PPI) network [41]. Here, we have used STRING to construct the PPI network for the DEGs that are shared between SARS-CoV-2 infection and BCG vaccination. STRING provides a knowledge base about known and estimated PPIs that comprises both physical and functional interactions, where nodes represent genes and edges indicates interconnection between them. At present, this database includes 24,584,628 proteins from 5090 organisms [42]. The medium confidence score 0.40 was set to generate this network. Proteins with different network characteristics such as having high-degree of interactions, may have a significant role in the cellular responses to a special physiological stimulus. We identified such highly interconnected nodes of the network, known as hub genes, using cytoHubba plugin of Cytoscape software [43] with the Degree topological algorithm [44]. These hub genes produce a highly dense module inside the interactome that could be of importance in effective drug discovery. We have extracted such highly concentrated modules by analysing the PPI network by another Cytoscape plugin, namely Molecular Complex Detection (MCODE) [45].

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) for a set of genes identifies their significant involvement in a certain molecular pathway or functional category to yield knowledge about the biological corollary, position on the chromosome, or regulation they share [46]. Such functional categories are defined by gene ontology (GO) terms that are further categorised as a biological process (BP), cellular component (CC) and molecular functions (MF) [47]. Similarly, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database provides functional knowledge regarding the cellular processes for analysing signalling pathway [48]. To have a better understanding of the metabolic pathways that are active in SARS-CoV-2, we have performed the GSEA for the shared DEGs, using a web-based graphical tool, namely ShinyGO V0.61 considering GO (BP, CC, and MF), and KEGG databases. ShinyGO provides a comprehensive analysis of a given set of genes for graphical representation of related molecular pathways and functional categories by incorporating GO and other data sources [49]. For statistical significance, the cut-off limit of adjusted -value 0.05 was set for the assessment of the enrichment results.

Regulatory analysis

Transcription factors (TFs) and micro-RNAs (miRNAs) usually regulate the expression pattern of a target protein at their transcription and post-transcriptional level, and thus have an impact on the biological processes [50]. We performed the gene regulatory networks (GRNs) analysis to obtain the regulatory factors that might influence the consequences of COVID-19 for not being vaccinated with BCG. For this, we analysed the common DEGs using NetworkAnalyst 3.0 web platform to obtain the TF-gene and gene-miRNA interactions. NetworkAnalyst 3.0 provides a free online platform to facilitate expression profiling, interactome analysis, and meta-analysis using transcriptomic data [51]. We have identified the TF-gene interactions using JASPAR database that offers open access to annotated and high-quality matrix‐based profiles of TF binding site [52]. For gene-miRNA interactome analysis, we considered the miRTarBase database, since it maintains a collection of manually curated and empirically validated miRNA targets [53].

Candidate drug and chemical agent identification

We performed the protein-drug interactions (PDI) and protein-chemical interactions (PCI) analysis using the overlapping DEGs for the potential drug target and chemical agents identification. For PDI, we incorporated Enrichr platform to explore the disease signature database DSigBD (http://dsigdb.tanlab.org/DSigDBv1.0/). EnrichR integrates a range of pre-compiled geneset libraries to facilitate enrichment analyses for a gene list of interest [54]. We have assessed the significance of the gained enrichment results by considering the adjusted p-value 0.05 for statistical significance. Again, the DSigDB is a collection of 22,527 gene sets related to the drug and small molecules considering the dysregulation in gene expression due to drug/compounds [55]. We have carried out PCI analysis using the NetworkAnalyst framework to exploit the Comparative Toxicogenomics Database (CTD) databases, which illuminates the effect of chemicals on diseases by providing manually curated information regarding protein-chemical and chemical-disease association [56].

Disease comorbidity assessment

To gain further insights into what implications COVID-19 may have on the overall health conditions, especially on those who have not been BCG vaccinated, we carried out the gene-disease association analysis for the shared 62 DEGs using the DisGeNET dataset through Enrichr [57]. DisGeNET is a publicly available database of gene-disease associations that comprises 21,671 genes with 30,170 human diseases. We obtained enrichment results by considering gene enrichment and adjusted -value . The obtained gene-disease associations (GDA) are then represented graphically as a bipartite network constructed with Cytoscape v3.8 software.

Results

Differential expression analyses found common DEGs between SARS-CoV-2 infection and BCG vaccination

PBMC samples from COVID-19 patients and healthy controls were compared to obtain the differential expression results of SARS-CoV-2 infection, which yielded a total of 1289 significant (892 up- and 397 down-regulated) DEGs. Similarly, blood samples from subjects without having BCG vaccination (case) were compared with the individuals having BCG vaccination (control). Thus, the resulted DEGs indicate genes to showcase altered expression if a person has not been vaccinated with m. bovis BCG. Cross comparison resulted in 11 down-regulated and 51 up-regulated DEGs that were common between SARS-CoV-2 and BCG vaccination datasets. The quantities of obtained DEGs overlapping between SARS-CoV-2 and BCG vaccination datasets are depicted in the Venn diagrams shown in Fig. 2 D and Fig. 2 E for up- and down-regulated DEGs, respectively. The volcano plots in Fig. 2(A, B and C) presented the expression pattern of the genes in three datasets. Note, all the downstream analyses were carried out considering these 62 common DEGs shown in Fig. 2 F.

Fig. 2

Differential gene expression and common DEGs. Volcano plots depict the genes expression in A) SARS-CoV-2 infected PBMCs, and two datasets for BCG vaccinations B) GSE90748 and C) GSE108363. Venn Diagram for finding common D) up- and E) down-regulated DEGs among three dataset, COVID-19, GSE90748 and GSE108363. F) The bubble plot shows the common DEGs between BCG vaccination and SARS-CoV-2 in PBMCs.

PPI network analysis for hub gene and module analysis

Next, we have queried the common DEGs in STRING for their PPI network, which is comprised of 54 nodes, each representing proteins and 166 edges indicating interactions among them (Fig. 3 A). Next, a set of highly connected modules (i.e. PPI sub-networks) were identified by MCODE, a Cytoscape plugin as shown in Fig. 3 B. The most densely connected module contains 10 genes including ADORA3, CCL4, CCR2, CXCL1, CXCL5, CXCL8, IFNG, OPRL1, PTGS2 and TLR5. Moreover, a set of hub genes were determined using cytoHubba, another Cytoscape plugin resulted in 10 hub genes including ITGB2, CXCL8, CXCL1, CCR2, IFNG, CCL4, PTGS2, ADORA3, TLR5 and CD33, which is shown in Fig. 3C. Next, this set of important genes (i.e. hub genes from cytoHubba and genes constituting the top module) was investigated in the literature, which revealed their pathogenic mechanism and associated disorders as tabulated in Table 1 .

Fig. 3

A) The protein-protein interaction network for the common DEGs between COVID-19 and BCG vaccination. B) The blue colored nodes indicate the top module of the network. C) The nodes having color from red, orange and yellow are the top significant hub genes.

Table 1

Particulars for the hub genes and the genes in the top module of PPI network.

Gene symbol	Name	Pattern	Pathogenetic mechanism	Associated Disorders	Ref.
ADORA3	Adenosine A3 receptor	Up	ADORA3 is highly regulated, most plentiful in the brain and several endocrine cells. G proteins mediate this receptor to inhibit adenylyl cyclase.	Ischemia and Ataxia, Sensory, 1, Autosomal Dominant.	[58,59]
CCL4	C–C Motif Chemokine Ligand 4	Down	CCL4 encodes mitogen-inducible monokine protein. It is one of the primary factors that the CD8+ T-cells produce and are suppressed in HIV. The protein expresses inflammatory and chemokine related processes.	Bacterial meningitis and Human Immunodeficiency Virus Infectious Disease.	[60]
CCR2	C–C Motif Chemokine Receptor 2	Up	CCR2 is a chemokine that mediates monocyte chemotaxis. This is responsible for infiltrating monocyte in inflammatory disorders such as rheumatoid arthritis and in the inflammatory reaction related to tumours.	Human Immunodeficiency Virus Type 1 and idiopathic Anterior Uveitis.	[61]
CD33	CD33 Molecule	Up	CD33 belongs to the sialic-acid-binding immunoglobulin-like lectin (Siglec) family that mediates cell-cell interactions and maintains rest for the immune cells	Alzheimer's Disease, Acute Leukemia and Acute Promyelocytic Leukemia.	[62]
CXCL1	C-X-C Motif Chemokine Ligand 1	Down	CXCL1 encodes CXC receptor 2, which is involved in inflammation and chemoattraction for neutrophils. Irregular expression of this protein plays a role to grow and develop certain tumours.	Alzheimer's Disease and Bacterial Meningitis	[63]
CXCL5	C-X-C Motif Chemokine Ligand 5	Down	The protein encoded by CXCL5 is a member CXC subfamily of chemokines that recruit leukocytes. It also participates to activate neutrophils.	pulmonary sarcoidosis, rheumatoid arthritis	[64]
CXCL8	C-X-C Motif Chemokine Ligand 8	Down	CXCL8 acts as a chemotactic element that activates neutrophils. It acts as basophils, and T-cells attractant, but not for monocytes. various cells release it as inflammatory responses.	Melanoma, bronchiolitis	[65]
IFNG	Interferon Gamma	Down	IFNG encodes cytokine that both the adaptive and natural immune system cells secret. Mutations in this gene are lined with an increase in vulnerability to the infections of viruses, bacteria and parasites as well as many autoimmune diseases.	Hepatitis C Virus, Tuberous Sclerosis 2.	[66]
ITGB2	Integrin Subunit Beta 2	Down	ITGB2 encoded proteins activate the immune response and leukocyte adhesion deficiency is resulted due to its defect. It also participates in the transmigration of leukocytes that includes T-cells and neutrophils.	leukocyte adhesion deficiency type i	[67]
OPRL1	Opioid Related Nociceptin Receptor 1	Up	OPRL1 encodes G-protein-coupled receptors belonging to the opioid family including kappa, delta and mu receptors. This receptor-ligand system regulates various biological processes and neuro-functioning, that include response to stress and anxious activities, memory and learning, locomotor action, and immune and inflammatory responses.	Drug dependence	[68]
PTGS2	Prostaglandin-Endoperoxide Synthase 2	Down	PTGS2 encodes isozymes that are inducible. Various stimulatory actions modulate this indicating its involvement in the prostanoid biosynthesis associated with mitogenesis and inflammation.	gastric ulcer, familial adenomatous polyposis	[69]
TLR5	Toll Like Receptor 5	Up	TLR5 identifies individual pathogen-related molecular models that are expressed in infections. It encodes proteins that can recognise bacterial flagellin which is a virulence component and the prime factor of bacterial flagella.	melioidosis, legionnaire disease	[70]

GSEA analyses identified pathways shared by SARS-CoV-2 and M. Bovis BCG

The GSEA was performed for the common DEGs considering BP, CC and MF for GO annotation as well as KEGG pathways to have better insight into their biological functions. The top 30 most significant GO terms and KEGG pathways based on FDR adjusted -values are shown in Fig. 4 . As depicted in the hierarchical clustering trees of enriched pathways/terms, inflammatory response and myeloid leukocyte activation were the top enriched GO terms of the biological process category. The cellular component part showed equally high enhancement of cytoplasmic vesicle and intracellular vesicle terms. Again, the molecular function subsection identified signalling receptor binding as the most significant whereas cytokine and chemokine related GO terms were found to be involved with the DEGs, as expected. On the other hand, most DEGs were found to be associated with lysosome and rheumatoid arthritis molecular pathways along with cytokine-cytokine receptor interaction, chemokine signalling pathway, IL-17 signalling pathway and TNF signalling pathway that could be of great biological interest.

Fig. 4

The hierarchical clustering representation of the 30 most significant a) GO-Biological process, b) GO-Cellular component, c) GO-Molecular functions, and d) KEGG pathways based on FDR adjusted p-value. Clustering is performed using the number of genes on the pathway and bigger dots indicate more significant P-values.

Transcription factors and miRNAs

From the gene regulatory network analysis, we obtained 56 transcription factors (TFs) as regulating the expression level of the common DEGs that could be involved in SARS-CoV-2 infection. As shown in Fig. 5 , the top 2 dominating TFs were GATA2 and FOXC1, regulating 37 and 34 DEGs, respectively. Whereas, CD36 was found to exhibit the highest regulation among all through 25 TFs. Moreover, we found 137 miRNAs that interacted with 40 DEGs among all as depicted in Fig. 6 . Among these identified miRNAs, hsa-mir-335–5p and hsa-mir-26b-5p shared the maximum 12 and 10 interactions, respectively, indicating their most influential role in gene regulation. On the other hand, ITGB8 and LIPA were the top 2 interacted DEGs having 38 and 35 interactions, respectively.

Fig. 5

TF-gene interactome where lime circles are the TFs while pink hexagons indicate the shared DEGs between SARS-CoV-2 infection and BCG vaccination.

Fig. 6

Gene-miRNA interaction network where lime hexagons represent shared DEGs and pink circles indicate miRNAs.

Potential drug target and chemical agents

After filtration of the obtained PDI data, we found 265 drug molecules to be significantly related to the common DEGs. The top 20 significant drug components are tabulated in Table 2 . The PCI analysis estimated 56 chemical agents to be associated with 44 DEGs, where both CXCL8 and PTGS2 interacted with maximum of 48 chemical compounds. The resulted network is depicted in Fig. 7 .

Table 2

The top 10 significant drug candidates obtained for the shared DEGs by SARS-CoV-2 and BCG vaccination.

Drug/Small molecule	Adj.p-val	Genes
Sodium dichromate	2.64E-13	CXCL8, RNASE6, RNASE1, RCBTB2, PYCARD, IFNG, NPC2, ADORA3, ALDH3B1, STAB1, CAT, CCL4, CD36, VSIG4, KCTD12
NICKEL SULFATE	3.49E-07	CXCL8, RNASE6, CXCL1, PTGS2, CXCL5, MS4A6A, ADORA2A, IFNG, GPNMB, CAT, CCL4, TIMP2, CD36, CCR2, C1QC
Lycorine	8.35E-07	CEBPA, GRN, SYK, IDH1, ITGB2, CYBRD1, CD1D, LIPA, RCBTB2, PYCARD, SCPEP1, MS4A6A, GPNMB, NPC2, ALDH3B1, KCTD12, CCR2
Medroxyprogesterone acetate	2.24E-06	MS4A6A, GPR34, IFNG, GPNMB, ADORA3, IDH1, CAT, VSIG4, C1ORF162, PTGS2, RNASE1, C1QC
Phorbol 12-myristate 13-acetate	2.00E-05	MSR1, GRN, CXCL8, IFNG, ITGB2, CAT, CCL4, TIMP2, ITGB8, CD36, PTGS2, CXCL5
1-chloro-2,4-dinitrobenzene	3.73E-05	CXCL8, GPNMB, IDH1, CAT, CCL4, CXCL1, CD36, PTGS2, CCR2
Anisomycin	3.85E-05	CEBPA, GRN, SYK, IDH1, ITGB2, CYBRD1, CD1D, LIPA, PYCARD, MS4A6A, ALDH3B1, PPT1, VSIG4, KCTD12, CD33, CCR2
Aspirin	6.39E-05	CXCL8, OAS1, ADORA2A, IFNG, ADORA3, ITGB2, TIMP2, RNASE6, CXCL1, CD36, PTGS2, RNASE1
RUTIN	7.47E-05	CXCL8, IFNG, GAA, CAT, PTGS2
Acetovanillone	9.67E-05	MSR1, CXCL8, CXCL1, CD36, PTGS2

Fig. 7

Protein-chemical interaction network for common DEGs, where olive octagons are the DEGs and pink circles indicate chemical agents. The network was constructed using degree and betweenness .

Comorbid diseases of COVID-19

After manual curation of the GDA, we obtained 30 highly significant diseases related to 53 distinct DEGs with up to 21 shared DEGs. The resulted GDA network is depicted in Fig. 8 . Diabetes mellitus was found to be associated with the highest number (21) of DEGs where both obesity and Alzheimer's disease were related to 20 DEGs. On the contrary, both IFNG and PTGS2 were linked with all the 30 diseases. We hypothesize that these DEGs and diseases association could yield further research to investigate the COVID-19 progression and its association with BCG vaccination.

Fig. 8

Gene-disease association network. In this bipartite network, circular nodes (blue) represent the shared DEGs while octagonal nodes indicate COVID-19 (red) and different diseases (lime).

Discussion

This study primarily evaluates the potentiality of the BCG vaccination to minimise the severity and/or mortality of COVID-19 disease at the molecular level by adopting a series of bioinformatics strategies. For this, we have compared the gene expression profiles of PBMC in COVID-19 patients with healthy individuals as well as subjects without and with BCG vaccination. Cross-comparison identified 62 genes exhibiting similar alteration of expression patterns indicating their significant protective contribution on COVID-19 as a result of BCG vaccination. Subsequently, we explored their biological functionalities by employing PPI analysis, shared GO terms and KEGG pathway identification, gene regulatory network analysis followed by drug target and chemical agent identification to estimate the influence of BCG vaccination over COVID-19.

The PPI network analysis, being an integral part of this study, identified the most densely connected sub-network containing ADORA3, CCL4, CCR2, CXCL1, CXCL5, CXCL8, IFNG, OPRL1, PTGS2 and TLR5 via MCODE algorithm. Hub gene identification algorithm also obtained 10 hub genes (ITGB2, CXCL8, CXCL1, CCR2, IFNG, CCL4, PTGS2, ADORA3, TLR5 and CD33) being the most interacting genes with each other. Many of these genes have already been implicated with COVID-19. ADORA3 interacts with Adenosine to moderate the anti-inflammatory mechanism by reducing the production and release of pro-inflammatory cytokines [71]. The elevated level of cytokine circulation, known as cytokine storm, has been reported to be highly associated with the COVID-19 pathogenesis resulting in lung damage [72]. Yong et al. reported increased expression of pro-inflammatory cytokines including CXCL1 and CXCL8 as well as chemokine CCL4 and CCR2 in SARS-CoV-2 infected patients [73]. Higher expression of CD33 is reported to be associated with the increased severity of COVID-19 [74]. An enhanced ratio of IFNG has been found in severe SARS-CoV-2 infected patients, which could exacerbate the cytokine storm [75]. ITGB2 is found to be co-expressed between angiotensin-converting enzyme II (ACE2) and leukocyte mediated immunity, playing a crucial role in the immune responses [76]. It also accelerates the lung repair and restoration after injury by negatively regulating WNT signalling in the lung [77]. PTGS2 plays a vital role in SARS-CoV-2 infection by encoding cytochrome c oxidase subunit II (COX2) [78]. SARS-CoV-2 activates COX2 and excites COX2 inflammatory cascades to cause lung inflammation [79]. TLR5 has been suggested as a putative therapeutic target to fight against COVID-19 as it stimulates early signalling for innate immunity generation [80]. Thus, the close association of these identified genes with the COVID-19 pathogenesis hints at the prospective influence of BCG vaccination against the severity of SARS-CoV-2 infection.

Functional enrichment analysis identified several significant GO terms and KEGG pathways that COVID-19 shares with BCG vaccination. Among the top 30 enriched GO terms in each category notable terms include immune response, inflammatory response, neutrophil activation, cytokine activation and chemokine activation. On the other hand, prospective shared molecular pathways include lysosome, IL-17 signalling pathway, TNF-signalling pathway, chemokine signalling pathway and cytokine-cytokine receptor interaction. Lysosome serves as the animal cell's primary digestive chamber and the drugs targeting the lysosomes are considered to be the prospective therapy against COVID-19 [81]. Again, Interleukin 17 (IL-17) plays a vital role in recruiting immune cells to the infected site as redemption and also promoting the reduced flow of chemokines and cytokines [82]. Thus, IL-17 activation is supposed to be implicated with SARS-CoV and MERS-CoV infections [83]. This pathway also induces pro-inflammatory cytokines and thus corresponds to SARS-CoV-2 infection [84]. Consequently, IL-17 has been suggested as a plausible target in developing effective therapies to treat severe COVID-19 [85]. Meanwhile, Yabo, et al. reported TNF-signalling pathway to be enhanced in intense SARS-CoV-2 infection [86]. Overall, these shared DEGs between BCG vaccination and SARS-CoV-2 infection and their immune response-related pathway activities suggest that BCG vaccination may potentially contribute to induce or boost the host response against SARS-CoV-2 infection and hence may reduce COVID-19 mortality rate, which is also currently reported to be evident from experimental data [87].

We have also investigated the association between COVID-19 and BCG vaccination from the perspective of regulatory mechanisms, e.g., TF-gene, gene-miRNA, protein-drug and protein-chemical interactions. PDI analyses revealed several potential drug candidates that could be further investigated with chemical experiments at a larger scale for verification. Among these, lycorine is a potential candidate to develop medicine against SARS-CoV [88]. This phytochemical also exhibits strong inhibitory effects against SARS-CoV-2 infection [89]. According to a recent study, anisomycin can impede the inflammation in macrophages, which eventually may obstruct the cytokine storm [90]. Initially, various in vitro results suggested prospective therapeutic possibilities in several drugs and chemicals to inhibit COVID-19. For instance, remdesivir and chloroquine had been proposed to fight against SARS-CoV-2 infection [91]. Among them, remdesivir was the first drug that the US Food and Drug Administration (FDA) has approved for use clinically to treat hospitalized COVID-19 patients [92]. Besides this, favipiravir evidenced promising effects to protect from COVID-19 [93]. Furthermore, hydroxychloroquine and azithromycin combination therapy have reported significant protection against COVID-19 [94]. But, at later stages randomized clinical trials (RCTs) at a larger global scale have evidenced insignificant certainty for the efficacy of most of these components. For example, no definite evidence has been found to endorse the effectiveness of chloroquine or hydroxychloroquine with or without azithromycin in the COVID-19 treatment [95]. Therefore, the urge of continuous endeavour for finding putative potential therapeutic candidates has stretched even more. However, the identified chemical agents through the PCI yield potential considerations against COVID-19. For example, Valproic acid shows antiviral effects that demonstrate its promising prospect against SARS-CoV-2 [96]. Furthermore, tretinoin is predicted as a repurposable drug target for COVID-19 [97]. Besides this, cyclosporine impedes the replication of SARS-CoV and MERS-CoV [98], hence it is speculated that cyclosporine treatment could be beneficial in COVID-19 [99]. Altogether, these findings are of great biological and clinical interest that may enhance our insight into SARS-CoV-2 infection and eventually promote the therapeutic strategy development to counter the COVID-19 pandemic. Some limitations of this study would be noted as samples of SARS-CoV-2 infection and BCG vaccination were taken at different time frames, the number of samples is small, different technologies (microarray and RNA-seq) were used for data extraction, and the lack of clinical validation of the identified signatures. Thus cautions have to be taken while interpreting the findings of the study. In future studies, datasets from the same patient samples with BCG vaccination status and their SARS-CoV-2 infection status along with their categorical distributions of the factors like age, sex, and comorbidities could be investigated and external validation of the findings could be administered.

Conclusion

In this study, we aimed to reveal whether BCG vaccination could boost the immunity against COVID-19 by employing a series of bioinformatics approaches. We found several hub genes that have a protective effect against COVID-19 and its severity. Some signature genes are involved in reducing the production and release of pro-inflammatory cytokines, and hence, modulate the anti-inflammatory mechanism. Enrichment analysis indicated that the BCG vaccine has immunomodulatory activities which are necessary to reduce the fatality of COVID-19 patients. We hope, these interesting findings would potentially open up further research direction on this hypothesis with thorough pathological investigations.

Contributions

All authors contributed to the manuscript. U.N.C., M.B.I., S.A., W.S. and M.A.M. conceived and designed the study. U.N.C., M.O.F. and M.M. analysed the data and wrote the R programming code for the development of the pipeline. U.N.C., M.O.F., M.B.I. and A.K.M.A. conducted all other bioinformatics analyses and wrote the manuscript. All authors read and approved the final version of the manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.