CeDR Atlas: a knowledgebase of cellular drug response.

Drug response to many diseases varies dramatically due to the complex genomics and functional features and contexts. Cellular diversity of human tissues, especially tumors, is one of the major contributing factors to the different drug response in different samples. With the accumulation of single-cell RNA sequencing (scRNA-seq) data, it is now possible to study the drug response to different treatments at the single cell resolution. Here, we present CeDR Atlas (available at https://ngdc.cncb.ac.cn/cedr), a knowledgebase reporting computational inference of cellular drug response for hundreds of cell types from various tissues. We took advantage of the high-throughput profiling of drug-induced gene expression available through the Connectivity Map resource (CMap) as well as hundreds of scRNA-seq data covering cells from a wide variety of organs/tissues, diseases, and conditions. Currently, CeDR maintains the results for more than 582 single cell data objects for human, mouse and cell lines, including about 140 phenotypes and 1250 tissue-cell combination types. All the results can be explored and searched by keywords for drugs, cell types, tissues, diseases, and signature genes. Overall, CeDR fine maps drug response at cellular resolution and sheds lights on the design of combinatorial treatments, drug resistance and even drug side effects.

INTRODUCTION

Drug response to many diseases varies dramatically due to the complex genomics and functional features and contexts (1–3). As the majority of cellular systems are heterogeneous (4), the bulk RNA sequencing (RNA-seq) technologies may only measure the average expression level that includes a diverse collection of cells concealing cell type specific signatures with different levels of response to treatment (5). Much of the diversity in individual drug response can be attributed to cellular heterogeneity including cell types, proportions, cell-cell communications, and temporal-spatial distributions and compositions of cells in tissues that are critical to diseases (6,7). The inability to efficiently measure transcriptional sensitivities across diverse cell contexts has limited our understanding of how drug response differs across genomic and molecular cell states, which could be critical for predicting the therapeutic outcome of patient tumors (8). The recent advance in single-cell RNA sequencing (scRNA-seq) technologies has made it possible to measure the transcriptomes at the single cell resolution (9,10). Over the past few years, nearly a thousand studies have been conducted using scRNA-seq technologies to study the cell types and related features in bulk tissues. However, a gap has been existing to both experimentally and computationally infer cellular drug response, which would provide references for drug repurposing, drug combination design, and new therapeutic development (11–16).

So far, many efforts have been made to develop computational approaches to predict drug response. The pioneer work led by the Connectivity Map (CMap) project has examined 1309 drugs with different doses in five cell lines to quantitatively measure drug induced gene expression profiles (17). With such a rich resource, many studies implemented the concept of ‘anti-correlation between drug and disease signatures’ for various drug related researches such as drug repurposing and combination (18–20). The underlying assumption is that, for a query transcriptome, candidate drugs can be prioritized if such drugs have drug-induced perturbations (as measured by gene expression) in the opposite direction as the query transcriptome. With the rapid accumulation of transcriptomics data, this method has been adapted and applied in various cases. For example, in a study searching for candidate drugs that can be repurposed for schizophrenia, So et al. searched for drugs that induced ‘anti-correlation’ gene signatures compared to the GWAS-imputed gene expression and interestingly, the set of candidate drugs they identified was significantly enriched with antipsychotics (21). Such an approach has also been implemented to identify candidate drugs in cancer. However, due to the substantial intra- and inter-tumor heterogeneity, advanced procedures are in pressing need to incorporate changes in cell types, cell states, cell–cell communications and proportions and accurately infer drug responses in the precision medicine practice.

To this end, we present CeDR Atlas, a knowledgebase reporting computational inference of cellular drug response for hundreds of cell types from the major organs and tissues in human and mouse. We took advantage of the high-throughput profiling of drug-induced gene expression available through the CMap project as well as hundreds of scRNA-seq data representing cells from a wide variety of organs/tissues, diseases, and conditions. For each cell type, we conducted anti-correlation screening and identified drugs that induced gene expression with the opposite direction compared to the cell transcriptome. CeDR is augmented with user-friendly interfaces and functions, allowing users to easily explore the cell types with associated drugs as well as the corresponding gene signatures. CeDR provides direct references for cellular drug response profiles including not only disease cell types but also normal cell types. The comprehensive cell-type specific drug sets can be used for design of combinatory treatments and identification of drug resistance and even drug side effects.

MATERIALS AND METHODS

Data collection

We conducted a comprehensive literature-mining and collected single-cell transcriptome data from hundreds of individual studies (Figure 1A). All the datasets we used in this study were derived from public resources and the detail information for each dataset can be found in Supplementary Data. In the current version, CeDR has included the following data.

Figure 1.

Schematic overview of CeDR. (A) Data curation and collection from different sources. (B) Data quality control and preprocessing. (C) Illustration of the algorithm to infer cellular drug response.

Human cell landscape

The raw expression matrix was downloaded from http://bis.zju.edu.cn/HCL/ (22) consisting of more than 720 000 single cells with cell type annotations from about 50 human tissues. This dataset presents a comprehensive annotation for major human normal organs/tissues.

Gene expression nebulas (GEN)

We collected the human scRNA-seq datasets from the GEN database (https://ngdc.cncb.ac.cn/gen/) (23) generated using the 10× Genomics platform. As a result, a total of 42 projects consisting of 294 high-quality datasets were derived and curated. All of these datasets were processed by using a one-stop analysis pipeline implemented by the CellRanger (24) (v3.1.0) software for quality control, sample de-multiplexing, barcode processing, and generation of feature-barcode matrices. Cell type annotation at the single cell level was performed by using the R package SingleR (25) based on five transcriptomic datasets as built-in references.

Manually curated human scRNA-seq data

We collected publicly available human scRNA-seq data as well as their metadata and cell type annotations by searching PubMed (26) and public databases, such as Gene Expression Omnibus (GEO) (27), Single Cell Portal of Broad Institute (https://singlecell.broadinstitute.org/single_cell) and ArrayExpress Archive of Functional Genomics Data (https://www.ebi.ac.uk/arrayexpress/) (28). As a result, 51 projects with 2 421 320 cells were collected containing data for 55 tissues and 1313 tissue-cell types.

Mouse cell atlas

The datasets in Mouse Cell Atlas (29) were downloaded from GEO with accession ID GSE108097. It contains more than 1 200 000 single cells from about 28 tissues and 170 tissue-cell type clusters.

Tabula muris data

The Tabula Muris (30) mouse atlas expression data and cell annotation files were downloaded from (https://tabula-muris.ds.czbiohub.org/). This dataset contains nearly 100 000 cells from 20 organs and tissues generated by either the FACS-based or droplet-based scRNA-seq technologies.

Manually curated mouse scRNA-seq data

We also manually curated 16 mouse projects with 655 950 cells from GEO, which finally resulted in 25 tissues and 745 tissue-cell types.

Cancer cell line single-cell dataset

The cancer cell line scRNA-seq datasets were downloaded from Single Cell Portal of Broad Institute with the accession ID SCP542 (31). It contains 19 cancer cell lines, including 195 cancer-cell types and 53 298 cells.

Quality control and preprocessing

In the current version of our database, we mainly collected the expression data that have well-annotated cell types. Datasets with unavailable annotations were discarded from further analyses. All datasets were preprocessed using the Python package Scanpy (32) and recorded as an AnnData object in Python. Genes detected in less than five cells were filtered out. Cells with a mitochondrial gene proportion (mtDNA%, the fraction of mitochondrial transcript counts of the total transcript counts) greater than 5% were further filtered out due to the low-quality. Raw count expression matrices were subsequently normalized using the normalize_per_total function and log transformed. For each project, we split it into sub-objects based on the phenotype or disease status. For each single-cell expression sub-object, we kept only the cell types that contained at least 20 cells and further removed the cell types that had less than 500 expressed genes (Figure 1B). In summary, we processed about 120 projects, resulting in 582 data objects, 140 phenotypes and 1250 tissue-specific cell types for human, mouse and cell lines. More detailed information for all studies can be found in Supplementary Table SI.

Drug induced gene expression data

To systematically analyze the cellular drug sensitivities across different tissues and perturbations, we assembled drug induced gene signatures from the CMap database (version: build 02) (17) for 1309 drugs with different doses originally generated using the Affymetrix Gene Chips. We downloaded the rank matrix based on the differentially expressed gene analyses of drug treatments (drug-treated versus no drug-treated). The Affymetrix (33) probe set identifiers were mapped with gene symbols for the following analyses.

Cell type-drug response analysis (CeDR)

To predict cellular drug response, we implemented the anti-correlation screening procedure. Specifically, for each query transcriptome (which can be from bulk tissues or single cells), candidate drugs can be prioritized if they have drug-induced perturbations (as measured by gene expression) in the opposite direction as the query transcriptome. To this end, for each cell type, we first rank the expressed genes (defined as those with non-zero average expression values) according to their average expression across all cells for the same type. The k most highly expressed genes and the k most lowly expressed genes were then selected as the signature gene set for the cell type (Figure 1C, ‘within cell type’). Next, we conducted differentially expressed gene (DEG) analysis for each cell type as compared to other cell types, which is a standard procedure in scRNA-seq analysis. DEGs identified in this way are deemed as cell-type specific genes. We similarly selected the k most highly specifically expressed genes and the k most lowly specifically expressed genes as the signature gene set for the ‘across cell type’ information (Figure 1C). The ‘within cell type’ gene set and the ‘across cell type’ gene set were combined to define the cellular signature genes. The parameter k can be defined manually. In this work, we used k = 300 and the cellular signature genes contain a total of 1200 genes.

For each compound, we also defined a signature gene set based on their rank as provided by the CMap data, i.e. the top k and bottom k genes. In this case, we selected k = 600 to match the size of cellular signature genes. Following the concept of ‘anti-correlation’, we next constructed contingency tables using the down-regulated signature genes of each drug and the cellular signature genes that were highly (or highly specifically) expressed, followed by the chi-square test for drug and cell type association test (denoted as P-value 1). Similarly, we constructed contingency tables using the up-regulated signature genes of each drug and the cellular signature genes that were lowly expressed or lowly specifically expressed, followed by the chi-square test (denoted as P-value 2). The two contingency tables should be constructed separately to ensure that the chi-square test was conducted purposely to identify the anti-correlation relationship. Moreover, we required that the expression of the overlapping genes from the two signature sets of the drug and the cell type should also be significantly anti-correlated, which is examined by Spearman correlation coefficient.

For each data object, to provide a relatively comprehensive resource, we provide in CeDR the drug and cell type associations with nominal significance (P-value 1 < 0.05 and P-value 2 < 0.05). In addition, we present the enrichment and anti-correlation P-values for each associated pair, as well as the signature genes. The functional enrichment analysis for the corresponding signature genes were performed by GSEA (34).

Description of the website and tools

CeDR aims to provide references for cell type and drug associations across different tissues and species. Table 1 summarized the data sets and results deposited in CeDR. Currently, CeDR organized all data into 109 projects, including 93 from human, 15 from mouse, and one from cell lines. These projects contained a total of 582 datasets (460 for human, 102 for mouse, and 20 for cell lines), providing information for 140 phenotypes (human: 93, mouse: 27, cell line: 20) and 1250 cell types (human: 684, mouse: 370, cell line: 196). As aforementioned, we selected associations with nominal significance based on both the chi-square test and the Spearman correlation analysis. This resulted in a total of 241 116 associations in CeDR.

Table 1.

Data summary of CeDR

Source	Projects	Subjects	Phenotype	Tissue-cell type	Association
Human	93	460	93	684	188 157
Mouse	15	102	27	370	42 660
CCLE	1	20	20	196	10 299
Total	109	582	140	1250	241 116

Web interface

As shown in Figure 2, the CeDR interface allows users to intuitively browse and search for any data in the database. The home page provides an overview of phenotypes across different tissues and a quick search function for users to query the database for species, tissues, cell types, or phenotypes conveniently (Figure 2A and B). Datasets have been classified by tissues where users can select relevant dataset in a pop-up window to retrieve the corresponding phenotypes.

Figure 2.

Screenshots of web pages for CeDR. (A) The web images in the home page allow users to switch sources and tissues with detailed data summary information and phenotypes. (B) Quick search functions in the home page. (C) The browse page allows users to browse and search all the project, tissues, phenotypes and other detailed information. (D) Summary of the selected dataset and cell type-drug association network. (E) Detailed results for each cell type with associated drugs and corresponding signature genes.

Users can also navigate the whole database through ‘Browse’ functions. The ‘Browse’ page displays a general table of all sub-objects and provides data advanced search functions (Figure 2C). The visible part of the table row contains the major meta information of the data, including Source, Project ID, Tissue, Phenotype, Cell type, Drug and other information specific to this sub-object. Users can select a sub dataset of interest and clicking the corresponding ‘Dataset ID’ button will enable the browse of detailed results for tissue specific cell type-drug associations. In particular, cell types-drug associations with significant P-values will be returned. The result page contains an intuitive table with more detailed information about the sub-object. Moreover, the visualization of cell types, fraction and predicted association network will be displayed (Figure 2D). Moving the mouse to a point in the UMap diagram, and the coordinate of the cell-type information will be displayed. Clicking the ‘Detail’ button for each association, the gene signature with GSEA prerank result in drug and single cell expression data will be further displayed (Figure 2E).

We also provide multiple search function pages for users to identify cell type, drug, disease or tissue of interest. For a user-input query string, we will search both the short names and the full names of interest. In addition, users can download all data via the ‘Download’ page. We also built a detailed tutorial for the usage in the documentation page to briefly describe the data collection, preprocessing, and analysis (https://ngdc.cncb.ac.cn/cedr/documentation).

Application of CeDR: pancreas as an example

CeDR provides direct references for cellular drug response profiles including not only disease cell types but also normal cell types which can be further used for design of combinatory treatments, drug resistance and even drug side effects. As shown in Figure 3A, treatment with a single drug on normal cells may lead to side effects; on the other hand, treatment on disease cells may not be always as efficient as expected due to the sample heterogeneity. Investigating drug response at single-cell level can provide insight into combination treatment which can maximize the efficiency by killing disease cells and minimize the side effect to the normal cells.

Figure 3.

Drug response at the single cell Level: Pancreas as an example. (A) Tissue subpopulations maintain diverse response to different drugs. (B) UMap visualization for major cell types across different datasets. (C) Cell type-drug associated network. The circular node represents cell types and the hexagon denotes associated drugs. Different interactions refer to treatments or side effects and the shape of the edge denotes the corresponding phenotypes. Here, only the top 10 drugs with corresponding cell types were shown in (C). (D) Functional enrichment analysis for signature genes referring to acinar-dapsone association. (E) Matrix plot of signature genes in pancreas cancer scRNA-seq dataset referring to buetulin-macrophage association.

Here, we select human pancreas as an example to indicate the application of CeDR in drug treatment. Pancreas is part of the human endocrine system and mainly produces hormones and digestive enzymes. To study the major cell type in pancreas, we integrated nine scRNA-seq data sets from different studies using our recently developed method INSCT which based on batch-aware triplet neural networks (35,36). These data sets included four normal pancreas tissues, three type I diabetes, one type II diabetes and one pancreatic cancer dataset and were expected with strong batch effects. By applying INSCT, we successfully corrected the potential batch effects and reported the primary cell types in pancreas such as acinar, duct, alpha, beta, delta, gamma, and epsilon cells (37,38) (Figure 3B). With CeDR, we screened the compounds in CMap across all the cell types based on the signature genes defined above. Among our results, some of the top ranked drugs have already been reported sensitive to pancreas. To better visualize the results, we construct a cell type-drug association network (Figure 3C).

In normal cells, we particularly examined the reversed activity. Our results showed that the drug dapsone had adverse effects with normal beta cells (P-value = 4.81 × 10–14), macrophage (P-value = 1.13 × 10–10), alpha cells (P-value = 1.15 × 10–12), acinar cells (P-value = 1.26 × 10–9) and fibroblast cells (P-value = 6.55 × 10–11). Dapsone is an antibiotics and anti-inflammatory medication typically used for skin disorders. Previous reports have indicated that dapsone can induce pancreatitis, especially acute pancreatitis (AP) (39), and beta cell dysfunction is involved in the early stages in pancreatitis (40). In addition, it is well known that macrophage cells can secrete proinflammatory cytokines and further expedite the formation of pancreatic fibrosis in chronic pancreatitis (CP) (41). Moreover, premature activation of digestive enzymes in acinar cells may lead to the onset of AP. Thus, the cell types we identified well supported previous studies. Importantly, the signature genes associated with dapsone and cell types are also significantly enriched in GTEx (42) pancreas and adipose tissues as well as cytoplasmic translation pathways (Figure 3D). Although the functions of pancreas associated adipose tissue (PAT) are still unknown, recent studies have indicated that PAT released adipokines may protect the pancreas against the dysfunction of metabolism (43). These results provided insights into understanding the potential mechanism of the dapsone induced side effect in pancreatitis.

On the contrary, by comparing the results in normal and tumor cells, CeDR results can be used to infer the potential drugs that can be repurposed individually for the treatment of pancreas cancer. For example, numerous reports have indicated that B cells can promote pancreatic tumorigenesis and have been investigated as a potential target in cancer treatment (44). In our results, the compound mefloquine, previously used to treat malaria, was suggested to be effective for B cells (P-value = 3.65 × 10–12). This is consistent with previous reports that malaria drugs may improve the effectiveness in cancer therapies. Moreover, tumor-associated macrophages (TAM) as versatile immune cells can lead to a variety of malignant changes in pancreatic cancer (45). Our results suggested that betulin and naringenin are potential drugs associated with TAM for pancreatic cancer treatment. Indeed, betulin has already been suggested as a potent anticancer agent based on a comprehensive review and naringenin was reported to decrease invasiveness and metastasis (46,47). Among the signature genes in pancreatic cancer scRNA-seq, there were well-known genes highly expressed in pancreatic cancer, such as OLR1, LSP1, FGL2 (Figure 3E) (48). Interestingly, our results can also be used to infer drug combinations in the treatment of cancer. For instance, the combined treatment with naringenin and mefloquine may potentially inhibit the growth and metastasis of pancreatic cancer cells without damaging the normal cells. Notably, naringenin has already been identified as a combination therapy with hesperetin for pancreatic cancer (49). Collectively, these results demonstrated that CeDR can be used as a valuable resource for better investigation of tissue heterogeneity and drug response at the single cell level.

DATABASE DESIGN AND UPDATES

CeDR is hosted by a local server with a Centos Linux 7.4 environment. It is constructed using the Java Spring Boot (https://spring.io/projects/spring-boot) as a back-end RESTful web service framework. The cell type and drug association results are hosted by the MySQL Database service (https://www.mysql.com). Front-end user interfaces are designed using the React (https://reactjs.org) and Umi (https://umijs.org) frameworks, which are scalable enterprise-class front-end application frameworks allowing flexible maintenance and extension in future. Ant Design (https://ant.design) is used as the UI library which contains a set of high-quality components and demos. Furthermore, we used the charting library available through HighCharts (https://www.highcharts.com/) to implement interactive charting and data visualization. The CeDR resource will be updated twice per year depending on the number of newly published scRNA-seq, especially those with cell type annotations.

CONCLUDING REMARKS AND FUTURE DEVELOPMENT

A large number of scRNA-seq datasets have been generated to decode the cell type compositions and expression heterogeneity across different tissues. With these available data, CeDR provides direct references for cellular drug response profiles and has implications in various future applications. First, many cells deposited in our database were identified from disease tissues, such as tumor tissues from cancer patients, post-mortem human brain samples from Alzheimer's disease, and pancreatic tissues from patients with type 2 diabetes. We identified candidate drugs that were consistent with previous reports. Second, the comprehensive cell-type specific drug sets can be used for design of combinatory treatments. In cancer samples, cellular heterogeneity plays important roles in shaping the actual drug response and eventually leads to drug resistance. Identification of drugs, or drug combinations, that can kill cancer cells but not normal cells would be of high interest in precision medicine. Third, CeDR includes not only disease cell types but also normal cell types. The cellular drug response to the normal cell types would provide new insights for drug side effects.

Due to the limitation of well-known tumor-related scRNA-seq as well as the cell annotation information, the drug response analysis for tumors did not include all tissues. However, with the rapidly evolving technologies in next-generation sequencing, we expect more data at the single cell level to be generated in the near future, especially tumor related datasets. We will continuously collect and curate the emerging single-cell transcriptome data and conduct cellular drug response analyses for different tissues and phenotypes. In addition to scRNA-seq, other single cell techniques, such as scDNA-seq, scRNA-seq and scATAC-seq, are also valuable resources to investigate drug response. We will incorporate more drug treatment application tools for various types of data and develop CeDR as a comprehensive database for cellular drug response. Future database updates will thus include and integrate the genomic and transcriptome profiles of individual cells. We believe that CeDR will be a useful resource for the single-cell and drug design research community.

Click here for additional data file.