GeneFriends: An online co-expression analysis tool to identify novel gene targets for aging and complex diseases

GeneFriends: an online tool for the research community

The aim of the project was to create a user-friendly tool, which can take a list of genes related to a given disease or process and quickly identify new candidate genes. Using co-expression profiling the genes are given in rank order to help prioritize candidates for experimental follow up. Underlying GeneFriends is a Mus musculus co-expression map created from 1,678 microarray datasets, comprising over 20,000 individual samples from previously published experiments. To create the co-expression map we employed a vote counting method. The co-expression map contains ≈ 427.5 million gene pairs (20,676 x 20,676) arranged in a matrix and given a score based on how often they are co-expressed across all microarrays (see Materials and Methods).

We feel this output will help researchers in various fields identify interesting genes for follow up studies in a quick and intuitive manner. To test if this novel tool can be used to derive biologically-relevant predictions we tested gene sets related to aging, cancer or mitochondrial complex I disease. Furthermore, we tested two predicted novel candidates experimentally, as detailed below.

Testing the co-expression map

The biological significance of the co-expression map was verified using nine genes known to be involved in three biological processes: cell division cycle, immune response and fatty acid metabolism. The top 5% of genes identified by GeneFriends as co-expressed were analysed by DAVID [17] to detect enriched biological processes and functions. The results confirmed that genes co-expressed were functionally related to either cell cycle (Enrichment score: 56.5, FDR: 3.0x10^-75), immune response (Enrichment score: 14.7, FDR: 7.8x10^-13) or fatty acid metabolism (Enrichment score: 20.98, FDR: 1.2x10^-23); detailed results are included in the supplementary data (Additional file 1). Furthermore, Figure 1 shows the clustering of the co-expression map’s network, demonstrating that co-expressed genes tend to be involved in the same biological processes.

GeneFriends uses a vote-counting method to rank co-expression. We compared GeneFriends to COXPRESSdb [9], which utilizes the more commonly used correlation value (Pearson or Fisher). To do so we selected 3 genes with known functions and retrieved output from both tools and used DAVID to determine enriched categories. The results show similar categories and scores although the overlap in the specific genes can vary (Additional file 2). When comparing the numbers of transcription factors present in the top 1000 co-expressed genes from GeneFriends and COXPRESSdb the results are similar. This demonstrates that our approach is similar to using Pearson correlation to create a co-expression map.

Candidate gene prediction from process/disease gene lists

We used GeneFriends to identify novel candidate genes associated with specific processes or diseases. The results show the numbers of times each of the 20,676 genes in the co-expression map were "friends" with genes in the disease gene seed and corresponding p-values indicate the statistical significance of the co-expression (see Materials and Methods). The p-value is calculated based on the number of seed genes a given gene is co-expressed with and the total number of genes it is co-expressed with (Materials and Methods). DAVID was used to interpret the broader biological significance of the results. All genes with a co-expression p-value <10^-6 were classified using the default settings in DAVID. This is a stringent cut-off since using a Bonferroni correction for multiple testing results in: 0.05/20,677 > 10_.^-6

Aging-related gene prediction and putative transcriptional mechanisms

The most significantly over-represented functional clusters were inflammatory response (enrichment score (ES) = 24.13, FDR = 1.97x10^-18), vasculature development (ES = 10.18, FDR = 2.31x10^-8) and lysosome (ES = 9.00, FDR = 2.25x10^-8). Since most of the genes in the seed list were classified in the categories related to the immune system, it was unsurprising to find similar results for the co-expressed genes.

Since these TFs are co-expressed with the aging-related genes it was expected that these genes, at least in part, would be regulated by the co-expressed TFs. Therefore, they would share transcription factor binding sites (TFBS) for these TFs. To identify over-represented binding motifs in the genes co-expressed with the aging genes (p-value < 10^-6) we employed factorY [23]. For the aging gene set this revealed Nfkb as the most significant result (Additional file 4). Some of the TFBS identified have co-expressed genes with the aging seed list such as NFKB1 (p_TFBS =1.48 x10^-5, p_Coexpress = 4.44 x10^-9), the C/ebp (p_TFBS = 6.95 x 10^-3, p_Coexpress =7.9 x 10^-34) genes and IRF1 (p_TFBS = 5.8 x 10^-4, p_Coexpress =8.17x10^-23) genes (Additional file 4). However, TFBS for Isre, Nfkb2 (p65) and Sp1 were identified as over-represented but not co-expressed and many co-expressed TFs did not have over-represented binding sites.

Cancer-related gene prediction

From the original seed list only 6 genes pass the p <10^-6 threshold co-expressed with the entire set of cancer genes, which could be due to the heterogeneity of cancer etiology. However, there were several significantly co-expressed genes, not included in the seed list, that have previously been associated with cancer. For example, Cdc25a, Cdc25b and Cdc25c, members of the Cdc25 family, are significantly co-expressed (p < 10^-6) with the cancer-related genes.

There were a high number of significantly co-expressed centromeric proteins co-expressed with the cancer seed list. These proteins play a role in chromosome segregation and incorrect segregation of chromosomes during the cell cycle can lead to cancer [24]. Cep152 is involved in centriole duplication [25]. Cenpp, as well as CenpN/F/H/J/I/C1/T/K/M/E/Q/A/L, are all co-expressed significantly with the cancer seed list (Additional file 5) and are part of the Cenp-A NAC complex. This complex is required for proper kinetochore function and mitotic progression and its disruption can lead to incorrect chromosome alignment and segregation that preclude cell survival despite continued centromere-derived mitotic checkpoint signalling [26, 27]. Plk1, Aurka, Aurkb and Cdca8 are in the top 50 co-expressed genes, these play an important role in cancer formation [28, 29].

Several un-annotated genes (Bc055324, E130303B06Rik, 4930547N16Rik, F730047E07Rik, 1110034A24Rik, and 4632434I11Rik) were co-expressed with the cancer-related genes, suggesting these genes might play a role in occurrence or pathophysiology of cancer. One of these un-annotated genes, Bc055324, is a predicted protein coding gene, which has a high co-expression ratio of more than 0.7 with the cancer genes Rad51 and Ccdc6 [30], indicating this gene is increased in expression in >70% of the cases when Rad51 is increased in expression. Many other cancer-related genes such as Brca1 and Brca2 also show a strong co-expression with the Bc055324 gene (Additional file 6). All genes co-expressed with Bc055324 show enrichment for the cell cycle (ES = 52, FDR = 1.7x10^-74). A BLAST analysis of the protein sequence shows no significant homology to other Mus musculus proteins. Similar sequences, however, are found in a large number of different multi-cellular species such as Gallus gallus, Bos taurus and Homo sapiens and there also is a significantly similar gene present in Arabidopsis thaliana, suggesting it is conserved in plants as well, indicating it plays a functional role.

Validating the role of C1ORF112 and C12ORF48 in growth of cancer cells

To test our predictions, we employed siRNAs to knock down the human homologs of Bc055324 (C1ORF112) and 4930547N16Rik (C12ORF48) in the widely-used HeLa cancer cell line. These two genes were selected for validation because they are co-expressed with genes involved in the cell cycle (Additional file 6), so knockdown should lead to a measurable phenotype. Furthermore, validated siRNAs are available (see Materials and Methods) for these genes. The results show that the growth rate of the cancer cells is significantly lower when either C1ORF112 or C12ORF48 are knocked down (Figure 2). These results support our predictions and demonstrate that C1ORF112 and C12ORF48 are functional.

Mitochondrial complex I disease-related gene prediction

The co-expressed genes include several mitochondrial complex I genes (not in seed list), multiple cytochrome c proteins and genes involved in the ATP synthase complex. Furthermore, there are approximately 50 poorly/un-annotated genes co-expressed. A pseudo gene, 3000002C10Rik, shows a co-expression ratio of >0.5 with 512 genes. Classification of these 512 genes using DAVID results in an enrichment score of 53.7 (FDR = 2.9x10^-70) for mitochondrial genes. Therefore, 3000002C10Rik may play a biologically relevant role in mitochondrial processes.

Predicting functions of un-annotated genes

GeneFriends: a genetics and genomics tool for the research community

Genefriends is freely available online (http://genefriends.org/) and is an intuitive tool, which can be used to identify the genes co-expressed from a user supplied gene list. This simple, yet powerful new tool can be a valuable resource for genome interpretation, annotation, mouse genetics, functional genomics and transcriptional regulation. It may also be useful to develop network analyses of mouse genes in a variety of studies.

We tested GeneFriends to determine whether it can give biologically relevant data. We also demonstrated how GeneFriends can be used to quickly identify interesting gene targets for follow up studies. Furthermore, we experimentally validated two un-annotated genes co-expressed with a cancer seed list. Below we discuss the findings we have obtained from our example analyses and the biological relevance of our results.

Validation of the co-expression map

The main assumption made when constructing our mouse co-expression map is that co-expressed genes tend to be involved in the same biological processes. Our results clearly support this, showing that genes known to play a role in certain biological processes such as immune system, cell cycle and lipid metabolism are co-expressed and form clusters in the co-expression network (Figure 1). Different clusters show different strengths of co-expression, therefore the effectiveness of the gene set analyses may vary. However, there is a high degree of functional coherence between co-expressed genes showing that our co-expression map can be used to obtain biologically-relevant information.

Given the intrinsic noisy nature of microarray data, we used a vote counting approach, which is a standard meta-analysis technique, to build our co-expression map [31]. The aim of this approach is to emphasize sensitivity. This is particularly important when combining large and diverse datasets and has been shown to increase sensitivity when studying aging, for example [18]. From a technical point of view, we aimed to include as much data as possible and because of the technical implications accompanied by this the vote counting method is more efficient and has faster run times than other methods.

Co-expression analysis of genes associated with aging

Using a guilt-by-association method we identified candidate genes related to seed lists of genes associated with diseases or processes. Our study not only identified genes relevant to current theories of aging, e.g. inflammation, but it also identified novel candidates for further research. C/ebp transcription factors showed the strongest co-expression and are therefore candidate activators of the altered expression patterns with age. TFBS for C/ebp genes were identified in the aging genes and there is some evidence of a transcriptional cascade via SP1 [32]. The two proteins encoded by the C/ebpβ gene are liver activating protein (LAP) and liver inhibiting protein(LIP), which have opposing effects [33, 34]. The LIP protein is also capable of inhibiting other C/ebp proteins. This could explain why C/ebp transcription factors themselves are not found to be increased/decreased in expression with age. This could also be due to the fact that TFs are sometimes expressed at low levels not detected by microarrays.

Replacement of the C/ebpα gene with C/ebpβ increases lifespan by 20% [35, 36], and may alter the rate of aging [37], indicating that altering the isoform expression of these genes can affect lifespan. Moreover, the life-extending drug rapamycin may affect isoform ratios of C/ebpβ. Rapamycin has been shown to increase lifespan via the suppression of Mtor [38] which in turn controls the isoform ratios of C/ebpβ [39]. Therefore, we speculate that rapamycin may in part exert its life extending effect through C/ebpβ.

Co-expression analysis of cancer genes and experimental validation of candidates

We used GeneFriends to identify new candidate genes for a role in cancer. Many of the cancer genes in the initial seed list were not present in the results, indicating they are not co-expressed with each other. This may be due to the fact that this set of cancer genes includes both oncogenes and tumour suppressor genes which are not expected to be co-expressed. Also cancer can arise through different mechanisms. Therefore the genes identified as co-expressed in this study are likely involved in common pathways leading to cancer, or are at least triggered by transformation.

Genes that are co-expressed with several oncogenes may prove to be useful targets in countering the proliferating effect of these genes in tumours. Examples of such genes that are already being studied are Cdc7 and Cdc25 genes, both of which were identified as co-expressed in our study. Cdc25 has been suggested as a therapeutic cancer target and on-going studies in this direction have shown some level of success [40–43]. Two compounds that target Cdc7 are currently in phase I clinical trials [44]. The fact that candidate genes identified by our method have already been suggested as potential drug targets shows that GeneFriends can be useful for the identification of candidate targets for cancer studies.

Bc055324 is one of the poorly annotated genes that is strongly co-expressed with a large number of cancer genes. Knock-down of the human homolog, C1ORF112, in HeLa cells diminishes cell growth, which adding the fact that Bc055324 knockout mice are not viable [45] (http://www.europhenome.org/), demonstrates that this gene is functional. Further studies of this gene in the context of cell cycle regulation, development and cancer are warranted. These results show that GeneFriends can indeed be used to identify novel targets for particular diseases. In addition it confirms that the functional enrichment of co-expressed genes can give indications about an un-annotated gene's function. The other un-annotated gene co-expressed with cancer we tested was 4930547NRik (C12ORF48). C12ORF48 was recently shown to be over-expressed in pancreatic ductal adenocarcinoma cells [46] and in other aggressive and therapy-resistant malignancies [46]. In line with our findings in HeLa cells, knock down of C12ORF48 significantly suppressed PDAC cell growth [46].

Co-expression analysis of mitochondrial I complex disease genes

Mitochondrial complex I diseases include isolated complex I deficiency, which is the most common enzymatic defect of the oxidative phosphorylation disorders and can cause a wide range of clinical disorders [47, 48]. These include macrocephaly with progressive leukodystrophy, nonspecific encephalopathy, cardiomyopathy, myopathy, liver disease, Leigh syndrome, Leber hereditary optic neuropathy, and some forms of Parkinson’s disease [49–51]. Mutations in the nuclear encoded mitochondrial genes have been associated with several pathologies [52, 53]. However, half of the patients with mitochondrial complex I (CI) deficiencies lack mutations in any known CI subunit suggesting that yet unidentified genes crucial for maturation, assembly, or stability of CI may be involved in these diseases [53]. We identified several un-annotated genes that show a strong co-expression with the mitochondrial disease gene set. As most of the other co-expressed genes encode mitochondrial proteins, these un-annotated genes most likely also encode mitochondrial proteins. This is further supported by the fact that a number of these un-annotated genes have been shown to be active in the mitochondria in another large-scale study [54]. Some of these genes could be responsible for the CI deficiency phenotype and are therefore promising candidates for further studies.

Data selection

To create the co-expression map, normalized microarray data obtained from the GEO database was used [55]. GEO files GSE1 to GSE18120 were downloaded containing 16,916 datasets in total. From these, 3,850 Mus musculus datasets containing 64,849 microarrays and the corresponding annotation files were extracted. As mouse experiments are generally better controlled than human studies and there is less variation caused by genotypic and environmental factors in the mouse, Mus musculus data was chosen over Homo sapiens to reduce noise. Using mouse data also allows more datasets, coming from a more diverse set of experiments [18]. This also potentially allows for the investigation of target genes in the different mouse models of aging and complex diseases.

All datasets containing annotation files that did not include gene symbols for at least 90% of the probes present in the data were removed. All microarray datasets containing values higher than 25 were log transformed, under the assumption this data was non-log-transformed data. To remove poor signal, low quality or nonsense values up to 10⁹⁹ datasets containing no values above 2log(5,000) or one or more values over 2log(20,000,000) were removed. Datasets with no reference to any annotation file were removed. Even though it is not feasible for us to perform a comprehensive evaluation of the quality of the data in each experiment, a meta-analysis is in its essence a technique to eliminate poor quality data and hence we are confident that there are no systematic errors in our analysis that artificially originate false results. After these steps 1,678 datasets containing 8,417 different conditions and 21,744 individual samples remained. The probe IDs were converted into gene symbols. If multiple probe IDs mapped to the same gene symbol, they were averaged. Within each dataset the experimental conditions were manually determined. Microarrays from individuals under the same conditions were averaged; in other words, duplicates were averaged. Missing values were ignored as long as there were duplicates; if there were no duplicates this gene symbol was removed.

Constructing the co-expression map

To create GeneFriends we first constructed a genome-wide co-expression map, using normalized Mus musculus microarray data from the GEO database. This describes which genes are related based on how often they are co-expressed. In total, 1,678 mouse datasets containing 8,417 different conditions and 21,744 individual samples met our data selection criteria. To construct our expression map the different conditions within each dataset were compared to each other. Since different datasets contain different probes mapping to different gene symbols a selection was made. Only those gene symbols that are present in gene platform file GPL1261 (Affymetrix GeneChip Mouse Genome 430 2.0 Array) were used. This platform contains 20,676 gene symbols and is the most common platform used for microarrays amongst those included in this work. All of these gene symbols were present in over 850 datasets.

Where x is the total number of comparisons and N the number of times two genes are differentially expressed (in the same direction) simultaneously.

For simplification the terms “UP” and “DOWN” are used in this formula. The actual directionality of the change is irrelevant as it is dependent on the direction of the comparison (i.e. whether one compares group 1 with group 2 or the other way around). Either way the results will be the same.

Where x is the total number of comparisons and i describes the current comparison between the different conditions.

From the values N and Q the co-expression ratio was deducted. The genes were then ranked based on their N:Q ratio. A ratio of 0.50 would indicate that if gene 1 is increased or decreased in expression in 50% of the cases gene 2 is also increased or decreased in expression. Each gene pair is present in at least 850 datasets and so the ratio is based on a large number of measurements.

Testing the co-expression map

To investigate the capacity of the co-expression map to provide biologically-significant results nine genes known to play a role in specific biological processes were investigated: three genes that are known to be active in fatty acid metabolism: Ppara, Acaa2 and Acadm; three genes known to be involved in immune response: Cd4, Cd8 and Il10; and three cell cycle genes: Cdc6, Cdc7 and Cdc8. These were selected before the analyses to reduce biases.

It was expected that genes co-expressed with these nine genes would be involved in the same biological processes; to test this assumption we inspected if any specific categories were over-represented by these groups of genes (see below for functional enrichment analysis).

Prediction of novel candidate genes in aging and complex diseases

In order to identify genes co-expressed with known disease genes, three disease related gene sets were included. The first of these was an aging gene set. It consisted of genes over-expressed with age obtained from a meta-analysis of aging microarray studies in mice, rats and humans that revealed several conserved genes increasing or decreasing in expression with age [18] (Additional file 9). The second gene set included was a set of cancer-related genes [56] (Additional file 9). This is a manually curated cancer set that includes only heritable cancer genes with strong evidence that mutations in these genes are causative for cancer. The third gene set added included genes known to cause diseases through mitochondrial complex I deficiencies. The genes in this set contain the nuclear mitochondrial complex I deficiency genes in the OMIM database (Additional file 9). Gene symbols that were not present in the co-expression map were not included in the analysis.

Where p is in effect the chance that a particular gene occurs in the top 5% of a random gene.

Where Pr(K > = k) is the probability that a gene would be “friends” with k or more genes in the disease gene set; k is the number disease gene “friends”; n is the number of genes in the gene set. When calculating p the number of occurrences of a gene in the top 5% of all genes was included. This is necessary since some genes tend to be co-expressed more often in general than other genes.

Experimental validation of cancer-related genes Bc055324 and 4930547N16Rik

To test the predictions from the analyses using GeneFriends, we took un-annotated genes that were the most co-expressed with the cancer disease gene list. Validated siRNAs were available from Qiagen for two the human homologs of the top un-annotated genes: Bc055324 (C1ORF112) and 4930547N16Rik (C12ORF48). The experiment was conducted in human HeLa cells using standard culture conditions. A negative and a positive control were also included (Qiagen). The positive control contained a mix of several apoptosis inducing siRNAs, demonstrating that the transfection was successful through elevated cell death. The negative control consisted of siRNAs targeting non-mammalian genes. The full protocol followed for this experiment is described in Additional file 10.

Gene function enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) [17] was used to identify enriched functional groups within these genes. The default settings were used in this analysis. The results were ranked based on p-value and genes with a p-value <10^-6 were selected. This is a stringent cut-off since using a Bonferroni correction for multiple testing results in: 0.05/20,677 > 10^-6_. In addition, a randomization test showed this is a very stringent cut-off in which no false positives are expected by chance. The test entailed the construction of several random sets of genes which were then used to find co-expressed genes. The p-values found for all results were >10^-5, indicating that genes with a p-value <10^-5 are very unlikely false positives.

Next, to understand the significance of the DAVID enrichment score, 1000 genes were randomly selected and used as an input for DAVID. This resulted in an enrichment score of 2.2 with an FDR score of 0.7 for the most significant category found. The same was done for smaller numbers resulting in similar scores. This indicates that the Enrichment scores of >10 and FDR <10^-10 are very significant and cannot be found by random chance.

Benchmarking using our co-expression map and COXPRESSdb revealed similar results (Additional file 2), suggesting our co-expression map is not inferior to those built using correlation measures. Therefore, our work demonstrates that vote counting is a viable method to build co-expression maps.

One concern is that some of the genes in functional categories were assigned based on their expression pattern; if this would be the case it could lead to circular reasoning. However, an analysis conducted by the DAVID team shows that less than 1% (116/20676) of the genes is grouped based solely on their expression pattern (see Additional file 11 for a list of these genes).

BLAST

A BLAST search was conducted on the Bc055324 gene. The protein sequence for this gene was recovered from GenBank version: Bc055324.1. This sequence was PSI-blasted against the non-redundant protein sequence database. Then the search was iterated twice including all sequences recovered in the initial search.