Prediction of severe adverse events, modes of action and drug treatments for COVID-19’s complications

Following SARS-CoV-2 infection, some COVID-19 patients experience severe host driven adverse events. To treat these complications, their underlying etiology and drug treatments must be identified. Thus, a novel AI methodology MOATAI-VIR, which predicts disease-protein-pathway relationships and repurposed FDA-approved drugs to treat COVID-19’s clinical manifestations was developed. SARS-CoV-2 interacting human proteins and GWAS identified respiratory failure genes provide the input from which the mode-of-action (MOA) proteins/pathways of the resulting disease comorbidities are predicted. These comorbidities are then mapped to their clinical manifestations. To assess each manifestation’s molecular basis, their prioritized shared proteins were subject to global pathway analysis. Next, the molecular features associated with hallmark COVID-19 phenotypes, e.g. unusual neurological symptoms, cytokine storms, and blood clots were explored. In practice, 24/26 of the major clinical manifestations are successfully predicted. Three major uncharacterized manifestation categories including neoplasms are also found. The prevalence of neoplasms suggests that SARS-CoV-2 might be an oncovirus due to shared molecular mechanisms between oncogenesis and viral replication. Then, repurposed FDA-approved drugs that might treat COVID-19’s clinical manifestations are predicted by virtual ligand screening of the most frequent comorbid protein targets. These drugs might help treat both COVID-19’s severe adverse events and lesser ones such as loss of taste/smell.

Large scale benchmarking of MEDICASCY MOA predictions. The first step uses MEDICASCY 12 to predict MOA protein targets of all diseases. For large scale benchmarking for MOA prediction, we map all drugs in our indication library to DrugBank drugs 13 (v5.09) and obtain their respective human protein targets. These are combined with those from the Therapeutic Target Database 14 . Using drug-indication relationships in our training library, we compiled indication-protein target relationships of 145,722 pairs for 3539 indications (with an average ~ 41 proteins/indication) for benchmarking. In benchmarking, any drug in the training library having a Tanimoto coefficient, Tc 15 ≥ 0.8 to the given drug whose indications are predicted is excluded. We define a MOA prediction for an indication when its p-value < 0.05 using the upper tailed null hypothesis. 43.7% of indications have correctly predicted protein-indications. If a Tc = 1 cutoff in training is used, this increases to 65.9%. However, the incompleteness of the known MOA targets suggests that these are lower bounds. Table S0, large scale benchmarking of LeMeDISCO on different clinical data sets shows high comorbid disease coverage and accuracy compared to alternative methods [16][17][18][19] . It is also superior to alternatives that rely solely on symptom data, which lack molecular mechanism-based associations and cannot provide information of confidence-ranked putative protein targets. In practice, LeMeDISCO's recall rate is close to 70% for a representative set of 2630 disease pairs.

COVID-19's clinical manifestations.
We first predicted COVID-19 comorbidities using the 332 high confidence human proteins that interact with SARS-CoV-2 6 . There are 916 significant comorbidities (with a p-value cutoff < 0.05), of which 458 map to a COVID-19 clinical manifestation group. The top two disease comorbidities ranked by their p-value are shown in Table 1, with an expanded list in Table S1. Also provided are comorbidity enriched protein targets. Without extrinsic information or training, MOATAI-VIR recapitulates many key COVID-19 phenotypes such as myelosuppression, immunodeficiency, neurotoxicity, blood indications, myocardial infarctions, stroke, and cytokine storm symptoms 2,6,20 . The ICD-10 code of the comorbid diseases was used to map them to the 30 complications in 11 . In practice, indications are mapped to 21/30 COVID-19 complications (see Table 1). Our library of 3608 indications 10 does not have these 4/30 complications: Dialysis initiation, Intracranial hemorrhage, Hypertensive crisis, Cardiogenic shock. Thus, the complication recall rate is 21/26 ~ 81%. These mapped indications are then used to prioritize MOA proteins, pathways, and predict drugs for each complication.
The 6 human genes near the 3p21.31 locus of the human genome identified in a GWAS study as strongly associated with respiratory failure in COVID-19 patients 5 (odds ratio 1.77) were next used to predict comorbidities. We determined 598 significant comorbidities having a p-value < 0.05, of which 360 map to a clinical manifestation group. As shown in Table 2, many severe clinical complications associated with COVID-19 are predicted including respiratory complications, myocardial infarction, and cytokine storms. Table S2 provides an expanded list including myocardial infarction, stroke, neurological manifestations, hearing disorders, hypoxemia, lung, cardiovascular and diabetic risk factors 2,20,21 . Excluding 4 indications not in our library, with GWAS risk gene input, the recall rate of the COVID-19 complications is also 21/26 ~ 81%. Since the 332 proteins of the SARS-CoV-2 human interactome and the 6 GWAS COVID- 19  www.nature.com/scientificreports/ olfactory receptors may be an alternative SARS-CoV-2 entry into the local host cells, which may lead to its spread into the central nervous system 23 .
Neurological symptoms. There have been a number of unusual COVID-19-related neurological symptoms such as stroke, confusion, and as previously mentioned, loss of sense of smell and taste 24 . The most frequent clinical manifestation group for the interactome results was neurologic whose top pathway is glutathione conjugation associated with facilitating xenobiotic metabolism. Dysregulation of glutathione plays a role in many diseases including neurodegenerative diseases and cancer 25,26 . Decreased glutathione levels can lead to oxidative stress, resulting in Parkinson's and Alzheimer's disease. Moreover, an imbalance in glutathione levels can impact the immune system 26 . Decreased glutathione concentration is highly associated with serious manifestations causing increased COVID-19 mortality, possibly from increased susceptibility to uncontrolled viral replication 27 . The most frequent clinical manifestation group from the GWAS results was also neurological. The second top pathway was Na + /Cl− dependent neurotransmitter sodium symporters, which use sodium and chloride electrochemical gradients to import/export several substrates. They are associated with Parkinson's disease, orthostatic intolerance, and depression 28 . A meta-analysis found that low blood sodium increases the risk and www.nature.com/scientificreports/ severity of COVID-19 29 . Thus, neurotransmitter transporters that depend on Na + could be dysregulated due to decreased blood sodium levels.
Blood clotting. Increased blood clots are seen in COVID-19 infected individuals which causes increased mortality 30 . The third most frequent clinical manifestation group for the interactome results was hematologic. RAB geranylgeranylation, the third most frequent pathway across all clinical manifestations from the interactome, is a post-translational modification that allows RABs to connect with intracellular membranes where they  Cytokine storms. Several COVID-19 patients face respiratory complications, which may be due to antiinflammatory cytokines. Using the GWAS input, we find the ADOR2AB mediated anti-inflammatory cytokine production pathway associated with the respiratory manifestation. Over-secretion of interferons can yield uncontrolled systemic inflammation 32 . The ADOR2AB mediated anti-inflammatory cytokine production pathway is involved in 7/21 clinical manifestations. GWAS input also yields the interleukin-1 signaling pathway associated with the CHF clinical manifestation. Interleukin-1 is up-regulated in CHF patients and is a target for treating heart-related diseases 33 .
The second most frequent clinical manifestation group was cardiovascular. One significant pathways is G alpha (s) signaling events, involving 7/21 GWAS clinical manifestation groups. This pathway activates adenylate cyclase producing cAMP. G-protein receptors are associated with heart disease. Among the proteins in this pathway, C5aR1 is a G-protein-coupled receptor. The C5a-C5aR1 complex is involved in COVID-19 progression and is part of a potential therapeutic strategy 34 .This complex is associated with the innate immune response, with C5 a key driver in complement-mediated inflammation 34 .
Examples of immune-related pathways from the interactome clinical manifestation results include regulation of actin dynamics for phagocytic cup formation, the NLRP1 inflammasome, and RUNX3 regulation of immune response and cell migration.
Uncharacterized manifestations. 480 significant comorbidities were not mapped to a known COVID-19 clinical manifestation group from the interactome results. To further understand their effects, we grouped them by their main ICD-10 classification and performed CoPathway analysis. The top three uncharacterized groups were neoplasms, Congenital malformations/deformations/chromosomal abnormalities, and digestive system diseases. The comorbidities, comorbidity enriched MOA proteins and pathways for the interactome results are shown in SI, Table S3.
285 diseases were not mapped to a clinical manifestation group from the GWAS results. Their comorbidities, comorbidity enriched MOA proteins and pathways for the GWAS results are shown in Table S4. The top three uncharacterized groups involve mental and behavioral disorders, diseases of the digestive system, and neoplasms.
Neoplasms. The interactome input resulted in 90 pathways with a p-value < 0.05 involving neoplasms, with many involving hormonal regulation. The top pathway, activation of AMPK downstream of NMDARs, is associated with the neuronal system. AMPK is an enzyme that regulates cellular energy and homeostasis via activating catabolic pathways while switching off cellular growth and proliferation 35 . AMPK has been targeted for cancer treatment because its activation can reduce cancer incidence. NMDARs control synaptic plasticity and memory. Increased expression of NMDARs occurs in a variety of cancers such as neuroblastoma, breast, small-cell lung, and ovarian cancer 36 . Anti-NMDAR encephalitis, characterized by abnormal neurological and behavioral symptoms, has been reported in both COVID-19 and herpes simplex virus 2 (HSV-2) 37,38 . Notably, HSV-2 can lead to an increased cervical cancer risk 39,40 .
GWAS provides 29 neoplasm related pathways with a p-value < 0.05. The top pathway is the nuclear receptor transcription pathway. Nuclear receptors are DNA-binding transcription factors capable of binding hormones, vitamins, small molecules, and other ligands. A number of underlying disease mechanisms associated with dysregulation of nuclear receptors that can result in cancer, diabetes, and hormone-related conditions. Nuclear receptors have been targeted by cancer therapeutics as they are key players in gene regulatory networks 41 . There has not been substantial research on the relationship between nuclear receptors and COVID-19; but we note that some viruses target nuclear receptors as part of their replication process 42 .
Is SARS-CoV-2 an oncovirus? As indicated above, there were many disease comorbidities associated with neoplasms from both the interactome and GWAS results. Neoplasms cause abnormal tissue growth, a significant cancer characteristic. Perhaps, SARS-CoV-2 hijacks the human host replication machinery or proliferation pathways 6 . Indeed, viruses can initiate signal transduction pathways leading to cytokine and chemokine expression. They also dysregulate signaling pathways to promote viral infection and cellular transformations 43 that elicit a proinflammatory response similar to cancer 44 . A salient example is Human Papillomavirus (HPV). Most cervical cancers 45 are caused by the cytokine flux associated with inflammation post-HPV infection 40 . Furthermore, the second most significant COVID-19 comorbid disease, T-cell leukemia, is linked to the human T-cell lymphotropic virus (HTLV-I), an RNA retrovirus. More generally, a number of oncoviruses cause cancer 46 . Certain viruses transform human cells causing loss of ability to regulate cell division.
Although we do not yet know the long-term consequences post-COVID-19 infection, these results raise the distressing possibility that SARS-CoV-2 is an oncovirus. To assess this potential relationship, we screened our comorbidity enriched MOA proteins associated with neoplasms from both the interactome and GWAS sets against the COSMIC 47 database gene set containing 723 oncogenes. There were 1488 and 97 neoplasm comorbidity enriched MOA proteins using the interactome and GWAS as input, respectively (comorbidity weighted frequency > 0.1). From this, 12.3% (n = 89) and 0.97% (n = 7) of the interactome and GWAS comorbidity enriched MOA proteins, respectively, are oncogenes in the COSMIC database (42); see SI Tables S6, S7. The overlap between the neoplasm comorbidity enriched MOA proteins from the interactome as input and the COSMIC database oncogenes resulted in a significant p-value of 2.5 × 10 -5 . The overlap between the GWAS neoplasm comorbidity enriched MOA proteins and COSMIC yielded an insignificant p-value of 0.08. www.nature.com/scientificreports/ As further substantiation of the conjecture that SARS-CoV-2 is an oncovirus, we compared the differential gene expression analysis of COVID-19 patients (n = 1918 differentially expressed genes with an adjusted p-value < 0.05), to the COSMIC database 47 . 11% (n = 82) of the genes overlap with oncogenes in the COSMIC database (see SI Table S8). Although the overlap p-value is not significant, this indicates that there are some overlapping oncogenes possibly associated with COVID-19. We next performed a 3-way merge between the interactome/GWAS neoplasm comorbidity enriched MOA proteins, the SARS-CoV-2 differentially expressed genes, and the COSMIC database oncogenes. We found 11 and 1 overlapping gene(s) from the interactome and GWAS 3-way merge, respectively. The pathway analysis on the overlapping COVID-19 differentially expressed genes and the COSMIC database indicates that viral replication and oncogenesis employ similar biochemical mechanisms. Indeed, a number of the identified pathways such as interferon-gamma signaling 48,49 , immunoregulatory interactions between lymphoid and non-lymphoid cell 50 , and antigen processing-cross presentation 51,52 are related to viral replication and oncogenesis. Clearly, additional investigation is needed to explore the possibility that SARS-CoV-2 might be an oncovirus. If this conjecture were true, it would provide an even greater incentive for people to get vaccinated.

Comparative study of viruses for oncogenic propensity.
To further support our conjecture that SARS-CoV-2 might be an oncogenic virus, we compiled the virus-human host interacting proteins of 13 viruses from the literature [53][54][55][56][57][58][59][60][61][62][63] . After applying the same LeMeDISCO procedure as done for SARS-CoV-2 to prioritize these interactomes, we examined the top 100 proteins and their overlap with the 723 COSMIC census putative cancer drivers 47 . Here, LeMeDISCO uses only the top 100 comorbid indications associated with Neoplasms (mapped according to their ICD-10 main codes) to prioritize the MOAs of these Neoplasm indications. We use the p-value (calculated using Fisher's exact test 64 ) of the overlapped proteins as the propensity that ranks the virus's likelihood of being oncogenic. In practice, if the p-value < 0.05, then we consider the virus is oncogenic.
The results are summarized in Table S9. For the 9 known oncoviruses, we failed for 3: Ad5, KSHV and HTLV. However, for the 4 viruses that are not oncoviruses, none has a p-value < 0.05. This indicates for these viruses, the false positive rate is 0. SARS-CoV-2 has a p-value of 0.016 that is close to those (0.04) of HIV and PyV oncoviruses. Thus, SARS-CoV-2 is closer to oncoviruses than to non-oncoviruses. We should point out that it may take over a decade or longer for cancer to emerge post-infection 65 . Thus, increased cancer rates might be a long term consequence of COVID-19. At present, there is a report lung cancer metastases have increased during the pandemic 66 as well as a conjecture that SARS-CoV-2 can induce glioma tumorigenesis 67 .

Predictions of possible repurposed drug to treat COVID-19's comorbid diseases.
To identify potential repurposed drugs from DrugBank 13 that might treat a given SARS-CoV-2 complication, MOATAI-VIR utilizes the LeMeDISCO disease profiles for identifying comorbidity-based treatments. For COVID-19's clinical manifestations, possible treatments are identified by CoMOAdrug and CoVLS. CoMOAdrug identifies drugs that might treat the comorbid diseases to the given disease via indication-based virtual ligand screening using MEDICASCY 12 . A drug is ranked by the fraction of comorbid diseases with that indication weighted by the drug-indication predicted precision inferred from MEDICASCY benchmarking 12 . This yields a rank ordered list, LIST CoMOAdrug , for drugs common to the comorbid diseases. CoVLS identifies efficacious drugs by FINDSITE comb2.068 virtual ligand screening of the comorbidity frequency weighted MOA proteins of the given adverse response. A drug is ranked by the product of the predicted molecule's binding precision times the p-value weighted frequency that the protein is a MOA protein of a comorbid disease divided by the summed binding precision of all its human targets. The resulting list is LIST CoVLS . A highly ranked drug often has multiple targets with high p-value weighted comorbidity frequency; for additional details, see SI.
Benchmarking of the LeMeDISCO approach to identify drugs. To prove that comorbidity-based drug ranking works in principle, using MEDICASCY we successfully predicted novel, anti-proliferative small molecules in 79.4% of the top ranked 20 of 1597 molecules from the NCI diversity set in 10 different NCI-60 cancer cell lines 69 . If protein targets are selected based on their p-value weighted comorbidity frequency and binding precision, then the success rate is 85.7% for the top ranked 20 of 1597 molecules.
To further benchmark the ability to map drugs to their predicted indications that also shows the generality of our methodology, we utilized 2,059 training drugs with 123,146 drug-indications pairs in a modified jackknife test. For each drug, when predicting its indications, we used models trained from other drugs having a Tanimoto Coefficient Tc < 0.8 to the given drug. We then evaluated the top 20 drug predictions for each indication using the above two strategies and compared the results to the single indication based MEDICASCY predictions (see Table S5). Since all indications of a given drug are not known, the estimated precision is a lower bound to the true precision 70 . CoVLS has the best lower bound precision of 72.6%, MEDICASCY has the smallest lower bound precision of 58.0%, and CoMOAdrug is in-between, having a lower bound precision of 64.2%. The mean enrichment factors within the top 20 by CoMOAdrug and CoVLS are 7.05 and 7.97, respectively. Thus, CoMOAdrug and CoVLS show significant enhancement in their ability to select efficacious drugs over random. Tables 1 and 2 also present the top 2 predicted repurposed drugs to treat COVID-19's major adverse complications using both CoMOAdrug and CoVLS ligand ranking approaches based on the SARS-CoV-2 human interactome and GWAS, respectively. The full list of high-ranking repurposed drugs for the 24 mapped complications for inputs of human interactome and GWAS risk genes are given in Tables S1 and S2,  www.nature.com/scientificreports/ Next, ClinicalTrials.gov data (data obtained October 2020) was mined to determine how many of the drugs predicted for the COVID-19 complications are currently undergoing clinical trials for their respective COVID-19 complication. For the 24 predicted complications, the average number of FDA-approved drugs undergoing clinical trials is 4.1 (~ 21%) within the top 20 predicted drugs (4.3 for CoMOAdrug and 3.9 for CoVLS). In a sense, this is a soft validation of MOATAI-VIR. Drugs undergoing clinical trials for their respective COVID-19 complication are marked with "*" in SI Tables S1 and S2. Furthermore, drugs undergoing clinical trials for generic treatment of COVID-19 are marked with "$" in SI Tables S1 and S2.

Potential candidate repurposed drugs to treat COVID-19 complications.
To further establish the plausibility of the drug predictions, a literature search of several drugs from the comorbidity-associated complications was done. For example, Prednisolone, typically used for treating allergies and infections 71 , is the top drug predicted using CoMOAdrug for the respiratory clinical manifestation from the interactome results. Prednisolone is a glucocorticoid with anti-inflammatory, immunosuppressive, antineoplastic, and vasoconstrictive MOA. Corticosteroids reduce adverse events attributed to ARDS 72 . Prednisolone was also the top drug for the asthma exacerbation clinical manifestation from GWAS, whose top pathway was olfactory signaling. MEDICASCY predicts that Prednisolone has 142 indications which overlap 13 of COVID-19's complications or their comorbid indications.
The top drug predicted using CoVLS for the respiratory clinical manifestation from GWAS was Betamethasone, another corticosteroid with immunosuppressive and anti-inflammatory properties. Vindesine, an inhibitor of mitosis and a chemotherapeutic 73 , was among the top 5 drugs found from both the CoVLS and CoMOAdrug methods for the interactome hematological clinical manifestation. Vindesine can lower the number of platelets in the blood, thus, preventing blood clots 74 .
The immunomodulatory drug, Pomalidomide, an FDA-approved thalidomide derivative for treating multiple myeloma 75 was among the top predicted drugs by CoMOAdrug for the neurologic clinical manifestation from the interactome input. Pomalidomide has been investigated for use in neurological conditions such as Parkinson's disease. It can improve age-related neurological impairment/motor disability 76 and can reduce ischemic brain injury in an in vivo study 77 .
Temazepam was in the top 20 predicted drugs using the CoMOAdrug method for acute myocardial infarction/unstable angina, hematologic, neurologic, and asthma exacerbation clinical manifestations from the GWAS results. Temazepam is used to treat insomnia and is suggested as a treatment for sleep-related disturbances associated with COVID-19 78 .
Systematic validation of predicted drugs to treat COVID-19 complications. The above results are promising. However, drugs in clinical trials might not work in practice. Nevertheless, to have an idea of how well our method can do compared to similar drug prediction methods as in 79 for predicting drugs that directly treat COVID-19, we performed a similar assessment by using the clinical trial drugs as true positives and others as true negatives to calculate the area under the ROC curve (AUROC). For the interactome input, the mean AUROCs of CoMOAdrug and CoVLS are 0.70 and 0.73, respectively. For the GWAS input, they are 0.71 and 0.73, respectively. These are comparable to the non-AI-based methods of 79 . Our method does not have any AIbased method trained using SARS-CoV-2 related information whereas AI-net of 79 does.
To systematically validate our predictions, we utilized our training drug-indication dataset in MEDICASCY that was curated by experts. Since all predictions are carried out by training models on drugs having Tc < 1 to the drug whose efficacy is predicted, they are true predictions and not just memorization of known indicationdrug relationships. We use an enrichment factor within the top 20 of the total 2095 screened drugs (~ top 1%). Since the AUROC depends on true positives in the middle of the ranking, we instead calculate the AUPRC (area under precision-recall curve) which depends on true positives ranked at the very top. The results are shown in Table 5. The mean enrichment factors of CoMOAdrug and CoVLS drugs for interactome input are 3.91 and 14.6, respectively. For GWAS input, they are 3.36 and 7.99, respectively. These are consistent with a benchmarking test (Table S5) and are far better than random selection. The mean AUPRCs (around 0.09 to 0.184) are consistent with an earlier MEDICASCY benchmark 12 and better than those of the indication prediction method of Himmelstein et al. 80 (~ 0.005 to 0.1) where a systematic integration of biomedical knowledge was used for computing drug features and a logistic regression machine learning was employed for learning and prediction.
Predicted repurposed drugs to treat loss of sense of smell. Using CoMOAdrug and CoVLS, we predicted repurposed, FDA-approved drugs that possibly target the olfactory signaling pathway as prospective treatments for loss of the sense of smell. The top 20 drugs provided by each approach are shown in Table 6 Theophylline, a drug known to treat respiratory diseases, such as COPD and asthma was predicted by CoVLS and is currently undergoing clinical trials for anosmia (loss of sense of smell) 78 .

Discussion
MOATAI-VIR can identify possible molecular mechanisms responsible for COVID-19's severe adverse consequences. Not only are most of COVID-19's severe symptoms successfully predicted, but Tables 1 and 2 suggest a list of possible repurposed and mostly untested drug treatments for these complications. They could also be combined with antiviral drugs that directly target SARS-CoV-2 proteins to kill the virus. The goal is to mitigate both COVID-19 infection and its subsequent adverse complications to improve clinical outcome. Thus, MOATAI-VIR provides a series of logical, systematic suggested treatments for COVID-19's adverse reactions. Equally important, MOATAI-VIR is a general methodology for antiviral drug repurposing that can be applied to new outbreaks of other novel viral infections as they emerge.
One possible limitation of the current method as with all in silico drug predictions is that the predictions are not 100% accurate. Even though we have an enrichment factor much better than random, there are still many www.nature.com/scientificreports/ www.nature.com/scientificreports/ false positives. Further studies that could improve drug predictions that we will undertake in the near future are as follows: For MEDICASCY, we plan to include the more accurate protein structure models from AlphaFold 2 81 . Then, to improve its ligand virtual screening component, we shall employ the better quality structure models of AlphaFold 2 and the better screening method of FRAGSITE 82 as well as its consensus ligand binding proteins identified in combination with FINDSITE comb2.068 . Another limitation is that the current method is population based, i.e., comorbidities and drugs are predicted for the whole population, not specific individuals. For specific individuals, some comorbidities may not occur, and the predicted drugs may not work. These problems could be addressed by including personal genetic information and gene expression profiles of the appropriate tissues (e.g. the lungs) of SARS-CoV-2 patients to identify the specific proteins in that that are affected by SARS-CoV-2.

Materials and methods
Here, we give a brief description of the methods used in MOATAI-VIR, with additional details provided in SI. We start by performing large-scale prediction of MOAs for 3,608 indications using MEDICASCY 12 . MEDICASCY predicts MOAs for given indication by combining its drug-indication predictions and drug whole human genome protein target predictions from FINDSITE comb2.068 . A p-value is derived for each indication-human protein pair using Fisher's exact test 64 . Using a p-value cutoff of 0.05, for each indication, we define a list of putative MOAs.
To eliminate false positive predictions, the human atlas protein expression data is utilized to exclude proteins "not detected" in relevant tissues of a given indication 83 . Next, we apply LeMeDISCO (see SI for details) to examine the overlap proteins of the input interactome and GWAS sets to the MOA proteins of each of the 3,608 library indications to obtain their co-morbid indications. LeMeDISCO calculates a Jaccard index for ranking and the corresponding p-value using Fisher's exact test to determine the comorbidity of two sets of proteins. With a p-value cutoff of 0.05 for comorbidity, we obtained comorbid indications for the interactome and GWAS inputs. Then, comorbid indications are mapped to their respective COVID-19 complications. We subsequently employ CoPathway to determine significant pathways associated with the most frequent comorbidity enriched MOA proteins. We assess the frequency of MOA proteins across the comorbidities for a desired group and then processes the top ranked (comorbidity ranked p-value weighted frequency > 0.1) MOA proteins through the Reactome 84 for global pathway analysis. Pathways with a p-value < 0.05 are deemed significant . Lastly, two methods, CoMOAdrug and CoVLS, for drug discovery of the mapped indications were performed. CoMOAdrug screens FDA-approved drugs by combining the screening from MEDICASCY on the top 100 comorbid indications of a given indication. CoVLS screens drugs for efficacy by screening all comorbid frequency ranked MOAs of a given indication against FDA-approved drugs using FINDSITE comb2.0 .

Data availability
In addition to the Supplementary information associated with this paper, all Supplementary Tables are available on our website at https:// sites. gatech. edu/ cssb/ moatai-vir/ https:// sites. gatech. edu/ cssb/ moatai-vir/.