A searchable cross-platform gene expression database reveals connections between drug treatments and disease

Drug treatment based profile SPIED queries

The CMAP contains expression change profiles as ranked array probe IDs for 6,100 individual treatments corresponding to 1,309 distinct drug-like compounds. Statistically filtered response profiles can be defined for 1,218 of the drugs as these have at least three instances in the database. The profiles can be mapped onto a non-redundant gene list by uniquely associating one probe ID to a given gene and dropping the other probe ID for this gene with less robust expression changes over the database. This is the same methodology underlying the SPIED database. We took the responder profiles for the 1,218 drugs and searched the SPIED for maximally (anti-)correlated expression change profiles. The objective is to see to what extent the CMAP transcriptional signatures correlate with transcriptional responses assimilated within our platform independent database of over 100,000 microarrays deposited by a very large number of groups to the public domain.

The CMAP is well populated with drugs that target the same or different steps in the PI3K-mTOR signalling cascade. In this context the results for LY-294002 (PI3K inhibitor) (61 instances), rapamycin (mTOR antagonist) (44 instances) and wortmannin (PI3K/mTOR antagonist) (18 instances) showed a high degree of overlap (rapamycin v LY-294002 r = 0.90 N = 3565, rapamycin v wotmannin r = 0.91 N = 1849, LY-294002 v wotmannin r = 0.88 N = 2217, where r is the Pearson correlation coefficient and N is the number of shared genes with significant fold values), see additional file 1 for the full fold change data. It is a straightforward matter to query the SPIED with these drug expression profiles. This is done by calculating the regression scores against the individual SPIED entries and retaining the top ~100 correlations, see Methods for details. For simplicity and uniformity of treatment, unless otherwise stated, we query SPIED with expression profiles containing 500 genes with the largest fold values passing the p < 0.05 significance threshold. It should be noted that results will be largely insensitive to the size of the query profile. The top SPIED correlate for all three drugs was the Pan-PI3K inhibitor GDC-0941 treated T47D breast cancer cells and the regression scores for the tree query signatures against all 6 samples in the series are shown in Figure 1A. The high degree of correlation is illustrated by regression plots for the three query profiles against the pooled GDC-0941 profile, see Figure 1B, C, D. All three inhibitor queries also pick out mTOR antagonist studies [16], but a more interesting correlation is with a glucocorticoid (dexamethasone) treatment of acute lymphoblastic leukaemia (ALL) cells [17], the rapamycin scores are shown in Figure 2A. The correlation increases with the length of drug treatment, being higher at 24 hours, Figure 2B, C. This result reveals another connection between mTOR antagonism and the corticosteroid mechanism as it has been shown that corticosteroid resistance in ALL can be overcome by mTOR antagonism [18]. Chronic myeloid Leukaemia (CML) and some instances of ALL are the result of the ABL tyrosine kinase translocation and fusion to BCR, the BCR-ABL fusion event [19]. This pathology has been targeted with rapamycin and our results support this approach based on the high degree of anti-correlation of the CMAP rapamycin profile with a transcriptional profile of BCR fusion construct transformed chord blood cells. The correlation scores are shown in Figure 3A. There is a clear anti-correlation of rapamycin profile with the BCR-ABL profiles pointing to a possible reversal of the phenotype, Figure 3B. Also, there is a high anti-correlation with the BCR-FGFR1 profile indicating a possible therapeutic role of rapamycin, Figure 3C.

Next we consider profiles derived from disease states. For brevity we focus on two unrelated pathologies: cancer and neurodegeneration.

Querying SPIED with cancer derived profiles

The class of diseases with the most extensive repository of expression data is cancer and therefore a cancer disease profile search of SPIED will be an ideal testing ground for the methodology. The original CMAP disease application implicated mTOR inhibition as a target for imparting sensitivity to dexamethasone treatment resistant ALL [18]. We searched the SPIED database with the dexamethasone resistant v sensitive profile to see if there are common features in published transcriptional studies. The query profile consisted of the 500 most highly regulated genes that passed the lowest significance test of p < 0.05, see additional file 1. As with the SPIED profiles the query profile also consists of a non-redundant gene list. Not surprisingly, the highest correlation scores came from the experiments from which the query profile was generated, see additional file 2 file. In addition, we found a high correlation to an independent later study of ALL sensitivity to corticosteroid (prednisolone) treatment [22]. This study generated transcriptional profiles of ALL patient leukaemia cells with the objective of uncovering a gene signature that can predict the sensitivity to prednisolone treatment. Combining the 27 infant and non-infant corticosteroid sensitive samples and the 25 resistant samples we can define a statistically filtered sensitivity profile to make a direct comparison with the query profile and we find a high degree of correlation (r = 0.94), see additional file 2. When the high scoring sample belongs to a relatively large sample series and the phenotype is binary we can perform a non-parametric significance test to measure the extent of enrichment of the given phenotype for high or low correlation scores. For example in the last case there were 25 resistant and 27 sensitive samples. Ranking the samples according to their correlation with the resistant versus sensitive query profile we find 20 resistant samples in the top 25 and 22 sensitive samples in the bottom 27. This is highly significant and can be quantified with a simple Fisher exact test. Explicitly, the probability p of 20 or more resistant samples in the top 25 correlations is less than 9 × 10^-7. The K-S significance score can be calculated by counting the number of times a random rearrangement of the samples gives a better enrichment, we find p < 3 × 10^-6. The enrichment plot is given in Figure 4A. As expected the top scoring correlations were dominated by samples from blood derived cells, for simplicity we restricted our analysis to the top 100 most significantly correlating samples. However, two studies in unrelated tissue pathologies were highly correlated with the corticosteroid resistant profile. These were a comparison of lung epithelia with cancer in smokers [23] and a differential expression between healthy and cancerous pancreatic tissue [24]. The smoking study consisted of non-diseased lung epithelia from 187 individual smokers 97 of whom were diagnosed with lung cancer. Ranking the samples according to query correlation score we find that in the top 97 there are 64 cancer cases and in the bottom 90 there are 57 non-cancer cases, with a significance score of p < 5 × 10^-5. The K-S significance is p < 2 × 10^-4. The enrichment for positive correlations with the corticosteroid resistance profile in the cancer cases is shown in Figure 4B. Interestingly, it has been shown that there is a down regulation of the glucocorticoid receptor (GR) in small cell lung cancer [25, 26] and reversing this promotes cancer cell apoptosis [27]. The pancreatic cancer study sought to establish a transcriptional signature of tumour versus normal pancreatic tissue by laser capture of cancerous and normal tissue from the same pancreas. In total 39 sample pairs were published and we find a high positive correlation with the corticosteroid resistance profile, p < 2 × 10^-6 and a K-S significance score of p < 3 × 10^-7. The enrichment curve is shown in Figure 4C. In this context it has been reported that loss of GR expression has been seen in pancreatic carcinoma relative to normal tissue [28] and elevating GR expression has been shown to inhibit pancreatic tumour growth in a hamster model [29]. The query results are given in the additional file 2.

Neurodegenerative Disease

The analysis of gene expression changes associated with neurodegenerative disease has been hampered by the difficulty of extracting high quality RNA from post-mortem tissue [30, 31]. One way of validating a disease associated gene expression profile is to show that it shares significant features with profiles derived from independent experiments on related pathologies. A positive result would validate the query profile and furthermore lead to a more robust core response profile based on multiple experiments. To this end we constructed three separate query profiles based on transcriptional profiles from the brains of patients with three degrees of severity of Alzheimer's disease (AD) [32], see additional file 1. The number of significant changes increases with severity of disease and we queried the SPIED with these three profiles, see additional file 2. Not surprisingly, the high scoring correlations are those from which the query profiles were derived. In addition to these the query returned correlations with other AD studies and various neurodegenerative diseases. The high scoring AD expression series was an extensive study of 161 samples from various brain regions of AD patients and age-matched controls [33]. Ignoring brain regions for now, there are 87 AD samples and 74 controls. Ranking the samples according to correlation score against the severe AD query profile we find a very significant enrichment of positive correlations with AD samples (p < 10^-8, based on the Fisher test as above). Pooling the samples from the different brain regions results in significant correlations for 5 out of the 6 brain regions, see Figure 5.

The PD correlation was with a study of 94 samples from three different regions of diseased and normal brains [35]. Pooling samples according to brain region we find that the severe AD profile had a high correlation with all three regions studied: superior frontal gyrus (SFG) r = 0.88; lateral substatia nigra (LSN) r = 0.77; medial substantia nigra (MSN) r = 0.82, see Table 3.

The chromosome 21 trisomy underlying DS leads to the development of many of the characteristics of AD pathology [36, 37]. Therefore, it is not surprising to find a high correlation in SPIED form a transcriptional profiling of DS brains. This study comprised 8 healthy and 7 DS individual brains [38]. Combining the expression data into a thresholded fold change profile we find that there is a significant but small positive correlation with the severe AD profile, with r = 0.58. Interestingly, the correlation is higher with the moderate AD profile, with r = 0.68, see Table 3.

The first transcriptional profiling of BD brains pointed to the down regulation of synaptic and mitochondrial proteins in the orbital frontal cortex [39]. This synaptic pathology picture of BD is further strengthened by our analysis of the AD profile correlates within SPIED. Pooling the 10 BD samples and 11 controls we find a high regression score with the severe AD profile, see Table 3[39]. It is important to note that this correlation is with a subset of the BD signature as it consists of genes that are also altered in AD. However, it is outside the scope of the present paper to combine profiles into disease specific queries.

Compiling the data

Microarray sample files, GSM files, were downloaded form the NCBI GEO database. Individual GSM files were assigned to GSE series and log scaled values scaled to linear and low level responders dropped. EF profiles were then generated based on ratio of individual condition to the average across the series. Expression data from five Affymetrix GeneChip platforms corresponding to three species were collected. These were: all samples from two Human array platforms corresponding to Human Genome U133 Array Set HG-U133A GPL96 (27,337 samples) and U133 Plus 2.0 Array GPL570 (55,196 samples); all samples from the Mouse Genome 430 2.0 Array GPL1261 (21,219 samples) chip; all samples from two Rat chips corresponding to Rat Genome U34 Array GPL85 (2,827 samples) and Rat Genome 230 2.0 Array GPL1355 (7,476 samples). The database thus totals 106,101 samples. Of course, this can always be extended to include more platforms from the same species and/or other species.

Non-redundant database

The individual GSM sample file expression values were transformed into EF values corresponding to the expression relative to the series mean. Expression values that have been logarithmically transformed are transformed back to a linear scale and low expression values dropped, that is are set to zero and don't contribute to the fold profile. We found that the results were relatively insensitive to the cut-off value and we set this to be 10% of the average expression value. All sample expression profiles within a series were scaled to the same average. The fold vales are defined as$f_k=2\frac{s_k-\bar{s}}{s_k+\bar{s}}$[13], where s _k is the expression level of the k ^th probe set and $\bar{s}=\frac{1}{n}\sum _{k=1}^n{s}_k$is the average over the series. For the database to be searchable with cross platform response profiles and gene lists it has to be rewritten as a database of expression profiles over non-redundant gene lists. The EF profiles across the probe sets were therefore mapped onto expression profiles for a non-redundant gene list. In general each gene is represented by multiple probe sets. For each platform we generated the EF statistics for each probe set across the totality of samples. The probe set with the most robust response across the samples was chosen to represent the gene. Explicitly, the probe set with the highest root mean square deviation form zero was chosen to represent the given gene. The number of genes defined on each platform were as follows: GPL96 11,807, GPL570 15,983 genes, GPL1261 13,202 genes, GPL85 chip with 3,844 genes, GPL1355 chip with 6,341 genes. The database totals 106,101 samples and is searchable on a reasonably fast desktop PC in ~10 minutes per query.

Searching the database

The query profile is a statistically thresholded non-redundant list of genes and associated fold values. Statistical significance is assigned to a fold change based on a simple Student's t-test between multiple control and treatment sample expression values. This is compared to each profile in the database by means of a simple Pearson regression analysis, with a correlation coefficient r. The experiments are ranked according to the significance. The significance is measured by scaling the correlation to the normal by a Fisher transformation and measuring the number of standard deviations from the mean. The Fisher transformation is $r^{\prime }=\frac{1}{2}ln\left(\frac{1+r}{1-r}\right)$ and the standard deviation is$\frac{1}{\sqrt{N-3}}$, where r is the Pearson correlation coefficient and N is the number of genes making up the correlation. The final ranking score is $s=\frac{\sqrt{N-3}}{2}\left|ln\left(\frac{1+r}{1-r}\right)\right|$.

CMAP combined profiles

The CMAP contains ranked lists of probes for 6,100 separate perturbagen treatments of four different human cell lines, with the ranking based on response level relative to control. The treatments are various multiples of 1,306 different drug-like compounds. To generate responder sets that can be used to search SPIED we combined rankings for each separate compound treatment and converted these into pseudo-fold values with associated statistics. The pseudo-fold value is defined by$f_i=1-2\frac{r_i-min}{max-min}$, where r _i is the rank of the i ^th gene and min/max are the minimal/maximal ranks. Remembering that the highest rank corresponds to the most up-regulated gene. The SPIED was searched with CMAP profiles corresponding to folds with a p < 0.05 threshold and with at least three replicates. This left 1,218 separate perturbagen probes. We sought to cluster the perturbagens based on predicted target and response profile similarity. The profiles are given in the additional file 1 file.

Availability of SPIED

The SPIED database and associated executables are available for download from ftp://ftp.hostedftp.com/~GWftpFILES/SPIED/. The download consists of the SPIED database together with executables for searching SPIED. Source code files to generate the database and perform query searches are provided together with the executables. Documentation on the database, the executables and source code files is also included.