{'fraction_non_alphanumeric': 0.031724754244861486, 'fraction_numerical': 0.017828418230563, 'mean_word_length': 5.316963025684448, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Drug genetic associations with COVID-19 manifestations: a data mining and network biology approach

Available drugs have been used as an urgent attempt through clinical trials to minimize severe cases of hospitalizations with Coronavirus disease (COVID-19), however, there are limited data on common pharmacogenomics affecting concomitant medications response in patients with comorbidities. To identify the genomic determinants that influence COVID-19 susceptibility, we use a computational, statistical, and network biology approach to analyze relationships of ineffective concomitant medication with an adverse effect on patients. We statistically construct a pharmacogenetic/biomarker network with significant drug-gene interactions originating from gene-disease associations. Investigation of the predicted pharmacogenes encompassing the gene-disease-gene pharmacogenomics (PGx) network suggests that these genes could play a significant role in COVID-19 clinical manifestation due to their association with autoimmune, metabolic, neurological, cardiovascular, and degenerative disorders, some of which have been reported to be crucial comorbidities in a COVID-19 patient.


INTRODUCTION
Since the first outbreak of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) infection in Wuhan, Hubei, China in December 2019, more than 263 million people have been infected by SARS-CoV-2, leading to the COVID-19 pandemic with a total of 5.2 million deaths globally [1]. COVID-19 has been shown to be a multifactorial disease, with the scientific interest focused on the host genetic variants associated with critical illness aiming at identifying targets for efficient therapeutics development.
Drug repurposing is a cost-effective approach to quickly investigate whether current drugs could be used to potentially improve the clinical outcome of a patient [2][3][4][5]. To our knowledge, there are no approved drugs till today which have been shown to be safe and effective in terms of reducing the number of deaths and preventing severe COVID-19 [6][7][8].
Due to the urgency, scientists investigate the effect of different drug mixtures to propose an alternative, perhaps more effective, treatment, often at any cost [9]. However, initiating ineffective therapies and trying drug cocktails can be crucial, especially to patients with comorbidities on top of acute respiratory distress syndrome that require mechanical ventilation [10]. In this sense, co-medications can cause phenoconversion, a case where there is a discrepancy between drug-metabolizing genotype and phenotype status of an individual [11]. Inflammation and infection processes [12] in COVID-19 could constitute the causative ground for therapeutic failures that may be due to genetic interactions, yet to be explored [13,14]. To this end, available reported drugdrug interactions (DDI) accompanied by a clinical outcome and information about possible adverse effects may help health professionals in treatment strategy decisions when comorbidities exist and co-administration of drugs is advised [15,16].
Pharmacogenomics (PGx) approach uses an individual's genetic information to predict drug response and guide optimal drug dosage for safer, and cost-effective treatments [17]. The recent advances in pharmacogenomics in combination with the rapid development of next-generation sequencing technology have led to remarkable findings such as the identification of various drugrelated genetic loci variations associated with complex diseases [18]. By knowing the genetic interactions which affect the comedication response, common genes could be considered as candidate biomarkers for optimum treatment and clinical management strategies. On the other hand, there are a limited number of methods to correctly estimate all possible drug-drug interactions (DDIs) by in silico approaches [19] and predict the successfulness of a COVID-19 treatment by just considering disease progression.
In this study, we use a computational approach to evaluate all possible pharmacogenes, with a recorded clinical outcome, which influences co-medication response in order to uncover the significant pharmacogenetic determinants that affect the drug response of COVID-19 patients [13]. We identify a cluster of genes, known risk factors for other diseases, and propose them as pharmacogenomic biomarkers of relevance in COVID-19 treatment and perhaps long-term outcomes. Our findings could help in clinical decision-making to combat COVID-19 pandemic with medication [13].

Drug-drug interactions (DDI)
All DDIs with recorded clinical adverse effects have been recorded in an excel file along with its pharmacogenomic evaluation through PharmGKB [27]. Using the annotation files from PharmaGKB along with the list of 375 drugs (from Liverpool datasets and driven by a data-mining algorithm), a drug-gene association list was created (Supplementary File 1, sheet 2). Importantly, each drug from various databases was tested for synonyms in the DrugBank database [23]. Drug-gene interactions extracted from PharmGKB for both drugs participating in a DDI were linked (Supplementary File 1, sheet 3). Drugs without PGx information were discarded from the analysis (Supplementary File 1, sheet 3).

Statistical assessment of DDIs
For each DDI, common (intersection) genes, the union of genes, and different genes were recorded. A statistical method was designed using the STATA package [28] to determine the significance of the overrepresentation of common genes (between two interacting drugs) that affect co-medication response. Briefly, a hypergeometric test (Fisher-EXACT) [29] was performed for each DDI using as background the total number of genes (998) associated with at least one drug. The test returned a score as the probability of randomly selecting successful common DDI pharmacogenes (PGx genes) in the set of union genes of a certain DDI compared with the finite population of the background genes (Supplementary File 1, sheet 4). The most significant DDIs are those enriched for common genes with a p-value < 0.05 and used as input for further analysis.

Functional enrichment analysis of the highly significant DDI PGx genes
The statistically significantly enriched genes (44) which commonly associated with both drug responses of a DDI were listed and used as input for pathway analysis (KEGG, Reactome) and functional (Gene Ontology) enrichment analysis (Supplementary File 1). The latter was done with the use of the FLAME web-tool [30]. Further functional enrichment analysis of the 44 PGx genes was performed with the use of STRING database at a protein-protein interaction level and by only selecting interactions which were curated and experimentally validated.

Gene-disease network
All 44 genes of statistically significant DDIs were further investigated in GWAS [24], OMIM [25], and GAD [26], to identify pathophysiological diseases (phenotypes) that are linked with their genetic variations using a recently developed bipartite network analysis methodology [31]. For this task, Cytoscape [32] network visualization tool was utilized.

Shared-disease gene-gene network
A shared-disease gene-gene network was constructed by directly connecting the 44 identified PGx genes (so-called PGx-biomarkers from now on) associated with a certain disease. Similarly to before, the network visualization was performed with Cytoscape [32].

Creation of the expanded PGx-biomarkers network
An expanded version of genetic interactions identified through the shareddisease gene-gene network and through PPI interactions (STRING) for the 44 (statistically significant genes of all DDIs) PGx-biomarkers was constructed by merging and overlaying both networks using the Dynet software [33] to visually evaluate common and different edges between nodes.

Drug-drug interactions (DDI)
The Liverpool dataset [20][21][22] was used to retrieve all Drug-Drug Interactions (DDIs), from COVID-19 treatments where clinical adverse effects were reported. 375 drugs (or medication schemes) were found in total, 26 of which were repurposed drugs for COVID-19 treatment. In total 1989 DDIs were recovered according to the data mining workflow described below (see Material & methods along with Fig. 1A). DDIs are divided into four categories according to their clinical significance: (a) drugs that should not be co-administered (score 3), (b) drugs with potential interaction that may require a dose adjustment or close monitoring (score 2) and (c) drugs with potential interaction likely to be of weak intensity (score 1) (Supplementary File 1, sheet 1). Drugs with no clinically significant interaction were not considered in our analysis. From the significant interactions (DDIs), 15% belong to the first category, 70% to the second, and 15% to the third category.
Drug associated genes Data mining in PharmaGKB using the 375 drugs as input, reported 998 pharmacogenes affecting the response of at least one drug (Fig. 1B). According to the PharmGKB level of significance, patients who carry variants in these genes should be prescribed with caution regarding drug safety and efficacy. However, nine out of the 26 repurposed COVID-19 treatments do not have PGx information in PharmaGKB (Supplementary File 1, sheet 3). Thus, we propose that more research is needed to identify PGx associations with such drugs.

Identification of PGx biomarkers
From the 1989 DDIs, only 1413 contained PGx information for both drugs. Hypergeometric tests were performed on each of the 1413 DDIs and resulted in 571 statistically significant DDIs. Interestingly, only seven interactions were found for COVID-19 Adjunct Therapies drugs, 127 for COVID-19 Immunotherapy drugs, and 443 DDIs were found for antiviral drugs (Table 1). A great deal of all these interactions fell into score category 2 i.e., drugs with potential interaction that may require a dose adjustment or close monitoring. Common genes of the two drugs in each DDI pair were considered if p-values < 0.05 according to the Fisher-Exact test. The sum of the common genes was 44 and presented the most significant association with drugs in DDIs (Supplementary File 1, sheet 4, Fig. 1C).

Protein-Protein interactions of the PGx biomarkers
When these 44 genes were analyzed for their curated PPI interactions (co-expression, experiments, co-occurrence, databases, interaction) in the STRING database, a low-density interaction network was formed with four main interacting groups, such as xenobiotics metabolism partially overlapping with arachidonic-linoleic metabolism, cholesterol metabolism, immune response, and signal transducers some of which participate in tyrosine kinase pathways (Fig. 1D).

3D visualization of the association of DDIs with PPIs
To visualize the proportion of COVID-19 drugs and other disease drugs that are associated with gene variants regulating drug response, a three-layer interaction network was created by Arena 3D web [34]. This 3D network connects the DDIs from Liverpool data with the PPIs (gene-gene interactions) from STRING with the aid of drug-gene interactions according to PharmGKB. The top layer shows the 571 drug-drug interactions (drugs combined in COVID-19 treatments) that have significant common PGx associations (Fisher Exact test, p-value < 0.05). The bottom layer represents the network between all 44 pharmacogenes significantly associated with toxicity due to drugs administered for COVID-19 treatment or coadministered for comorbidities. Edges in each layer represent the significant drug-drug interactions with a recorded clinical effect and the public protein interactions found in the STRING database (Fig. 2). In this way, we see which interacting drugs are associated with genes that they, themselves, interact with others. Upon enrichment, we then investigated whether certain drug categories participate in one or multiple DDIs. As shown in Table 2, there is no specific trend for any drug category to continuously present one type of adverse effect, apart from the fact that most DDIs in general (70%) are in score 2. Importantly, antiviral therapies present adverse effects mainly due to genes participating in xenobiotic/drug metabolic pathways (Table 2). Adverse effects from adjunct therapy drugs are attributed to genes responsible for Autoimmune, metabolic, and neurological diseases as well as to genes important for Drug Metabolism and related blood disorders.
Gene-disease association network Next, a gene-disease network of the 44 pharmacogenes significantly associated with DDIs (Fig. 3, Supplementary File 3) was constructed depicting genes associated with diseases according to GWAS, OMIM, and GAD [31]. Visualization was performed with Cytoscape. Three major clusters of diseases became prominent: (a) Autoimmune, metabolic, and neurological diseases, (b) Cardiovascular and other degenerative diseases, (c) Drug Metabolism and related blood disorders. Additionally, a small number of genedisease associations became prominent, however, they are disconnected from the major gene-disease clusters.
Shared-disease gene-gene network A shared-disease gene-gene network was then constructed by directly connecting the genes of the gene-disease network (associated with the same disease), with one edge (Fig. 4, and .05 to determine the statistical significance of having overrepresented common genes in both drugs of the same co-medication with a clinical interaction. D Significant DDIs (p-value < 0.05) are selected, and their gene list has been used as input in the STRING database to find curated interactions between drug-interacting genes of each DDI. E Gene-disease associations of the above identified genes, based on GWAS, OMIM and GAD datasets, were retrieved and a new, sheared-disease gene-gene network was constructed with edges depicting genes associated with a common disease. F Networks of (D) and (E) are combined together providing an expanded PGx biomarker network that shows statistically significantly associated genes with COVID-19 treatment adverse effects. . Remarkably, the extra information from the PPIs resulted in expanding each gene cluster rather than connecting the clusters, suggesting that the text-mining approach (from Liverpool data) we followed was of high specificity. The resulting expanded network expands, in a statistically significant manner, the number of genes that can be associated with a cluster of related diseases.

DISCUSSION
Bioinformatics has been very helpful in identifying potential molecular mechanisms and key genes involved in COVID-19 [35]. Due to the limited data on common PGx affecting concomitant medications response, to our knowledge, this is the first attempt   which tries to combine bioinformatics, data mining, and network biology approaches to mine and analyze drug interactions of repurposed COVID-19 treatments. This effort aims at predicting the most significant genes which minimize the effectiveness and/ or safety of drugs and could be possible early prognostic markers in relation to COVID-19. The predicted PGx network could initiate the development of clinically applicable diagnostic and prognostic tests for patients at high risk of COVID-19 and those who cannot be vaccinated due to health conditions and allergies [36]. Extended network analysis of these markers' genetic interactions indicates that patients who carry variants of these genes may have long-term cardiovascular, immune, or neurological effects after infection with SARS-CoV2 as they are found to be known risk factors of heart, neuron, liver, and metabolic diseases. Several predicted genes in the network are currently used as prognostic biomarkers (Table 3), however, some others (14 genes) are not yet verified in their relation to COVID-19, and thus, through this work, they could be of great benefit for the prevention and assessment of disease severity and could play a role for the disease outcome. The receptor neurotrophic tyrosine kinase (NTRK1) is involved in the development and the maturation of the central and peripheral nervous system, and it is the only unverified gene in the first  (Fig. 4, Supplementary File 3). Gray edges are common edges between the two networks. The line width indicates the number of common diseases between two genes as shown in Fig. 4. The size of a node indicates the degree connectivity. cluster of the network. NTRK1 is associated with insensitivity to pain and thyroid carcinoma but it doesn't share common disease with VEGFA. However, there is a curated interaction between NTRK1 and VEGFA gene suggesting that both are participating in the same signaling pathway. A recent study suggests that SARS-CoV-2 spike protein interferes with endogenous ligand VEGFA, promotes a signal arriving at the central nervous system, and stimulates specialized sensory receptors in the peripheral nervous system inducing analgesia [37,38]. Thus, small-molecule inhibitors of this signaling for the treatment of neuropathic pain and cancer are being tested and could have added potential of inhibiting SARS-CoV-2 virus entry [39].
Network analysis reveals some genes that are central and tend to connect different clusters of diseases (bottlenecks) that may play a significant role in disease outcome in patients with comorbidities. Such genes in the first cluster are the ITGA2/ITGB3 and IRS1 that connect pathways of the neurological system with cardiovascular respectively. Moreover, IRS1 has a central node that interacts with NTRK1 and ITGB3 via curation suggesting that it is an immediate connector of genes related to cardiovascular and degenerative diseases. All genes of the second cluster seem to be central with dense connectivity with apolipoproteins APOC1, APOE, APOC3 which are early prognostic biomarkers for severe COVID-19 and share common diseases of cholesterol and mental disorders with genes such as TOMM40, the mitochondrial import receptor subunit. TOMM40 is an unverified gene in relation to COVID-19 but it is known as a risk factor for metabolic or neurological and mental disorders.
Interestingly, dementia was among the common comorbidities and was associated with higher mortality due to the APOE homozygous genotyping in hospitalized COVID-19 patients [40][41][42]. These findings strongly highlight the possible cognitive  -DRB1  HLA-DQA1  CTLA4  VEGFA  ITGB3  ITGA2  IRS1  MT-ND3  NOS3  ABCB2 HLA is associated with risk for severe COVID-19. CTLA4 expression levels were correlated with viral levels. VEGFA increases endothelial dysfunction and correlates with COVID-19 disease severity, stimulates sensory receptors in central and peripheral nervous systems, cooperates with sars-cov2 spike protein to induce analgesia. CTLA4: HLA: VEGFA are highly expressed and associated with progression, immune regulation, and symptoms in COVID patients. ITGB2/3 SARS-CoV-2 infection revealed changes in genes related to coagulation (ITGB3). IRS1 Sars-cov-2 viral infection activates stress response, which induce IRS-1 phosphorylation and insulin resistance. MT-MD3 mitochondrial gene downregulated in sars-cov-2. NOS3 constitutes an important endothelial protection mechanism. However, ARDS diffuse inflammatory process triggers a vasodilatation cascade in non-ventilated parts of the lungs (deregulating hypoxic vasoconstrictive mechanisms) and vasoconstriction in ventilated areas. The imbalance between vasoconstricting and vasodilating pathways leads to endothelium dysfunction. ABCB contributes to drug resistance, lysosomal accumulation of COVID-19 treatments and related to virus replications. impact of COVID-19 due to the APOE/APOC cluster with TOMM40 and their interactions [43][44][45] and could be also potential genetic biomarkers of COVID-19 severity. In addition, network analysis shows SORCS2 as a bottleneck that connects the 2nd and the 3rd cluster of drug metabolic genes through PTGS1 and shares common diseases of myeloid leukemia and cholesterol. PTGS1 according to literature has an unknown role in COVID-19, but its inhibition with COVID-19 treatments triggers upregulation of IL10 gene expression and represses platelet aggregation [46]. A distant connected sub-cluster of SORCS2:B4GALT2:GPX1 genes is associated with blood protein levels and intelligence. SORCS2 and B4GALT2 genes were not found as verified in relation to COVID-19 however, they are found as risk factors of diseases in cholesterol and protein levels, and they may indirectly increase the susceptibility of patients to SARS-CoV-2 and increase the risk of death. The NR1I2 gene, a member of the nuclear receptor subfamily 1 Group, is associated with disease of C-reactive protein, maybe another significant node as in COVID-19 the elevated levels of CRP might be linked to the overproduction of inflammatory cytokines in severe patients [47]. We suggest that elevated levels of C-reactive protein may be an early marker to predict risk for severity of COVID-19 and NR1I2 may be a prognostic biomarker as well as its neighbors [47]. In the third cluster, several CYP genes are known as PGx biomarkers affecting the metabolism of repurposed drugs and known risk of liver disease. Our analysis has predicted multiple bilirubin-associated genes including ABCC2, SLCO1B1, G6PD, UGT1A1/3/6/7, and other essential hypertension and blood pressure-associated genes PEAR1, G6PD, SLO1B1, ABCC2, ENOSF1 densely interacting with CYPs. Future work is needed to identify PGx associations with repurposed COVID-19 treatments that do not have PGx information in PharmaGKB and include more interactions in the network. A potential solution to this problem would be to include variants and genes having predicted effects derived from existing publicly available computational prediction methods [ref. [1][2]. However, such an approach has many problems related to the large number of potential "targets" that need to be evaluated and added to the model. Nevertheless, it is a possible extension of our approach that needs to be investigated in future studies. In addition, evaluation of population-specific PGx landscape genetic analysis of the genes proposed as biomarkers is needed that might aid in a better understanding of the inconsistency in therapy response [48,49].

DATA AVAILABILITY
All data generated or analyzed during this study are included in this published article and its Supplementary Information files.