{'fraction_non_alphanumeric': 0.032661276084884006, 'fraction_numerical': 0.027258924473005895, 'mean_word_length': 5.242450958783337, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bioinformatics and machine learning approach identifies potential drug targets and pathways in COVID-19

Abstract Current coronavirus disease-2019 (COVID-19) pandemic has caused massive loss of lives. Clinical trials of vaccines and drugs are currently being conducted around the world; however, till now no effective drug is available for COVID-19. Identification of key genes and perturbed pathways in COVID-19 may uncover potential drug targets and biomarkers. We aimed to identify key gene modules and hub targets involved in COVID-19. We have analyzed SARS-CoV-2 infected peripheral blood mononuclear cell (PBMC) transcriptomic data through gene coexpression analysis. We identified 1520 and 1733 differentially expressed genes (DEGs) from the GSE152418 and CRA002390 PBMC datasets, respectively (FDR < 0.05). We found four key gene modules and hub gene signature based on module membership (MMhub) statistics and protein–protein interaction (PPI) networks (PPIhub). Functional annotation by enrichment analysis of the genes of these modules demonstrated immune and inflammatory response biological processes enriched by the DEGs. The pathway analysis revealed the hub genes were enriched with the IL-17 signaling pathway, cytokine–cytokine receptor interaction pathways. Then, we demonstrated the classification performance of hub genes (PLK1, AURKB, AURKA, CDK1, CDC20, KIF11, CCNB1, KIF2C, DTL and CDC6) with accuracy >0.90 suggesting the biomarker potential of the hub genes. The regulatory network analysis showed transcription factors and microRNAs that target these hub genes. Finally, drug–gene interactions analysis suggests amsacrine, BRD-K68548958, naproxol, palbociclib and teniposide as the top-scored repurposed drugs. The identified biomarkers and pathways might be therapeutic targets to the COVID-19.


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a newly evolved virus, first identified in Wuhan, China in December 2019. Pneumonia caused by the SARS-CoV-2 was referred to as coronavirus disease , which was declared as the COVID-19 as a pandemic by the World Health Organization (WHO) [1]. The main symptoms of COVID-19 are fever, cough, pneumonia and shortness of breath [2]. The SARS-CoV-2 has infected almost 65 943 003 peoples with more than 1 519 137 deaths globally (as of 5 December 2020) [3].
The spread of the virus infection may be controlled through early detection of COVID-19 patients. However, the currently used methods including real-time polymerase chain reaction (RT-PCR) are subject to limited sensitivity and specificity as well as time-consuming. Moreover, careful sample collections and preparations, and skilled manpower are required, which are a tremendous drawback for developing countries. Thus, transcriptomic analysis of SARS-CoV-2 PBMC may provide candidate biomarkers. Past studies have conducted a transcriptomic analysis of various organs including lung epithelial cell [4][5][6], PBMC [7,8]. Most of the previous reports detected the hub genes in COVID-19 from each module either via the PPI network or module membership criterion [5,[9][10][11][12][13]. Extensive gene expression analysis to identify differentially expressed genes (DEGs) and associated gene ontologies have been proposed by previous reports [5,[9][10][11][12][13]. Detection of specific gene modules has not been performed in COVID-19, and identification of key gene module hubs and targeting those critical genes for drugs repurposing is crucial to combat COVID-19. However, despite important findings from those studies, integrative analysis is needed to detect novel dysregulated genes and pathways for the pathogenesis of COVID-19.
The weighted gene coexpression network analysis (WGCNA) identifies significant modules (clusters) of highly correlated genes [14]. It explains the correlation patterns within genes and samples and biologically interprets the function of gene modules. The hub genes of differentially coexpressed modules provide more significance than the usual hub DEGs. WGCNA is widely used for biomarkers identification in various diseases [15][16][17], it has a great prospect in COVID-19. Identification of repurposable drug candidates for COVID-19 may reverse the DEGs of COVID-19, we have decided to utilize the anti-signature approach [5]. Identification of this in silico based approach offer opportunities to identify potential candidate drugs that might be considered for drug repositioning for COVID-19 treatment.
In this study, we implemented a system biology approach to key gene modules (identified via WGCNA) that were DEGs in COVID-19 PBMC. The hub genes were then identified from the key gene modules based on gene module membership (MMhub) and protein-protein interaction network (PPIhub), respectively. Then, we employed machine learning methods to determine the validity of these hub genes. Finally, we identified several candidate drugs considering these hub genes as therapeutic targets.
Our results may provide novel insights into the pathogenesis of COVID-19 and the potential molecular targets.

RNA-sequencing datasets
In this study, we used two RNA-Sequencing PBMC datasets of SARS-CoV-2 (COVID-19). One of the COVID-19 gene expression raw counts dataset was obtained with the accession number GSE152418 under the platform GPL24676 from the NCBI Gene Expression Omnibus (GEO) [18]. Recently this dataset was deposited by Arunachalam et al. [7] that contained 34 samples (17 COVID-19 samples and 17 healthy control samples). The other dataset was obtained from the Chinese Academy of Science with the accession number CRA002390 that contained PBMC samples from three COVID-19 infected patients and three healthy donors [6]. In this study, the GSE152418 discovery dataset was used to analyze WGCNA, and the CRA002390 dataset was used for independent validation.

Data preprocessing and differential gene expression analysis
The transcriptomics dataset, GSE152418, of COVID-19 contained a large number of genes (60 683 genes). For data preprocessing, the low expressed genes (sum of the gene counts for all samples <100) were excluded from this dataset. Then the differential gene expression analysis of the dataset was carried out through the DESeq2 package in the R language [19]. For the CRA002390 normalized count dataset, we used the limma package in R [20] to identify the DEGs. We considered adjusted P-value with Benjamini-Hochberg FDR correction techniques (FDR < 0.05) and log 2 (FC) ≥ 1 statistical threshold parameters for DEGs identification.

Weighted gene coexpression networks construction
The gene coexpression network construction was executed after removing the outlier samples (if there exist). The cluster dendrogram of the samples was constructed to check the outliers through the hclust function in R. We used the WGCNA package in R [25] to construct the weighted gene coexpression network. For finding numerous soft-thresholding powers β over R 2 , we used the pickSoftThreshold function. Then we picked the value of β for which the value of R 2 maximum. The adjacency matrix and Topological Overlap Matrix (TOM) were then constructed using this soft threshold power with the transformed gene expression matrix. Then the dissimilarity of TOM (dissTOM) was computed to construct a network heatmap plot and for further analyses.
For the detection of the module, the dendrogram of genes was constructed with a dissTOM matrix using hclust function with different colors. The Automerged technique was used to get modules using the parameters: deepSplit = 2 and minClus-terSize = 30 for avoiding the generation of small or large modules. Here MEDissThres = 0.25 was taken for merging similar modules [15].

Preservation analysis for the key module selection
To find the key modules, we used module preservation analysis. The modulePreservation function [21] was used to evaluate each of the modules whether it was robust and reproducible across datasets or not. If preservation statistics Z summary>10, then the module is considered as preserved [15]. It is apparent that the module preservation and preservation statistics-median Rank are negatively correlated and there is present a positive correlation between Z summary statistic and module preservation.

Gene ontology and pathway enrichment analysis
The high connectivity of genes inside the coexpression modules may represent crucial information about the similar biological roles within the same module. The functional enrichment analysis of the genes was studied in each selected key modules via Gene Ontology (GO) and pathway analysis [22,23]. The GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was executed via the Database for Annotation, Visualization, and Integrated Discovery and visualized (DAVID) tools [24] and the result was visualized through cluster-Profiler package in R [25]. A statistical threshold criterion with an adjusted P-value <0.05 was used to select significant GO terms and KEGG pathways.

Identification and validation of hub genes
In a gene module, a series of genes with the greatest degree of connectivity was detected as hub genes that investigate the characteristics of a module. Also, a module's connectivity defining the hub genes was calculated using the absolute value of the Pearson's correlation (|cor.geneModuleMembership| > 0.8) [15]. Furthermore, we uploaded all genes of the key modules into the STRING database, choosing confidence score cutoff >900 to construct protein-protein interaction (PPI). In the PPI network, genes with a connectivity degree of ≥8 were also defined as hub genes [26]. The present study further analyzed the CRA002390 validation dataset to confirm the role of these hub genes as molecular signature genes for COVID-19. The common genes among the MMhub, PPIhub and the DEGs of CRA002390 will be considered for further network analyses.

Performance evaluation of the hub genes with classification algorithms
To check the validity of the identified common hub genes in two different analyses results (MMhub and PPIhub), the five popular classification algorithms, support vector machine (SVM) [27] radial basis kernel function, random forest (RF) [28], Poisson linear discriminant analysis (PLDA) [29], negative binomial linear discriminant analysis (NBLDA) [30] and voom-based diagonal linear discriminant analysis (voomDLDA) [31] were conducted through MLSeq package in R [32]. For the SVM and RF classifiers, we used DESeq normalization and VST transformation on the count dataset. Here, we considered the COVID-19 GSE152418 transcriptomic dataset for the classification analysis separately with MMhub genes and PPIhub genes, respectively. We computed four performance measures namely, accuracy, area under the ROC curve (AUC), sensitivity and specificity based on the data with MMhub and PPIhub genes, respectively.

PPI network analysis
The common genes among MMhub, PPIhub and DEGs of CRA002390 were considered for further network analysis. The PPI network for these genes was constructed via the STRING web tools [33]. STRING provides the PPI network that shows how the identified hub genes (proteins) interrelate functionally and physically with each other through encoding the gene list as input. Through the STRING information, the PPI network was constructed via Cytoscape [34] which is an open-source platform. The hub genes from this network were chosen based on degree connectivity via Cytohubba in Cytoscape. These hub genes were considered for final biomarkers of COVID-19.

Transcription factor and miRNAs identification
The significant transcription factors (TFs) were identified through a freely accessible database of TFs repository-JASPAR [35] by executing the interaction of the TFs-target genes via Net-workAnalyst [36]. The significant miRNAs were identified from miRNAs-target gene interaction analysis through the Tarbase [37] and mirTarbase [38] database via NetworkAnalyst [36]. These networks were visualized with Cytoscape and the significant hub TFs and miRNAs were selected via the CytoHubba plugin in Cytoscape based on the degree connectivity.

Drug-gene interaction analysis
To predict potential drugs for the treatment of COVID-19, we performed a transcriptomic anti-signature approach [5,39] using the L1000FWD web-based tool [40], which measures the similarity score between input DEGs and expression signature and LINCS-L1000 data to detect drugs that may reverse the input gene signature. LINCS-L1000 contained drugs induced gene signature of about 50 human cell lines in response to 20 000 compounds. The significant drugs were chosen based on the criterion, q-value < 0.05.

Proposed bioinformatics pipeline
The workflow of the proposed bioinformatics methods to identify significant pathway and drug targets were illustrated in Figure 1 with followings: 1. RNA-sequencing datasets: Differential gene expression analysis was achieved using RNA-seq data. One dataset was used to analyze the weighted gene coexpression network and the other independent dataset was used for validation. 2. Coexpression network reconstruction: The weighted gene coexpression network was constructed through the WGCNA package in R. 3. Key module selection: The Module Preservation function was used to identify robust and reproducible modules. 4. Enrichment analysis: Biological insights of key module genes were determined via the DAVID tools. 5. Hub genes identification: Hub gene signature was identified based on module membership (MMhub) statistics which calculated using the absolute value of the Pearson's correlation and PPI networks (PPIhub) by using degree metrics. The common genes as hub signature were selected among the MMhub, PPIhub and the DEGs of CRA002390.. 6. Validation analysis with machine learning methods: For validation of hub signatures, the five popular classification algorithms, SVM-radial basis kernel function, RF, PLDA, NBLDA and voomDLDA were conducted through MLSeq package in R. 7. Performance evaluation of hub genes: The performance measures including accuracy, AUC, sensitivity and specificity of MMhub and PPIhub genes were calculated by using an independent dataset. 8. Transcriptional regulator identification: Potential regulators (i.e.: TF and miRNA) of hub gene signature were determined using JASPAR, Tarbase and mirTarbase databases via NetworkAnalyst. 9. Candidate drug identification: LINCS-L1000 data were used to detect significant putative drugs that may reverse the input hub gene signature.

DEGs identification
The discovery dataset (GSE152418) contained 60 683 genes with 17 COVID-19 and 17 health control samples. After excluding the lowly expressed genes, we selected 20 251 genes for DEGs identification. We identified a total of 1520 DEGs (1299 upregulated and 221 downregulated) from this dataset (FDR < 0.05) that were considered for further weighted gene coexpression network analysis. We identified 1733 DEGs (1139 upregulated and 594 downregulated) from the validation dataset (CRA002390) in the COVID-19 samples. Table 1 presented the information regarding the datasets used in this study.

Weighted gene coexpression networks and module preservation analysis
The gene coexpression network analyses were performed with identified 1520 DEGs of 16 COVID-19 samples in the GSE152418 dataset. The cluster dendrogram of these samples was visualized in Figure 2A that revealed, no outlying samples presented in this dataset. We selected the optimized soft-thresholding power, β = 6 with R 2 = 0.90 as the scale-free topology criteria ( Figure 2B). The coexpression networks were then constructed with this soft threshold power β = 6 and obtained 10 coexpressed modules namely, black, blue, brown, green, gray, magenta, pink, red, turquoise and yellow. The cluster dendrogram of these 10 modules presented in Figure 2D. We observed 63, 273, 259, 95, 46, 58, 77, 426 and 198 genes for black, blue, brown, green, magenta, pink, red, turquoise and yellow modules, respectively. The gray module contained 25 genes that were tied up as noncoexpressed.   Figure 2C showed the network heatmap of all genes with these nine modules. The interactions among these coexpressed modules were presented with the module eigengene dendrogram and eigengene network heatmap ( Figure 2E).
In the module preservation analysis, we used the GSE152418 dataset with 17 health control samples as the test dataset. We identified turquoise, blue, brown and pink modules as the most stable through preservation analysis. The remaining modules were considered nonstable since their Z summary statistic <10 ( Figure 2G). Figure 2F presented that the turquoise, blue, brown and pink modules were the best-preserved among all modules since their medianRank statistic were minimum than other modules.

Biological insights of the four-module genes
To obtain further biological insight into the genes of the selected four modules, the GO and KEGG pathway analysis was conducted in this study. The significant biological process (BP) mainly enriched in the immune response, division and fission related BP ( Figure 3A). The significant molecular function (MF) mainly enriched in the binding related functions ( Figure 3B). The most significant cellular components (CC) for the four modules are enriched in several cell compartments ( Figure 3C). The KEGG pathways for the genes of four modules significantly enriched in several pathways such as infection-related pathways (i.e: herpes simplex virus 1 infection, human papillomavirus infection), autoimmune diseases related pathways (i.e.: systemic lupus erythematosus, rheumatoid arthritis, type I diabetes mellitus), ECM-receptor interaction, IL-17 signaling pathway and p53 signaling pathway ( Figure 3D and Table 2). Interestingly, alcoholism and systemic lupus erythematosus were significantly enriched for the genes of these four modules.

Identification and validation of hub genes
We identified 422 hub genes (MMhub) in turquoise, blue, brown and pink modules with high connectivity using the module connectivity threshold criterion |cor.geneModuleMembership | > 0.8. Additionally, we identified 240 hub genes (PPIhub) for the four modules from the PPI through the STRING database with a connectivity degree ≥8. The present study further analyzed the CRA002390 validation set to confirm the role of these hub genes as candidate biomarker genes for COVID-19. We identified 52 common genes among MMhub, PPIhub and the DEGs of CRA002390 ( Figure 5A). The summary of these 52 genes was described in Table 4. The expression values of these common genes over the COVID-19 and health control samples were presented in Figure 5C. The heatmap showed the two main clusters of these 52 hub genes in Figure 5B in terms of COVID-19 infected samples and control samples.

Performance evaluation of the hub genes with a classification algorithm
To investigate the validity of the identified hub DEGs in two different analyses (i.e.: MMhub and PPIhub), the five popular classification algorithms were executed in this study. The   performance measures have been computed based on the datasets with MMhub and PPIhub genes, respectively. We executed these calculations 20 times using 5-fold crossvalidation and the average performance measurement values were computed and summarized in Table 3. The boxplot of the five machine learning classifiers based on test accuracies and AUC were presented in Figure 4. In Table 3 and Figure 4

PPI network analysis with identified common genes
The PPI networks for the 52 common genes were constructed via the STRING in Cytoscape. Figure 6A presented the network interaction among these genes and identified 10 hub genes (PLK1, AURKB, AURKA, CDK1, CDC20, KIF11, CCNB1, KIF2C, DTL and CDC6) based on a higher degree of connectivity.

Drug repositioning based on drug-gene overrepresentation analysis
We have identified 50 candidate drugs by reversal gene signature-based approach with q-value<0.05. Among them, the top 20 drugs have been summarized in Table 5. According to our analysis, 'amsacrine', 'BRD-K68548958', 'naproxol', 'palbociclib' and 'teniposide' were the top significantly identified drugs among others. These candidate drugs or components may be used for therapeutic applications in COVID-19.

Discussion
The COVID-19 is affecting severely millions of people and taking thousands of precious lives every day over the globe due to its pandemic behavior. Though there are several candidate drugs and vaccines that were studied and proposed to treat the disease, no solid cure is available yet. The current study employed a gene coexpression network analysis to decode the critical genes and pathways of COVID-19. We identified 1520 and 1733 DEGs for the GSE152418 and CRA002390 RNA-sequencing PBMC based datasets, respectively. The four key modules were identified in GSE152418 COVID-19 data via WGCNA and module preservation analysis. The GO and pathway enrichment analysis were conducted for these key modules. We checked the validity of these identified hub genes (MMhub and PPIhub) with machine learning classifiers. We found 52 common genes in MMhub, PPIhub and the DEGs of CRA002390. The PPI networks, transcriptional regulatory networks of the common hub genes were constructed. We found 10 hub genes from the PPI of the identified common genes, and those were considered as the final candidate molecular blood signatures of COVID-19. These findings may provide new insights into the COVID-19 pathogenesis. Our employed approach is significantly different from previous bioinformatics reports in COVID-19 research [5,9,12,13,41,42] which relied on the identification of genes by differential analysis. However, in order to provide systems biology insights, we have implemented methodologies particularly gene coexpression module analysis that provides key gene modules rather than finding DEGs. Firstly it clustered genes in a specific condition followed by detection of hub genes which are highly connected nodes in each module network. There is a possibility of detecting some inappropriate genes as the hub if we use only one approach. To detect the hub genes from each module more accurately, we performed joint-analysis which revealed common hub genes by module membership significance analysis and PPI analysis, which is a robust approach in selecting hub genes in COVID-19. Our analysis focused on PBMC gene expression analysis, to obtain further insights regarding the potential utilization of the identified hubs in diagnostic development for COVID-19, we decided to perform classification by widely used machine learning classifiers. The performance of state-of-theart methods was evaluated the classification performances in which obtained a higher score in accuracy, AUC, specificity and sensitivity.
To elucidate the roles of the identified DEGs, the GO and KEGG pathway analysis were executed in this study. Among the identified GO terms, the immune response, response to cytokines, cytokine-mediated signaling and response to external stimuli play crucial roles in restricting viral infections. Among the identified pathways of four modules the cytokine-cytokine receptor interaction, the IL-17 signaling pathway was highly enriched in COVID-19. Viral infectious diseases such as herpes simplex virus 1 infection, human papillomavirus infection and viral carcinogenesis pathways were detected. Human papillomavirus infection pathways are significantly related to colorectal cancer. An autoimmune disease like rheumatoid arthritis and type 1 diabetes mellitus were also expressively enriched. These pathways show massive significance to drug repurposing chances in COVID-19.
The application of machine learning classifiers have been widely used in different bioinformatics tasks [43][44][45]. We executed the machine learning classification algorithms based on MMhub and PPIhub genes data respectively to check their validity. We observed the satisfactory sample classification performance (accuracy >0.90; except PLDA) between COVID-19 and health control samples for both datasets (MMhub, PPIhub). The results indicate the validation of our identified hub genes through module membership and PPI networks.
We identified PLK1, AURKB, AURKA, CDK1, CDC20, KIF11, CCNB1, KIF2C, DTL and CDC6 hub genes representing that, they had a high association with clinical trait along with vital BPs and some of them were detected in COVID-19. Among them, the Polo-like kinase 1 (PLK1) gene was detected as down-regulates in Parainfluenza Virus 5 [46]. The Aurora Kinase B (AURKB) and Aurora Kinase A (AURKA) were found SARS-CoV-2 as DEGs in Caco-2 cells [47]. The Cyclin-Dependent Kinase 1 (CDK1) genes interact with the thrombocytopenia syndrome virus which initiates the cells into the M phase [48]. The Cell Division Cycle 20 (CDC20) and CDK1 were also found as potential biomarkers for hepatocellular carcinoma [49]. The identified hub TFs and miRNAs are also significantly associated with viral infectious diseases.
Finally, we detected the candidate drugs using the reversal gene signature-based approach [5,9]. Among them, Naproxen is a nonsteroidal anti-inflammatory drug that was studied to use for the treatment of critically COVID-19 infected patients and to limit the spread of the virus [50,51]. The drug Teniposide was suggested to evaluate the treatment of SARS-CoV-2 infected patients [52]. We proposed to send these candidate drugs for biological and clinical experimentation for the possible use in COVID-19 treatment.  Our analysis focused on PBMC gene expression analysis, to obtain further insights regarding the potential utilization of the identified hubs in diagnostic development for COVID-19, we decided to perform classification by widely used machine learning classifiers. Our analysis showed a good level of specificity in classification performances.
However, several limitations of the study should be noted as findings of this study relied on bioinformatics analysis without functional studies in wet-lab, thus caution should be taken in interpreting the results. Moreover, the transcriptomic analysis and candidate drugs were identified by reversal of PBMC gene expression in COVID-19 but the primary affected organ by SARS-CoV-2 is lung tissues, thus further research is now proposed to explore biological insights in COVID-19.

Conclusions
The present research aimed to identify key genes and molecular pathways altered in response to SARS-CoV-2 in blood cells compared to normal blood cells. We detected four key modules through module preservation analysis. The 52 common genes were identified from resultant 422 and 240 hub genes based on module membership statistics and PPI networks and from 1733 detected DEGs of CRA002390. The 10 hub genes (PLK1, AURKB, AURKA, CDK1, CDC20, KIF11, CCNB1, KIF2C, DTL and CDC6) were identified from the PPI networks of these 52 genes. The TFs (FOXC1, GATA2, YY1, E2F1, NFIC, FOXL1 and SRF) were also found as potential regulators of the hub genes. The naproxol, teniposide, amsacrine, BRD-K68548958, palbociclib were identified as the top-scored repurposed drugs for COVID-19 pathogenesis. The identified drugs should be judged with wet-lab experiments before clinical studies. Our results may provide novel insights into the pathogenesis of COVID-19 and the potential molecular targets for novel interventional approaches.

Key Points
• COVID-19 pandemic has emerged as a massive threat to humankind limited by the unavailability of effective drugs.
• This study has performed a comprehensive clinical bioinformatics and systems biology analysis of all available SARS-CoV-2 infected peripheral blood mononuclear cell (PBMC) transcriptomic datasets to identify gene modules by gene coexpression analysis.
• A robust four key gene modules and hub gene signature were detected based on gene module membership statistics and protein-protein interaction networks and machine learning methods. • Functional annotation by enrichment analysis of the genes of these modules demonstrated immune and inflammatory response biological processes enriched by the gene signature.
• Several potential candidate drugs based on the reversal of transcriptomic signature were also detected that may be effective treatment candidate for COVID-19.