key: cord-1030303-5qh0ohzn
authors: Pantaleón García, Jezreel; Kulkarni, Vikram V; Reese, Tanner C; Wali, Shradha; Wase, Saima J; Zhang, Jiexin; Singh, Ratnakar; Caetano, Mauricio S; Kadara, Humam; Moghaddam, Seyed Javad; Johnson, Faye M; Wang, Jing; Wang, Yongxing; Evans, Scott E
title: OBIF: an omics-based interaction framework to reveal molecular drivers of synergy
date: 2022-04-05
journal: NAR Genom Bioinform
DOI: 10.1093/nargab/lqac028
sha: 64001eb45dd3476c27f75d18d2cfbee9dbd4f68c
doc_id: 1030303
cord_uid: 5qh0ohzn

Bioactive molecule library screening may empirically identify effective combination therapies, but molecular mechanisms underlying favorable drug–drug interactions often remain unclear, precluding further rational design. In the absence of an accepted systems theory to interrogate synergistic responses, we introduce Omics-Based Interaction Framework (OBIF) to reveal molecular drivers of synergy through integration of statistical and biological interactions in synergistic biological responses. OBIF performs full factorial analysis of feature expression data from single versus dual exposures to identify molecular clusters that reveal synergy-mediating pathways, functions and regulators. As a practical demonstration, OBIF analyzed transcriptomic and proteomic data of a dyad of immunostimulatory molecules that induces synergistic protection against influenza A and revealed unanticipated NF-κB/AP-1 cooperation that is required for antiviral protection. To demonstrate generalizability, OBIF analyzed data from a diverse array of Omics platforms and experimental conditions, successfully identifying the molecular clusters driving their synergistic responses. Hence, unlike existing synergy quantification and prediction methods, OBIF is a phenotype-driven systems model that supports multiplatform interrogation of synergy mechanisms.

Superior treatment outcomes are achieved for many disease states when more than one therapeutic agent is administered (1) (2) (3) (4) . Indeed, there are many well documented instances when the therapeutic benefit of two agents administered together substantially exceeds the benefit that would be predicted by the additive effects of the agents administered individually. These clinical observations parallel the many instances of synergy observed in nature (1) (2) (3) (4) (5) . Widespread availability of high throughput technologies has allowed multi-level study of complex biological responses from genome to phenome (6) (7) (8) . Yet, there remains lack of consensus regarding the appropriate analysis of statistical and biological interactions found in antagonistic or synergistic responses (5, 9, 10) . Moreover, previously proposed strategies to analyze these interactions frequently lack sufficient generalizability to study these processes outside of their home Omics platforms (1, (10) (11) (12) (13) (14) . Thus, while synergistic therapeutic combinations may be empirically derived from fortuitous clinical experiences or through screening of bioactive small molecule libraries, the absence of established means to investigate these favorable drug-drug interactions ultimately precludes understanding of their underlying mechanisms. Similarly, quantification of biologic synergy or antagonistic interactions can be accomplished with existing tools, but the mechanisms of their interactions will remain elusive. Consequently, development of a methodology to integrate the statistical and biological components of synergistic interactions in diverse Omics settings can advance the rational design of combination therapies while affording deeper understanding of molecular interactions in a broad range of biological contexts.

Pneumonia is a major worldwide cause of death (15) (16) (17) (18) . Further, the ongoing COVID-19 pandemic highlights the need to develop novel anti-pneumonia interventions with activity against both established and emerging pathogens (19) (20) (21) . We have previously reported that a therapeutic dyad of immunostimulatory small molecules induces robust protection against a broad range of pneumonia-causing pathogens (22) (23) (24) (25) , including coronaviruses (26) , while the individual components confer only modest protection. This combination (hereafter, 'Pam2-ODN') Comprises a Tolllike receptor (TLR) 2/6 agonist, Pam2CSK4 ('Pam2'), and a TLR 9 agonist, ODN M362 ('ODN'), which stimulate protective responses from lung epithelial cells (22) . This biological response, termed inducible epithelial resistance, promotes survival benefits and microbicidal effects that significantly exceed the additive effects of the individual ligands (23, 27) . Thus, understanding the molecular mechanisms underlying this unanticipated synergy may allow optimized manipulation of epithelial antimicrobial responses and support new generations of host-based therapeutics against infections.

In the absence of a systems theory to interrogate synergistic mechanisms (5), we introduce Omics-Based Interaction Framework (OBIF) to identify molecular drivers through integration of statistical and biological interactions in synergistic biological responses. Unlike existing statistical models that are principally designed to either predict synergistic drug combinations (11) (12) (13) (14) (28) (29) (30) (31) (32) (33) or to quantify synergy between drugs (1, (34) (35) (36) , OBIF is a phenotype-driven (6) synergy interrogation model that is applicable to any combination of exposure types (e.g. drugs, genotypes, environmental conditions, etc.). To explain the molecular mechanisms underlying empirically observed instances of synergy, OBIF performs full factorial analysis (37) (38) (39) of feature expression data from single vs. dual factor exposures to identify molecular clusters that reveal synergy-mediating pathways, functions and regulators. To demonstrate the utility of OBIF, we applied this strategy to multi-Omics experimental data from epithelial cells exposed to Pam2-ODN to identify biologically relevant, unanticipated cooperative signaling events that we subsequently confirmed to be required for the synergistic pneumonia protection. Then, to demonstrate generalizability, OBIF was applied to datasets from diverse types of Omics platforms and experimental models, successfully identifying molecular clusters driving their combined responses.

All mouse experiments were performed with 6-to 10-weekold C57BL/6J mice of a single sex. Immortalized human lung epithelial (HBEC3-KT) cells were cultured in keratinocyte serum-free media (KSFM) supplemented with human epidermal growth factor and bovine pituitary extract. Mouse lung epithelial (MLE-15) cells were cultured in RPMI supplemented with 2% fetal bovine serum. Cultures were maintained in the presence of penicillin and streptomycin. The cell lines used were authenticated by the MD Anderson Characterized Cell Line Core Facility. All mouse experiments were performed in accordance with the Institutional Animal Care and Use of Committee of The University of Texas MD Anderson Cancer Center, protocol 00000907-RN01.

S-[2,3-bis(palmitoyloxy)-propyl]-(R)-cysteinyl-(lysyl ) 3lysine (Pam2 CSK4) and ODN M362 were reconstituted in endotoxin-free water, then diluted to the desired concentration in sterile PBS. For in vivo experiments, as previously described (24, 25) , the indicated ligands were placed in an Aerotech II nebulizer driven by 10 L/min air supplemented with 5% CO 2 for 30 min. The nebulizer was connected by polyethylene tubing to a polyethylene exposure chamber. Twenty-four hours prior to infections, 10 ml of PBS, Pam2 (24, 25) .

To simultaneously evaluate the expression of 161 regulatory proteins and phospho-proteins in HBEC3-KT cells after exposure to either PBS, Pam2, ODN or Pam2-ODN, a targeted high-throughput screening proteomic assay was performed by the Reverse Phase Protein Array Core Facility at The University of Texas MD Anderson Cancer Center (40, 41) . The RPPA included 4 biological replicates per treatment condition, and data is available at GitHub (www.github.com/evanslaboratory/Datasource).

For in vivo infections, frozen stock (2.8 × 10 7 50% tissue culture infective doses [TCID50] ml − 1) of influenza A H3N2 (Wide et al, 1977), virus was diluted 1:250 in 0.05% gelatin in Eagle's minimal essential medium and delivered by aerosolization for 30 min to achieve a 90% lethal dose (LD90) to LD100 (∼100 TCID50 per mouse). Mouse health was followed for 21 d post infection, n = 15 mice per condition. Animals were weighed daily and sacrificed if they met euthanasia criteria, including signs of distress or loss of 20% pre-infection body weight. For in vitro infections, IAV (multiplicity of infection [MOI] of 1.0) with or without NF-B inhibitor IMD-0354 at 25 ng/L was added to cells in submerged monolayer and viral burden was assessed 24 h post infection.

To measure transcript levels of IAV nucleoprotein (NP) gene, samples were harvested in RNAlater and RNA was extracted using the RNeasy extraction kit. About 500 ng total RNA was reverse transcribed to cDNA by using an iScript cDNA synthesis kit and submitted to quantitative reverse transcription-PCR (RT-PCR) analysis with SYBR green PCR master mix on an Bio-Rad CFX Connect Real-Time PCR Detection System. Host 18S rRNA was similarly probed to determine relative expression of viral transcripts (24, 42) .

To perform analysis of interaction effects across both the Omics and feature level, OBIF performs a sequence of transformations to the unscaled original data matrix m that are tailored to meet the statistical assumptions needed for two-way ANOVA analysis of interaction terms across (43, 44) different Omics platforms: (i) background correction is applied if it has not been corrected for abundant negative or low signal intensities (45) (46) (47) ; (ii) a log 2transformation of count-(i.e. raw counts from sequencingbased platforms) or continuous-based (normalized counts, signal intensity and peak area from sequencing-, array-and spectrometry-based platforms respectively) data are applied to provide a Gaussian-like data distribution, if the expression values are not normally distributed (48, 49) ; and (iii) quantile normalization is applied to all datasets to equalize inter-sample variances (50, 51) . These transformations allow OBIF to have an analysis-ready (non-negative, nonzero values) (52, 53) data matrix m ( Figure 1D ) with expression values and of dimensions f x n, where f is the number of features as rows and n is the number of samples S as columns. The appropriate sample order in dimensions n of m is:

The subscripts denote the condition of the samples: exposed to neither factor (0,0), exposed to factor A alone (1,0), exposed to factor B alone (0,1) or exposed to both factors A and B (1,1). The superscripts represent the sample replicates from 1 to i within each of the four conditions.

To evaluate non-additive phenotypes of combined exposures as synergistic, we used an effect-based definition of synergy using the response additivity metric (43, 44, 54) . First, the observed effects of single and dual factor exposures are calculated from their main differences relative to the unexposed (sham) group: 0) ) for individual exposure to Factor A, E B = (S (0,1) -S (0,0) ) for individual exposure to Factor B, E AB = (S (1,1) -S (0,0) ) for combined exposure to Factor A and Factor B.

Hence, the observed combination effect (E AB ) is compared to the expected non-interacting additive effect of individual exposures (E ADD = E A + E B ) where E AB > E ADD indicates synergism, while E AB = E ADD indicates additivity, and E AB < E ADD indicates antagonism. Then, a significant interaction effect (P-value < 0.05) is systematically evaluated (1, 5) to discriminate true synergistic or antagonistic combination effects from a non-interacting cooperation between factors when E AB = E ADD (43, 44, 54) . Significant interaction effects were calculated from a multivariate Cox model with an interaction term for survival data, a two-way ANOVA for viral NP expression data and using full factorial analysis for feature expression data.

To evaluate significant interaction terms between factors at the Omics level within a dataset, OBIF performs a multiple linear regression across its feature expression values:

where the Omics-based interaction (I O ) regression formula is equivalent to a two-way ANOVA analysis where the intercept is referenced to the control samples (0) and returns a statistical summary of terms for the individual factor A (F A ), factor B (F B ) and their interaction (F A · F B ). Goodness of fit is calculated from the adjusted R 2 values, and overall significance is determined by the P-values of F-statistics that result from the I O regression.

OBIF performs full factorial analysis across the expression values of m to calculate the factorial effects of both expression and contrast analysis. The fixed-effect model fitted to the expression levels of features is:

where the expression levels of features (L f ) requires an intercept parameter to reference the control samples (β 0 ) and estimates ␤ coefficients for the main effects of individual factor A (F A ) and factor B (F B ) and their combination (F AB ). Subsequent model fitting for each treatment condition (Factor A, B and AB) involves ß 0 and their respective beta coefficient (ß 1-3 ). Therefore, calculation of main effects for each condition requires subtraction of ß 0 from their respective fitted model (ß 0 + ß 1 for Factor A, ß 0 + ß 2 for Factor B, or ß 0 + ß 3 for Factor AB) to accurately represent the group comparisons used for differential expression for selection of differentially expressed molecules (DEMs):

To adjust for multiple testing, the resulting P-values are controlled for false discovery rate (FDR) using the Benjamini-Hochberg procedure to select F A DEMs, F B DEMs and F AB DEMs, respectively. The Bonferroni procedure and Tukey's range test can be applied to datasets that require control for family-wise error rate. Similarly for detection of interactive DEMs (iDEMs), the remaining multifactor effects are estimated using contrast analysis for the simple main effects (SME) and their interaction effect:

Inter acti on e f f ect Bioinformatics, 2022, Vol. 4, No. 2 A significance threshold of P-value < 0.05 is used to determine if F AB DEMs are susceptible to any SME. iDEMs are defined as F AB DEMs with a statistically significant interaction term (P-value is < 0.05) during full factorial analysis that results in their synergistic or antagonist expression after combinatory dual exposures. Hence, iDEMs are more representative of synergistic or antagonistic statistical interactions than of molecular interactions between DEMs.

To evaluate non-additive expression levels of individual features, OBIF uses an Interaction Score (IS) where the absolute ratio of the Log 2 fold change with combined factors (FC AB ) over the linear sum of the Log 2 fold change of individual factors (FC A + FC B ) can identify antagonistic, synergistic or additive features using fold change principles where additivity is represented when FC AB = FC A + FC B . A subsequent log 2 -transformation can then quantify the size of the interaction effect by re-scaling the additivity threshold around 0 for additive features, with IS > 0 for synergistic and IS < 0 for antagonistic features.

Antagonistic and synergistic features (iDEMs) are differentiated from additive features with a significant interaction effect (P-value < 0.05) (44, 55) calculated during full factorial analysis.

All differentially expressed molecules (DEMs) were represented in heatmaps after hierarchically clustering using Ward's minimum variance method with Euclidean distances of Log 2 FC values to compute dissimilarity by rows (features) and by columns (samples). Column dendrograms were plotted to represent the distance between samples, vertical side bar colors summarize DEMs according to their and horizontal side bars colors represent sample types by factors. Color scale keys indicate the levels of feature expression with upregulation in red and downregulation in green. DEMs with F AB (Pam2-ODN) were clustered by principal component analysis based on the mean linear fold change difference to reveal the expression patterns biologically present across all factors F A , F B and F AB . Principal components 1 and 2 were used for plotting DEMs with F AB and the variability between features is marked in each axis. Expression profiles (EPs) were identified in the clusters for each individual feature according to Figure 3 C.

We compared OBIF's capacity to detect interaction effects at the level of individual features against a mixed-effect model (56) :

where the expression level of features in a mixed-effect model (L f-Mix ) is a function of the estimated ␤ coefficients for the fixed effects of individual factor A (F A ) and factor B (F B ) and their interaction (F A · F B ) with a random effect (1|S) for all sample conditions (S (0,0) , S (1,0) , S (0,1) , S (1, 1) ). After fitting the mixed-effect model, empirical Bayesian shrinkage is used to estimate the model fitness from the t-and F-statistics of the L f-Mix regression formula. The significance threshold of the interaction term uses a Pvalue < 0.05 to determine interaction between F A and F B .

To compare the performance of OBIF against other models and bioinformatic tools, we applied a beta-uniform mixture (BUM) model using the ClassComparison R package to distributions of P-values in order to compute the 'ground truth'. The BUM model estimates a ground truth reference for a large set of P-values by calculating the occurrence of true positive (TP), false positive (FP), true negative (TN) and false negative (FN) detections under the fundamental properties that such P-values are distributed uniformly and can be expressed as mixture of false null hypothesis and true null hypothesis (57, 58) . This approach has been validated to address the multiple-testing errors often found in large sets of P-values calculated from two-way ANOVA models (as OBIF), correlation analyses, Kruskal-Wallis tests and others (57, 59) . Using these estimates, we systematically compared their discrimination ability by calculating the following parameters:

Applying BUM estimates, we compared the model fitness of OBIF against other models or bioinformatic tools by calculating the following parameters: Then, a receiver operating characteristic area under the curve (ROC AUC) is calculated for each model or tool by plotting their sensitivity against their 1-specificity.

To provide biological interpretation of the full factorial analysis and classification of features, enrichment analysis was integrated in the pipeline (60) to determine candidate effectors and regulators, biological pathways and functional processes. OBIF allows export of feature lists and expression values of DEMs, EPs and iDEMs. These are formatted to be used with common bioinformatic enrichment tools, including GSEA (61,62), EnrichR (63, 64) and Ingenuity Pathway Analysis (IPA) (65) . Data presented here were initially analyzed with IPA (QIAGEN Inc.

https://www.qiagenbioinformatics.com/products/ ingenuitypathway-analysis) using core analysis with the expression values of DEMs, EPs, and iDEMs. Both gene and chemical Ingenuity Knowledge Base modules are used as enrichment reference, considering only experimentally observed confidence levels for identification of direct and indirect relationships. The thresholds of significance for canonical pathways, upstream analysis, diseases and functions, regulator effects and network analysis were selected from those with a z-score ≥ 2 for activation or ≤ -2 for inhibition, and < 5% false discovery rates for all predictions (65) . Validation of enrichment results of iDEMs from IPA was accomplished with EnrichR using the TRRUST Transcription Factor 2019 reference database and the top enriched terms were selected and clustered using their combined score (P-value + z-score) (63, 64) .

HBEC3-KT were grown to 80-100% confluence in 24-well plates and treated with PBS, Pam2, ODN or Pam2-ODN for the indicated durations. Measurements of DNA-binding activity of NF-B and AP-1 transcription factor subunits were performed from whole cell lysates were made using their respective TransAM Kit according to product directions. For signal detection, samples were read immediately for absorbance at 450 nm with reference wavelength at 655 nm on a microplate reader. Experiments were repeated in triplicate and statistical analysis was performed with unpaired Student's t test using GraphPad Prism 8.0 with a significance threshold of P-value < 0.05.

HBEC3-KT were grown to 80-100% confluence in 100 mm dishes and treated for 30 min with PBS, Pam2, ODN or Pam2-ODN with or without pretreatment with NF-B inhibitor IMD-0354 at 25 ng/l for 16 h. Cells were detached from the plate with a 5 min incubation at 37ºC degrees with 3 ml of Accutase to prevent additional activation of transcriptional activity. Cells were pelleted in individual 15 ml tubes at 500 g for 5 min and suspended in 500 l of eBioscience FOXP3 fixation/permeabilization buffer for 15 min at room temperature. Cells were stained with a LIVE/DEAD Fixable Near IR Dead Cell Dye and with a 1:1000 dilution of NF-B/p50 (E-10) Alexa Fluor 647, NF-B/RelA (F-6) Alexa Fluor 488, AP-1/cJun (G-4) Alexa Fluor 594 and AP-1/cFos (D-1) Alexa Fluor 546 conjugated antibodies for 1 h on ice and protected from light. After incubation, cells were pelleted and washed with 200 l of sterile PBS 4 times, then resuspended in 100 l sterile PBS. After the last wash, cells were pelleted and resuspended in 50 l of sterile PBS and nuclear DAPI staining at 0.5 g/ml was performed just prior to data acquisition on ImageStreamX MII. HBEC3-KT images were acquired using INSPIRE software on the ImagestreamX Mark II imaging flow cytometer (Amnis Corporation) at 60 × magnification, with lasers 405 nm (50.00 mW), 488 nm (200.00 mW), 561 nm (200.00 mW), 642 nm (200.00 mW) and side scatter (782 nm) (2.25 mW). About 10 000-25 000 images per sample acquired include a brightfield image (Channel 1 and 9), p65 Alexa Fluor 488 (Channel 2), c-Fos Alexa Fluor 546 (Channel 3), c-Fos Alexa Fluor 594 (Channel 4), side scatter (Channel 6), DAPI (Channel 7) and p50 Alexa Fluor 647 (Channel 11). The laser outputs prevented saturation of pixels in the relevant detection channels as monitored by the corresponding Raw Max Pixel features during acquisition. For image compensation, single color controls were stained with all fluorochromes and 500 events were recorded with each laser for individual controls. Fluorescent images were taken in all channels with brightfield LEDs and scatter lasers turned off to accurately capture fluorescence. Individual single-color control file was then merged to generate a compensation matrix and all sample files were processed with this matrix applied. After compensation for spectral overlap based on single color controls, analysis was performed, and individual cell images were created using IDEAS® software version 6.1. Cell populations were hierarchically gated first by single cells, then cells in focus, then negative selected for live cells, and finally as double positive for both DAPI and the transcription factor subunit of interest (Supplementary File S3D). The spatial relationship between the transcription factors and nuclear images was measured using the 'Similarity' feature in the IDEAS software to quantitate the mean similarity score in the cell populations per sample. A similarity score > 1 represents nuclear translocation, and the shift in distribution of nuclear translocation between two samples was calculated using the Fisher's Discriminant ratio (Rd value) (66) .

Statistical analyses were performed using Prism 8 (Graph-Pad, San Diego, CA, USA) and R. Analysis of survival was performed by calculation of Kaplan-Meier curves and multivariate Cox model was used for group comparisons relative to PBS. Analysis of viral NP expression was performed using a two-way ANOVA with post hoc Tukey analysis for paired comparisons that was adjusted for multiple testing. Analysis of DNA-binding activity in vitro was performed using a Student's t test for comparisons between two groups or using one-way ANOVA for comparison between multiple groups. Grouped data are shown as means ± standard error of the mean. To verify the statistical assumptions for each test, Gaussian distribution was evaluated with Saphiro-Wilk test, and equal variance between two samples was evaluated with F-tests, or for more than two samples with Barlett's or Levene's test. Simultaneous multiple outlier detection was performed using the robust regression and outlier removal (ROUT) method with a q value of 5% (maximum FDR). Treatment allocation of animals was randomized in the experiments, though assessment could not be blinded. A pre-specified minimum requirement of three biological replicates for in vitro studies and 10 for in vivo studies.

Data preprocessing and analysis with OBIF was performed in RStudio on both Mac-and Windows-based computers. For Mac-based analysis, R version 3.6.3 (2020-02-29) with x86 64-apple-darwin15.6.0 platform was used on a personal computer equipped with 2.7 GHz Quad-Core Intel Core i7 with 16GB of RAM under macOS 11.6 operating system. For Windows-based analysis, R version 4.1.0 (2021-05-18) with x86 64-w64-mingw32 platform was used on a desktop 6 NAR Genomics and Bioinformatics, 2022, Vol. 4, No. 2 computer equipped with 3.4 GHz Intel Core i-6700 CPU with 16GB RAM under Windows 10 Enterprise operating system. To provide a processing power reference for OBIF users in other computers, we use the benchmarkme R package (https://github.com/csgillespie/benchmarkme) to assess the calculation speed of matrix benchmark functions. We used the bm matrix cal lm() function to model a linear regression over a 5000 × 5000 matrix and obtained an elapsed time of 0.857 and 0.870 s in a Mac-and Windows-based computer, respectively. We used the bm matrix cal manip() function to create, transpose and deform a 5000 × 5000 matrix and obtained an elapsed time of 0.458 and 0.630 s in a Mac-and Windows-based computer, respectively.

Our laboratory's interest in synergistic interactions arises from our experience investigating single versus dual immunostimulatory treatments to prevent pneumonia (9, (22) (23) (24) (25) 27, 67) . As a demonstrative example, data are presented here from influenza A virus (IAV) challenges of different models following pretreatment with Pam2 alone, ODN alone or the Pam2-ODN combination (i.e. both treatments delivered together). When mice are challenged with IAV 24 h after the indicated inhaled treatments, we observed little increase in survival after the individual treatments compared to sham-treated control mice, whereas mice treated with the Pam2-ODN combination demonstrated profound antiviral protection ( Figure 1A ). Similarly, when isolated mouse lung epithelial (MLE-15) cells were challenged with IAV 4 h after pretreatment with the individual ligands, we observed no significant reductions in the viral burden relative to PBS-treated cells. However, cells pretreated with Pam2-ODN showed a significant reduction in viral nucleoprotein (NP) gene expression as assessed by qPCR relative to host 18s gene ( Figure 1B ). Comparing the effect of combination ligand treatment (E AB ) to the expected response additivity of the individual ligand treatments (E ADD = E A + E B ) (44) reveals synergistic effects on both in vivo survival benefits and in vitro viral clearance ( Figure 1C ). To better understand the molecular mechanisms driving such unanticipated synergy, we developed OBIF as a phenotype-driven model (6) to understand the mechanisms underlying observed favorable combination exposure interactions that mediate synergistic responses and outcomes.

To formally test whether the effect of combined factors (F AB : Pam2-ODN) is greater than the expected linear sum of its individual factors (F A : Pam2; F B : ODN), an initial 2-level 2-factor (22) factorial design is required (43, 44) . Our strategy adapts the traditional analysis of variance (ANOVA) approach into a model that links the empirical analysis of synergy (43, 44) with the high-throughput capacity and high-dimensionality of Omics datasets (68,69) (Fig-ure 1D ). As summarized in Figure 1E , OBIF integrates statistical and biological interactions in Omics data matrices from single versus dual factor exposures, leveraging Omics screening to promote discovery of the molecular drivers of synergy, and facilitating the biological validation of such regulators. The analytical pipeline is freely available as an R package at GitHub (www.github.com/evanslaboratory/ OBIF) (Table 1) . Naturally, the experimental approach to validating computational findings must be tailored to the individual tools and characteristics of the biological responses being studied.

OBIF allows sequential transformation of an unscaled original data matrix with background correction, Log 2transformation, and quantile normalization (Figure 2A ) in order to improve performance of an Omics-based interaction analysis between the studied factors ( Figure 2B ). This allows better adjusted R 2 and F-statistic P-values of the two-way ANOVA with improved detection of interaction terms (Supplementary File S1A) across multiple platform types. We further evaluated OBIF's overall model fitting using multiple standard test performance metrics calculated from 40 publicly available (4, 27, 55, 56, (70) (71) (72) (73) (74) (75) (76) (77) (78) (79) (80) Omics datasets (Supplementary File S5) that contain a broad range of samples (n = 4-30) and data sizes (161 to 217 507 features). Here, we found that OBIF's overall performance is best when the total sample size is n ≥ 12 (i.e. ≥3 biological replicates/condition) because that sample size adjusted R 2 values ≥ 0.6 (Supplementary File S1B), area under the receiver-operating characteristic (ROC) curve (AUC) analysis ≥ 0.8 (Supplementary File S1D) and diagnostic odds ratios (DOR) ≥10 (Supplementary File S1F).

Analysis of factorial effects in a data matrix from single versus combined factor exposures can statistically differentiate whether stochastic feature expression in a combination is correlated with the effect of an individual factor (simple main effect, SME) or their influence on each other (interaction effect) (81) (82) (83) . Based on this principle, OBIF performs full factorial analysis through paired comparisons of calculated ␤ coefficients in each condition (84) to discover main effects during expression analysis and multi-factor effects (SMEs and interaction effect) from contrast analysis ( Figure 2C) (85) . These statistical effects can be readily visualized ( Figure 2D Group comparisons to PBS used a multivariate Cox model for survival data, and a two-way ANOVA with post hoc Tukey analysis for viral NP expression data: * P < 0.05, ** < 0.01, *** P < 0.001, # P < 0.05 of interaction effect. Grouped data are shown as means ± 95% confidence interval for survival data and standard error of the mean for viral NP expression data. S1D). Similarly, OBIF outperforms other profiling tools commonly used to analyze multiple Omics data types, as demonstrated by superior values for area under the receiver operating characteristic (ROC) curve (AUC) and better false discovery rates, precision and recall (Supplementary File S1E) to select differentially expressed features from combination treatments with an improving overall performance in its diagnostic odds ratio as the sample size increases (Supplementary File S1F). Together, incorporation of full factorial analysis to dissect both main effects and hidden multi-factor effects allow enhanced analytical capacities over other available Omics profiling and processing tools ( Table 2) .

To investigate the mechanisms underlying Pam2-ODN synergy, we used OBIF to re-analyze previously published (27) lung homogenate transcriptomic data from mice inhalationally treated with single versus combination ligands (GSE28994). Despite unsurprising overlap of detected differentially expressed features, OBIF outperformed commonly used methods of differential expression analysis with increased R2 values and more significant P values in all treatment conditions (Supplementary File S2A) . Once feature expression is fitted ( Figure 2D ), this analysis identifies 3456 features as differentially expressed molecules (DEMs) Figure 3A ). Despite the fact that 52% (1617/3138) of DEMs were shared by Pam2-ODN and the individual ligands, enrichment analysis using IPA software revealed an over-representation of 11 canonical cellular immune response pathways that were activated by Pam2-ODN ( Figure  3B ). Of these 11 immune and host defense pathways, most were only enriched by Pam2-ODN treatment, with a few enriched to a lesser degree by Pam2, and none by ODN.

Rather than relying on potentially redundant DEM clusters, OBIF classifies F AB DEMs into eight expression profiles (EPs) that characterize cooperative and competitive biological interactions of individual factors ( Figure 3C ). EPs are defined by expression directionality (up-or downregulation) of individual features and are not biased by the expression analysis of DEMs. Cooperative EPs have accordant expression directionality induced by F A and F B (i.e. both individual exposures induce increased feature expression or both induce decreased expression), while competi-tive EPs have opposite directionalities induced by F A and F B . Among the cooperative EPs, concordant profiles result when F AB directionality corresponds with the single factor effects (EPs I and II), and discordant profiles occur when F AB directionality opposes the single factors (III and IV). Alternatively, among the competitive EPs, factor-dominant profiles are defined by F AB directionality correspondence with one factor (F A -dominant EPs V and VI; F B -dominant EPs VII and VIII). Principal component analysis ( Figure  3D ) demonstrates that concordant EPs I and II were the most abundant in our dataset, followed by Pam2-dominant EPs V and VI. This abundance of EPs I and II better emphasizes the cooperative effects of both factors than does conventional DEM clustering alone (Supplementary File S2B) .

Notably, enrichment analysis reveals that molecular effectors clustered by EPs correspond with Pam2-ODNinduced functions (Figure 3E) , suggesting a biological basis for both the synergy and feature selection. Specifically, we found that features in EPs I and V contributed to host survival functions. Considered from an organizational perspective, induction of resistance to infection at the organismal level correlated with features in concordant EP I, at the cellular level with concordant EP II, and by leukocytes with Pam2-dominant EP V.

Main effects determined significant DEMs per condition, while multi-factor effects explained whether Pam2-ODN DEMs and EPs resulted from SMEs and/or an interaction of individual ligands ( Figure 3F ). This analysis showed that most features in concordant EPs I and II are influenced by at least one multi-factor effect, while all features in discordant EPs III and IV are influenced by all multi-factor effects simultaneously. Not surprisingly, Pam2-dominant EPs V and VI and ODN-dominant EPs VII and VIII expression mainly results from their respective SMEs. This analysis also revealed that 67% (2116/3138) of Pam2-ODN DEMs are driven by the interaction effect of Pam2 and ODN as interactive DEMs (iDEMs). Thus, OBIF reconciled the biological interactions from EPs with the statistical interactions from multi-factor effects of Pam2-ODN.

Synergistic or antagonistic responses result from strong interaction effects between two factors in a combination (43, 54) . Hence, iDEMs integrate this principle during feature selection based on significant interaction effects between factors, allowing quantification of feature expression in a narrower set of DEMs. Using an interactions score (IS) that quantifies the effect size of the interaction effect relative to the additivity threshold, we can separate non-interactive DEMs (non-iDEMs) from true antagonistic (IS < 0) and synergistic (IS > 0) iDEMs ( Figure 3G ). This allows more focused enrichment analysis (Supplementary File S3A), in this case supporting NF-B/RelA and AP-1/cJun as key transcriptional upstream regulators of Pam2-ODN's interaction effect and synergistic expression ( Figure 3H) S3A,B). To visually represent the DEMs, EPs and iDEMs with the integrative analysis of biological and statistical interactions, OBIF generates a custom Circos plot that summarizes the main characteristics of these synergy-mediating feature clusters ( Figure 3I ).

Downstream analyses of DEMs, EPs and iDEMs have the capacity to discern the contributing roles of individual factors to a combination treatment regardless of the Omics platform being studied. To demonstrate the cross-Omics function of OBIF, we analyzed our reverse-phase protein array (RPPA) data from single-or dual-treated human lung epithelial cells (Supplementary File S2C) and were able to visually compare the prevalence of DEMs, EPs and iDEMs in transcriptomics and proteomics datasets (Supplementary File S2D). Interestingly, the use of EPs improved the analysis of hierarchical clustering given that the latter approach alone did not reveal distinctive gene clusters to explain the synergistic response while the abundance of concordant EPs I and II highlighted the contribution of ODN to the protective combination that might otherwise be overlooked, as enrichment analysis showed far fewer signaling pathways ( Figure 3B ) and had a greater clustering distance from Pam2-ODN samples (Supplementary File S2B) from the microarray. Similar to our previous findings with iDEMs ( Figure 3H ), transcription factors from many path-ways were involved, and NF-B and AP-1 transcription factors subunits remained central elements of the regulatory network analysis from our microarray results while our RPPA data identified the top phospho-signaling DEMs and cross-validated STAT3, NF-B/RelA and AP-1/cJun as transcriptional regulators involved in the Pam2-ODN signaling network (Supplementary File S2E).

Prompted by the foregoing results, we tested whether NF-B and AP-1 subunits were biologically relevant synergy regulators of Pam2-ODN-induced epithelial resistance. DNA-binding activity of a panel of NF-B and AP-1 family subunits in human lung epithelial cells (HBEC3-KT) after stimulation with Pam2-ODN confirmed that NF-B/RelA and AP-1/cJun subunit activation was strongly increased after 15 min of treatment without significant contribution of other subunits ( Figure 4A and Supplementary File S3C). Indeed, NF-B/RelA and AP-1/cJun subunits exhibited surprisingly similar activation kinetics after Pam2-ODN treatment, further supporting cooperation or coordination ( Figure 4B ). Investigating this co-activation of non-redundant transcriptional families, single-cell nuclear translocation of canonical NF-B (p50/RelA) and AP-1 (cFos/cJun) complexes was assessed in HBEC3-KT cells by imaging flow cytometry. We found that all transcriptional subunits exhibited an increased nuclear translocation (sim- 

Tests for interaction effects ilarity score > 2) after 15 min of Pam2-ODN treatment relative to the PBS-treated cells ( Figure 4C ). However, neither Pam2 nor ODN alone induced the same magnitude of nuclear translocation, whether assessed by change in similarity scores (Rd value) or by the percentage of translocated cells ( Figure 4D ) relative to PBS-treated cells.

To differentiate transcriptional cooperation from coincidental transcriptional activation after Pam2-ODN treatment, we assessed the Pam2-ODN-induced nuclear cotranslocation of NF-B (p50/RelA) and AP-1 (cFos/cJun) complexes in the presence or absence of NF-B inhibitor IMD-0354 (IMD). As expected, pre-treatment with IMD alone reduced the Rd Value and percentage of translocated cells for NF-B/RelA and NF-B/p50 subunits without significantly modifying the percentage of translocation for AP-1/cJun and AP-1/cFos subunits. However, NF-B inhibition with IMD also unexpectedly reduced the Pam2-ODN-induced similarity score shifts and nuclear translocation of AP-1 subunits, particularly of AP-1/cFos ( Figure  4E ). This indicates that NF-B inhibition impaired Pam2-ODN-induced AP-1 nuclear translocation, supporting the cooperative regulation of these two non-overlapping signaling pathways. Representative images shown in Figure  4F demonstrate that inhibition with IMD reduced Pam2-ODN-induced heterodimerization and nuclear translocation of NF-B (p50/RelA) and AP-1 (cFos/cJun) complexes. Further, we confirmed that disruption of this transcriptional cooperation was sufficient to impair the in- (ii) EPs, where inner sectors represent individual profiles (I to VIII) followed by their statistically significant F A -SME (green), F B -SME (orange) or F A ·F B interaction effect (pink); and (iii) iDEMs, represented by their synergistic (purple) or antagonistic (yellow) interaction scores. FC, fold change; DEMs, differentially expressed molecules; EPs, expression profiles; iDEMs, interactive DEMs; SME, simple main effect. n.s., non-significant threshold of Z-score (dashed line). ducible viral burden reduction seen with Pam2-ODN (Figure 4G ).

To demonstrate its generalizability, we used OBIF to interrogate the synergistic mechanisms of different datasets using a single-( Figure 5A ) or a multi-Omics approach (Figure 5B) . In our first case study, OBIF analyzed transcriptome data from untargeted total RNA-seq of an in vivo mouse model of K-ras mutant lung adenocarcinoma (CC-LR) (56) . In this model, the exposure factors are sex ( (G) Virus burden of mouse lung epithelial cells challenged with influenza A with or without NF-B inhibition. n = 4 samples/condition. Comparison to baseline activity used Student's t test for individual TF subunits, and one-way ANOVA for multiple TF timepoints: *, P < 0.05; **, P < 0.01; ***, P < 0.001; n.s., non-significant. Grouped data are shown as means ± standard error of the mean. (56) , but OBIF provides the first mechanistic support for that hypothesis. Additionally, network analysis of iDEMs revealed four previously unrecognized regulators (CCL28, CAGA, CAGB and CSF3R) that help explain the immunosuppressive molecular mechanisms that are permissive for the increased tumor burden in Male + Stat3 K.O. CC-LR mice.

In our second case study, OBIF analyzed both transcriptome data from an untargeted microarray and proteomic data from a targeted RPPA of an in vitro cell culture model of a chemotherapy-resistant adenosquamous lung carci-noma. In this model, cancer cells are exposed to a drug (F A : Volasertib, a polo-like kinase 1 inhibitor) and/or a cytokine (F B : TGF-ß). It has been shown that tumor growth is inhibited only by combined exposure (F AB : Volasertib + TGFß) (72) . Here, iDEMs upstream regulator analysis and hierarchical clustering from both platforms co-identified transcription factor cJun from antagonistic iDEMs as a central regulator of the observed exposure synergy. Further, building on the investigators' report that the cMet/FAK/SRC axis regulates tumor growth inhibition in this model (72) , network analysis of iDEMs yields downstream mechanistic insights into tumor growth regulation and epithelial to mesenchymal transformation through cJun that were not previously hypothesized. OBIF was similarly applied to additional datasets derived from Omics platforms including microarray (27) 

Synergistic and antagonistic interactions are common in nature and frequently promote efficacy of therapeutic interventions (1) (2) (3) (4) (5) . While synergy quantification from doseresponse data, combinatorial screening of molecule libraries, and other synergy prediction may suggest poten-tially beneficial conditions or treatments, they do not provide substantive insights into the underlying molecular mechanisms (1). Thus, synergy-mediating pathways cannot be strategically targeted in rational drug development only with synergy quantification or prediction models.

Our interest in synergy arose from our observations of the strikingly favorable interactions of one such empirically derived combination, Pam2-ODN. While we could easily quantify the superiority of protection conferred by the dual treatment, in the absence of a systems theory to interrogate synergistic mechanisms (1,5), we were limited in our capacity to use available Omics datasets to deduce the mechanisms mediating the synergy. This is important because, although this lack of mechanistic understanding does not limit the utility of the current combination, it precludes development of next generation interventions that more precisely (perhaps, more efficaciously) target the synergydriving pathways with fewer off-target (potentially toxic) effects. In contrast to prediction and quantification models, OBIF was developed as a synergy interrogation model with the explicit intent to investigate observed synergistic events instead of combined exposures with additive-only effects or phenotypes (a potential misuse) for which generalized additive models are already available (86, 87) . As such, it is inherently a phenotype-driven model that performs full factorial analysis on feature expression data from single vs. dual factor exposures to identify molecular clusters that reveal synergy-mediating pathways, functions and regulators.

Current synergy quantification methods and tools (34) (35) (36) use accepted statistical approaches to detect drug-drug interactions (34) (e.g. Loewe additivity, Bliss additivity, Chou-Talalay Method, response surface methods, isobologram analysis, etc.). These approaches use dose-effect metrics to identify unestablished interactions (screening) that cannot be easily applied to Omics-based datasets because of the large number of observations and resources needed, while OBIF was designed to investigate mechanisms of observed instances of synergy (phenotype-driven). Although OBIF is not designed to predict toxicities, as do some of these synergy quantification tools, enrichment analysis of cytotoxic pathways or functions may also provide insights into the safety and tolerability of certain combined exposures. Further, these approaches require a series of observations derived from dose-response data (continuous variables) to conform to their required principles for synergy detection, face validity and formal validity (34, 43, 44, 54, (88) (89) (90) , whereas OBIF can be applied to beyond continuous variables (i.e. dichotomous or categorical factors). This principle expands the potential applications of OBIF to study multiple factor classes (i.e. cytokines, cell lines, drugs, environmental factors, genetic alterations, gender differences, time), system levels (i.e. in vitro, in vivo, ex vivo, clinical) and both targeted and untargeted platform formats (i.e. array-based, sequencing-based, mass spectrometry-based) (Supplementary File S5).

Other Omics-based models rely on molecular interactions (30), pathway- (28, 32) or network-based analysis (12, 28, 31) , and drug response similarity scores (29, 33) to predict synergistic drug combinations. While these predictive strategies have facilitated the discovery of some novel drug combinations (1, 7, 10) , they largely depend on modeling of drug targets, drug actions, and organism-or diseasespecific data inputs that limit their capacity to explore potential mechanisms and hinder their transferability to other Omics platforms or research fields. In contrast, OBIF is capable of revealing mechanistically relevant molecular clusters (DEMs, EPs and iDEMs) directly from single-or multi-Omics expression data following combined exposures and can be easily transferred for downstream analysis with a wider span of bioinformatic tools. Thus, while established synergy prediction and quantification tools can offer important contributions to the discovery of drug-drug interactions, the objectives of a synergy interrogation tool are entirely different. In contrast to these types of tools, OBIF is applied only after synergy has been observed in any of a wide range of biological exposures and contexts, and focuses exclusively on the mechanistic interrogation of the drivers of synergistic interactions rather than on detecting undiscovered/unanticipated synergistic combinations de novo, and is capable of achieving this via independent single Omics dataset analysis with high correspondence when OBIF is serially applied to both cross-platform or multi-Omics analysis.

Using Pam2-ODN datasets as demonstrative examples, OBIF identified unanticipated transcriptional cooperation between non-redundant transcription factors, NF-B (p50/RelA) and AP-1 (cFos/cJun) complexes, as a molecular mechanism of inducible synergistic protection against IAV. Thus, by facilitating understanding of combined factor exposures in terms of the individual components, a computational discovery allowed focused experimental validation of a discrete, novel mediator of a synergistic biological re-sponse. Perhaps as importantly, the computational analyses were accomplished by integration of data from different Omics platforms, different specimen types, and even different host species ( Figure 5 and Supplementary File S4). In all three detailed examples, OBIF revealed molecular mechanisms that were congruent with previously published findings and suggested novel mechanistic hypothesis for experimental validation. These applications already include single ( Figure 3A -H) and multi-Omics (Supplementary File S2B-E) tested scenarios for Pam2-ODN's specific RelA and cJUN interaction ( Figure 4A -F) and broad antiviral response ( Figure 4G ). These applications also include single ( Figure 5A ) and multi-Omics ( Figure 5B ) untested scenarios where specific ( Figure 5B ) and broad ( Figure 5A ) molecular drivers and mechanisms have been discovered using OBIF. Thus, OBIF can be used to rationally focus investigation into phenotypically observed synergy in a broad array of contexts, including currently unexplored synergistic gene mutation or copy number alterations from mapped reads in DNA copy number variation sequencing (CNV-seq) or from intensity signals obtained from array-based comparative genomic hybridization (aCGH) (91) .

Unlike most 2 2 designs, OBIF dissects factorial effects of combined factor exposures through full factorial analysis of feature expression data in a single processing step better than conventional differential expression pipelines (i.e. LIMMA, DESeq2, MetaboAnalyst) that are limited to a few platform-specific Omics dataset. This allows simultaneous identification of DEMs directly from main effects of single or combined factors, overcoming pairwise comparisons to control and repetitive analysis of each condition. While this simultaneous identification of DEMs could alternately be performed with a mixed-effect model, we showed how the mixed-effect model is suboptimal to detect interaction effects at the level of individual features and to facilitate iDEM selection when compared with full factorial analysis. Additionally, clustering by DEMs, EPs and iDEMs improves the specificity of enrichment analysis to disentangle the signaling pathways, functions and regulators of this synergistic combination and to capture their specific driving features. Further, quantification of multi-factor effects (SMEs and interaction effects) reveals whether particular features, molecular clusters or functions enriched by combined exposures are the result of individual factors or their crosstalk.

These statistical relationships have biological analogues that are integrated during feature discovery of DEMs, EP and iDEMs. Because each study of combined exposures using Omics datasets is unique, fixed cut-offs are not appropriate for every application. These statistical effects are visualized in volcano and Q-Q plots during model fitting to tailor each application. In fact, concordant EPs I and II rescued the underrepresentation of ODN observed in distance-based clustering and enrichment analysis. Further, iDEMs derived from features with significant interaction effects allow focusing discovery on molecular regulators and the calculation of interaction scores allows quantification of their synergistic or antagonistic expression. Thus, unlike most systems models of synergy, OBIF facilitates integrative analyses of biological and statistical interactions that are easily discoverable and interpretable through molecular

Systematic synergy modeling: understanding drug synergy from a systems biology perspective

Emerging concepts for immune checkpoint blockade-based combination therapies

Combination therapy is the new gene therapy?

Comparative metabolomics reveals key pathways associated with the synergistic killing of colistin and sulbactam combination against multidrug-resistant acinetobacter baumannii

Development of fangjiomics for systems elucidation of synergistic mechanism underlying combination therapy

Multi-omics approaches to disease

Empowering biologists with multi-omics data: colorectal cancer as a paradigm

A mechanistic pan-cancer pathway model informed by multi-omics data interprets stochastic cell fate responses to drugs and mitogens

Safety, tolerability, and biomarkers of the treatment of mice with aerosolized Toll-like receptor ligands

Network-based drug discovery by integrating systems biology and computational technologies

DeepSynergy: predicting anti-cancer drug synergy with deep learning

Driver network as a biomarker: systematic integration and network modeling of multi-omics data to derive driver signaling pathways for drug combination prediction

CDA: combinatorial drug discovery using transcriptional response modules

Target inhibition networks: predicting selective combinations of druggable targets to block cancer survival pathways

Estimates of the global, regional, and national morbidity, mortality, and aetiologies of lower respiratory tract infections in 195 countries: a systematic analysis for the global burden of disease study

Global, regional, and national estimates of pneumonia morbidity and mortality in children younger than 5 years between 2000 and 2015: a systematic analysis

Diagnosis and treatment of adults with Community-acquired pneumonia. An official clinical practice guideline of the american thoracic society and infectious diseases society of america

Global and regional burden of hospital admissions for pneumonia in older adults: a systematic review and meta-analysis

Covid-19 -Navigating the uncharted

Conducting clinical trials in outbreak settings: points to consider

Clinical characteristics of coronavirus disease 2019 in china

Lung epithelial cells are essential effectors of inducible resistance to pneumonia

Synergistic interactions of TLR2/6 and TLR9 induce a high level of resistance to lung infection in mice

Inducible lung epithelial resistance requires multisource reactive oxygen species generation to protect against viral infections

Inducible lung epithelial resistance requires multisource reactive oxygen species generation to protect against bacterial infections

Inducible epithelial resistance against coronavirus pneumonia in mice

Synergistic TLR2/6 and TLR9 activation protects mice against lethal influenza pneumonia

Predict effective drug combination by deep belief network and ontology fingerprints

A simple gene set-based method accurately predicts the synergy of drug pairs

In-silico prediction of synergistic anti-cancer drug combinations using Multi-omics data

Predicting drug synergy for precision medicine using network biology and machine learning

Synergy from gene expression and network mining (SynGeNet) method predicts synergistic drug combinations for diverse melanoma genomic subtypes

Stratification and prediction of drug synergy based on target functional similarity

SynergyFinder 2.0: visual analytics of multi-drug combination synergies

SynergyFinder: a web application for analyzing drug combination dose-response matrix data

Combenefit: an interactive platform for the analysis and visualization of drug combinations

Full factorial analysis of mammalian and avian influenza polymerase subunits suggests a role of an efficient polymerase for virus adaptation

Design of Experiments for Engineers and Scientists

Computational Phytochemistry

Reverse phase protein array: validation of a novel proteomic technology and utility for analysis of primary leukemia specimens and hematopoietic stem cells

A technical assessment of the utility of reverse phase protein arrays for the study of the functional proteome in Non-microdissected human breast cancers

The ultimate qPCR experiment: producing publication quality, reproducible data the first time

The statistics of synergism

Analysis of drug combinations: current methodological landscape

Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis

Compositional uncertainty should not be ignored in high-throughput sequencing data analysis

Compositional Data Analysis: Theory and Applications

To transform or not to transform: using generalized linear mixed models to analyse reaction time data

A phylogenetic transform enhances analysis of compositional microbiota data

A protocol to evaluate RNA sequencing normalization methods

Comparison of normalization and differential expression analyses using RNA-Seq data from 726 individual drosophila melanogaster

Harmonization of quality metrics and power calculation in multi-omic studies

A field guide for the compositional analysis of any-omics data

Understanding synergy

Synergistic gene expression during the acute phase response is characterized by transcription factor assisted loading

Sex specific function of epithelial STAT3 signaling in pathogenesis of K-ras mutant lung cancer

Estimating the occurrence of false positives and false negatives in microarray studies by approximating and partitioning the empirical distribution of P-values

Applications of beta-mixture models in bioinformatics

Sources of variation in false discovery rate estimation include sample size, correlation, and inherent differences between groups

Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, cytoscape and enrichmentmap

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

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

Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool

Enrichr: a comprehensive gene set enrichment analysis web server 2016 update

Causal analysis approaches in ingenuity pathway analysis

Simultaneous assessment of NF-B/p65 phosphorylation and nuclear localization using imaging flow cytometry

Inhaled innate immune ligands to prevent pneumonia

The role of analytical frameworks for systemic research design, explained in the analysis of drivers and dynamics of historic land-use changes

Multi-level and hybrid modelling approaches for systems biology

Genome-Wide CRISPR-Cas9 screening identifies NF-kappaB/E2F6 responsible for EGFRvIII-Associated temozolomide resistance in glioblastoma

Integrated single-cell and bulk gene expression and ATAC-seq reveals heterogeneity and early changes in pathways associated with resistance to cetuximab in HNSCC-sensitive cell lines

Non-canonical cMet regulation by vimentin mediates plk1 inhibitor-induced apoptosis

Pneumonia recovery reprograms the alveolar macrophage pool

Unique genetic responses revealed in RNA-seq of the spleen of chickens stimulated with lipopolysaccharide and short-term heat

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

Synergistic gene expression signature observed in TK6 cells upon co-exposure to UVC-Irradiation and protein kinase C-Activating tumor promoters

Untargeted metabolomics analysis reveals key pathways responsible for the synergistic killing of colistin and doripenem combination against acinetobacter baumannii

Metabolic analyses revealed time-dependent synergistic killing by colistin and aztreonam combination against multidrug-resistant acinetobacter baumannii

Drought and plant litter chemistry alter microbial gene expression and metabolite production

C9orf72 suppresses systemic and neural inflammation induced by gut bacteria

Gene expression patterns combined with bioinformatics analysis identify genes associated with cholangiocarcinoma

Combination of gene expression patterns in whole blood discriminate between tuberculosis infection states

Use of a combined gene expression profile in implementing a drug sensitivity predictive model for breast cancer

Beyond the one-way ANOVA for 'omics data

A Comprehensive Guide to Factorial Two-Level Experimentation. SpringerNY

Trajectory-based differential expression analysis for single-cell sequencing data

Negative binomial additive model for RNA-Seq data analysis

Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors

Searching for drug synergy in complex dose-response landscapes using an interaction potency model

Drug combination studies and their synergy quantification using the chou-talalay method

CNV-seq, a new method to detect copy number variation using high-throughput sequencing

The opinions expressed are those of the authors and do not represent views of these funding agencies. J 

clusters representing the complex dynamics of dual exposures.OBIF is available as an open-source R package with a semi-automated pipeline to facilitate its broad application to unscaled original data from various targeted and nontargeted Omics platforms, factor classes and biological systems. Compared to other commonly used analysis tools and modeling methods, OBIF demonstrates superior differential expression analysis for combined factor exposures and detection of interaction effects at the whole Omics and feature level. We have shown that OBIF can be fitted to perform full factorial analysis and that it adequately identifies DEMs, EPs, iDEMs and their attendant values and scores to promote discovery of molecular drivers of synergy in multiple, diverse datasets.In summary, OBIF provides a phenotype-driven systems biology model that allows multiplatform dissection of molecular drivers of synergy. And, we encourage the application of OBIF to provide holistic understanding in research fields where greater-than-additive beneficial combinations remain understudied.

The codes generated during the current study are available in the EvansLaboratory GitHub repository for the OBIF R Package (www.github.com/evanslaboratory/ OBIF) and OBIF R Code Extensions (www.github.com/ evanslaboratory/Extensions). The datasets generated during this study are available in the EvansLaboratory GitHub repository for the Reverse-phase protein array data (www.github.com/evanslaboratory/Datasource). Further information and reasonable requests for resources, data and reagents used in this work should be directed to and will be fulfilled by the corresponding author, Scott E. Evans (seevans@mdanderson.org).The