key: cord-0842947-r7dqnoq9
authors: Peng, Lihong; Shen, Ling; Xu, Junlin; Tian, Xiongfei; Liu, Fuxing; Wang, Juanjuan; Tian, Geng; Yang, Jialiang; Zhou, Liqian
title: Prioritizing antiviral drugs against SARS-CoV-2 by integrating viral complete genome sequences and drug chemical structures
date: 2021-03-18
journal: Sci Rep
DOI: 10.1038/s41598-021-83737-5
sha: e920fd8b20d773e93d2ec4af09568bf10e5d01ae
doc_id: 842947
cord_uid: r7dqnoq9

The outbreak of a novel febrile respiratory disease called COVID-19, caused by a newfound coronavirus SARS-CoV-2, has brought a worldwide attention. Prioritizing approved drugs is critical for quick clinical trials against COVID-19. In this study, we first manually curated three Virus-Drug Association (VDA) datasets. By incorporating VDAs with the similarity between drugs and that between viruses, we constructed a heterogeneous Virus-Drug network. A novel Random Walk with Restart method (VDA-RWR) was then developed to identify possible VDAs related to SARS-CoV-2. We compared VDA-RWR with three state-of-the-art association prediction models based on fivefold cross-validations (CVs) on viruses, drugs and virus-drug associations on three datasets. VDA-RWR obtained the best AUCs for the three fivefold CVs, significantly outperforming other methods. We found two small molecules coming together on the three datasets, that is, remdesivir and ribavirin. These two chemical agents have higher molecular binding energies of − 7.0 kcal/mol and − 6.59 kcal/mol with the domain bound structure of the human receptor angiotensin converting enzyme 2 (ACE2) and the SARS-CoV-2 spike protein, respectively. Interestingly, for the first time, experimental results suggested that navitoclax could be potentially applied to stop SARS-CoV-2 and remains to further validation.

www.nature.com/scientificreports/ LRLSHMDA 16 . These three methods were applied to biological association prediction in other application areas and obtained better prediction performance. We found that remdesivir and ribavirin come together on three datasets. Molecular docking is a key bioinformatics modeling tool for drug discovery and used to predict the "best-fit" intermolecular binding between a small molecule and a target or two proteins at the atomic level. It characterizes the behavior of ligands in the binding sites of target proteins as well as elucidates fundamental biochemical processes 17 . The docking process comprises two basic steps: predicting conformation, position, and orientation of ligands within the binding sites and ranking these conformations based on the binding affinity 18 . We used AutoDock 19 , a molecular docking software, to measure the molecular activities of the predicted two compounds (remdesivir and ribavirin) binding to the SARS-CoV-2 spike protein/human receptor angiotensin converting enzyme 2 (ACE2). The docking showed that remdesivir and ribavirin have higher binding energies of − 7.0 kcal/ mol and − 6.59 kcal/mol with the structure of the spike protein receptor-binding domain bound to the ACE2 receptor, respectively.

Experimental settings. In this section, we conducted extensive experiments to investigate the performance of our proposed VDA-RWR method. For the VDA matrix Y n×m from n viruses and m drugs, fivefold Cross-Validations (CVs) were performed under the following three different experimental settings.

• Fivefold Cross Validation 1 (CV1): CV on viruses, that is, random rows in Y (i.e., viruses) were selected for testing. • Fivefold Cross Validation 2 (CV2): CV on drugs, that is, random columns in Y (i.e., drugs) were selected for testing. • Fivefold Cross Validation 3 (CV3): CV on virus-drug pairs, that is, random entries in Y (i.e., virus-drug pairs)

were selected for testing.

Under CV1, in each round, 80% of rows in Y were used as training set and the remaining 20% of rows were used as test set. Under CV2, in each round, 80% of columns in Y were used as training set and the remaining 20% of columns were used as test set. Under CV3, in each round, 80% of entries in Y were used as training set and the remaining 20% of entries were test set. These three settings CV1, CV2, and CV3 specially refer to potential VDA identification for (1) new viruses (especially for SARS-CoV-2), (2) new drugs, and (3) new virus-drug pairs, respectively.

Parameters r , µ , and α denote the global restart rate, the transition probability, and the weight between the virus network and the drug network, respectively. For these three parameters, we performed cross validations on the training set to find the optimal values. In addition, the iteration stopped when p t+1 − p t 2 ≤ 1e − 11 .

SMiR-NBI need not set the parameters. For the parameters in NGRHMDA and LRLSHMDA, we conducted grid search to find the optimal values. The detailed settings are shown on Table 1 .

Evaluation metrics. Sensitivity, specificity, F1 score, accuracy and AUC were widely applied to evaluate the proposed methods. Sensitivity denotes the ratio of correctly predicted positive VDAs to all positive VDAs. Specificity is the ratio of correctly predicted negative VDAs to all negative VDAs (all the unknown associations were labeled as negative). F1 score denotes the harmonic mean of recall and precision. Accuracy represents the ratio of correctly predicted positive and negative VDAs to all positive and negative VDAs. We used these five metrics to evaluate the performance of VDA-RWR. They were defined as follows: (4) (5) .

For these five evaluation metrics, higher values represent better performance.

Performance evaluation under three fivefold CVs. We compared VDA-RWR with NGRHMDA 14 , SMiR-NBI 15 and LRLSHMDA 16 . NGRHMDA was presented to find potential microbe-disease associations by integrating neighbor-based collaborative filtering and graph-based scoring 14 . SMiR-NBI can comprehensively identify new pharmacogenomic biomarkers by constructing a heterogeneous network connecting genes, drugs, and miRNAs 15 . LRLSHMDA was applied to predict human microbe-disease associations based on Laplacian regularized least squares 16 . These three state-of-the-art approaches obtained good performance in their corresponding applications. We performed these four methods for 100 times on three different fivefold CV settings on three datasets. The final performance was averaged over the five rounds for 100 times. The results are shown in Tables 2, 3, and4. The best results were shown in bold in each column.

On dataset 1 and dataset 3, VDA-RWR outperformed other three methods in terms of specificity, accuracy, F1 score and AUC under three CVs. On dataset 2, although the sensitivity of VDA-RWR was slightly lower, VDA-RWR computed better specificity, accuracy, F1 score and AUC under majority of conditions. The slight difference can be produced by different data structures. AUC is one more important evaluation metric compared to other four measurements. AUC = 0.5 represents random performance and AUC = 1 shows perfect performance. www.nature.com/scientificreports/ VDA-RWR obtained the best AUCs under majority of conditions. In general, VDA-RWR is proper to discover potential VDAs. In addition, under CV1, VDA-RWR computed better specificity, accuracy, F1 score and AUC on the three datasets. This result showed that VDA-RWR can effectively find possible antiviral drugs for new viruses (for example, SARS-CoV-2). Under CV2, VDA-RWR outperformed other three methods in terms of specificity, accuracy, F1 score and AUC on dataset 1 and dataset 3. Although the sensitivity, specificity and accuracy values of VDA-RWR were slightly lower than other individual methods on dataset 2, it obtained the best F1 score and AUC. Thus AUC can identify potential viruses associated with new drugs. Under CV3, VDA-RWR calculated the best specificity, F1 score and accuracy. Figure 1 showed the AUC values of four methods under CV1, CV2, and CV3. The results demonstrated that VDA-RWR obtained relatively higher AUCs under three different CVs. It suggested that VDA-RWR could be used to infer potential VDAs.

Case study. In this section, we want to find possible drugs for SARS-CoV-2 after verifying the performance of our proposed VDA-RWR method. We predicted the top 10 drugs with the highest association scores with SARS-CoV-2 on three datasets. The results were shown in Tables 5, 6, and 7, respectively. Among the predicted top 10 small molecules associated with SARS-CoV-2 on dataset 1, all drugs were supported by recent works. Among the predicted top 10 chemical agents related to SARS-CoV-2 on dataset 2, there were 9 VDAs validated by current literatures, that is, 90% chemical agents were reported. Among the predicted top 10 antiviral drugs against SARS-CoV-2 on dataset 3, all compounds were validated by recent publications.

The results on Tables 5, 6, and 7 showed that there were two FDA-approved drugs coming together on three datasets, that is, remdesivir and ribavirin. Remdesivir is a small molecular compound undergoing a clinical trial and shows superior antiviral activity against many RNA viruses including orthocoronavirinae, filoviridae, paramyxoviridae, and pneumoviridae [20] [21] [22] . Sheahan et al. 17 presented that it can improve pulmonary function and reduce severe lung pathology in mice. Similar to SARS-CoV-2, both Ebola virus (EBOV) and MERS-CoV may result in severe acute respiratory diseases. And remdesivir has been used as inhibitors of EBOV and MERS-CoV 20,21 . More importantly, an array of works have reported that remdesivir is highly effective in controlling SARS-CoV-2 infection and has been directly applied to the treatment of COVID-19 6, 7, 9, [23] [24] [25] [26] [27] [28] . Specially, on October 22, 2020, FDA approved remdesivir for use in adults, pediatric patients with age of 12 years, and older and weighing at least 40 kg 29 . All these results showed that remdesivir may be the best anti-SARS-CoV-2 drug.

Ribavirin is identified as another anti-SARS-CoV-2 drug with a higher association score. Huang et al. 5 found that 28 of 38 patients treated by ribavirin have been discharged. Zhang et al. 30 reported that a patient has been treated with antiviral drugs including ribavirin. Therefore, ribavirin may be applied to treat COVID-19 caused by SARS-CoV-2. Interestingly, for the first time, experimental results suggested that navitoclax could be potentially applied to stop SARS-COV-2. Navitoclax has been applied to boost the treatment and basic science of chronic lymphoid leukemia, hematological malignancies, non-Hodgkin's lymphoma, solid tumors, and EGFR activating mutation. Molecular docking. The molecular docking between the above two antiviral drugs (remdesivir and ribavirin) and the spike protein and ACE2 are described in Table 8 . The results showed that remdesivir and ribavirin have higher binding energies of − 7.0 kcal/mol and − 6.59 kcal/mol with the structure of the spike protein receptor-binding domain bound to the ACE2 receptor, respectively. The subfigure in each circle denotes the residues at the binding site of the spike protein/ACE2 and their corresponding orientations. For example, the amino acids K68 and Q493 were predicted to be the key residues for remdesivir binding to the SARS-CoV-2 spike protein/ACE2 while K353, R403, Q493 and G496 were predicted as the key residues for ribavirin binding to these two target proteins.

In Table 8 , green denotes the structure of ACE2 and cyan denotes the SARS-CoV-2 spike protein in the figures of molecular docking. 

Finding possible antiviral drugs against SARS-CoV-2 is extremely urgent with the rapid spread of COVID-19. However, it seems very difficult to design a novel drug for COVID-19 within a very short time. One of efficient ways is to identify new clues of the treatment from FDA-approved drugs. In our proposed VDA-RWR method, we computed the association scores for each virus-drug pair to predict potential antiviral drugs against SARS-CoV-2 based on random walk with restart and biological information www.nature.com/scientificreports/ of viruses and drugs. The originality of our proposed VDA-RWR method remains, constructing three small datasets and inferring possible antiviral chemical agents against SARS-CoV-2 from FDA-approved drugs. The comparative experiments showed better performance of the VDA-RWR method. Higher AUC values under three fivefold CVs on three datasets and molecular binding energies indicated that the selected small molecules are likely to be used to stop the transmission of COVID-19. VDA-RWR can obtain superior performance under the three fivefold CVs on three datasets. This observation may be attributed to random walk with restart, a state-of-the-art model that can randomly walk on the heterogeneous virus-drug network and effectively compute association score for each virus-drug pair. More importantly, VDA-RWR integrated various biological information including the complete genome sequences of viruses and chemical structures of chemical agents.

The proposed VDA-RWR method is also helpful in design and interpretation of pharmacological experiment related to COVID-19. More importantly, VDA-RWR can be further applied to predict antiviral drugs against novel viruses without any associated chemical agents.

Virus-drug association data. Dataset 1. Virus data. We considered 11 viruses similar to SARS-CoV-2.

These viruses include influenza A viruses including A-H1N1 32 , A-H5N1 33 , and A-H7N9 34 , chronic hepatitis C virus (HCV) 35 Drug data. We manually curated drugs associated with these 11 viruses from the DrugBank 45 and NCBI 43 databases and published literatures reported by the PubMed database 46 and collected 78 small molecules after removing macromolecules. Based on the assumption that two drugs are more similar if they share more chemi- We obtained 770 VDAs from 69 viruses and 128 drugs after removing the viruses whose complete genome sequences are unknown from the database. The chemical structure of drugs and the complete genome sequences of viruses were downloaded from the DrugBank database and the NCBI database, respectively. Similar to dataset 1, we used RDKit and MAFFT to calculate drug similarity and virus similarity.

Dataset 3. We retrieved 407 VDAs from 34 viruses and 203 drugs by searching documents related to viruses and drugs based on text mining techniques. Similar to dataset 1, we computed drug similarity matrix and virus similarity matrix. The details of three datasets are shown in Table 9 .

In this study, the set of known VDAs was considered as the 'gold standard' dataset and was applied to evaluate the performance of our proposed VDA-RWR method. We described the known VDAs as a matrix Y: The VDA-RWR method. Inspired by the method provided by Valdeolivas et al. 50 , we developed a VDA prediction method based on Random Walk with Restart on the heterogeneous network (VDA-RWR). The proposed VDA-RWR method comprised two steps. First, a random walk-based model integrating various biological data was learned to explain the constructed 'gold standard' dataset. Second, this model was used to find potential VDAs for viruses and drugs absent from the 'gold standard' dataset. The details are shown Fig. 2 . We first considered virus-virus similarity graph G v , drug-drug similarity graph G d , and VDA graph G a , which formed a heterogeneous network. We defined S v(n×n) , S d(m×m) , and Y (n×m) as their corresponding adjacency matrices. The adjacency matrix of the heterogeneous network can be denoted as: The transition probability from virus v i to drug d j was defined as www.nature.com/scientificreports/ The transition probability from drug d i to drug d j was defined as

The transition probability from drug d i to virus v j was defined as For a given virus/drug, the particle can either jump between graphs or stay in the current graph with a defined probability r ∈ (0, 1) . Therefore, we finally defined the random walk with a restart probability r as:

where p t represented the computed association probability at the t-th step random walk. We defined the initial probability as:

, where u 0 and t 0 denoted the initial probability on the drug network and the virus network, respectively. If we tend to identify possible drugs for a given virus v i , it is considered as a seed node in the virus network. Here, v i was assigned as 1 and other nodes as 0, constructing the initial probability of the virus network t 0 . All nodes in the drug network u 0 were assigned as equal probabilities with the sum of 1.

For example, to find potential antiviral drugs against SARS-CoV-2, we set SARS-CoV-2 as a seed node, and all drugs in the drug network were assigned as the same probabilities with the values of 1/m . The parameter α was used to control the weight of the virus network and the drug network. In addition, a virus is new if it does not associate with any drugs, and a drug is new if it is not applied to any viruses.

abilities between the predicted anti-SARS-CoV-2 drugs and the SARS-CoV-2 spike protein/human ACE2. The chemical structures of drugs were downloaded in the form of the PDB format files from the DrugBank database. We used AutoDockTools to convert these PDB files into pdbqt files needed by AutoDock4. The structures of SARS-CoV-2 spike receptor-binding domain bound with ACE2 (PDB ID: 6M0J) were downloaded from the RCSB Protein Data Bank 51 . The spike protein and ACE2 were used as receptors, and the predicted anti-SARS-CoV-2 drugs were used as ligands for the molecular docking. We first removed solvent and organic compounds and preprocessed the receptor proteins based on PyMOL 31 (https ://githu b.com/schro dinge r/pymol -open-sourc e, release 2.4.0, open source license: BSD-like). The receptors' atoms were assigned the AD4 type and Gasteiger charges were considered before docking. Molecular docking software, AutoDock 19 (http://autod ock.scrip ps.edu/, AutoDock 4.2.6, open source license: GPL), was then used to conduct molecular docking. The binding pocket was defined by AutoGrid4, the grid size was set to 82 × 154 × 84 with a spacing of 0.375 Å, and the grid center was placed at the area of SARS-CoV-2 spike receptor-binding domain bounding with ACE2 (x = − 36.884, y = 29.245, z = − 0.005). The LGA (Lamarckian genetic algorithm) with default parameter provided by AutoDock4 was used as the search method. The docking contained two main processes: computation of conformation, position, and orientation of ligands within the binding sites and ranking of these conformations based on the binding affinities 18 .

To find potential antiviral drugs, in this study, we integrated the complete genome sequences of viruses, the chemical structures of drugs, and the VDA network. We then developed a VDA prediction method based on random walk with restart on the heterogeneous network. The results suggested that remdesivir and ribavirin may be applied to the treatment of COVID-19. In the emergency situation, this study focused more on finding antiviral drugs. In the future, we will further integrate more biological data and design more powerful models to improve the accuracy of VDA identification. We hope that our proposed VDA-RWR method could help the screening of drugs for preventing COVID-19.

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health: The latest 2019 novel coronavirus outbreak in Wuhan, China

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

WHO. Coronavirus disease 2019 (COVID-19) Situation Report-COVID-19 Weekly Epidemiological Update

US $675 Million Needed for New Coronavirus Preparedness and Response Global Plan [Internet]. Geneva: World Health Organization

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Drug treatment options for the 2019-new coronavirus (2019-nCoV)

Potential natural compounds for preventing SARS-CoV-2 (2019-nCoV) infection

Therapeutic drugs targeting 2019-nCoV main protease by high-throughput screening

Drug and vaccine design against novel coronavirus (2019-nCoV) spike protein through computational approach

Deep learning based drug screening for novel coronavirus 2019-nCov

Predicting commercially available antiviral drugs that may act on the novel coronavirus (SARS-CoV-2) through a drug-target interaction deep learning model

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Prediction of microbe-disease association from the integration of neighbor and graph with collaborative recommendation model

Network-based identification of microRNAs as potential pharmacogenomic biomarkers for anticancer drugs

LRLSHMDA: Laplacian regularized least squares for human microbe-disease association prediction

The performance of current methods in ligand-protein docking

Molecular docking: A powerful approach for structure-based drug discovery

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Mechanism of inhibition of ebola virus RNA-dependent RNA polymerase by remdesivir

Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV

Broad spectrum antiviral remdesivir inhibits human endemic and zoonotic deltacoronaviruses with a highly divergent RNA dependent RNA polymerase

Learning from the past: Possible urgent prevention and treatment options for severe acute respiratory infections caused by 2019-nCoV

Coronavirus infections: more than just the common cold

Prediction of the 2019-nCoV 3C-like protease (3CLpro) structure: Virtual screening reveals velpatasvir, ledipasvir, and other drug repurposing candidates

Reducing mortality from 2019-nCoV: Host-directed therapies should be an option

What to do next to control the 2019-nCoV epidemic

Emerging novel coronavirus (2019-nCoV)-current scenario, evolutionary perspective based on genome analysis and recent developments

FDA's approval of Veklury (remdesivir) for the treatment of COVID-19-The Science of Safety and Effectiveness

Clinical features and treatment of 2019-nCov pneumonia patients in Wuhan: Report of a couple cases

The PyMOL molecular graphics system

Pandemic potential of a strain of influenza A (H1N1): EARLY FINDINGS

Human influenza A H5N1 virus related to a highly pathogenic avian influenza virus

Human infection with a novel avian-origin influenza A (H7N9) virus

Peginterferon Alfa-2a plus ribavirin for chronic hepatitis C virus infection

Three-dimensional structure of aspartyl protease from human immunodeficiency virus HIV-1

Primary human immunodeficiency virus type 2 (HIV-2) isolates, like HIV-1 isolates, frequently use CCR5 but show promiscuity in coreceptor usage

Isolation of hendra virus from pteropid bats: A natural reservoir of Hendra virus

Analysis of the protein-coding content of the sequence of human cytomegalovirus strain AD169

Commentary: Middle east respiratory syndrome coronavirus (MERS-CoV): Announcement of the coronavirus study group

The respiratory syncytial virus vaccine landscape: Lessons from the graveyard and promising candidates

Severe acute respiratory syndrome coronavirus (SARS-CoV) infection inhibition using spike protein heptad repeat-derived peptides

Database resources of the national center for biotechnology information

MAFFT online service: Multiple sequence alignment, interactive sequence choice and visualization

DrugBank 5.0: a major update to the DrugBank database for 2018

PubMed: the bibliographic database

Extended-connectivity fingerprints

Open-Source Cheminformatics Software

Discovery and development of safe-in-man broad-spectrum antiviral agents

Random walk with restart on multiplex and heterogeneous biological networks

The protein data bank

Source codes and datasets are freely available for download at https ://githu b.com/plhhn u/VDA-RWR/.Received: 18 July 2020; Accepted: 18 December 2020

We are thankful for help from Guangyi Liu, Ming Kuang, and Longjie Liao from Hunan University of Technology, and Ruyi Dong, Lebin Liang, Qinqin Lu, and Jidong Lang from Geneis (Beijing) Co. Ltd. We would like to thank all authors of the cited references.

Authors Geng Tian and Jialiang Yang were employed by the company Geneis (Beijing) Co. Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

The online version contains supplementary material available at https ://doi. org/10.1038/s4159 8-021-83737 -5.Correspondence and requests for materials should be addressed to J.Y. or L.Z.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.