Principled multi-omic analysis reveals gene regulatory mechanisms of phenotype variation

Abstract

Recent studies have analyzed large-scale data sets of gene expression to identify genes associated with interindividual variation in phenotypes ranging from cancer subtypes to drug sensitivity, promising new avenues of research in personalized medicine. However, gene expression data alone is limited in its ability to reveal cis-regulatory mechanisms underlying phenotypic differences. In this study, we develop a new probabilistic model, called pGENMi, that integrates multi-omic data to investigate the transcriptional regulatory mechanisms underlying interindividual variation of a specific phenotype-that of cell line response to cytotoxic treatment. In particular, pGENMi simultaneously analyzes genotype, DNA methylation, gene expression, and transcription factor (TF)-DNA binding data, along with phenotypic measurements, to identify TFs regulating the phenotype. It does so by combining statistical information about expression quantitative trait loci (eQTLs) and expression-correlated methylation marks (eQTMs) located within TF binding sites, as well as observed correlations between gene expression and phenotype variation. Application of pGENMi to data from a panel of lymphoblastoid cell lines treated with 24 drugs, in conjunction with ENCODE TF ChIP data, yielded a number of known as well as novel (TF, Drug) associations. Experimental validations by TF knockdown confirmed 41% of the predicted and tested associations, compared to a 12% confirmation rate of tested nonassociations (controls). An extensive literature survey also corroborated 62% of the predicted associations above a stringent threshold. Moreover, associations predicted only when combining eQTL and eQTM data showed higher precision compared to an eQTL-only or eQTM-only analysis using pGENMi, further demonstrating the value of multi-omic integrative analysis.

Figure 1.

Plate diagram of pGENMi model. A latent variable z_g represents whether a gene g mediates the influence of a TF on phenotype, and its enclosing rectangle denotes G such genes. p_g denotes the TWAS P-value between the gene's expression and phenotype variation for each of the G genes. If z_g = 1, we expect an enrichment for significant TWAS P-values, and p_g is modeled by a Beta distribution parameterized by α; otherwise, p_g is modeled as being distributed uniformly in [0,1]. We also observe one or more lines of evidence supporting the TF's influence on the expression of each gene g, such as the existence of a cis-eQTL within a ChIP peak of the TF near that gene. These “regulatory evidences” are denoted by the binary variables r_gm, and there exist M such types of evidence (m = 1 … M), with relative weights w_m. These evidences combine in a logistic function to determine Pr(z_g = 1). The weights w_m are learned over all genes, and as such, are shown outside of the rectangle enclosing G. The H variable indicates whether w_m is free or restricted to zero (null model) for hypothesis testing (see “Probabilistic graphical model” section in Methods).

Figure 2.

The 90 (TF, Drug) associations predicted by eQTL + M with LLR ≥ 3, colored by LLR. Labels indicate whether eQTL-only or eQTM-only analysis corroborated the prediction at LLR ≥ 1.74 (labels “S” and “M,” respectively). An “X” indicates both eQTL and eQTM analysis supported the prediction, while “–” indicates the prediction is unique to eQTL + M.

Figure 3.

All 25 experimentally tested (TF, Drug) dosage-response curves. Red outlines indicate significant shifts in cytotoxicity between siRNA negative and siRNA TF conditions. Curves with a gray background are eQTL + M predictions, while those with a black background are negative controls. We confirmed seven out of 17 predictions and only one out of eight negative controls.