Genome-wide modeling of complex phenotypes in Caenorhabditis elegans and Drosophila melanogaster

Abstract

Background

The genetic and molecular basis for many intermediate and end stage phenotypes in model systems such as C. elegans and D. melanogaster has long been known to involve pleiotropic effects and complex multigenic interactions. Gene sets are groups of genes that contribute to multiple biological or molecular phenomena. They have been used in the analysis of large molecular datasets such as microarray data, Next Generation sequencing, and other genomic datasets to reveal pleiotropic and multigenic contributions to phenotypic outcomes. Many model systems lack species specific organized phenotype based gene sets to enable high throughput analysis of large molecular datasets.

Results and discussion

Here, we describe two novel collections of gene sets in C. elegans and D. melanogaster that are based exclusively on genetically determined phenotypes and use a controlled phenotypic ontology. We use these collections to build genome-wide models of thousands of defined phenotypes in both model species. In addition, we demonstrate the utility of these gene sets in systems analysis and in analysis of gene expression-based molecular datasets and show how they are useful in analysis of genomic datasets connecting multigenic gene inputs to complex phenotypes.

Conclusions

Phenotypic based gene sets in both C. elegans and D. melanogaster are developed, characterized, and shown to be useful in the analysis of large scale species-specific genomic datasets. These phenotypic gene set collections will contribute to the understanding of complex phenotypic outcomes in these model systems.

Background

Traditional experimentation in animal model systems such as the worm Caenorhabditis elegans and the fly Drosophila melanogaster often results in complex molecular and phenotypic outcomes. Frequently a targeted deletion or ectopic expression of a single gene product results in pleiotropic phenotypes. Similarly, broad high-throughput multiplex experimental strategies such as microarray based gene expression, RNA interference (RNAi) screens, or next-generation DNA and RNA sequencing, analyzing phenomena such as development, behavior, mating, diet, and life span, typically produce large datasets requiring complex analytical approaches.

Gene sets are collections of keyword terms with annotated genes derived from multiple sources of a priori information. They have been used in computational analysis of gene expression data [1–3] with the goal of identifying higher order relationships beyond simple gene list results, as well as in analysis of population based GWAS in humans [4, 5]. The most commonly used gene sets include those derived from GO annotations [6], biological pathways from KEGG [7] or BioCarta, expression modules, DNA binding sites, or other sources of molecular information [1, 3, 8]. Each collection of gene sets has its own unique qualities and features which are useful in different ways. For instance, KEGG emphasizes metabolic and biochemical pathways; GO annotations, while having some phenotypic content, emphasizes molecular function, cellular component, and biological processes, while MSigDB [8] emphasizes gene expression signatures. This information is often closely related, or “proximal” to gene and molecular function, rather than more “distal” information regarding phenotypic outcomes and disease susceptibility. Recently, phenotype based gene sets have been derived exclusively from genetically determined phenotypic associations for mouse phenotypes and common human disease [9, 10], resulting in gene sets for specific phenotypes, organized by a structured systematic ontology.

Here, we present gene sets for worm and fly, which use the structured ontology found in the Worm Phenotype Ontology from the C. elegans database - WormBase [11] and phenotypic descriptions for D. melanogaster found in FlyBase [12]. These gene sets are derived from information on gene-phenotype relationships based on genetically determined phenotypes. We use these collections in large scale phenotypic modeling in worms and flies and demonstrate their utility in complex analysis in multiple ways, including analysis of gene expression datasets representing complex phenotypic and biological phenomena in both C. elegans and D. melanogaster. In this way, we integrate large scale genome analysis with large scale phenotypic analysis in these two model systems.

Results

Derivation of worm gene sets

The worm gene sets presented here are derived from two lists of genes and assigned phenotypes provided by Gary Schindelman and Paul Sternberg as a component of the Worm Phenotype Ontology [13]. These two lists originated from information curated from RNAi experiments and genetic variations (VAR) as archived in WormBase [14].

Two worm gene set files (CE- RNAi-GS and CE-VAR-GS) were produced by parsing each gene list separately into non-redundant lists of unique phenotypic terms with all genes assigned to their corresponding phenotypic terms. This produced two non-redundant gene set files containing 850 and 1109 gene sets for RNAi and VAR, respectively. In addition, we developed a master worm file by combining the original RNAi and VAR gene lists into a combined file (CE-Combined-GS) containing 1,385 non-redundant phenotypes and their associated gene sets.

Derivation of fly gene sets

The Drosophila gene sets described here are derived from phenotypic data provided in FlyBase (see Methods). A file containing 259,162 phenotypic descriptions with assigned Drosophila genes was collapsed and parsed resulting in a non-redundant gene set file of 11,999 unique phenotypic terms with annotated genes. This file named DM-narrow-GS was used for systems biology and gene expression analysis.

Table 1 shows representative examples of individual gene sets from the C. elegans and D. melanogaster gene set files. Official gene symbols are shown where available, locus tags (C. elegans) where gene symbols are not available. As in other gene set collections, as the number of genes in any given gene set decreases, the phenotypes progress from broad categories to more specific phenotypic descriptors. The full gene set lists consist of a wide range of developmental, structural, metabolic and behavioral phenotypes, representing a large majority of the experimentally determined phenotypes found in worms and flies. They range from broad phenotype categories such as “sterile”, “slow_growth”, or “larval_arrest” in worms and “viable”, “lethal” and “fertile” in flies; to narrow phenotypic descriptors such as “flaccid”, “DNA_synthesis_variant” or “no_posterior_pharynx” in worms and “ejaculatory_bulb”, “dorsal_vessel_primordium”, or “dense_body” in flies. In addition, there is often overlap of the genes found in related gene sets in both species, emphasizing the contributions of the same genes to multiple phenotypic traits. The complete C. elegans (Additional file 1: Table S1: Additional file 2: Table S2: Additional file 3: Table S3) and D. melanogaster (Additional file 4: Table S4) gene set files are available at this address http://www.grc.nia.nih.gov/branches/rrb/dna/index/Worm-fly_gene_sets_5-9-12.html.

Table 1 Selected Phenotype gene sets

General uses of phenotype based gene sets in both worm and fly

As described here, a single gene set is essentially a single phenotypic term followed by a single row of genes that have been associated with that phenotype. A collection of gene sets consists of a list of phenotypic terms with their corresponding gene sets. Gene sets can be used individually, as a collection, or compared across collections in a number of ways including network analysis, genome-wide model representations, hierarchical clustering, gene set analysis (GSA) of microarray data, and principal component analysis (PCA) of gene set values; among others. A property of this collection of gene sets is that they describe complex intermediate and end stage phenotypes as opposed to molecular function or lists of coordinately regulated genes. They can be used in a variety of bioinformatics applications to reveal higher order or emergent biological and phenotypic relationships and to provide insight into the biological relevance of complex molecular datasets.

Network analysis

Each individual gene set can be used to build networks to determine transcriptional regulation or protein-protein interactions. Figure 1 is a representative network of six networks showing regulatory relationships analyzed by Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, http://www.ingenuity.com) from a single 169 gene, C. elegans gene set, “life span variant”, found in the worm CE-Combined-GS 7-28-2011 file. This analysis identifies members of the gene set (shaded) as well as regulatory or transcriptional partners not found (unshaded) in the original gene set. This network highlights the central role of insulin, ERK family members, and PI3 Kinase as important contributors to longevity in worms.

Figure 1

Network Diagram of the single C. elegans gene set “life_span_variant” as produced by Ingenuity Pathway Analysis (IPA) shows members of the geneset (shaded) as well as regulatory or transcriptional partners not found (unshaded) in the original geneset. The network highlights the central role of insulin, ERK family members, and PI3 Kinase as important contributors to longevity in worms.

An example of a network showing regulatory relationships from a single 82 gene “long_lived” gene set, found in the fly gene set file (DM-narrow-GS 9-7-2011), is also shown in Figure 2. Like in the worm, insulin is central in this fly network, as well as ERKs, AKT, and histones, demonstrating significant overlap in age related biochemical pathways between worms and flies. Each individual gene set (one phenotype with one row of annotated genes) produces multiple network diagrams showing the transcriptional neighbors and protein-protein partners of the core genes, while the entire collection of thousands of gene sets would produce many thousands of individual networks relative to phenotypic descriptions.

Figure 2

Network Diagram of the single D. melanogaster gene set “long_lived” as produced by Ingenuity Pathway Analysis (IPA) as in the worm, insulin is central in this fly network, as well as ERKs, AKT, and histones, demonstrating significant overlap in age related biochemical pathways between worms and flies.

Genome-wide phenotypic modeling in worms and flies

In addition to analysis of a single gene set, a collection of phenotypic gene sets can be compared to itself to reveal biological relationships between all members of the collection. Figure 3 shows a dendrogram of the combined C. elegans file (CE-Combined-GS), using gene sets, having three or more genes, compared to each other based on the degree of gene sharing between individual gene sets. The overall worm tree (Figure 3) is composed of eleven large branches enriched for related biological functions. Moreover, local relationships within a specific branch suggest functional relationships between closely spaced individual gene lists. For instance, in branch 2 (Additional file 5: Figure S1) cell cycle phenotypes such as “cell cycle timing”, “cell cycle delayed” and “cell cycle variant” are closely positioned in space and close to spindle assembly phenotypes. Likewise, in branch 6 (Additional file 6: Figure S2) Dauer phenotypes are closely aligned with multiple lifespan phenotypes based on individual gene sharing within their respective gene sets. Close apposition of related phenotypes as determined by gene sharing between gene sets is a pervasive feature of these dendrogram displays and represents overlap of related phenotypes being influenced by shared genes.

Figure 3

Genome-wide modeling of 931 C. elegans gene sets, each with a minimum of 3 genes shows eleven large branches enriched for related biological functions. Local relationships within a specific branch suggest functional relationships between closely spaced individual gene lists.

The Drosophila gene set collection also produced a similar complex dendrogram of phenotypic functional groups based on gene sharing between gene sets (Figure 4). Like the worm dendrogram, individual branches of the fly dendrogram display a functional relatedness within subregions in each branch. For example, chromosome related phenotypes are grouped in branch 2 (Additional file 7: Figure S3) with mitotic and meiotic phenotypes, including meiotic telophase phenotypes, being closely aligned to each other, as well as spermatid and spermatocyte phenotypes. Behavioral, neuronal, and sensory response phenotypes are shown closely aligned in branch 11 of Figure 4 (Additional file 8: Figure S4), demonstrating overlapping genetic control of related complex phenotypes.

Figure 4

Genome-wide modeling of 1,503 D. melanogaster gene sets, each with a minimum of 5 genes showing broad phenotypic groups represented by branches and a functional relatedness within subregions in each branch.

Phenotype Gene Set Analysis (GSA) of microarray data and Principal Components Analysis (PCA) of gene sets

C. elegans: In addition to comparisons of gene sets either individually or collectively to themselves, these phenotype gene sets are useful in analysis of microarray based gene expression datasets in worm and fly. Figure 5a illustrates statistically significant gene sets resulting from Gene Set Analysis (GSA) of a single 4 day old larva versus 15 day old whole genome gene expression comparison in a C. elegans aging microarray dataset [15]. This dataset (GEO # GSE21784) represents a 15 day time course with incremental stages of infection with P. areuginosa. Statistically significant up-regulated gene sets include germ cell gene groups, as well as meiosis and cell division gene sets, among others. Down-regulated gene sets include gene groups involved in body vacuoles, as well as alae and cuticle formation. Figure 5b is a heat map of the significant changes across the entire time course.

Figure 5

a Phenotype gene set analysis of whole genome microarray expression values from 4 day larva versus 15 day C. elegans (GEO# GSE21784). This was performed using the combined C. elegans gene set file on the Disease/Phenotype WEB-PAGE GSA web tool found here: http://dpwebpage.nia.nih.gov/ which shows age-related changes in the worm. b Heatmap analysis of whole genome microarray expression values from 4 day larva to 15 day C. elegans (GEO# GSE21784). This was performed using the combined C. elegans gene set file on the Disease/Phenotype WEB-PAGE GSA web found here: http://dpwebpage.nia.nih.gov/ which shows age-related changes in the worm.

Figure 6 shows changes in selected gene sets from a different aging time course in C. elegans over 24 days [16] (GEO # GSE12290). Aging related increases (Figure 6a) or decreases (Figure 6b) in gene groups related to locomotion, energy metabolism, and life span are highlighted.

Figure 6

Selected changes in gene sets over a C. elegans aging time course. (GEO# GSE12290) showing increases (a) and decreases (b) in gene set scores over time. The X axis shows days 4,8,12,14,16,20,24 and the X axis shows Z-score for each gene set shown.

In addition to GSA of microarray data the gene set values derived from gene expression data can be further analyzed by principal components analysis (PCA) using the Z-score values of the original gene set data output. This is in contrast to more commonly described PCA resulting from individual gene expression values. Figure 7 shows tight grouping of individual biological samples within three groups; larvae, adult day 6, and adult day 15, and dramatic separation of time points within the experiment, based solely on PCA analysis of the gene sets values from the previous gene set analysis. This demonstrates that there is useful biological information content in the aggregate gene set results, in addition to that found in any individual gene set, which can discriminate between discrete biological states.

Figure 7

Two dimensional Principal Components Analysis (PCA) of gene set results from whole genome microarray expression values from C. elegans larva, 6 day adult worms, and 15 day Adult C. elegans (GEO# GSE21784) showing age as the principal variation component.

D. melanogaster: In a similar fashion to the worm (above), microarray data from young versus aging flies was analyzed with the Drosophila gene set file DM-narrow-GS containing 11,999 gene sets. Gene set analysis was performed using the WEB-PAGE gene set analysis tool [10] on a dataset of gene expression values from young versus old flies [17] (GEO# GSE22437). The top 100 statistically significant enriched gene sets using Z ratios of the expression values from day 10 versus day 40 fly heads is shown in Figure 8. Over enriched gene sets include minute phenotypes, life span, as well as developmental growth rate phenotypes, among others. The discriminative ability of PCA using gene set Z-scores (as opposed to individual gene values) is illustrated using the individual samples of day 10 versus day 40 fly heads in Figure 8.

Figure 8

a Phenotype gene set analysis of average whole genome microarray expression values from 10 day versus 40 day D. melanogaster heads (GEO# GSE22437). This was performed using the D. melanogaster gene set file on the WEB-PAGE GSA web tool showing age-related changes in the fly. b Principal Components Analysis of gene set results from of whole genome microarray. Gene expression values from D. melanogaster 10 day and 40 day Drosophila heads (GEO# GSE22437) showing age as the principal variation component.

Conclusion

Here we describe genome-wide phenotypic modeling using gene sets based on gene-phenotypic assignments in C. elegans and D. melanogaster. Unlike previous gene set collections such as KEGG, GO, MSigDB, in these and other species, every gene in every gene set described here is based on genetic evidence contributing to each specific phenotype. Although very useful, these gene sets should be considered a first generation. They may not be complete. Some may describe certain phenotypes in different developmental contexts, or in particular applications and not in others. In addition, many subtleties and details were not included in deriving these gene sets including penetrance of different alleles, strain differences, and environmental modifiers. Moreover, these gene sets may produce different results depending on the statistical algorithms used in complex analysis.

However, we have demonstrated these gene sets can be used to identify complex higher order biological and genetic relationships through network analysis, whole genome phenotypic modeling, and analysis of complex molecular datasets. They will help elucidate complex multigenic relationships between genes and phenotypes in worms and flies in many experimental and biological contexts and will provide a bridge for phenotypic comparisons between model and intermediate species.

Methods

Derivation of phenotypic gene sets

Worm

Phenotype-gene lists obtained from WormBase on 4/24/11 were titled RNAi and VAR. RNAi, consisted of 34,433 gene phenotype pairs having 7,289 unique genes and 850 unique phenotypes. These phenotypes were the results of observations of phenotypes from knockdown of the gene products (RNAi experiments). The list VAR contained 8,440 records, having 2,165 unique genes, and 1,109 unique phenotypes and was the result of observations of phenotypes from genetic mutations as deposited in WormBase. The overlap between each file consists of 1,410 genes and 237 phenotypes.

Phenotype gene set files were created by parsing the original gene lists into non-redundant phenotype lists with annotated genes using a custom Perl script as previously described [9]. This was done for RNAi and VAR independently, as well as combined to create the gene set files; CE-RNAi-GS 7-26-11, CE-VAR-GS 7-26-11, and CE-Combined-GS 7-28-11. The resultant individual Phenotype Gene set names are identical to the Phenotype descriptors found in the original WormBase Phenotype file. These files can be downloaded here: http://www.grc.nia.nih.gov/branches/rrb/dna/index/Worm-fly_gene_sets_5-9-12.html.

Fly

Phenotypes and gene assignments were obtained from FlyBase on 9-11-11 at this web address: http://FlyBase.org/static_pages/downloads/FB2011_07/alleles/allele_phenotypic_data_fb_2011_07.tsv.gz. This file began with 259,162 phenotypic descriptions with assigned Drosophila genes. Redundant phenotype-gene combinations were removed resulting in a list of 154,428 unique phenotype-single gene pairs. Parsing of this file resulted in a non-redundant gene set file of 11,999 unique phenotypic terms with annotated genes. The resultant individual Phenotype Gene set names are identical to the Phenotype descriptors found in the original FlyBase Phenotype file. This Phenotype gene set file named DM-narrow-GS 9-7-2011, can be downloaded here: http://www.grc.nia.nih.gov/branches/rrb/dna/data/worm-fly/DM-narrow-GS_9-7-2011.txt.

Gene set nomenclature

It should be noted that nomenclature of many phenotype gene sets in both worm and fly often have a directionality in the name which may or may not be relevant to any given microarray or other analysis. Please see Additional file 9: S7 for an explanation of directionality in gene set nomenclature and interpretation in their use.

Network analysis

Networks for C. elegans and D. melanogaster were produced using Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, http://www.ingenuity.com). Using the “life_span_variant” gene set in C. elegans generated on 7-26-2011, and the “long_lived” gene set in D. melanogaster generated on 12-07-2011. The input and output files can be found here for C. elegans (Additional file 10: Table S5) and D. melanogaster (Additional file 11: Table S6).

Genome-wide phenotypic modeling

Genome-wide dendrograms were produced by a unique method similar to phylogenetic classification as previously described [9]. Briefly, the distance between each phenotypic gene set was calculated by pairwise comparison of every gene set pair by finding the number of common genes between each pair and dividing that number by the number of genes in the smallest group of the pair, resulting in a correlation value between 1 and 0 for each pair. This was done for all gene sets to produce a distance matrix. This number was then subtracted from 1because if two lists are identical (100 % match) then the resultant distance should be 0. This is represented as:

$$d_{i,j}=1-\frac{N\left(C_i\cap C_j\right)}{min\left[N\left(C_i\right),N\left(C_j\right)\right]}$$

when i ≠ j. If i = j then d = 0

Where: C _k : Genes in each disease set (where k = i,j) ; N(C _k ): Number of genes in each disease set (where k = i,j) ; d_ij is the pairwise distance ; i,j: index of genes in each disease set where; i = 1,2,3,………,n ; j = 1,2,3,………,m.

The gene set relationships were calculated from the distance matrix using the Fitch program [18]. It calculates the relationships based on the Fitch and Margoliash method of constructing the phylogenetic trees using the following formula (from the Phylip manual):

$$\mathit{Sum}\_\mathit{of}\_\mathit{squares}={\displaystyle \sum _i{\displaystyle \sum _j\frac{n_{\mathit{ij}}{\left(D_{\mathit{ij}}-d_{\mathit{ij}}\right)}^2}{D_{\mathit{ij}}{}^P}}}$$

where D is the observed distance between gene sets i and j and d is the expected distance, computed as the sum of the lengths of the segments of the tree from gene set i to gene set j. The quantity n is the number of times each distance has been replicated. In simple cases n is taken to be one. If n is chosen more than 1, the distance is then assumed to be a mean of those replicates. The power P is what distinguished between the Fitch and Neighbor-Joining methods. For the Fitch- Margoliash method P is 2.0 and for Neighbor-Joining method it is 0.0. The resulting coefficient matrix file was displayed using the Phylodraw graphics program [19].

Gene set analysis

This analysis used the Disease/Phenotype WEB-PAGE GSA web tool[10] using the PAGE algorithm [2] with the CE-Combined-GS gene set file excluding gene sets containing over 500 and less than 3 genes. Briefly, for each gene set a Z score was computed as, $Z_{\mathit{phenotype}}\left(i\right)=\frac{\sqrt{n_i-1}•\mathit{dif}{f}_i}{{\sigma }_a}$ In which the phenotype index i = 1,2,…,K; where K is the total number of the disease phenotypes we included in our data set; n_I is the number of genes in the sub-group of phenotype i in the current sample array; σ_A is the standard deviation of the current gene expression changes of the sample. Diff(i): is the difference between the mean value of gene expression changes in the subgroup disease phenotype (i) (GC_I ) and the mean value of the gene expression changes on the whole sample (GC_A) i.e. $\mathit{dif}{f}_i=\overline{G{C}_i}-{\overline{\mathrm{G}C}}_{\mathrm{a}}$. The empirical p-value of the disease phenotype i changes is described by: $p\left(i\right)=2\left[1-\Phi \left(\frac{\mathit{dif}{f}_i}{\sigma \left(\mathit{dif}{f}_i\right)}\right)\right]$ in which Φ(x) is the standard normal distribution function with the variable as X = DIFF_I/σ(DIFFF_I). σ(DIFF_I) is the standard deviation of the difference for gene expression changes between phenotype subgroup (i) and the whole array $\sigma \left(\mathit{dif}{f}_i\right)=\sqrt{\frac{{\sigma }_i{}^2}{n_i}+\frac{{\sigma }_a{}^2}{n_a}}$ σ_I is the standard deviation of the average gene expression changes in the disease phenotype (i). N_A is the total number of genes in the whole sample set. The plots were drawn with R-statistical programming language (R Development Core Team 2005) using either calculated or absolute z-score values.

Principal components analysis

Principal components analysis was performed on the gene set Z values using DIANE 8.0 a JMP based software package (http://www.grc.nia.nih.gov/branches/rrb/dna/diane_software.pdf) based on the Singular Value Decomposition (SVD) function in JMP 9.0. In short, the data was organized as m × n matrix where m is the different samples (columns) and n is gene set Z-values (rows), mean of each row was subtracted and SVD was calculated using JMP’s in-built SVD function as illustrated in this document: http://www.cs.princeton.edu/picasso/mats/PCA-Tutorial-Intuition_jp.pdf and also used in this script: http://abs.cit.nih.gov/MSCLtoolbox.

Data access

The complete C. elegans and D. melanogaster gene set files are available at this address: http://www.grc.nia.nih.gov/branches/rrb/dna/index/Worm-fly_gene_sets_5-9-12.html.