Bayesian approach to transforming public gene expression repositories into disease diagnosis databases

Results

Construction of a Disease Diagnosis Database.

Our analysis scheme is sketched in Fig. 1. The first phase concerns data preprocessing and it involves two tasks: (i) standardizing the expression data to remove platform or laboratory differences, and (ii) transforming the heterogeneous phenotype information, embedded in texts, into a workable format.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Major steps of the disease diagnosis system: (1) Preprocess the public microarray repositories to build the diagnosis database with standardized expression and phenotype data. (2) Diagnose a query profile via a two-stage Bayesian approach: at the first stage, we build Bayesian classifiers for each UMLS concept; at the second stage, we integrate the individual predictions with a Bayesian network model to allow collaborative error-correction over all classes in the hierarchy (red nodes represent diagnosed disease concepts).

For the first task, we performed the standardization by ranking the expression values within each array profile, and then taking the logarithm of the expression rank ratio between two arrays, one of a disease and another of normal conditions, within the same dataset (14). The resulting log-rank-ratio vector is termed a standardized profile; its components reflect the level and direction of differential expression in disease-related genes. Given two standardized profiles, we measure their similarity by Pearson correlation. This measure is favored due to its sensitivity to large absolute values, which in our case correspond to highly differentially expressed genes. Our standardization is in principle similar to the Gene Set Enrichment Analysis (13). The rationale is that differentially expressed genes are likely to carry the most stable information on disease characteristics regardless of differences in platform or lab. With our method, researchers need not choose a threshold of significance for differentially expressed genes. Extensive testing has demonstrated the effectiveness of our method (SI Text). Note that a typical GEO dataset consists of several subsets, each subset having a group of replicated samples. We term the set of standardized profiles derived from all replicated arrays in a pair of disease and normal subsets, a standardized dataset. Naturally, the standardized profiles within a standardized dataset are considered replicates. For this initial study, we focus on those GEO datasets that contain at least one pair of disease and normal subsets. We note that the requirement of normal subsets poses a potential limitation for our standardization approach. Further discussions on our approach and its possible alternatives are in SI Text.

To standardize the disease information provided with microarray datasets, we use the Unified Medical Language System (UMLS) (15). Its associated text-mining tools (MetaMap) (16) allow for the automatic extraction of phenotype concepts from text annotations (17 and 18). Note that each microarray dataset has phenotype descriptions at both the dataset level and the subset level. In our Bayesian diagnosis analysis, we give the subset-level annotations a higher credence due to their stronger association with individual samples.

This preprocessing phase results in a new database consisting of standardized expression profiles and disease concepts (i.e., UMLS concepts). Hereafter, we refer to this collection as the “disease diagnosis database.” More details on constructing this database are in Methods and SI Text.

A Bayesian Framework for Automated Disease Diagnosis.

We formulate the task of automated disease diagnosis as a hierarchical multilabel classification problem. The goal is to place a query profile into one or multiple UMLS disease classes, following the hierarchical UMLS disease taxonomy. We introduce a two-stage Bayesian learning approach: first a classifier is built for each disease class, and then the predictions of individual classifiers are combined to allow collaborative error correction across classes in the hierarchy. We summarize the key steps of this method below; a complete description is in SI Text.

Given a query array to be diagnosed and a control array profiling a normal sample (of the same lab and platform), we derive the log-rank-ratio vector x; this is the standardized query profile. Let Q_x,k = 1 when x is diagnostic of the UMLS concept U_k, and Q_x,k = 0 otherwise. That is, Q_x,k is a binary label indicating membership of x in the disease class k. To build a Bayesian classifier for disease class k, it is equivalent to derive the posterior distribution of Q_x,k given the information including: (i) s = {s_x,i,i = 1,…,M}, where M is the number of standardized datasets in our database, and s_x,i is a list of similarity scores (i.e., Pearson correlation coefficients) quantifying the similarities between x and the profiles in the ith standardized dataset. (ii) e = {e_i,k,i = 1,…,M}, where e_i,k tells whether the UMLS concept U_k is an annotation of the ith standardized dataset (see Methods). We define e_i,k = 2 if U_k occurs at the subset level (see Methods), e_i,k = 1 if it occurs at the dataset level, and e_i,k = 0 if the annotation does not occur.

The relationship between Q_x,k, e and s is illustrated in Fig. 2. It is obvious that e provides no useful information to infer the value of Q_x,k unless the similarity scores s are also known. That is P(Q_x,k|e) = P(Q_x,k). In these terms, the target Bayesian posterior can be expressed as [1]The prior P(Q_x,k) can be empirically estimated from the database. The computation of P(s|Q_x,k,e) is more involved due to the complex properties of e and s. For instance, the UMLS annotations e include items at the dataset and subset levels, and may also suffer from text-mining errors. To take into account such complexities and facilitate the modeling of P(s|Q_x,k,e), we introduced a set of latent binary random variables T = {T_i,k} whose values are not observable but can be inferred from e: T_i,k = 1 when the ith standardized dataset is related to the UMLS concept U_k, and T_i,k = 0 otherwise. (The principles for defining P(T|e) are in Methods.) Accordingly, we can express the target posterior in Eq. 1 as [2]We now need to determine P(s|Q_x,k,T). By assuming independence among the standardized datasets, we can decompose P(s|Q_x,k,T) into the computation of individual terms P(s_x,i|Q_x,k,T_i,k), where s_x,i is a vector denoting the similarity scores between x and the profiles in the ith standardized dataset(i = 1,…,M). Due to the difficulty in modeling and deriving P(s_x,i|Q_x,k,T_i,k) (SI Text), we alternatively considered the ratio P₁(s_x,i)/P₀(s_x,i), referred to simply as “P₁/P₀” hereafter. P₁ denotes the distribution of s_x,i when the query x and the ith standardized dataset are associated with a common disease, and P₀ the distribution when they are not. We modeled P₁/P₀ by a log-linear regression: log(P₁/P₀) = λ₀ + λ₁ × Mean(s_x,i) + ϵ_i, where λ₀ and λ₁ are estimated independently of the query profile and ϵ_i is a Gaussian error term. Further details are in Methods and SI Text.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Information chart for the posterior inference of Q_x,k. We wish to estimate the probability that a query profile x is diagnosed with the UMLS concept U_k, given e_i,k and s_x,i with i = 1,…,M.

Putting all the details together, we can infer Q_x,k (SI Text). This process is repeated for all UMLS disease classes. Those UMLS concepts with the posterior probability above a predefined threshold are considered significant, and generate an initial set of diagnostic annotations for the query profile x. Note that in practice, the query data may contain replicates of disease and normal arrays. In that situation, we derive the log-rank ratio vectors from all possible pairs of query disease and normal samples, and x represents a list of replicated standardized query profiles. s_x,i will then include all the similarity scores between every replicated standardized query profile and every profile in the ith standardized dataset. The inclusion of replicates can enhance the robustness of the P₁(s_x,i)/P₀(s_x,i) estimation, but all other procedures remain the same. We also note that in Bayesian analysis, the effects of the prior distribution and the significance threshold selection tend to vanish as the data accumulate. We have also demonstrated that our Bayesian classifiers, built by carefully modeling the specific properties of noisy data such as e and s, outperformed Support Vector Machine, the most commonly used classification method (SI Text).

Next, we use the UMLS hierarchical disease taxonomy to leverage the predictions made for individual UMLS concepts. In particular, we exploit a Bayesian network model defined on the UMLS hierarchy to resolve inconsistencies in the initial set of diagnostic predictions (19) (details are in SI Text). Given the high level of information exchange and integration among the UMLS concepts along the disease hierarchy, this procedure is expected to improve the accuracy and robustness of the diagnosis.

Performance Assessment of the Automated Disease Diagnosis System.

To validate our framework, we used an initial set of GEO microarray datasets containing at least one disease subset and at least one normal subset. In total, we collected 9,169 microarray experiments and constructed 110 disease classes, each containing from 3 to 62 standardized datasets. The 110 classes covered a wide spectrum of diseases: cardiovascular disease, neoplasms, CNS disorders, skin disorders, and metabolic diseases, to name a few.

Using the leave-one-out cross-validation approach detailed in Methods, our diagnoses achieved an overall accuracy of 95% (precision 82% and recall 20%). The recall rate of 20% is comparable to that observed in an analogous hierarchical multilabel classification problem for predicting gene functions. In the mouse model, the best performance in such an application was achieved with a recall rate of 20% and a precision of 41%.(20). Varying the threshold for classification, we plot the precision and recall curves in Fig. 3A. Not surprisingly, the performance of our method is significantly enhanced after applying collaborative error correction along the disease hierarchy. An example diagnosis in the context of the UMLS hierarchy is depicted in Fig. 3B. The prediction results for a subset of prevalent diseases are listed in Table 1.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Validation results and case examples. (A) Precision-recall plots by pooled disease classes. The blue curve shows the performance after Stage I diagnosis, and the red curve shows the final performance after Stage II refinement. (B) An example illustrating the error correction by the Stage II refinement. The query profile studies uterine leiomyomas obtained from fibroid afflicted patients (GDS484). The profile is annotated with four concepts by UMLS text mapping: Connective/Soft Tissue Neoplasm, Muscle tissue neoplasm, fibroid tumor, and uterine fibroids. The Stage I diagnosis predicted four concepts (red nodes) with one false positive (lymphoblastic leukemia), and one false negative (uterine fibroids). The false positive prediction is later corrected by Stage II refinement. (C) The figure presents the 110 disease classes and their hierarchical relationships. The red nodes represent diagnosed disease concepts for GDS563: (1) Nervous system disorder (2) Neuromuscular diseases (3) Myopathy (4) Musculoskeletal diseases (5) Congenital, Hereditary, and Neonatal diseases and abnormalities (CHNDA) (6) Genetic diseases, inborn (7) Genetic diseases, x-linked (8) Muscular disorders, atrophic (9) Muscular dystrophies (10) Muscular Dystrophy, Duchenne. (D) The prediction performance decreases with the data reduction.

We further exemplify the performance of our approach using the NCBI GEO dataset GDS563. This dataset was produced to identify modifying factors and pathogenic pathways involved in Duchenne Muscular Dystrophy (DMD). It consists of 24 microarrays from two subsets: 12 from DMD patients and 12 from unaffected control patients, which form a standardized dataset with 12 × 12 = 144 replicated standardized profiles. Masking the known phenotypic annotations and querying this standardized dataset against our database, we predicted 10 UMLS concepts (see Fig. 3C). All of the 10 predictions, positioned coherently along the UMLS hierarchy, agree 100% with the known annotations. Another example with less perfect yet more typical prediction performance is shown in SI Text.

A closer examination of the results shows further interesting features of our method. One example comes from the result for a query profiling the T-cells of HIV patients (GDS2649). Even though HIV is not included in the 110 disease classes of our diagnosis database due to the lack of sufficient training data, we obtain the relevant concept RNA virus infection that can describe the characteristic of the HIV disease. This implies that our system can not only diagnose known diseases, but may also identify important features of understudied or unknown diseases.

In general, the prediction performance on individual disease classes increases with the number of datasets in the class. For example, among disease classes containing only three datasets, the best precision achieved was 41% with a recall of 23%. For disease classes containing seven datasets, the same precision (41%) was achieved when the recall was 43%. In fact, for classes with seven datasets, the best achieved precision was 97% with 33% recall. To further confirm the important role of class size, we randomly reduced the number of datasets in the disease diagnosis database by 20%, 40%, or 80%. Fig. 3D demonstrates that both precision and recall significantly increase with the number of datasets. This behavior highlights the advantage of multiple dataset integration and demonstrates that the power of our approach can increase significantly with the continued growth of public gene expression repositories.

Construction of a Disease-Drug Connectivity Map.

Our approach can be generalized from disease diagnosis to building disease/phenotype connections when the query’s phenotype is known. The connections between drugs and diseases are of special interest. These can be discovered by applying our diagnosis system to queries involving drug (or small molecule) treatments. The diagnosis results would link the queries (involving drug-treatment effects) into one or more well-known disease classes, facilitating the establishment of unique links between diseases and drugs.

We constructed a connectivity map using 1,248 queries characterizing the phenotypic differences between “drug treated” and “untreated” subjects (these 1,248 queries were not in our established disease diagnosis database). We only kept predicted diseases that were not already included in the query’s annotations. A new drug-disease connection was considered significant if queries concerning a given drug treatment were predominantly classified into the same disease class (see Methods and SI Text for details). In total, we found 234 significant drug-disease links, unique to the queries’ phenotype description, connecting 99 drug concepts to 43 disease classes (Fig. 4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Disease-drug connectivity map. The map contains 234 significant connections between 99 drug concepts (pink nodes) and 43 disease concepts (blue nodes). (A) The network structure of the connectivity map. (B) Close-up of the Doxorubicin subnetwork. (C) Close-up of the obesity subnetwork.

To verify that the linked diseases and drug treatments truly share common molecular mechanisms, we used an independent resource: the Online Mendelian Inheritance in Man database (OMIM) (21). For each disease or drug, we compiled a list of associated genes based on the knowledge in OMIM (SI Text). We then assessed the statistical significance of the “link” between a disease and a drug based on the intersection of the disease-related genes and the drug-related genes using the hypergeometric test. Strikingly, 42.4% of the significant drug-disease links identified emerge also as statistically significant by this criterion with hypergeometric p value ≤ 0.05.

On close examination, Fig. 4A (see SI Text for a detailed version) reveals many known drug side effects as well as some unique associations. To take an interesting example, we found that the anticancer drug, doxorubicin, is linked to several diseases besides cancer/tumor (Fig. 4B). A predicted connection to rheumatoid arthritis confirms the drug’s newly proposed role as an antiarthritic agent (22); a connection to skin disorder reflects its known side effect of skin erupts (23). Most notably, the connection to cardiovascular diseases points to doxorubicin’s potentially fatal toxicity to the heart muscles, which is cumulative over the patient’s lifetime (24). This dosage-limiting cardiotoxicity, whose mechanism is not yet understood, has so far severely limited the usage of doxorubicin. Previous work has suggested that doxorubicin may break down the myofilament protein titin, leading to myocyte cell death (25). Our prediction, based on a global differential gene expression comparison, suggests a more fundamental transcriptional mechanism: doxorubicin may trigger cardiomyocyte-specific expression signatures. In fact, two transcription regulators, GATA4 and CARP, have already been implicated in such cardiomyocyte changes (26). Our result thus points out a unique direction in the design of cardioprotection strategies that would allow wider application of this effective oncologic agent.

Also of interest are the numerous links between “obesity” and anticancer drugs shown in Fig. 4C. Cancer and obesity have both become increasingly prevalent, and a physical basis for this relationship can be found in certain shared key molecules such as insulin, PDGF, VEGF, and cytokines. Accumulating evidence suggests that leptin and adiponectin, key hormones in regulating appetite and metabolism, may play important roles in regulating cancer cell growth and proliferation (27). In fact, weight gain is a common side effect in women undergoing adjuvant chemotherapy for breast cancer (28). Previous findings suggest that weight gain during adjuvant chemotherapy is associated with increased recurrence and poorer survival (29). Thus, understanding the molecular mechanism of this side effect is of therapeutic importance. The global similarity between obesity-specific expression and anticancer drug-treatment expression may lead to the discovery of shared mechanisms at the transcription level, allowing optimized treatments to be designed.