iDOMO: identification of drug combinations via multi-set operations for treating diseases

Abstract

Combination therapy has become increasingly important for treating complex diseases which often involve multiple pathways and targets. However, experimental screening of drug combinations is costly and time-consuming. The availability of large-scale transcriptomic datasets (e.g. CMap and LINCS) from in vitro drug treatment experiments makes it possible to computationally predict drug combinations with synergistic effects. Towards this end, we developed a computational approach, termed Identification of Drug Combinations via Multi-Set Operations (iDOMO), to predict drug synergy based on multi-set operations of drug and disease gene signatures. iDOMO quantifies the synergistic effect of a pair of drugs by taking into account the combination’s beneficial and detrimental effects on treating a disease. We evaluated iDOMO, in a DREAM Challenge dataset with the matched, pre- and post-treatment gene expression data and cell viability information. We further evaluated the performance of iDOMO by concordance index and Spearman correlation on predicting the Highest Single Agency (HSA) synergy scores for four most common cancer types in two large-scale drug combination databases, showing that iDOMO  significantly outperformed two existing popular drug combination approaches including the Therapeutic Score and the SynergySeq Orthogonality Score. Application of iDOMO to triple-negative breast cancer (TNBC) identified drug pairs with potential synergistic effects, with the combination of trifluridine and monobenzone being the most synergistic. Our in vitro experiments confirmed that the top predicted drug combination exerted a significant synergistic effect in inhibiting TNBC cell growth. In summary, iDOMO is an effective method for the in silico screening of synergistic drug combinations and will be a valuable tool for the development of novel therapeutics for complex diseases.

drug repurposing, drug combination, multiple set operations, synergistic effect, drug gene signature, disease gene signature

Introduction

Genetic and phenotypic heterogeneity among cancer patients causes substantial drug resistance [1–3]. Therapies with combinations of two or more drugs are designed to overcome drug resistance, under the assumption that cancer cells resistant to one drug might be responsive to another one. A study analysed many effective and approved combinations concluding that therapeutic benefit is likely to arise from independent drug action rather than drug additivity or synergy [4]. Currently, available data from clinical studies suggest the best treatment outcome is achieved by a combination of targeted drugs with conventional chemotherapy [5, 6]. However, the identification of successful drug cocktails is not a simple task and is hampered by the lack of standards in terminology, experimental protocols and models as well as data analysis. Given millions of possible drug combinations from thousands of drugs, it is crucial to accurately predict effective drug combinations using various functional genomic data in silico. This is necessary for the identification and successful implementation of combination therapy. Generally, both similar and distinct gene expression signatures can lead to the synergy of drug combinations [7]. Several approaches have been developed to predict drug synergy from gene expression-derived drug signatures [8–13], protein–protein interaction network-based [14, 15] or deep-learning/machine-learning [15–19]. However, those approaches did not consider disease signatures, a synergistic combination of drug signatures does not necessarily mean it is beneficial for treating a specific disease. For example, Therapeutic Score (TS) [12] was developed to predict synergy score based on disease signature between disease and control, and drug signatures of before and after treatment. However, TS doesn't take the drug signature size into consideration, while drugs inducing a large number of non-treatment-related gene expression changes imply significant side effects. Another study developed SynergySeq Orthogonality Score (OS) [20] from transcriptional consensus signatures (TCS) across multiple replicates and cell lines. However, the TCS was obtained from various model systems, which cannot exactly measure drug responses in each specific cell line. Two drugs are likely to have preferred synergistic effect if the combination of their signatures have enhanced beneficial effect and reduced side effect.

Transcriptomic data of drug perturbations such as the Connectivity Map (CMap) [21] and the Library of Integrated Network-based Cellular Signatures (LINCS) project [22] have been widely used to search for potential novel indications through matching signatures of diseases and drugs/compounds based on genome-wide gene expression phenotypes. Under the assumption that drugs that can reverse gene expression alterations in diseases would have potential therapeutic effect on restoring disease states to healthy conditions, we previously developed a drug repositioning algorithm EMUDRA [23] based on drug-induced gene signatures and disease gene signatures. On top of individual drugs predicted by methods like EMUDRA, this study develops a new approach, termed ‘Identification of Drug Combinations via Multi-Set Operations (iDOMO)’, to prioritize drug combinations using multi-set operations based quantification of a combination's potential synergistic effect using gene expression signatures of diseases and drugs. We systematically compared iDOMO with two existing drug combination methods, TS and SynergySeq OS, using two drug combination databases including DrugCombDB and DrugComb. The two databases contain drug combination tests from high-throughput screening studies, external databases and manual curations and in particular, they include experimentally tested data from over 1 million drug combinations across >100 cell lines/tissues.

Results

Development of the iDOMO approach

We developed a new approach iDOMO to identify drug combinations by modeling each drug pair’s beneficial and detrimental effects on diseased patients based on drug and disease gene expression signatures. Specifically, the up and down-regulated genes in a disease are determined by comparing gene expression data in a disease-relevant tissue from patients with the disease and those from healthy controls while the up and down-regulated genes for a drug are identified by comparing gene expression data in model systems from treatment experiments with the drug and those from control experiments. Next, these six up or down-regulated gene signatures from the disease and two drugs are further divided into 27 classes by multi-set operations including intersection and difference (Fig. 1A). Specifically, drug or drug-combination-induced gene expression changes that are opposite to the changes in the disease are defined as beneficial changes, while drug- or drug-combination-induced changes that are consistent with those in the disease are detrimental (Fig. 1B). In addition, gene expression changes that are induced by two drugs but in opposite directions are neutral. Then, the significance of the overlap among 3 gene sets of each class could be calculated using the Super Exact Test (SET) [24] |$P\left(x\ge \left|A\cap B\cap C\right|\right)=1-\sum_{x=0}^{x-1} \mathit{\Pr}\left(|A\cap B|=j;|A|=a,|B|=b\right)\ \mathit{\Pr}\left(|A\cap B\cap C|=x;|A\cap B|=j,|C|=c\right)$| (Fig. 1B). Scores of beneficial classes and detrimental classes (see Methods for details) are summed up as beneficial scores (BS) and harmful scores (HS), respectively. Similarly, for each drug in a combination, an individual drug score (IDS) can be calculated for each drug based on Fisher’s exact test (FET) to evaluate its overall beneficial effect. Lastly, the final score for the drug combination, termed total drug combination score (TDCS), is calculated as the difference between the drug combination scores (BS and HS);

$$ TDCS\left(d,{g}_A,{g}_B\right)= BS- HS $$

and a combination synergy score (CSS), is calculated as the difference between the TDCS and the maximum of the IDS against a disease signature:

$$ CSS\left(d,{g}_A,{g}_B\right)= BS- HS-\max \left({IDS}_A,{IDS}_B\right). $$

Figure 1

Evaluation of the effect of the combination of two drugs on a disease of interest. (A) Venn diagram of 27 subsets from the multi-set operations of the disease and individual drug signatures. Left panel shows the Venn diagram of down-regulated genes in a disease signature and two drug signatures of a combination. The right panel shows the Venn diagram of up-regulated genes in a disease signature and two drug signatures in combination. Red indicates up-regulated genes in a signature and blue indicates down-regulated genes. (B) Scoring 27 subsets of the disease signature |$d$| and opposite drug signatures |$-{g}_A$| and |$-{g}_B$|⁠. -1 indicates the gene subsets down-regulated in the disease signature or down-regulated in the drug signatures, while 1 indicates the gene subsets up-regulated in the disease signature or down-regulated in the drug signatures, and 0 indicates no significant change subsets in signatures. |$-{g}_A\ast d$| indicates the drug effect of drug a on the disease |$d$|⁠, and |$-{g}_B\ast d$| indicates the drug effect of drug b on the disease |$d$|⁠, while |$\left(-{g}_A\ast d\right)\ast \left(-{g}_B\ast d\right)$| indicates the effects of drugs a and b on disease |$d$|⁠. P indicates the p-value calculated by the super exact test. B indicates a beneficial overlap between the disease and drug signatures, while N indicates neutral changes and H indicates a detrimental (harmful) overlap between disease and drug signatures.

Performance evaluation through simulations

We further designed a procedure to evaluate the robustness of TDCS. First, we generated disease signatures by 1) combining up-regulated genes in a target drug pair as disease-down-regulated genes, and 2) combining down-regulated genes in the same drug pair as disease-up-regulated genes. Then, we added random genes as noise to the disease signature and then calculated TDCS for each drug combination. If the drug pair is ranked at the top with various noise levels, it means the approach can robustly identify the right drug combination for a disease signature. We performed 1000 simulations for each noise rate level. As shown in Fig. 2, the average rank is less than 10 with 35% or less of random genes in a signature. However, the average rank is significantly increased with more than 70% or more random genes added. The result indicates that the drug combination approach can robustly identify drug combinations.

Figure 2

Robustness of the proposed drug combination method iDOMO. Two drugs were randomly chosen as a target true combination. Disease signatures were generated by sampling genes from the signatures of the target drug pair and then adding a certain number of random genes that are not in the drug signatures. The mean ranks and standard deviations were based on 1000 repeats.

Performance evaluation of iDOMO and comparison with two existing approaches

We then evaluated the performance of iDOMO based on a DREAM challenge dataset in which cell viability was measured for 91 drug pairs of 14 individual drugs in a human B-cell lymphoma (BCL) cell line OCI-LY3 [9]. Gene expression data were also generated from the samples before and after treatment of each drug. We first identified the signatures of the forteen drugs from the gene expression data by comparing the samples after treatment to those before treatment. We further identified a BCL signature from a B-cell lymphoma gene expression dataset (GSE72026) [25]. The BCL signature is composed of 3816 up-regulated genes and 2378 down-regulated genes from the comparison of the BCL samples and the peripheral blood B cell control. Finally, we applied our SET approach to the BCL signature and the forteen drug signatures to calculate synergy scores for the 91 drug pairs.

As in the DREAM challenge study [9], we compared our prediction score of the 91 drug pairs and their matched EOB, which measures the drug pair's inhibition effect that excess additive effect in experiments (see Methods). A drug pair with EOB = 0 has an additive effect, whereas a drug pair with positive (or negative) EOB values has a synergistic (or antagonistic) effect. As shown in Fig. 3, TDCS has a significant positive correlation with the observed EOB (Spearman correlation ρ = 0.40, P-value = 8.32E-05). Furthermore, we also assessed the performance of iDOMO using a probabilistic concordance index (PC-index), which was used as the main performance measurement in the DREAM challenge study [9]. The iDOMO has a PC-index of 0.6207, which is better than the best approach, the Drug Induced Genomic Residual Effect (DIGRE) model (PC-index = 0.6130) which was tested in the DREAM challenge paper [9].

Figure 3

Performance evaluation of the iDOMO in the DREAM challenge dataset. The observed excess over bliss (EOB) data was downloaded from the DREAM challenge paper [9]. iDOMO scores were calculated based on the 91 pairs of the 14 individual drug signatures and the B-cell lymphoma signature identified from the GSE72026 [25].

We further evaluated the performance of the iDOMO in the publicly available drug combination screening databases, DrugCombDB [26] and DrugComb [27]. We first extracted the Highest Single Agency (HSA) model-based drug combination synergy score, for four widely used cell lines (MCF7, PC3, HT29, and A549), corresponding to the top four cancer types with the highest number of new cases in the US [38230766], i.e. breast carcinoma (BRCA), prostate adenocarcinoma (PRAD), colorectal adenocarcinoma (COAD) and lung adenocarcinoma (LUAD). For each of the four cell lines, we then matched the drug combinations to the drugs tested in the LINCS and extracted the combinations with both drugs tested in the LINCS. Next, to calculate the iDOMO drug combination scores for those drug combinations matched to LINCS, we identified the drug signatures between the drug-treated samples and vehicle control samples for the four cell lines. Furthermore, we identified cancer signatures for the corresponding four cancer types using the gene expression samples from the TCGA under FDR < 0.05 and fold change >2, which were summarized in Supplementary Table S1, which is available online at http://bib.oxfordjournals.org/. Finally, we calculated the iDOMO score using the cancer signature and drug signatures for the matched drug combinations for each of the four cancer types. We then used the concordance index (C-index) and Spearman correlation to assess the performance of the iDOMO scores on the prediction of the HSA synergy score. As shown in Table 1, the C-index is about 0.60 across the four cancer types in both the DrugCombDB and DrugComb databases. Furthermore, the iDOMO TDCS is significantly positively correlated with the HSA synergy score (Fig. 4). For instance, the C-index of the BRCA is 0.6093 for the prediction of the 550 drug combinations in the DrugCombDB, and the iDOMO TDCS is significantly correlated with the HSA synergy score (Spearman’s ρ = 0.3243, P = 6.20E-15) (Fig. 4A). Moreover, iDOMO performed better for BRCA in the DrugComb. We matched 597 drug combinations for the LINCS MCF7 in the DrugComb database. The C-index of the iDOMO TDCS is 0.6431 and Spearman’s ρ = 0.4200 (P < 2.20E-16) (Fig. 4B). We further tested iDOMO TDCS with the cancer signatures identified by FDR < 0.05 and fold change >1.5. The performance (Supplementary Table S2 available online at http://bib.oxfordjournals.org/">http://bib.oxfordjournals.org/ and Fig. S2) is similar to the cancer signatures based on FDR < 0.05 and fold change >2, suggesting the iDOMO TDCS is robust to the thresholds for signature identification. Furthermore, we tested the CSS in the same databases using the same signatures. The performance is slightly inferior to the TDCS but still good (Supplementary Table S3, which is available online at http://bib.oxfordjournals.org/">http://bib.oxfordjournals.org/).

Table 1

The performance of iDOMO total drug combination score on the prediction of experimental HSA synergy score.

Database	Signature*	Cell line	# Drug Combs	C-Index	Rho	P-value
DrugCombDB	BRCA	MCF7	550	0.6093	0.3243	6.20E-15
	PRAD	PC3	614	0.5578	0.1655	3.79E-05
	COAD	HT29	299	0.5897	0.2641	3.66E-06
	LUAD	A549	246	0.6052	0.3139	5.03E-07
DrugComb	BRCA	MCF7	597	0.6431	0.4200	< 2.20E-16
	PRAD	PC3	652	0.5681	0.2016	2.10E-07
	COAD	HT29	367	0.6187	0.3443	1.18E-11
	LUAD	A549	248	0.6056	0.3113	5.69E-07

^aCancer signatures were identified based on adjusted P-value <0.05 and folder change >2.

Figure 4

Performance evaluation of iDOMO in the DrugCombDB and DrugComb databases and in comparison with two existing approaches, therapeutic score (TS) and TCS based SynergySeq orthogonality score (TCS OS). Scatter plots of the total drug combination score (TDCS) and highest single agency (HSA) synergy score for four TCGA cancer types (BRCA, PRAD, COAD and LUAD) in the DrugCombDB (A) and DrugComb (B) databases. (C) Comparison of the correlation coefficients between HSA synergy scores and prediction scores of iDOMO TDCS, TS and TCS OS.

We also compared the performance of our iDOMO TDCS and two existing gene-signature-based approaches TS score [12] and SynergySeq Orthogonality Score (TCS OS) [20], in predicting the drug combination synergy score HSA in the DrugComb and DrugCombDB databases. We first matched the drug combination TS and TCS OS values of the four TCGA cancer types (BRCA, PRAD, COAD, and LUAD) to the HSA synergy scores of the corresponding cancer type cell lines (MCF7, PC3, HT29, and A549) of the same drug combinations in the two databases. We then calculated each cancer type’s C-index and Spearman correlation coefficient between TS or TCS OS prediction and HSA synergy score. Generally, TS and TCS OS showed worse performance than our iDOMO TDCS measured by the C-index of predicting the HSA synergy score and the correlation between the prediction scores and the corresponding HSA observations (Supplementary Table S4 available online at http://bib.oxfordjournals.org/, Fig. 4C and Fig. S3), except for TS prediction for the PRAD in the DrugComb database. Among the three approaches, TCS OS is the worst. For example, the TS values extracted from the BRCA signatures were not correlated with the HSA scores in the DrugCombDB when the same drug combinations were considered (Spearman’s ρ = 0.0876, P = 0.691, Fig. S3A and Table S4). The C-index was 0.5079 as a random prediction (Table S4). Notably, the C-indices of TCS OS were very close to 0.5 (Fig. 4C and Table S5) and Spearman’s correlation coefficients were close to 0 (Fig. S4) except for the prediction for COAD (Table S5). Therefore, the extensive performance comparison demonstrates iDOMO’s superior performance over two popular approaches for predicting drug combinations.

In summary, iDOMO outperforms the existing approaches for predicting drug combinations in both the DREAM challenge and public drug combination databases.

Application of iDOMO to triple negative breast cancer

We then applied iDOMO to predict drug combinations for treating triple negative breast cancer (TNBC). In our previous EMUDRA study [23], we predicted 20 drugs for TNBC using a TNBC signature of 4776 differentially expressed genes (DEGs) between the TNBC samples and the matched controls in the TCGA cohort [28]. We calculated the TDCS for the 190 combinations of the 20 drugs using the 4776 DEGs and drug signatures identified from the CMap [21]. As shown in Table 2, the combination of two compounds, trifluridine and monobenzone, was predicted to be most effective with a SET combination score of 123.63. Notably, trifluridine is a component in 7 of the predicted top 10 most synergistic drug pairs, while monobenzone is a component in 3 of the 10 drug pairs. As shown in Fig. 5, trifluridine and monobenzone significantly reversed the TNBC signature individually and in combination.

Table 2

The top 10 pairs of compounds that are predicted to be synergistic for treating TNBC.

Drug A	Drug B	BS-HS	|${IDS}_A$|	|${IDS}_B$|	DCS
Trifluridine	Monobenzone	123.63	40.94	26.74	55.95
Trifluridine	0173570–0000	113.16	40.94	21.1	51.12
Monobenzone	0173570–0000	95.97	26.74	21.1	48.13
Trifluridine	Ciclopirox	107.34	40.94	24.25	42.15
Trifluridine	Vidarabine	106.54	40.94	24.84	40.76
Trifluridine	Methotrexate	115.68	40.94	35.23	39.51
Resveratrol	Trifluridine	101.76	30.42	40.94	30.4
Etoposide	Monobenzone	110.79	64.53	26.74	19.52
Etoposide	Vidarabine	99.41	64.53	24.84	10.04
Etoposide	Trifluridine	110.98	64.53	40.94	5.51

Figure 5

SET analysis of the drug pair trifluridine (Trif) and monobenzone (mono). (A) Venn diagram of 27 subsets from the multi-set operations of the triple-negative breast cancer (TNBC) signature and Trif and mono signatures. The left panel shows the Venn diagram of down-regulated genes in the TNBC signature and the up- and down-regulated genes of the Trif and mono signatures. The right panel shows the Venn diagram of up-regulated genes in the TNBC signature and the up- and down-regulated genes of the Trif and mono signatures. Red indicates up-regulated genes in a signature and blue indicates down-regulated genes. (B) Scoring 27 subsets of the TNBC signature |$d$| and drug signatures trifluridine and monobenzone. The GO enrichment column shows the most significant enriched GO biological processes of the genes shared by the TNBC and drug signatures.

Experimental validation of the top predicted drug combination

As the combination of trifluridine and monobenzone was predicted by iDOMO to be most effective for treating TNBC, we further experimentally validated the drug pair in vitro. We first calculated the correlation coefficients of the gene expression levels of the untreated TNBC cell lines in CCLE [29] and the primary TNBC samples from TCGA [28]. Among the cell lines that showed significant correlation with most of the TCGA primary TNBC patient samples (Fig. S5), SUM159 is the only cell line available in our lab. So, we chose SUM159 for our in vitro validation experiment cell to validate the synergistic effects of trifluridine and monobenzone. Trifluridine is a fluorinated thymidine analog, which can incorporate into DNA and inhibit DNA replication leading to the death of cancer cells [30]. The exact mechanism of action of monobenzone is not fully understood, but it may inhibit cancer cell growth by inhibiting ribonucleotide reductase [31] and inducing oxidative stress and inflammatory cytokines [32]. Thus, the two drugs may have synergistic anti-tumor effects by targeting different pathways.

We treated SUM159 cells with trifluridine or monobenzone or a combination (ratio 1:2) at seven concentrations. Compared to the dimethyl sulfoxide (DMSO) control, the cell number was significantly decreased with trifluridine or monobenzone treatment at each of the seven concentrations (Fig. 6A and B). Moreover, the combination of the 2 drugs showed a synergistic effect at the lowest four doses (T25 + M50, T12.5 + M25, T6.25 + M12.5, T3.125 + M6.25) (Fig. 6C and D). We quantitatively analyzed the synergism of the combination of trifluridine and monobenzone using Chou-Talalay’s combination index method. The results indicated that the combination exerted a synergistic effect when cell growth inhibition was less than 50%, with the combination index (CI) ranging from 0.785 to 0.577 for the fraction affected (Fa) ranging from 0.694 to 0.896, as shown in the Fa-CI plot (Fig. 6C). Thus, this experiment validated our prediction that there is a synergistic effect in the combination of trifluridine and monobenzone.

Figure 6

Validation of the top predicted drug combination for treating TNBC. (A) the cell effect of trifluridine (multiple of 3.125 μM), monobenzone (multiple of 6.25 μM) and their combination (ratio 1:2) at 7 concentrations compared to DMSO baseline. The differences between individual drugs and combination were tested by t-test and shown by markers *** (P < 0.0001) and * (P < 0.01). (B) the dose-effect plot of trifluridine, monobenzone and their combination. The effect is shown by the fraction affected (FA). (C) the FA-CI (combination index) plot of the trifluridine and monobenzone combination at a 1:2 ratio. CI < 1 indicates synergy; CI = 1 indicates additive; CI >1 indicates antagonism. CI is shown in the log10 scale. (D) Dose-reduction index (DRI) (see methods) plot with different FAs, DRI is shown in the log10 scale.

Discussion

In this study, we developed a new method iDOMO to predict the synergistic effect of drug combinations based on disease and drug gene expression signatures. Based on the concept that reversing disease gene expression changes would reverse disease-related phenotypes, we modeled the beneficial and detrimental effects of a combination using the drug and disease gene signatures. Several tools have been developed to predict drug synergy using machine learning methods [33–35]. However, those tools require multiple experimental data points to calculate dose–response matrices. It is costly and time-consuming to generate dose–response matrices as it requires drug treatment experiments across multiple dosages and time points. iDOMO takes into account all possible combinations of the up- and down-regulated gene signatures of disease and two drugs for computing the final combination score. Thus, the synergy of drug combinations can be predicted for disease signatures derived from massive drug treatment data, e.g. the CMap and LINCS. In the CMap, drug combination can be calculated among over 1300 drugs in three cell lines, while in the LINCS dataset, our approach can be applied to over 12 000 FDA-approved and investigational compounds in over 70 cell lines. Nevertheless, gene expression change patterns are essential for drug synergy but not sufficient to detect a synergy in a specific disease [7]. In theory, it is ideal to target different parts of a disease network or pathway with multiple drugs to maximally normalize the dysregulated network or pathway [14]. On the other hand, concurrent inhibition of parallel subnetworks or pathways that drive proliferation and anti-apoptosis can acquire the best efficiency for cancer treatment [36, 37].

Application of iDOMO to TNBC identified the combination of trifluridine and monobenzone as a potential novel therapy for treating TNBC. Trifluridine is a thymidine-based nucleoside analog, which is incorporated into DNA via phosphorylation to inhibit DNA replication and cell proliferation. Trifluridine has been used as an antiviral agent to treat herpetic keratitis [38]. Currently, it is usually co-administrated with tipiracil, a thymidine phosphorylase inhibitor, to treat metastatic colorectal and gastric cancers [36, 39]. Monobenzone is a potent skin-bleaching agent used as a treatment for over-pigmentation in the clinic [40]. It was also reported that monobenzone is used to treat cutaneous metastatic melanoma patients in a phase 2 trial [37]. Recently, monobenzone was shown to be a potent inhibitor for lysine-specific demethylase 1 (KDM1A) [41], which is an oncogene that promotes the initiation and development of various cancers via diverse cellular signaling pathways [42, 43]. Monobenzone can inhibit the proliferation and migration of various cancer cells by competitively inhibiting KDM1A reversibly [41]. In our analysis, both trifluridine and monobenzone can down-regulate the gene expression of DNA replication, DNA repair and mitotic cell cycle. Furthermore, trifluridine can down-regulate the gene expression of ribonucleoprotein complex biogenesis, tRNA processing and rRNA processing, while monobenzone down-regulated more DNA replication and cell cycle-related genes as well as microtubule-based process, RNA methylation and protein ubiquitination genes. Besides, monobenzone up-regulated genes of response to cAMP and cytokine biosynthesis, which may increase the immune response to cancer cells. The above functional analysis of the genes induced by the trifluridine and monobenzone combination suggested that the synergy may be achieved by targeting the same pathways from different directions and different pathways essential to cancer. Additional pathways may be targeted by the combination of the two drugs. Nevertheless, drug synergy-induced gene expression changes should be measured to explore the synergistic mechanisms of the combination.

Previous studies have suggested that synergy is condition-specific [44–46]. Cell type, disease subtype, mutation and gene expression status may influence the synergy of drug combinations. In the drug combination DREAM Challenge, the results from the best-performed algorithms suggested that the similarity between the gene expression changes from monotherapies is associated with synergy [9]. A recent study showed that synergistic combinations could generate unique transcriptomic changes that were not induced by any of the monotherapies [7]. This study also demonstrated that the correlation of transcriptomic response of monotherapies was necessary but not sufficient for drug synergy [7]. In this study, we have shown that comparison of drug-gene signatures and different parts of a disease signature can be used to predict synergy with improved performance (PC-index) than the algorithms developed through the DREAM Challenge. Several factors may affect the performance of iDOMO. First, it is essential to identify the gene signature of a disease under study from the primary tissue(s) responsible for the disease. Second, selecting disease-relevant cell lines in LINCS is crucial, as different cell lines may exhibit varied molecular signatures in response to drug treatment. This selection can be achieved by comparing the expression profiles of LINCS cell lines with those from disease-causing tissues, allowing us to choose the most similar cell lines for drug repurposing. Third, identifying optimal drug combinations from a large pool of drugs is computationally intensive. Thus, optimizing iDOMO to quickly screen drug combinations is necessary. Lastly, the current performance evaluation is based only on a limited set of drug combinations available in both LINCS and the DrugCombDB or DrugComb databases. The drug repurposing research community should collaborate to generate larger datasets, including matched pre- and post-treatment expression data and viability information across many cell lines for future studies. With the pre- and post-treatment data from the CMap and LINCS phase I and phase II, our approach could be an efficient and cost-effective method for in silico screening of synergistic drug combinations for a variety of diseases such as cancers and neurodegenerative diseases.

Methods

Data preprocessing and normalization

The microarray data of drug and DMSO treatments were downloaded from the CMap [21]. Raw data were preprocessed as described in the EMUDRA paper [23]. Briefly, raw intensities of probe sets were preprocessed using the robust multi-array average (RMA) procedure [47] and summarized to a gene-level expression matrix of 12 328 genes. The expression matrix was then corrected for batch effects using a linear model [48]. For the LINCS dataset, the normalized level 3 data matrices of compound treated and control samples were downloaded from the GSE92742 (LINCS v1.1) of the Gene Expression Omnibus (GEO) database [49] and the LINCS data portal (LINCS beta) (https://clue.io/releases/data-dashboard). The LINCS data matrix includes the expression levels of directly measured landmark genes and the gene expression data inferred from those directly measured genes. In total, 12 328 genes were included in the analysis of this study. The data matrices were also corrected for batch effects. The level 3 RNA-seq data of TNBC samples and adjacent normal samples were downloaded from The Cancer Genome Atlas (TCGA) data portal [50].

The RNA-Seq gene expression data for the four cancer types (BRCA, PRAD, COAD, and LUAD) were downloaded from TCGA.

Signature identification

To identify a signature for each drug/compound of a cell line in the CMap, the gene expression profiles of experiments treated by various doses and treatment durations were combined and compared to DMSO control profiles of the same batch. P-values of differential expression were calculated using a linear model from the R package limma [37]. Limma uses an Empirical Bayes moderated t-statistic to calculate p-values with modified standard error and degrees of freedom. The standard error is moderated across genes with borrowed information from all genes to improve the inference of any single gene. The degrees of freedom are also adjusted using a factor representing a priori degrees of freedom for the model. To correct the multiple-testing problem, P-values were adjusted using the Benjamini–Hochberg (BH) procedure [51]. Drug signatures were then identified using the thresholds of adjusted P-value <0.05. Analogously, drug signatures of the LINCS dataset were also identified using the same procedure and thresholds for each cell line. As each cell line was treated with each drug in six different concentrations in more than a dozen batches, those profiles were treated as replicates and put together to identify one high-quality drug-induced signature for each cell line. Similarly, the TNBC signature was also identified using the limma package [48] by comparing the TNBC samples to the matched normal controls in TCGA. The signatures of the four cancer types were also identified by comparing the cancer and normal samples using the linear model in the limma package [48] and P-values were corrected with the BH procedure [51]. Significant differentially expressed genes (DEGs) of the signatures were filtered by FDR < 0.05 and FC > 2 or 1.5.

Drug combination scores based on detrimental and beneficial effects

We propose a SET-based drug combination score (DCS) to model the overall effect of a disease signature and drug signatures. This DCS evaluates how the different intersections of the drug and disease signatures exert a beneficial or detrimental effect on the disease under study. Specifically, the DCS is based on the following assumptions: (i) a change (up- or down-regulation) induced by one or both drugs which are opposite to that induced by a disease under study is beneficial; (ii) an expression change in a combination is detrimental if the change is in the same direction as that induced by the disease; (iii) a change that is in the opposite directions by two drugs is neutral (Fig. 1B).

Suppose there are three gene sets A, B, and C with a, b and c genes, respectively. Suppose j genes overlap between A and B, and x genes overlap between A, B, and C. The probability of the intersection is exact x can be calculated using the following equation:

$$ \begin{align*} \Pr &\left(\left|A\cap B\cap C\right|=\ x;\left|A\right|=a,\left|B\right|=b,\left|C\right|=c\right)\\ &=\ \Pr \left(\left|A\cap B\right|=j;\left|A\right|=a,\left|B\right|=b\right) \\ &\quad \times\Pr \left(\left|A\cap B\cap C\right|=x;\left|A\cap B\right|=j,\left|C\right|=c\right) \end{align*}$$

Note that

$$ \begin{equation} 0\le j\le \min \left(a-x,b-x,c\right) \end{equation} $$

(1)

Suppose that there is a disease signature |$d$| comprised of the up-regulated gene set |${d}_{up}$|⁠, down-regulated gene set |${d}_{dn}$| and gene set |${d}_0$| without significant change. Similarly, a drug signature |${g}_A$| (⁠|${g}_B$|⁠) includes the up-regulated gene set |${g_A}_{up}$| (⁠|${g_B}_{up}$|⁠), the down-regulated gene set |${g_A}_{dn}$| (⁠|${g_B}_{dn}$|⁠), and the unchanged gene set |${g_A}_0$| (⁠|${g_B}_0$|⁠). As shown in Fig. 1, there are 27 subsets from the multi-set operations on the aforementioned gene signatures. Given three gene sets, the probability of obtaining no less than a given number of overlapping genes can be calculated using the Equation 1 [24]. The P-values of the six subsets with beneficial effects can be calculated as follows:

$$ {p}_1=\Pr \big(x\ge |{d}_{dn}\cap{g_A}_{up}\cap{g_B}_{up}|\big) $$

$$ {p}_2=\Pr \big(x\ge |{d}_{dn}\cap{g_A}_{up}\cap{g_B}_0|\big) $$

$$ {p}_3=\Pr \big(x\ge |{d}_{dn}\cap{g_A}_0\cap{g_B}_{up}|\big) $$

$$ {p}_{25}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_0\cap{g_B}_{dn}|\big) $$

$$ {p}_{26}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_{dn}\cap{g_B}_0|\big) $$

$$ {p}_{27}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_{dn}\cap{g_B}_{dn}|\big) $$

The |${p}_1$|⁠, |${p}_2$|⁠, |${p}_3$|⁠, |${p}_{25}$|⁠, |${p}_{26}$| and |${p}_{27}$| can be summarized into a beneficial score (BS) in Equation 2 as follows:

$$ \begin{equation} BS=\sum_{k=1}^3{-\mathit{\log}}_{10}\left({p}_k\right)+\sum_{k=25}^{27}{-\mathit{\log}}_{10}\left({p}_k\right) \end{equation} $$

(2)

Similarly, the gene subsets with detrimental effects can be calculated using equation 1 as follows:

$$ {p}_7=\Pr \left(x\ge |{d}_{dn}\cap{g_A}_0\cap{g_B}_{dn}|\right) $$

$$ {p}_8=\Pr \left(x\ge |{d}_{dn}\cap{g_A}_{dn}\cap{g_B}_0|\right) $$

$$ {p}_9=\Pr \left(x\ge |{d}_{dn}\cap{g_A}_{dn}\cap{g_B}_{dn}|\right) $$

$$ {p}_{10}=\Pr \big(x\ge |{d}_0\cap{g_A}_{up}\cap{g_B}_{up}|\big) $$

$$ {p}_{11}=\Pr \big(x\ge |{d}_0\cap{g_A}_{up}\cap{g_B}_0|\big) $$

$$ {p}_{12}=\Pr \big(x\ge |{d}_0\cap{g_A}_0\cap{g_B}_{up}|\big) $$

$$ {p}_{16}=\Pr \left(x\ge |{d}_0\cap{g_A}_0\cap{g_B}_{dn}|\right) $$

$$ {p}_{17}=\Pr \left(x\ge |{d}_0\cap{g_A}_{dn}\cap{g_B}_0|\right) $$

$$ {p}_{18}=\Pr \left(x\ge |{d}_0\cap{g_A}_{dn}\cap{g_B}_{dn}|\right) $$

$$ {p}_{19}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_{up}\cap{g_B}_{up}|\big) $$

$$ {p}_{20}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_{up}\cap{g_B}_0|\big) $$

$$ {p}_{21}=\Pr \big(x\ge |{d}_{up}\cap{g_A}_0\cap{g_B}_{up}|\big) $$

Then detrimental score (HS) can be summarized in Equation 3 by the sum of |${p}_7$|⁠, |${p}_8$|⁠, |${p}_9$|⁠, |${p}_{10}$|⁠, |${p}_{11}$|⁠, |${p}_{12}$|⁠, |${p}_{16}$|⁠, |${p}_{17}$|⁠, |${p}_{18}$|⁠, |${p}_{19}$|⁠, |${p}_{20}$| and |${p}_{21}$| as follows:

$$ \begin{equation} HS=\sum_{k=7}^{12}{-\mathit{\log}}_{10}\left({p}_k\right)+\sum_{k=16}^{21}{-\mathit{\log}}_{10}\left({p}_k\right) \end{equation} $$

(3)

4, 5, 6, 13, 14, 15, 22, 23 and 24 are neutral (neither beneficial nor detrimental). Therefore, the BS and HS can be summarized into aTDCS in Equation 4:

$$ \begin{equation} TDCS\left(d,{g}_A,{g}_B\right)= BS- HS \end{equation} $$

(4)

Similarly, we defined a score for each drug signature based on the P-values the FET of the overlap between the drug signature and a disease signature. As shown in Supplemental Fig. S1, the significance of the overlap between individual drug signature and disease signature was log 10 transformed and added up in Equation 5 as follows to obtain IDS against a disease signature.

$$ \begin{align} {IDS}_i=&\ {-\mathit{\log}}_{10}\left({p}_{i1}\right){-\mathit{\log}}_{10}\left({p}_{i3}\right)+{\mathit{\log}}_{10}\left({p}_{i2}\right)+{\mathit{\log}}_{10}\left({p}_{i4}\right) \nonumber \\ &+{\mathit{\log}}_{10}\left({p}_{i6}\right)+{\mathit{\log}}_{10}\left({p}_{i8}\right) \end{align} $$

(5)

Finally, the BS, HS and IDSs are summarized into a TDCS and a combination synergy score (CSS) in Equation 6:

$$ \begin{equation} CSS\left(d,{g}_A,{g}_B\right)= BS- HS-\max \left({IDS}_A,{IDS}_B\right) \end{equation} $$

(6)

A combination is predicted to be synergistic if the TDCS and CSS are larger than 0.

Combinatorial therapeutic score and TCS based OS

Combinatorial therapeutic score (TS) scores of BRCA, PRAD, COAD, and LUAD were extracted from the Supplementary Table S10 of the TS paper [12]. Drug combinations with both drug names that matched the drug combinations in the databases DrugComb and DrugCombDB were further extracted to evaluate the performance of the TS in the 4 cancer types for the prediction of HSA.

Transcriptional Consensus Signature (TCS) based OS was calculated for the TCGA BRCA, PRAD, COAD and LUAD signatures using the OS algorithm [20] and TCS from the LINCS L1000 data, as described in the original SynergySeq study [20]. Specifically, TCS were calculated from the LINCS L1000 dataset, including the small molecule treatment transcriptional data of the MCF7, PC3, HA1E, HCC515, VCAP, A375, HEPG2, HT29, A549, ASC, NEU, and NPC cell lines. Aggregation was performed for each cell line first. By measuring gene expression data before and after drug perturbation, a gene was included in the aggregated expression profile of a cell line if it had a mean |z-score| > 1 across the samples using both the same small molecule and the same cell line. Next, the aggregated gene expression profiles were collapsed for each the small molecule across all cell lines to produce the TCS. A gene was included in the final TCS of a small molecule across cell lines if it was consistently up/down-regulated by a |z-score| > 1 in more than 30% of the cell lines that were treated by the same small molecule. The concordance ratio (CR) and the disease discordance ratio (DR) were then calculated for each drug pair [20]. Herein, CR is defined as the ratio of TCS genes of a compound/small molecule that in the same direction as the reference TCS gene signature to the ones in the opposite directions (Equation 7). DR is defined as the ratio of TCS genes induced by a drug in an opposite direction to the ones in the same direction against a disease signature, giving that the disease signature genes that are missing in the reference signature (Equation 8):

$$ \begin{align} \mathrm{CR}=\frac{\sum_{i=1}^{\#\mathrm{genes}}\left[{a}_i\right]}{\sum_{i=1}^{\#\mathrm{genes}}\left[{b}_i\right]}\ \ \mathrm{with}\ {a}_i=\left\{\begin{array}{@{}c}1,\quad if\ {z}_i\cdot{r}_i>0\ \\{}0,\quad if\ {z}_i\cdot{r}_i<0\end{array}\right. \nonumber \\ \mathrm{and}\ {b}_i=\left\{\begin{array}{@{}c}1, if\ {z}_i\cdot{r}_i<0\\{}0, if\ {z}_i\cdot{r}_i>0\end{array}\right. \end{align}$$

(7)

$$ \begin{align} \mathrm{DR}=\frac{\sum_{i=1}^{\#\mathrm{genes}}\left[{b}_i\cdot{c}_i\right]}{\sum_{i=1}^{\#\mathrm{genes}}\left[{a}_i\cdot{c}_i\right]}\ \ \mathrm{with}\ {a}_i=\left\{\begin{array}{@{}c@{}}1,\quad if\ {z}_i\cdot{d}_i>0\ \\{}0,\quad if\ {z}_i\cdot{d}_i<0\end{array}\right., \nonumber \\{b}_i=\left\{\begin{array}{@{}c}1,\quad if\ {z}_i\cdot{d}_i<0\\{}0,\quad if\ {z}_i\cdot{d}_i>0\end{array}\right.\ \mathrm{and}\ {c}_i=\left\{\begin{array}{@{}c}1,\quad if\ {r}_i=0\\{}0,\quad if\ {r}_i\ne 0\end{array}\right. \end{align} $$

(8)

where z, d, and r denote the TCS signatures of a small molecule/compound, the disease and a reference compound/small molecule, respectively. After scaling CR with DR into the range [0 and 1], the orthogonality score (OS) of each compound to a reference compound is defined as (Equation 9).

$$ \begin{equation} \mathrm{OS}=\sqrt{{\left(1-\mathrm{CR}\right)}^2+{\mathrm{DR}}^2} \end{equation} $$

(9)

Drug combination data of DrugCombDB and DrugComb databases

The drug combination data with experimentally tested synergy scores were downloaded from the DrugComb and DrugCombDB databases. Combinations with HSA synergy scores were extracted. For the TS drug combination prediction in the Table S10 of the TS paper [12], HSA synergy scores of corresponding combinations with both drugs matched to the combinations in each of the two databases were extracted to evaluate the performance of the TS. For the iDOMO TDCS, CSS and SynergySeq OS scores, HSA synergy scores of drug combinations with both drugs mapped to compound/small molecule in the LINCS data were extracted to evaluate their performance based on prediction scores.

Gene ontology enrichment analysis for drug-induced genes

Gene ontology (GO) analysis [52] was performed for each gene subset in Fig. 1 to obtain the top biological processes that were significantly enriched with the subset genes. The analysis was based on the GO.db and org.Hs.eg.db packages of the Bioconductor [53]. P-values of enrichment were calculated by the hypergeometric distribution using an R package GO-function [54] and corrected by the Benjamini–Yekutieli procedure [55]. Redundant biological processes were removed through the procedures in the GO-function [54].

Performance evaluation by simulation with drug signatures identified from the CMap

Target drug pairs were randomly generated from the top 20 drugs predicted for the TNBA signature using the EMUDRA [23]. For each pair, a fake disease signature was generated by replacing the genes in the drug signatures induced by the golden standard drug pair with a proportion of noise genes. Those noise genes were generated from the complementary gene set of the drug signatures that were not changed by the two drugs. Then, the fake disease signature was used to test iDOMO and obtain the rank for all the drug pairs. The rank of the golden standard drug pair was used to assess the performance. This procedure was repeated 1000 times under each noise rate (proportion). The average rank and standard deviation are reported in Fig. 2.

Testing in the DREAM challenge dataset

The signatures of the 14 drugs in the DREAM Challenge dataset [9] were identified using the same procedure for identifying the CMap signatures as described above. Briefly, the gene expression data of 282 samples of the drug and DMSO-treated B-Cell Lymphoma cell line OCI-LY3 were downloaded from the GEO accession GSE51068 [9]. The probes were normalized using the RMA procedure using the Bioconductor package in R to obtain log2 transformed probeset expression levels. Probesets were then annotated to gene symbols using the platform Affymetrix Human Genome U219 Array (HG-U219). Probesets annotated to none or multiple gene symbols were discarded. Probesets annotated to the same gene symbol were averaged to obtain gene-level expression data. Lowly expressed genes with the mean expression levels at the bottom 20% were removed from further analysis. In total, 16 747 genes were included in the analysis. For each drug, differential expression analysis was conducted by comparing profiles of various concentrations and treatment duration time of the same drug to the DMSO-treated profiles. P-values were calculated using the Empirical Bayes moderated t-statistic of the R package limma [48]. BH procedure [51] was applied to the P-values to correct the multiple testing problem.

The B-cell lymphoma signature was identified from the gene expression profiles of the GEO dataset GSE72026 [25]. Normalized and log2-transformed gene expression data matrix was downloaded from the GEO database. Probesets were annotated to gene symbols using the platform Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381 (GPL17077). Probesets annotated to lincRNA were removed from further analysis. As described above, probeset-level expression values were summarized to gene-level expression values using the same procedure as the HG-U219. In total, 23 675 genes were included in the following analysis. Differential expression analysis was conducted using the R package limma [48] and P-values were corrected by the BH procedure [51]. Under a 5% FDR, 4671 up-regulated and 3544 down-regulated genes were identified by comparing 4 intravascular large B-cell lymphoma samples to four peripheral blood B cell control samples. After removing the genes that are not on the HG-U219 platform, 3816 up-regulated and 2378 down-regulated genes were used in the drug synergy prediction.

We used the probabilistic concordance index (PC-index) to evaluate the performance of the prediction, which is modified from the concordance index (C-index). C-index computes the proportion of concordance between the predicted SET scores and observed Excess over Bliss (EOB) ranks of drug pairs to quantify the performance of the prediction. If drugs |${D}_x$| and |${D}_y$| have an additive effect, the expected fractional inhibition |${f}_{add}$| of their combination should be:

$$ {f}_{add}=1-\left(1-{f}_x\right)\times \left(1-{f}_y\right) $$

where |${f}_x$| and |${f}_y$| are experimentally measured fractional inhibitions of drug |${D}_x$| and |${D}_y$|⁠, respectively. EOB is determined by subtracting the expected fractional inhibition |${f}_{add}$| from the fractional inhibition induced by the drug combination (of |${D}_x$| and |${D}_y$|⁠) |${f}_{xy}$|⁠.

$$ EOB={f}_{xy}-{f}_{add} $$

Therefore, a drug pair with 0 EOB has an additive effect, while a drug pair with positive (or negative) EOB values has a synergistic (or antagonistic) effect. The experimentally measured EOB data of the 91 drug pairs and their standard error of the mean (S.E.M.) were downloaded from the DREAM Challenge dataset [9].

To calculate the C-index, the 91 drug pairs of the 14 drugs were first ranked by the experimentally determined EOB, from the most synergistic to the additive and then the most antagonistic. Suppose |${r}_i$| and |${r}_j$| are the observed ranks of drug pairs i and j, and |$p{r}_i$| and |$p{r}_j$| are the predicted ranks of the same drug pairs i and j. The concordance score |${C}_{ij}$| of the drug pairs i and j is 1 if the observed ranks and predicted ranks are concordant; otherwise |${C}_{ij}=0$|⁠.

$$ {C}_{ij}=\left\{\begin{array}{@{}c}1,\quad if\ \left({r}_i>{r}_j\ and\ p{r}_i>p{r}_j\ or\ {r}_i<{r}_j\ and\ p{r}_i<p{r}_j\right)\\{}0,\quad if\ \left({r}_i>{r}_j\ and\ p{r}_i<p{r}_j\ or\ {r}_i<{r}_j\ and\ p{r}_i>p{r}_j\right)\ \end{array}\right. $$

The C-index is then defined as follows:

$$ C. index=\frac{2}{90\times 91}\sum_{i=1}^{90}\sum_{j=i+1}^{91}{C}_{ij} $$

Considering the deviations of the drug inhibition experiments, the c-index is modified to the PC-index as in the DREAM Challenge paper [9].

Validation experiments in TNBC cell line SUM159

SUM159 cells were cultured in F12-HAM’s medium with 5% FBS and supplemented with 5 μg/ml insulin and 1 μg/ml hydrocortisone. SUM159 cells were seeded into 96 wells at a density of 4x10^3/well. The next day, different doses of trifluridine (200, 100, 50, 25, 12.5, 6.25, and 3.125 μM) and monobenzone (400, 200, 100, 50, 25, 12.5, and 6.25 μM) were added either individually or in combination (ratio 1:2). The cells were then treated for 72 hr before being fixed with ice-cold trichloroacetic acid for 1 hr at 4°C. The cells were then washed with water and air-dried. About 100 μl of 0.057% sulforhodamine B (w/v) was added to stain the cells for 30 min at RT. The cells were then rinsed quickly with 1% acetic acid for four times to remove the unbound dye. The plates were air-dried. About 200 μl of 10 mM Tris base solution (pH 10.5) was added to each well and shaken for 30 min at RT to ensure the dye was homogeneously dissolved.

Then OD510 was recorded with a plate reader [56]. The synergistic effect of the two drugs was calculated using Compusyn following the Chou-Talalay method [57]. The synergistic effect was shown as CI values. (CI < 1 synergy, CI = 1 additive, CI > 1 antagonism). The CompuSyn analysis involved determining the CI and dose-reduction index (DRI) values for each drug combination tested at a consistent dose ratio. CI values reflect the combined effect of multiple drugs (synergistic, additive, or antagonistic). In contrast, DRI values signify fold decrease in drug dose required in a combination to achieve the same efficacy (Fa) as the individual drug. Under this context, drug combinations with CI values < 1 are deemed synergistic, and DRI values ≥ 1 are considered favorable, mainly due to concerns about toxicity associated with combining multiple drugs. To evaluate drug combinations for cancer treatment, these criteria are especially crucial if met at high effect (Fa) levels, where the goal is to eliminate a significant proportion of tumor cells. Hence, we designated any drug combination with CI < 1 and DRI ≥ 1 (for both drugs in the combination) at Fa ranging from 0.694 to 0.896 as a promising combination.

Key Points

• Previous drug combination synergy prediction approaches did not take into account differences across diseases.

• A new computational approach iDOMO was developed to predict synergistic effects of drug combinations based on multi-set operations of drug and disease gene signatures.

• iDOMO significantly outperforms two popular methods for predicting drug combinations.

• iDOMO was applied to triple-negative breast cancer and the top predicted drug combination was experimentally validated.

Acknowledgements

This work was supported in part by the computational resources and staff expertise provided by Scientific Computing at the Icahn School of Medicine at Mount Sinai.

Conflict of interest: None declared.

Funding

This work was supported in part by grants from the National Institutes of Health (NIH)/National Institute on Aging [R01AG085182, RF1AG074010, U01AG046170, RF1AG057440, R01AG057907, U01AG052411, R01AG062355, U01AG058635, RF1AG054014, RO1AG068030, and R56AG058655], and NIH/National Institute of Diabetes and Digestive and Kidney Diseases [R01DK118243].

Data availability

The data underlying this article are available in the GitHub, at https://github.com/zxx59834/iDOMO