key: cord-0979566-fth8g8x4
authors: Le, Brian L.; Andreoletti, Gaia; Oskotsky, Tomiko; Vallejo-Gracia, Albert; Rosales, Romel; Yu, Katharine; Kosti, Idit; Leon, Kristoffer E.; Bunis, Daniel G.; Li, Christine; Kumar, G. Renuka; White, Kris M.; García-Sastre, Adolfo; Ott, Melanie; Sirota, Marina
title: Transcriptomics-based drug repositioning pipeline identifies therapeutic candidates for COVID-19
date: 2021-06-10
journal: Sci Rep
DOI: 10.1038/s41598-021-91625-1
sha: a5053852096d86b75c828ffb72c5f9f67627a8e3
doc_id: 979566
cord_uid: fth8g8x4

The novel SARS-CoV-2 virus emerged in December 2019 and has few effective treatments. We applied a computational drug repositioning pipeline to SARS-CoV-2 differential gene expression signatures derived from publicly available data. We utilized three independent published studies to acquire or generate lists of differentially expressed genes between control and SARS-CoV-2-infected samples. Using a rank-based pattern matching strategy based on the Kolmogorov–Smirnov Statistic, the signatures were queried against drug profiles from Connectivity Map (CMap). We validated 16 of our top predicted hits in live SARS-CoV-2 antiviral assays in either Calu-3 or 293T-ACE2 cells. Validation experiments in human cell lines showed that 11 of the 16 compounds tested to date (including clofazimine, haloperidol and others) had measurable antiviral activity against SARS-CoV-2. These initial results are encouraging as we continue to work towards a further analysis of these predicted drugs as potential therapeutics for the treatment of COVID-19.

SARS-CoV-2 has already claimed at least a million lives, has been detected in at least 40 million people, and has likely infected at least another 200 million. The spectrum of disease caused by the virus can be broad ranging from silent infection to lethal disease, with an estimated infection-fatality ratio around 1% 1 . SARS-CoV-2 infection has been shown to affect many organs of the body in addition to the lungs 2 . Three epidemiological factors increase the risk of disease severity: increasing age, decade-by-decade, after the age of 50 years; being male; and various underlying medical conditions 1 . However, even taking these factors into account, there is immense interindividual clinical variability in each demographic category considered 3 . Recently, researchers found that more than 10% of people who develop severe COVID-19 have misguided antibodies-autoantibodies-that attack the innate immune system. Another 3.5% or more of people who develop severe COVID-19 carry specific genetic mutations that impact innate immunity. Consequently, both groups lack effective innate immune responses that depend on type I interferon, demonstrating a crucial role for type I interferon in protecting cells and the body from COVID-19. Whether the type I interferon has been neutralized by autoantibodies or-because of a faulty gene-is produced in insufficient amounts or induced an inadequate antiviral response, the absence of type I IFN-mediated immune response appears to be a commonality among a subgroup of people who suffer from life-threatening COVID-19 pneumonia 3 .

While numerous efforts are underway to identify potential therapies targeting various aspects of the disease, there is a paucity of clinically proven treatments for COVID- 19 . There have been efforts to therapeutically target the hyperinflammation associated with severe COVID-19 4 , as well as to utilize previously identified antiviral www.nature.com/scientificreports/ medications 5, 6 . One of these antivirals, remdesivir, an intravenously administered RNA-dependent RNA polymerase inhibitor, showed positive preliminary results in patients with severe COVID-19 7 . In October 2020, the FDA approved remdesivir for the treatment of COVID-19 8 . Dexamethasone has also been shown to reduce the mortality rate in cases of severe COVID-19 9 . Nevertheless, the lack of treatments and the severity of the current health pandemic warrant the exploration of rapid identification methods of preventive and therapeutic strategies from every angle. The traditional paradigm of drug discovery is generally regarded as protracted and costly, taking approximately 15 years and over $1 billion to develop and bring a novel drug to market 10 . The repositioning of drugs already approved for human use mitigates the costs and risks associated with early stages of drug development, and offers shorter routes to approval for therapeutic indications. Successful examples of drug repositioning include the indication of thalidomide for severe erythema nodosum leprosum and retinoic acid for acute promyelocytic leukemia 11 . The development and availability of large-scale genomic, transcriptomic, and other molecular profiling technologies and publicly available databases, in combination with the deployment of the network concept of drug targets and the power of phenotypic screening, provide an unprecedented opportunity to advance rational drug design.

Drug repositioning is being extensively explored for COVID-19. High-throughput screening pipelines have been implemented in order to quickly test drug candidates as they are identified [12] [13] [14] [15] . In the past, our group has successfully applied a transcriptomics-based computational drug repositioning pipeline to identify novel therapeutic uses for existing drugs 16 . This pipeline leverages transcriptomic data to perform a pattern-matching search between diseases and drugs. The underlying hypothesis is that for a given disease signature consisting of a set of up and down-regulated genes, if there is a drug profile where those same sets of genes are instead downregulated and up-regulated, respectively, then that drug could be therapeutic for the disease. This method based on the Kolmogorov-Smirnov (KS) test statistic has shown promising results for a variety of different indications, including inflammatory bowel disease 17 , dermatomyositis 18 , cancer [19] [20] [21] , and preterm birth 22 .

In existing work from Xing et al. 23 , this pipeline has been used to identify potential drug hits from multiple input disease signatures derived from SARS-CoV or MERS-CoV data. The results were aggregated to obtain a consensus ranking, with 10 drugs selected for in vitro testing against SARS-CoV-2 in Vero E6 cell lines, with four drugs (bortezomib, dactolisib, alvocidib and methotrexate) showing viral inhibition 23 . However, this pipeline has not yet been applied specifically to SARS-CoV-2 infection.

A variety of different transcriptomic datasets related to SARS-CoV-2 were published in the spring of 2020. In May 2020, Blanco-Melo et al. studied the transcriptomic signature of SARS-CoV-2 in a variety of different systems, including human cell lines and a ferret model 24 . By infecting human adenocarcinomic alveolar basal epithelial cells with SARS-CoV-2 and comparing to controls, the authors generated a list of 120 differentially expressed genes. They observed two enriched pathways: one composed primarily of type-I interferon-stimulated genes (ISGs) involved in the cellular response to viral infection; and a second composed of chemokines, cytokines, and complement proteins involved in the humoral response. After infecting the cell lines, Blanco-Melo et al. did not detect either ACE2 or TMPRSS2, which are the SARS-CoV-2 receptor and SARS-CoV-2 protease, respectively 25 . However, supported viral replication was observed, thereby allowing the capture of some of the biological responses to SARS-CoV-2.

In May 2020, another study by Lamers et al. examined SARS-CoV-2 infection in human small intestinal organoids grown from primary gut epithelial stem cells 26 . The organoids were exposed to SARS-CoV-2 and grown in various conditions, including Wnt-high expansion media. Enterocytes were readily infected by the virus, and RNA sequencing revealed upregulation of cytokines and genes related to type I and III interferon responses.

A limited amount of transcriptomic data from human samples has also been published. One study detailed the transcriptional signature of bronchoalveolar lavage fluid (of which responding immune cells are often a primary component) of COVID-19 patients compared to controls 27 . Despite a limited number of samples, the results revealed inflammatory cytokine profiles in the COVID-19 cases, along with enrichments in the activation of apoptosis and the P53 signaling pathways.

On the drug side, data are available in the form of differential gene expression profiles from testing on human cells. Publicly-available versions include the Connectivity Map (CMap) 28 , which contains genome-wide testing on approximately 1300 drugs, wherein the differential profile for a drug was generated by comparing cultured cells treated with the drug to untreated control cultures.

Here, we applied our existing computational drug repositioning pipeline to identify drug profiles with significantly reversed differential gene expression compared to several diverse input signatures for SARS-CoV-2 effects on human cells. By taking into account a broader view of differentially expressed gene sets from both cell line and organoid disease models and human samples, the predictions are complementary to other drug discovery approaches. We identified 102 unique drug hits, from which 25 were identified in at least two of the signatures, several of which have been already investigated in clinical trials. We furthermore explore our findings in the context of other computational drug repurposing efforts for COVID-19. Finally, we tested 16 of our top predicted hits in live SARS-CoV-2 antiviral assays. Four of the top predicted inhibitors were tested for virus inhibition in a human lung cell line, Calu-3, infected with SARS-CoV-2 with quantitation of the secreted virus assessed by RT-qPCR assay. Thirteen predicted inhibitors (including one tested in Calu-3) were incubated with SARS-CoV-2 infected human embryonic kidney 293T cells overexpressing ACE2 (293T-ACE2) with viral replication determined using an immunofluorescence-based assay.

In this study, we applied our drug repositioning pipeline to SARS-CoV-2 differential gene expression signatures derived from publicly available RNA-seq data (Fig. 1) . The transcriptomic data were generated from distinct types of tissues, so rather than aggregating them together, we predicted therapeutics for each signature and then www.nature.com/scientificreports/ combined the results. We utilized three independent gene expression signatures (labelled "ALV", "EXP", and "BALF"), each of which consisted of lists of differentially expressed genes between SARS-CoV-2 samples and their respective controls. The ALV signature was generated from human adenocarcinomic alveolar basal epithelial cells by comparing SARS-CoV-2 infection to mock-infection conditions 24 . The EXP signature originated from a study where organoids, grown from human intestinal cells expanded in Wnt-high expansion media, were infected with SARS-CoV-2 and then compared to controls 26 . The BALF signature was from a contrast of primary human BALF samples from two COVID-19 patients versus three controls 27 . Each of these signatures was contrasted with drug profiles of differential gene expression from CMap. For each of the input signatures, we applied a significance threshold false discovery rate (FDR) < 0.05. We further applied minimum fold change thresholds in order to identify the driving genes. The ALV signature had only 120 genes, with 109 genes shared with the drug profiles; in order to maintain at least 100 genes for the pattern-matching algorithm to work with, we applied no fold-change threshold. For the EXP signature, we applied a |log 2 FC| > 2 cutoff, resulting in 125 genes for the expansion signature (108 shared with the drug profiles). For the BALF signature, we processed the raw read count data to calculate differential gene expression values. We applied a |log 2 FC| > 4 cutoff, with the BALF data yielding 1349 protein-coding genes for the lavage fluid signature (941 shared with the drug profiles). As the sample types across the three datasets were very different, we used transcriptomics-based bioinformatics approach for drug repositioning. We generated lists of statistically significant differentially expressed genes from the analysis of three published studies of SARS-CoV-2 and COVID-19. The drug repositioning computational pipeline compares the ranked differential expression of the COVID-19 disease signature with that of drug profiles from CMap. A reversal score based on the Kolmogorov-Smirnov statistic is generated for each disease-drug pair. If a drug profile significantly (FDR < 0.05) reverses the disease signature, then the drug could be therapeutic for the disease. Across all datasets, a total of 102 drugs have been identified as potentially therapeutic for COVID-19. Twenty-five drugs were identified in analyses of at least two of the three datasets. We further conducted pathways analyses and targeted analyses on the results, focusing on the 25 shared hits. Finally, we validated 16 of our top predicted hits in live SARS-CoV-2 antiviral assays. www.nature.com/scientificreports/ different fold change thresholds to identify the appropriate gene signatures to be used for drug repurposing. The gene lists for each of these signatures can be found in the supplement (Tables S1, S2, S3). We used GSEA (Gene Set Enrichment Analysis) 29, 30 to annotate enriched (FDR 0.05) Hallmark pathways from each of the input signatures ( Fig. 2A) . A number of pathways common to at least two signatures were found. Interferon alpha response and interferon beta response were upregulated in the ALV and EXP signatures. Adipogenesis and cholesterol homeostasis pathways were downregulated in the EXP and BALF signatures. KRAS signaling, and mTORC1 (mammalian target of rapamycin complex 1) signaling were enriched in all three signatures, but not in the same direction, showing the diversity of effects SARS-CoV-2 may have on human cells, and highlighting a need for utilization of diverse profiles as we do in the present study. When we look at the www.nature.com/scientificreports/ contributing genes within the three signatures ( Fig. 2B ), we found one overlapping upregulated gene-Dickkopf WNT Signaling Pathway Inhibitor 1 (DKK1). We used the publicly available single-cell RNAseq dataset GSE128033 31 composed of 13 patients (4 healthy, 3 presenting with mild COVID-19 symptoms, and 6 presenting with severe COVID-19 symptoms) to further characterize the expression of DKK1 ( Figure S1 ). Data were re-analyzed following the standard Seurat pipeline. From the analyses of the single-cell data, DKK1 is highly expressed in COVID-19 patients compared to controls, specifically in severe patients and it is expressed by epithelial cells.

After analyzing the input SARS-CoV-2 signatures, we utilized our repositioning pipeline to identify drugs with reversed profiles from CMap ( Fig. 1) . Significantly reversed drug profiles were identified for each of the signatures using a permutation approach: 30 hits from the ALV signature (Table S4) , 15 hits from the EXP signature (Table S5) , and 86 hits from the BALF signature (Table S6) . When visualizing the gene regulation of the input signatures and their respective top 15 drug hits, the overall reversal pattern can be observed ( Fig. 2C-E) . Interestingly, we found several drugs shared across datasets that significantly reversed the disease signature. For example, haloperidol, highlighted in purple in Fig. 2C -E, was shown to reverse the disease gene signature from three datasets, whereas levopropoxyphene, shown in green in Fig. 2D ,E, was observed to reverse the disease gene signature from two of the datasets. In total, our analysis identified 102 unique drug hits (Table S7) . Twenty-five drug hits reversed at least two signatures (p = 0.0334, random sampling), and four drug hits reversed three signatures (p = 0.0599, random sampling) ( Table 1 , Fig. 3A ).

We further characterized the common hits by examining their interactions with proteins in humans. We used known drug targets from DrugBank 32 and predicted additional targets using the similarity ensemble approach (SEA) 33 . We visualized the known interactions from DrugBank in a network. Figure 3B shows the connectivity across compounds highlighting both single drug genes (such as SIGMAR1 for haloperidol) and genes shared across drugs, such as ADRA2A and DRD1 for haloperidol and co-dergocrine mesilate. The proteins with the most known interactions with our list of 25 drugs included adrenergic receptors (particularly α2 adrenoreceptors), dopamine receptors, and serotonin receptors.

To confirm the validity of our approach, the inhibitory effects of 16 of our drug hits which significantly reversed multiple SARS-CoV-2 profiles were assessed in live pCC50 antiviral assays. Next, we wanted to test the potential of the predicted compounds to inhibit viral activity robustly using different human cell lines-Calu-3 Table 1 . Therapeutic hits reversing at least 2 of input SARS-CoV-2 signatures. A wide variety of drugs were identified by the analysis of multiple signatures. Drug reversal scores are normalized for each signature; drug entries marked "n.s. " were not significant for reversing that signature. www.nature.com/scientificreports/ and 293T-ACE2. The respective selection of 16 and 13 compounds for testing was based on side effect profiles and compound availability. The inhibitory effects of haloperidol, clofazimine, valproic acid, and fluticasone were evaluated in SARS-CoV-2 infected Calu-3 cells (human lung epithelial cell line), with remdesivir also tested as a positive control. From these five, remdesivir and haloperidol inhibited viral replication (Fig. 4A) , and the inhibitory effect was also observed by microscopy (Fig. 4B ). Toxicity assessments for haloperidol, clofazimine, valproic acid, and fluticasone were evaluated using viability assays (Alamarblue) in Calu-3 cells treated with each compound for 72 h (n = 1, with 2 technical replicates). No significant differences between controls and the biological replicates were detected using a non-parametric test (Kruskal-Wallis) ( Figure S6 ). Fluticasone and bacampicillin showed some toxicity in a dose-dependent manner at the highest doses tested. Haloperidol, clofazimine and valproic acid did not show significantly reduced cell viability ( Figure S6 ).

Additionally, 13 drugs (bacampicillin, ciclopirox, ciclosporin, clofazimine, dicycloverine, fludrocortisone, isoxicam, lansoprazole, metixene, myricetin, pentoxifylline, sirolimus, tretinoin) were independently assessed in a live SARS-CoV-2 antiviral assay. Remdesivir was again used as a positive control. This testing involved six serial dilutions of each drug to inhibit the replication of SARS-CoV-2 in 293T-ACE2 cells using an immunofluorescence-based antiviral assay 34 . All antiviral assays were paired with cytotoxicity assays using identical drug concentrations in uninfected human 293T-ACE2 cells. Positive control remdesivir and 10 of our predicted drugs (bacampicillin, ciclopirox, ciclosporin, clofazimine, dicycloverine, isoxicam, metixene, pentoxifylline, sirolimus, and tretinoin) showed antiviral efficacy against SARS-CoV-2, reducing viral infection by at least 50%, that was distinguishable from their cytotoxicity profile when tested in this cell line (Fig. 5) . Several inhibitors showed micromolar to sub-micromolar antiviral efficacy, including clofazimine, ciclosporin, ciclopirox, and metixene. These results not only confirm our predictive methods, but have also identified several clinically-approved drugs with potential for repurposing for the treatment of COVID-19.

Here, we used a transcriptomics-based drug repositioning pipeline to predict therapeutic drug hits for three different input SARS-CoV-2 signatures, each of which came from distinct human cell or tissue origins. We found significant overlap of the therapeutic predictions for these signatures. From 102 total drug hits, 25 drugs reversed at least two signatures (p = 0.0334) and 4 drugs reversed three signatures (p = 0.0599). The diversity of such signatures yet overlap of highlighted drugs underscores the utility of the current pipeline for identification of drugs which might be therapeutic for the diverse effects of SARS-CoV-2 infection. www.nature.com/scientificreports/ Twenty-five of our drug hits reversed at least two of the three input signatures (Fig. 3) . Notably, 14 of the 15 hits from the EXP signature were also hits for the BALF signature, despite being generated from different types of tissue. The EXP signature was generated from intestinal tissue, whereas the BALF signature was generated from constituents of the respiratory tract. Among the common hits reversing at least two of the signatures were two immunosuppressants (ciclosporin and sirolimus), an anti-inflammatory medication (isoxicam), and two steroids (fludrocortisone and fluticasone). Sirolimus (or rapamycin), an immunosuppressant and an mTOR inhibitor, is currently undergoing investigation in several clinical trials in COVID-19 patients (NCT04371640, NCT04341675, NCT04461340). Other hits currently in clinical trials for COVID-19 treatment include ciclosporin (NCT04412785, NCT04392531), niclosamide in combination with diltiazem (NCT04558021), and clofazimine in combination with interferon beta-1b (NCT04465695).

Among our four drug hits that reversed all three signatures, three drugs demonstrated in vitro antiviral efficacy-bacampicillin, clofazimine, and haloperidol with no toxicity effects ( Figure S6 ). Our group found haloperidol decreased viral growth in SARS-CoV-2 infected Calu-3 cells (Fig. 4B) in a dose-dependent manner (Fig. 4A ). Haloperidol is a psychiatric medication that is indicated for the treatment of psychotic disorders including schizophrenia and acute psychosis. By blocking dopamine (mainly D2) receptors in the brain, haloperidol eliminates dopamine neurotransmission which leads to improvement of psychotic symptoms 35 . Haloperidol can also bind to the sigma-1 and sigma-2 receptors, which are implicated in lipid remodeling and cell stress response 12 . As reported by Gordon et al. 12 , the SARS-CoV-2 proteins Nsp6 and ORF9c interact with the sigma-1 receptor and the sigma-2 receptor, respectively. Moreover, they found that haloperidol decreased viral replication in SARS-CoV-2-infected Vero E6 cells with, based on their reported pIC50 and pCC50 values, a calculated Selectivity Index (SI) of 53.7 12 . An SI greater than 10 is the generally accepted minimum cut-off for an antiviral worth pursuing 36, 37 . In another more recent study, Gordon et al. found in their analysis of a national electronic medical record database that fewer hospitalized COVID-19 patients who were newly prescribed haloperidol and other sigma-binding typical antipsychotic medications progressed to requiring mechanical ventilation compared to those who were newly prescribed atypical antipsychotic medications that do not bind to sigma receptors 14 .

Our testing of clofazimine demonstrated submicromolar antiviral effects of this drug in SARS-Co-V-2 infected 293T-ACE2 and Vero E6 cells (Figs. 4 and S3) . Clofazimine is an orally administered antimycobacterial drug used in the treatment of leprosy. By preferentially binding to mycobacterial DNA, clofazimine disrupts the cell cycle and eventually kills the bacterium 38 . In addition to being an antimycobacterial agent, clofazimine also possesses anti-inflammatory properties primarily by inhibiting T lymphocyte activation and proliferation 39 . Yuan et al. found that clofazimine inhibits SARS-CoV-2 replication by interfering with spike-mediated viral entry and viral RNA replication 40 . Their work also demonstrated that clofazimine has antiviral efficacy against SARS-CoV-2 www.nature.com/scientificreports/ in human embryonic stem cell-derived cardiomyocytes and in an ex vivo human lung culture system, as well as antiviral synergy with remdesivir demonstrating the potential of clofazimine as part of a combination treatment regimen for COVID-19 40 . www.nature.com/scientificreports/ Our group found bacampicillin to have micromolar antiviral efficacy in SARS-Co-V-2 infected 293T-ACE2 cells. Bacampicillin is an orally administered prodrug of ampicillin typically prescribed for treating bacterial infections 41 . As identified by SPOKE 42 , bacampicillin was found to downregulate the GDF15 gene and upregulate the NFKB2 (Nuclear Factor Kappa B Subunit 2) gene in studies by Cmap 28 and LINCS 43 . The GDF15 protein acts as a cytokine and is involved in stress response after cellular injury, and the NFKB2 is a central activator of genes involved with inflammation and immune function 44 . Circulating levels of GDF15 have been found to be significantly higher in COVID-19 patients who die 45 . Zhou et al. 's work revealed NF-kappa B signaling as one of the main pathways of coronavirus infections in humans. While the rapid conversion of bacampicillin to ampicillin in vivo makes this prodrug a less optimal therapeutic candidate for COVID-19, our findings nevertheless provide insights into the immunologic and inflammatory landscape from SARS-CoV-2 infection.

Overall, in testing of our drug hits across two human cell line assays, 11 of 16 exhibited inhibition of SARS-CoV-2 infection. In particular, three of our four consensus drug hits demonstrated antiviral efficacy, with haloperidol showing reproducible inhibition in Calu-3 cells, and bacampicillin and clofazimine inhibiting viral activity in 293T-ACE2 cells without cytotoxicity. Many of our tested drugs can be administered orally, and several are on the WHO Model List of Essential medications, including ciclosporin, clofazimine, and haloperidol 46 . These results suggest that our drug repositioning pipeline can rapidly identify readily available potential therapeutics in antiviral contexts.

There are several limitations of our approach that should be recognized. In general, drug repositioning pipelines are reliant upon data being encoded in a computable format. Considering variability of conditions under which experimental and clinical data are obtained, this also implies a particular set of limitations and biases to the ensuing results. The general basis of the pipeline operates under the hypothesis that reversing the set of transcriptomic changes will correlate with antiviral response. While it is hard to say whether reversing expression will stymie viral replication, we hypothesize these drugs could treat the symptoms and/or assist in withstanding the infection while the immune response fights off the virus. It is likely that secondary cellular responses influence a specific signature including pro-viral as well as antiviral response and should be further explored in future studies.

In our work, we applied the KS-based similarity metric on the CMap database. Rank-based methods, such as the KS statistic, may suffer from high false positive rates, as genes not differentially expressed can be ranked high and contribute to the similarity measurement for drug and disease signatures 47, 48 . Moreover, rank-based approaches can also miss many potential drugs as ranking captured just a small part of information in a gene expression profile. Alternative methods have been proposed, such as EMUDRA 48 and XSum 49 ; however, they have not been widely adopted by the community. Future work might include evaluating multiple similarity metrics on larger datasets 48 . The studies that we leveraged here are also limited because of their small sample size, which might explain the small gene overlap across the signatures. While the different sample types were chosen due to availability and to capture the heterogeneity of the response to viral infection, data generated from cell lines and organoids (the ALV and EXP signatures, respectively) might not accurately represent the biological changes and responses in human infection, especially with regards to the most severe respiratory infections caused by the SARS-CoV-2 virus. Furthermore, the data sources presented here offer no insight into the temporal aspects of response to viral infection. In vivo experiments are needed to better understand the biology of the disease and the effect of these drugs on it. Moreover, although the BALF signature was generated from fluid recovered from lavage of infected human tissues, this primary response data was aggregated from a very limited sample size (2 cases and 3 controls). As the sample types across the three datasets were very different, we used different fold change thresholds to identify the signature genes to be used for drug repurposing. Future work leveraging additional datasets can provide further insights into therapeutic predictions and drug repositioning opportunities. Gathering samples from a larger number of patients should generate a more robust gene expression signature and better inform therapeutic predictions. Furthermore, the drug profiles from CMap were generated from cell line data; drug data generated from more relevant tissue cultures (e.g. lung tissue) may generate more appropriate comparisons. Finally, our validation approaches focus on in vitro studies, which are limited. While a reduction of 50-70% in viral RNA levels meets the generally accepted minimum threshold change of 1 log, next steps would involve confirming the antiviral effects via established virological methods (e.g. plaque assay) in addition to further in vivo testing of the proposed compounds.

The drug development response for SARS-CoV-2/COVID-19 is rapidly developing. One drug, remdesivir, recently received FDA approval for the treatment of COVID-19, and numerous other drugs are being actively explored for possible therapeutic value in COVID-19 cases. Utilizing a diverse set of transcriptomic SARS-CoV-2 signatures, our drug repositioning pipeline identified 25 therapeutic candidates. Validation experiments revealed antiviral activity for 11 of 16 drug hits. Further clinical investigation into these drug hits, in vivo assays as well as potential combination therapies is warranted to further investigate both the antiviral as well as side effect profile of the drugs.

Study design. We have previously developed and used a transcriptomics based bioinformatics approach for drug repositioning in various contexts including inflammatory bowel disease, dermatomyositis, and spontaneous preterm birth. For a list of differentially expressed genes, the computational pipeline compares the ranked differential expression of a disease signature with that of a profile 16, 19, 28 . A reversal score based on the Kolmogorov-Smirnov statistic is generated for each disease-drug pair, with the idea that if the drug profile significantly reverses the disease signature, then the drug could be potentially therapeutic for the disease. www.nature.com/scientificreports/ propagated from Vero E6 cells (GSE147507, 67 samples) 24 . Due to the fast-moving nature of the research topic, we opted to use this cell line data in lieu of waiting for substantial patient-level data. This work identified 120 differentially expressed genes (DEGs)-100 upregulated and 20 downregulated. We used these 120 genes as the ALV signature for our computational pipeline (Table S1 ). Lamers et al. performed RNA-seq on their organoid samples, from which differentially expressed genes were calculated 26 . These samples were grown in a medium with a Wnt surrogate supplement and infected with SARS-CoV-2 propagated from Vero E6 cells (GSE149312, 22 samples). They detected 434 significant DEGs (FDR < 0.05). We additionally applied a fold-change cutoff (|log 2 FC| > 2), resulting in 125 genes used as the EXP signature (Table S2) .

Xiong et al. performed RNA-seq analysis of BALF samples from two COVID-19 patients (two samples per patient) and three healthy controls 27 . We processed their raw read counts in order to construct a differential signature (see below for details). FASTQ files were downloaded from the Genome Sequence Archive 50,51 under accession number CRA002390. Paired-end reads were mapped to the hg19 human reference genome using Salmon (v.1.2.0) and assigned Ensembl genes. After read quality control, we obtained quantifications for 55,640 genes in all samples. In order to identify genes differentially expressed between cases and controls for the BALF samples, we quantified gene expression as raw counts. Raw counts were used as inputs to DESeq2 (v.1.24.0 R package) to call differentially expressed genes (DEGs). After adjusting for the sequencing platform, the default settings of DESeq2 were used. Principal components were generated using the DESeq2 function ( Figure S2 ), and heat maps were generated using the Bioconductor package pheatmap (v.1.0.12) using the rlog-transformed counts ( Figure S3 ). Values shown are rlog-transformed and row-normalized. Volcano plots were generated using the Bioconductor package EnhancedVolcano (v.1.2.0) ( Figure S4 ). Retaining only protein-coding genes and applying both a significance threshold and a fold-change cutoff (FDR < 0.05, |log 2 FC| > 4), we obtained 1349 genes to be used as the BALF signature (Table S3) .

Functional enrichment gene-set analysis for GSEA (Gene Set Enrichment Analysis) was performed using fgsea (v.1.12.0 R package) and the input gene lists were ranked by log2 fold change. The 50 Hallmark Gene Sets used in the GSEA analysis were downloaded from MSigDB Signatures database 29, 52 . For GO (Gene Ontology) terms, identification of enriched biological themes was performed using the DAVID database 53 .

Drug gene expression profiles. Drug gene expression profiles were sourced from Connectivity Map (CMap), a publicly-available database of drugs tested on cancer cell lines 28 . CMap contains a set of differential gene expression profiles generated from treating cultured human cells with a variety of different drugs and experimental compounds. These profiles were generated using DNA microarrays to assay mRNA expression. These drug profiles are ranked genome-wide profiles (~ 22,000 genes) of the effects of the drugs on various cell lines. 6100 gene expression profiles are presented in CMap. A total of 1309 compounds were tested in up to 5 different cell lines. The overlap between the gene lists of CMap and the SARS-CoV-2 signature is 109 genes.

Computational gene expression reversal scoring. To compute reversal scores, we used a non-parametric rank-based method similar to the Kolmogorov-Smirnov test statistic. This analysis was originally suggested by the creators of the CMap database and has since been implemented in a variety of different settings [16] [17] [18] [19] 22, 28 . As also described by others, the drug signature is compared with the gene expression profiles. By splitting the gene signature into two lists containing only upregulated genes and downregulated genes, a so-called connectivity score is estimated via several auxiliary variables using a nonparametric rank-ordered Kolmogorov-Smirnov (KS) test) 47 . Similar to past works, we applied a pre-filtering step to the CMap profiles to maintain only drug profiles which were significantly correlated with another profile of the same drug. Drugs were assigned reversal scores based on their ranked differential gene expression profile relative to the SARS-CoV-2 ranked differential gene expression signature. A negative reversal score indicated that the drug had a profile which reversed the SARS-CoV-2 signature; that is, up-regulated genes in the SARS-CoV-2 signature were down-regulated in the drug profile and vice versa.

Statistical analysis. P-values were adjusted using the false discovery rate (FDR; Benjamini-Hochberg) procedure. P-values for individual drug hits were obtained by comparing reversal scores to a distribution of random scores. Negative reversal scores were considered significant if they met the criterion FDR < 0.05. For drugs tested multiple times (e.g. different cell lines), we used the most reversed profile (lowest negative score).

For significance values of the number of drugs reversing multiple signatures, we constructed distributions of the common reversal (reversing two of three signatures) and the consensus reversal (reversing three of three signatures) by randomly sampling the same number of drug profiles for each signature from CMap.

Seurat object for each sample was generated using Seurat v.3. While the data has been filtered by 10x's algorithm, we still needed to ensure the remaining cells are clean and devoid of artifacts. We calculated three confounders for the dataset: mitochondrial percentage, ribosomal percentage, and cell cycle state information. For each sample, cells were normalized for genes expressed per cell and per total expression, then multiplied by a scale factor of 10,000 and log-transformed. Low quality cells were excluded from our analyses-this was achieved by filtering out cells with greater than 5000 and fewer than 300 genes and cells with high percentage of mitochondrial and ribosomal genes (greater than 10% for mitochondrial genes, and 50% for ribosomal genes). SCTransform is a relatively new technique that uses "Pearson Residuals" (PR) to normalize the data. PRs are independent of www.nature.com/scientificreports/ sequencing depth 54 . We "regress out" the effects of mitochondrial and ribosomal genes, and the cell cycling state of each cell, so they do not dominate the downstream signal used for clustering and differential expression. We then performed a lineage auto-update disabled r dimensional reduction (RunPCA function). Then, each sample was merged together into one Seurat object. Data were then re-normalized and dimensionality reduction and significant principal components were used for downstream graph-based, semi-unsupervised clustering into distinct populations (FindClusters function) and uniform manifold approximation and projection (UMAP) dimensionality reduction was used. For clustering, the resolution parameter was approximated based on the number of cells according to Seurat guidelines; a vector of resolution parameters was passed to the FindClusters function and the optimal resolution of 0.8 that established discernible clusters with distinct marker gene expression was selected. We obtained a total of 21 clusters representing the major immune and epithelial cell populations. To identify marker genes driving each cluster, the clusters were compared pairwise for differential gene expression (FindAllMarkers function) using the Likelihood ratio test assuming an underlying negative binomial distribution (negbinom Infection experiments. Work involving live SARS-CoV-2 was performed in the BSL3 facility at the Gladstone Institutes with appropriate approvals. Calu-3 cells were seeded in 96-well plates for 24 h, infected with SARS-CoV-2 at a multiplicity of infection (MOI) of 0.05, and treated with compounds. 72 h post infection, supernatant was collected for RNA extraction and the RNA was analyzed using RT-qPCR to quantify viral genomes present in the supernatant. SARS-CoV-2 specific primers targeting the E gene region: 5ʹ-ACA GGT ACG TTA ATA GTT AAT AGC GT-3ʹ (Forward) and 5ʹ-ATA TTG CAG CAG TAC GCA CACA-3ʹ (Reverse) were used to quantify cDNA on the 7500 Fast Real-Time PCR system (Applied Biosystems). Cells were fixed with paraformaldehyde and used for immunofluorescence analysis with dsRNA antibody (SCICONS) and DAPI stain. Images were acquired and analyzed using ImageXpress Micro Confocal High-Content Imaging System.

In vitro microneutralization assay for SARS-CoV-2 serology and drug screening. For studies at Mount Sinai, SARS-CoV-2 was propagated in Vero E6 cells (ATCC CRL-1586) and 293T-ACE2 cells (ATCC CRL-3216), as previously described 12, 34 . Two thousand cells were seeded into 96-well plates in DMEM (10% FBS) and incubated for 24 h at 37 °C, 5% CO 2 . Then, 2 h before infection, the medium was replaced with 100 μl of DMEM (2% FBS) containing the compound of interest at concentrations 50% greater than those indicated, including a DMSO control. The Vero E6 cell line used in this study is a kidney cell line; therefore, we cannot exclude that lung cells yield different results for some inhibitors. Plates were then transferred into the Biosafety Level 3 (BSL3) facility and 100 PFU (MOI = 0.025) was added in 50 μl of DMEM (2% FBS), bringing the final compound concentration to those indicated. Plates were then incubated for 48 h at 37 °C. After infection, supernatants were removed and cells were fixed with 4% formaldehyde for 24 h before being removed from the BSL3 facility. The cells were then immunostained for the viral NP protein (an in-house mAb 1C7, provided by Dr. Thomas Moran) with a DAPI counterstain. Infected cells (488 nM) and total cells (DAPI) were quantified using the Celigo (Nexcelcom) imaging cytometer. Infectivity is measured by the accumulation of viral NP protein in the nucleus of the Vero E6 cells and 293T-ACE2 cells (fluorescence accumulation). Percentage infection was quantified as ((infected cells/total cells) − background) × 100 and the DMSO control was then set to 100% infection for analysis. The IC50 and IC90 for each experiment were determined using the Prism (GraphPad) software. Cytotoxicity was also performed using the MTT assay (Roche), according to the manufacturer's instructions. Cytotoxicity was performed in uninfected VeroE6 cells with same compound dilutions and concurrent with viral replication assay. All assays were performed in biologically independent triplicates.

The data used for the repositioning pipeline are all publicly available. The code for the drug repositioning pipeline was adapted from reference 19 and is available at https:// github. com/ brian lle/ sirota_ lab_ covid_ drug_ repos ition ing.

Predicted COVID-19 fatality rates based on age, sex, comorbidities and health system capacity

Extrapulmonary manifestations of COVID-19

Inborn errors of type I IFN immunity in patients with life-threatening COVID-19

COVID-19: Consider cytokine storm syndromes and immunosuppression

A protein interaction map identifies existing drugs targeting SARS-CoV-2

In silico discovery of candidate drugs against Covid-19

Remdesivir for the treatment of Covid-19-Preliminary report

FDA Approves First Treatment for COVID-19. FDA

Dexamethasone in hospitalized patients with covid-19-Preliminary report

The price of innovation: New estimates of drug development costs

Old drugs-new uses

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2

Comparative host-coronavirus protein interaction networks reveal pan-viral disease mechanisms

In silico studies on therapeutic agents for COVID-19: Drug repurposing approach

Discovery and preclinical validation of drug indications using compendia of public gene expression data

Computational repositioning of the anticonvulsant topiramate for inflammatory Bowel disease

Identification of alpha-adrenergic agonists as potential therapeutic agents for dermatomyositis through drug-repurposing using public expression datasets

Computational discovery of niclosamide ethanolamine, a repurposed drug candidate that reduces growth of hepatocellular carcinoma cells in vitro and in mice by inhibiting cell division cycle 37 signaling

Reversal of cancer gene expression correlates with drug efficacy and reveals therapeutic targets

A drug repositioning approach identifies tricyclic antidepressants as inhibitors of small cell lung cancer and other neuroendocrine tumors

Computational discovery of therapeutic candidates for preventing preterm birth

Analysis of infected host gene expression reveals repurposed drug candidates and time-dependent host response dynamics for COVID-19

Imbalanced host response to SARS-CoV-2 drives development of COVID-19

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

SARS-CoV-2 productively infects human gut enterocytes

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

The Connectivity Map: Using gene-expression signatures to connect small molecules, genes, and disease

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

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

Relating protein pharmacology by ligand chemistry

An in vitro microneutralization assay for SARS-CoV-2 serology and drug screening

More dopamine, more D2 receptors

Novel antiviral activity of mung bean sprouts against respiratory syncytial virus and herpes simplex virus-1: An in vitro study on virally infected Vero and MRC-5 cell lines

Antiviral activities of extremophilic actinomycetes extracts from Kazakhstan's unique ecosystems against influenza viruses and paramyxoviruses

Clofazimine: Current status and future prospects

Clofazimine broadly inhibits coronaviruses including SARS-CoV-2

Systematic integration of biomedical knowledge prioritizes drugs for repurposing

A Next Generation Connectivity Map: L1000 platform and the first 1,000,000 profiles

Entrez gene: Gene-centered information at NCBI

Circulating levels of GDF-15 and calprotectin for prediction of in-hospital mortality in COVID-19 patients: A case series

WHO | WHO Model Lists of Essential Medicines. WHO

A review of connectivity map and computational approaches in pharmacogenomics

EMUDRA: Ensemble of multiple drug repositioning approaches to improve prediction accuracy

Systematic evaluation of connectivity map for disease indications

Database Resources of the National Genomics Data Center in 2020

Genome sequence archive

Molecular signatures database (MSigDB) 3.0

Database for annotation, visualization, and integrated discovery

Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression

SARS-CoV-2 Orf6 hijacks Nup98 to block STAT nuclear import and antagonize interferon signaling

This work was funded in part by the Program for Breakthrough Biomedical Research (PBBR) Grant. This research was also partly funded by the Defense Advanced Research Projects Agency (HR0011-19-2-0020); by CRIP (Center for Research for Influenza Pathogenesis), a NIAID supported Center of Excellence for Influenza Research and Surveillance (CEIRS, contract # HHSN272201400008C); by supplements to NIAID grant U19AI135972 and DoD grant W81XWH-20-1-0270; by the generous support of the JPB Foundation and the Open Philanthropy Project (research grant 2020-215611 (5384)); and by anonymous donors to A. García-Sastre. M. Ott acknowledges support through a gift from the Roddenberry Foundation and by NIH 5DP1DA038043. We would like to thank members of the Sirota Lab for useful discussion, Dr. Boris Oskotsky for technical support, Edna Rodas for administrative support, Randy Albrecht for support with the BSL3 facility and procedures at the ISMMS, and Richard Cadagan for technical assistance.

B.L., T.O. and M.S. designed and coordinated the study. B.L. led the drug repurposing efforts. G.A., K.Y., I.K., and C.L. helped with data analyses. SARS-CoV-2 virus assays were led by A.V.G., G.R.K., K.L., R.R., K.W., A.G.S., and M.O. All the authors contributed to making figures, writing and editing the manuscript.

M.S. is on the advisory board of twoXAR. The García-Sastre Laboratory has received research support from Pfizer, Senhwa Biosciences and 7Hills Pharma. A.G.S. has consulting agreements for the following companies involving cash and/or stock: Vivaldi Biosciences, Contrafect, 7Hills Pharma, Avimex, Vaxalto, Accurius and Esperovax. Other authors declare no competing financial interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-91625-1.Correspondence and requests for materials should be addressed to M.S.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.