Properties of genes essential for mouse development

Essential genes are those that are critical for life. In the specific case of the mouse, they are the set of genes whose deletion means that a mouse is unable to survive after birth. As such, they are the key minimal set of genes needed for all the steps of development to produce an organism capable of life ex utero. We explored a wide range of sequence and functional features to characterise essential (lethal) and non-essential (viable) genes in mice. Experimental data curated manually identified 1301 essential genes and 3451 viable genes. Very many sequence features show highly significant differences between essential and viable mouse genes. Essential genes generally encode complex proteins, with multiple domains and many introns. These genes tend to be: long, highly expressed, old and evolutionarily conserved. These genes tend to encode ligases, transferases, phosphorylated proteins, intracellular proteins, nuclear proteins, and hubs in protein-protein interaction networks. They are involved with regulating protein-protein interactions, gene expression and metabolic processes, cell morphogenesis, cell division, cell proliferation, DNA replication, cell differentiation, DNA repair and transcription, cell differentiation and embryonic development. Viable genes tend to encode: membrane proteins or secreted proteins, and are associated with functions such as cellular communication, apoptosis, behaviour and immune response, as well as housekeeping and tissue specific functions. Viable genes are linked to transport, ion channels, signal transduction, calcium binding and lipid binding, consistent with their location in membranes and involvement with cell-cell communication. From the analysis of the composite features of essential and viable genes, we conclude that essential genes tend to be required for intracellular functions, and viable genes tend to be involved with extracellular functions and cell-cell communication. Knowledge of the features that are over-represented in essential genes allows for a deeper understanding of the functions and processes implemented during mammalian development.

Introduction Essential genes are those whose presence is imperative for the survival of an organism. However, the complete set of genes that are absolutely vital to sustain life are still unknown for most organisms [1]. In mammals, knowledge of essential genes is required to understand development, maintenance of major cellular processes and tissue-specific functions that are crucial for life. As such, essential genes are the key minimal set of genes needed for all steps of development. Genes that are not needed for development are termed non-essential or viable genes. Mammalian essential genes can be identified using experimental techniques [2], which include single gene knockouts [3][4][5], conditional knockouts [6,7], forward genetic screens [8], RNA interference [9,10], and transposon mutagenesis [11]. Though these experimental methods are the gold standard, they are time consuming and expensive. Nevertheless, major programs are currently underway to systematically knockout every mouse gene and characterise the resulting phenotypes [12]. An initial set of essential genes has been identified through these experimental approaches [13,14].
These recent data offer a valuable resource to help understand which processes are critical for mammalian development and to discover what makes a gene essential or viable. We hypothesised that essential and viable genes are distinguishable by various attributes. We explored a wide range of sequence and functional features of mouse genes in order to characterise essential and viable genes in mammals. We have discovered numerous gene and protein features that vary significantly between essential and viable genes in mouse, some of which were previously found to be associated with essentiality in E. coli [15,16], S. cerevisiae [17][18][19], mouse [20] and human [21]. These features thus reveal the key genetic functions required for development in mammals.

Datasets
The Mouse Genome Informatics (MGI) database [22] incorporates published gene data on mouse knockout phenotypes. We collected a total of 1,271 essential and 4,378 viable mouse genes from MGI (accessed on 1 November, 2013), based on phenotype annotations of null alleles of targeted deletion knockout mice. Mutant phenotypes generated from other experimental methods were not included in our dataset, since we could not exclude the possibility that essential genes might have hypomorphic alleles with viable phenotypes in gene trap, knockdown, or chemical mutagenesis experiments. We considered a gene as essential if it produced lethality in either the heterozygous or homozygous state, and did not differentiate between these two types of genes in the dataset. We defined essential genes as those that are required for an animal to survive past post-natal day 3. A total of 1,335 genes had both 'essential' and 'viable' annotations in the MGI database so were individually checked in the literature. We further manually checked each gene to ensure that our datasets contained only protein-coding genes, to allow for an analysis of features specific to protein function. This resulted in a total dataset of 1,301 essential and 3,451 viable mouse genes (S1 Data and S2 Data).
The proteins encoded by these essential and viable genes share significant levels of sequence identity. A protein sequence dataset is considered redundant if it includes a pair of proteins that are highly similar or homologous. The presence of redundancy is a barrier in using a dataset effectively as it increases the size of the dataset; also it could potentially create bias towards any conclusions drawn from the overall analysis using the dataset due to the over-representation of similar features. This problem can be overcome by removing redundant proteins from the dataset until all the proteins in the dataset share sequence similarity less than a predefined threshold. We therefore used Leaf [23] to remove redundant proteins from our datasets. We generated four culled or non-redundant essential and viable datasets from our original dataset, where the sequence similarity between all proteins is less than a threshold of 20%, 40%, 60% and 80%, respectively ( Table 1).

Analysis of genomic features
The functionality of a gene may rely on its inherent sequence features at the genomic level. Analysing these gene sequence based features may provide valuable insights into their contributions to gene essentiality.
GC content, gene length, and transcript diversity. We anticipated that genomic features such as gene length and GC content could be indicative of gene essentiality, and determined if these features differ between our essential and non-essential gene sets. We found that essential genes tend to be longer in length compared to viable genes ( Table 2; Fig 1A). Total gene length is comprised of individual exon and intron lengths. We therefore also measured these features in our datasets, finding that essential genes tend to have longer exons and introns than viable genes ( Table 2; Fig 1B and Fig 1C).
We also examined transcript diversity, finding that essential genes tend to have more transcripts than viable genes ( Table 2; Fig 1D). To quantify whether or not the number of exons could differentiate between essential and viable genes, the ranking of the number of exons from the longest transcript of each gene was analysed. We found that essential genes are likely to have more exons than viable genes ( Table 2; Fig 1E). Here, EN and VN refer to essential and viable genes in the non-culled dataset. Ex and Vx define essential and viable genes in the culled dataset where all coded proteins share sequence similarity less than x0%. In this box plot, the top and bottom of the box denote the upper and lower quartiles; the line inside the box denotes the median; and individual points denote the outliers. Top 5% essential and viable genes with longest gene length (A) and longest introns (F) were excluded from the datasets to make plots more readable.
When the distributions of GC content in essential and viable genes were examined, we observed that viable genes have a higher percentage of GC content only for the culled dataset where all coded proteins have a sequence identity < 60%; this observation was not statistically significant for other culled datasets (Table 2; Fig 1F).
Gene expression. Examining the temporal specificity of gene expression can identify genes that are active in a particular biological process. We therefore expected that expression could serve as an important indicator of essentiality, as developmentally essential genes should be expressed during embryonic development. We obtained mouse gene expression data for 1,301 essential and 3,409 viable genes from the UniGene database [24] covering 13 developmental stages. Essential genes are more highly expressed than viable genes at every stage of mouse development (Fig 2). However, the χ 2 tests with the Bonferroni correction analysis showed that these differences are not statistically significant at later stages of development (juvenile and adult), as nearly all genes are expressed at those stages (Table 3). Essential genes were found to be highly expressed, whereas viable genes are more likely to be found in the group of genes with zero transcripts present in developmental samples (Fig 3). Evolutionary age. The evolutionary age of a gene represents the time that has passed since the gene evolved from its ancestor, either by duplication or speciation. Studies in bacteria and yeast found essential genes to be evolutionarily more conserved than viable genes [5,15,25]. We therefore also expected that gene evolutionary age could be informative for distinguishing mammalian gene essentiality.
For mammalian genes that have been duplicated, the evolutionary age reported in millions of years ago (MYA) of the duplicate common ancestor (DCA) and the most recent duplication (MRD) event were collected from the Ensembl (release 75) gene trees. For mammalian genes without duplicates, the gene age was determined to be that of the singleton common ancestor (SCA). We observed 16 representative phylogenetic age groups for our mouse genes (Table 4A). We found ages for 1,276 (98.1%) essential and 3,358 (97.3%) viable genes. The oldest genes arose approximately 1215 MYA, whereas the youngest genes belong to the class Murinae arising approximately 25 MYA. We compared the enrichment of essential and viable genes in different age groups. We found that essential genes tend to be older than viable genes for both non-culled and culled datasets (Fig 4). We observed that a significantly greater percentage of essential genes have evolutionary origins of 1215 and 937 MYA, compared to viable genes in the non-culled dataset (Table 4B and Table 4C). The majority of the viable genes arose 400 MYA. Using MRD ages, we found that viable genes are more likely to have ages of 25 and 162 MYA (Table 4B). We further observed a significantly greater percentage of viable genes that have DCA ages arising at 296, 371, 414 and 535 MYA (Table 4C). We found similar trends for the culled datasets, which further confirms that genes essential for mouse development are more evolutionarily ancient.

Analysis of protein features
Prior research established that different physical, functional and evolutionary properties of proteins can facilitate the prediction of gene essentiality [15,18,20,26]. Here, we explore a number of protein properties, obtained from mouse protein sequence data, to test their efficacy at distinguishing essential genes from viable genes in mouse.
Simple sequence features. We found that essential proteins have significantly longer lengths than viable proteins (529aa versus 452aa (median length); p-value = 1.03×10 −21 Mann-Whitney U test). The distributions of protein lengths between essential and viable proteins within the non-culled and culled datasets are variable and discriminate between these classes ( Fig 5A). We also found variations in the frequencies of amino acids found in the proteins encoded by essential and viable genes (Table 5). Proteins encoded by essential genes in the non-culled dataset tend to have higher proportions of Ala, Asp, Glu, Lys, Gln and Ser. Distributions of Lys residues demonstrated the same trend for all culled datasets. Essential proteins in the 40%, 60% and 80% culled dataset also had more Asp, Glu and Gln compared to viable proteins (S3 Data). Viable proteins have more Leu, Cys, Phe, Val and Trp.
Protein average molecular weight, charge, isoelectric point and frequencies of different amino acid categories were computed using the tool Pepstats [27]. Proteins encoded by essential genes have a significantly higher average molecular weight (MW) compared to proteins encoded by viable genes (Table 6). Differences for charge, isoelectric point, tiny and small residues were not statistically significant. Essential proteins were found to have greater proportions of polar (Fig 5B), charged (Fig 5C), basic ( Fig 5D) and acidic ( Fig 5E) amino acids. In contrast, proteins encoded by viable genes have significantly higher proportions of aliphatic ( Fig 6A), aromatic ( Fig 6B) and non-polar residues ( Fig 6C). However, the Mann-Whitney U test showed that differences of aliphatic (p-value = 0.42) and aromatic (p-value = 0.19) residues between essential and viable proteins in the 20% culled datasets are not statistically significant (Table 6). Enzyme class. Almost all cellular processes are dependent on the presence of enzymes. Enzymatic function thereby could be another indicator of gene essentiality. We extracted the annotations of the six primary enzyme classes from UniProt [28] and counted the number of essential and viable proteins belonging to each of these classes (Table 7). In the non-culled datasets, 29.8% (388/1301) of the total number of essential proteins exhibit enzymatic activity compared to 27.7% (956/3451) of viable proteins, though this difference is not statistically significant. The culled datasets also show variations within each class in the percentage of proteins that function as enzymes ( Table 7). Analysis of the enzyme classifications shows that the proteins encoded by essential genes are rich in transferases and ligases as compared to those encoded by viable genes. Hydrolases were found to be strongly associated with viable proteins in the non-culled dataset. No statistically significant differences between the datasets were observed for oxidoreductases, lyases and isomerases.
Post-translational modifications and transcription. We investigated the frequency of annotations for four different post-translational modification keywords ('phosphoprotein', 'glycoprotein', 'acetylation' and 'transcription') as obtained from UniProt protein annotations. Protein phosphorylation plays crucial roles in regulating various cellular and metabolic processes, such as cell differentiation, cell division, survival etc. Around 30% of all eukaryotic proteins are estimated to be phosphorylated [29]. We found that essential proteins within the non-culled dataset are significantly more likely to be phosphorylated than viable proteins (51.42% versus 35.50%, p-value = 8.93×10 −15 ). We observed the same trend for culled datasets (Table 8). Greater than 50% of all proteins are glycosylated [30]. Glycoproteins are crucial for protein folding, solubility and localization [31]. A large number of them are secreted extracellular proteins, or are cell membrane proteins, and they therefore have roles in transport and cell-cell interactions. Viable protein are significantly more likely to be N-linked glycoproteins than essential proteins (Table 8).
Acetylated proteins in eukaryotes are those proteins that are post-translationally modified by the addition of an acetyl group at the N-terminus or on Lys side chains. The acetylation process is important for gene expression and metabolism. N-acetylated proteins also have vital roles in regulation of protein-protein interactions [32]. Proteins encoded by essential genes are more likely to have at least one acetyl group than proteins encoded by viable genes for all Here, EN and VN refer to essential and viable genes in the non-culled dataset. Ex and Vx define essential and viable genes in the culled dataset where all coded proteins share sequence similarity less than x0%. Top 2% longest proteins (A) were excluded from the datasets to make plots more readable. In this box plot, the top and bottom of the box denote the upper and lower quartiles; the line inside the box denotes the median; and individual points denote the outliers.
https://doi.org/10.1371/journal.pone.0178273.g005 datasets ( Table 8). Essential proteins in all datasets are thus more likely to be associated with regulating the transcription of genes, since acetylation is used to control gene expression.
Signal peptides. Signal peptides are short peptide sequences (usually 5-60 amino acids long) located at the N-terminus of a large number of newly synthesized proteins. They control Table 5. Differences in the frequency of usage of the 20 amino acids between essential and viable mouse proteins in the non-culled dataset. The p-value for the Bonferroni correction is 0.0025.  Here, EN and VN refer to essential and viable genes in the non-culled dataset. Ex and the targeting and translocation of secreted or cell membrane proteins. Signal peptides direct proteins to different cellular locations (e.g. nucleus, mitochondria, endoplasmic reticulum, endosome, Golgi apparatus). Signal peptide motifs (computed using UniProt annotation and SignalP servers [33] are significantly more frequent in proteins encoded by viable genes as compared to proteins encoded by essential genes (Table 9). Transmembrane domains. Transmembrane proteins extend through the lipid bilayer and span from the interior to the exterior of the cell. Transmembrane proteins usually adopt an α-helical structure while passing through the lipid bilayer once (single-pass proteins) or multiple times (multiple-pass proteins). Due to this structure, transmembrane proteins can mediate cellular functions both inside and outside of the cell. Transmembrane proteins are important for cell-cell communication, maintenance of cell structure, signalling, and ion transport. Many receptor proteins have a number of α-helical transmembrane domains spanning the cell membrane. Thus, the presence of transmembrane domains in protein encoded by essential and viable genes could be informative for functional annotation.

Amino acid
We found that the non-culled viable dataset is significantly enriched in transmembrane proteins (p-value = 1.86×10 −15 ). Approximately 20% of essential proteins are annotated as transmembrane proteins, whereas the corresponding percentage is 34% for viable proteins. A total of 10.5% essential proteins consist of a single transmembrane helix, whereas this number is 17% for viable proteins. Also, 2% of essential proteins have seven transmembrane helices, compared to 6.5% of viable proteins. Overall, a greater number of viable proteins have transmembrane domains, and viable proteins have significantly more transmembrane helices per protein than essential transmembrane proteins.
Vx define essential and viable genes in the culled dataset where all coded proteins share sequence similarity less than x0%. In this box plot, the top and bottom of the box denote the upper and lower quartiles; the line inside the box denotes the median; and individual points denote the outliers.
https://doi.org/10.1371/journal.pone.0178273.g006 Table 7. Differences in the frequencies of different enzyme class observed between essential and viable mouse proteins. The Bonferroni corrected p-valu for the Chi-squared tests is 0.0083. Gene Ontology (GO) [34] is the most widely used scheme for classifying gene functions. The GO consortium provides a set of controlled vocabularies (ontology) to annotate the functional properties of gene and gene products across all species. Gene functions are annotated by means of three aspects: (a) molecular function (b) cellular component and (c) biological process. Here, we test whether GO term distributions vary between essential and viable genes. Cellular component. Protein functions are closely related to the locations where they reside within a cell. Subcellular localisation has been shown to be important for predicting essential genes in prior studies [15,18,35]. As an example, eukaryotic proteins located in the nucleus carry out essential functions including DNA replication, mRNA synthesis and recombination. Subcellular localisation therefore should be useful in distinguishing mouse essential genes.

Datasets
GO terms were extracted from the DAVID v6.8 functional annotation tool [36] by submitting the Ensembl IDs of mouse essential and viable genes. A total of 225 cellular component GO terms for essential genes and 149 terms for viable genes were retrieved, of which 53 and 82 terms were found significant, utilising the Bonferroni corrected p-value 0.05 from the functional annotation output of DAVID. Tables 10 and 11 summarise these cellular component GO terms favoured for essential and viable genes, respectively. Lists of the 50 most enriched GO terms for each class are listed in S4-S9. A majority of essential genes are intracellular. Terms most frequently associated with essential genes include: "nucleus", "transcription factor complex", "nucleoplasm", "nucleolus", and "intracellular membrane-bounded organelle". Fifty-seven percent of total essential genes were found to be present in the nucleus. In contrast, as shown by our analysis of transmembrane helices, many viable genes are membrane bound. Viable genes were enriched for cellular component terms including "membrane", "plasma membrane", "cell surface", "extracellular region", "extracellular space" and "lysosome". A high percentage of essential (46%) and viable (41%) genes were also found with the annotation of cytoplasm. Notably, an individual protein can have more than one subcellular localisation annotation.
Subcellular locations were also analysed using the UniProt annotation and the WoLF PSORT tool [37]. Table 12 summarises the results of the UniProt analysis. We found that a significantly higher proportion of viable proteins are localised in plasma membrane (23%), membrane (15%) and extracellular region (14%), compared to essential proteins. A higher percentage of essential proteins are found within the nucleus (48%), as compared to viable proteins (23%). The same trend was also observed for culled datasets.
Subcellular location prediction results from WoLF PSORT are summarised in Table 13. While the absolute numbers are often different from UniProt, the trends in differences between essential and viable are similar. In this case, the most significant enrichment for subcellular localisation of essential proteins was the nucleus. We observed that 70% of total essential proteins are located in the nucleus compared to 49% of viable proteins. The analysis of WoLF PSORT prediction results further confirmed the tendency for viable genes to be membrane bound (36%) and extracellular (39%). Viable proteins were also enriched for localisation to endoplasmic reticulum (18%) and lysosome (11%) as compared to essential proteins. These analyses of cellular localisations indicate that proteins encoded by essential genes are commonly located in the nucleus, whereas viable proteins are more likely to be extracellular or membrane bound. Viable proteins are also more likely to be located in the lysosome.
Biological processes. A total of 1,575 biological process terms were retrieved for essential genes, with 1,777 terms for viable genes, of which 323 terms for essential and 315 terms for viable datasets were significant meeting the Bonferroni corrected p-value 0.05. Table 14 lists the top 20 biological process terms significantly favoured for essential genes. Essential genes are often involved in developmental processes, as expected (Table 15). Significant enrichment for processes related to "transcription", "cell proliferation", "cell differentiation", "organ morphogenesis", "cell division", and "DNA replication" is observed in the essential genes dataset. Biological process terms commonly annotated for viable genes include "inflammatory response", "signal transduction", "ion transport", "immune response", "response to drug", "response to stimulus", "behaviour", "transmembrane transport", "aging" and "regulation of apoptotic process" (Table 15). Molecular function. Analysing the molecular function output generated by DAVID, a total of 265 terms for essential genes and 105 terms for viable genes were retrieved, of which 75 and 81 terms were significant, respectively (Tables 16 and 17). Essential genes are more likely to be annotated as being involved in "DNA binding", "transcription factor activity", "transcription factor binding", and "transferase activity". Viable genes are more likely to have the annotations of "signal transducer activity", "ion channel activity", "hydrolase activity", "transporter activity", "calcium ion binding", "receptor binding", "SH3 domain binding", and "lipid binding". A higher percentage of essential and viable genes were also found to be annotated as being involved in "protein binding", "ATP binding", "protein kinase binding", and "protein kinase activity".

Protein domains
Protein domains are spatially distinct structural and/or functional units of a protein. They carry out particular functions or interactions, thereby contributing towards the overall functionality of a protein. We obtained domain data for essential and viable mouse proteins by analysing the functional annotation output of DAVID (Tables 18 and 19). We observed a total of 11 and 30 Pfam domains [38] that are significantly enriched in essential and viable proteins, respectively. Domains such as homeobox, helix-loop-helix DNA-binding domain, T-box, protein kinase domain, Zinc finger, and C4 type domain (many of which are found in transcription factors) showed enrichment in essential proteins. Domains including 7-transmembrane receptor, SH2, ion transport, Fibronectin type III domain (fn3), and SH3 (many of which are found in membrane proteins) were more frequently found in viable proteins. Although viable proteins were annotated with having protein kinase and zf-c4 domains, these domains were more frequently found within essential proteins.

Protein-protein interactions
Protein-protein interactions (PPI) are intrinsic to almost all biological processes. Since the majority of proteins interact with each other to expedite accurate functionality, knowledge about their interactions is crucial to understand the molecular mechanisms of cellular processes. A prior study found significant differences in PPI network properties between the essential and viable genes of S. cerevisiae and E. coli [39]. Network-based attributes were also found to be fundamental to elucidate proteins activities within the cell [21]. We therefore expected that the study of PPI networks could be an indicator of essentiality of mouse proteins. Mouse protein-protein interaction data was obtained from the I2D database [40]. The PPI data was examined with the intention of learning whether essential PPI networks differ in their network properties from their viable counterparts. We analysed both known and predicted mouse PPIs to ensure high quality PPIs. Two PPI networks namely Known (K) and Known-Predicted (KP) were constructed from all mouse PPIs. After removing self and duplicate interactions, the network of proteins encoded by essential genes (essential-K) contained 3,988 protein nodes and 8,074 interactions; the network of proteins encoded by viable genes  Our results demonstrated that essential proteins have more interactions (higher degrees) than viable proteins in both K and KP interaction networks (Fig 7, Table 20). The mean degree of essential proteins was higher than viable proteins for K (10.5 versus 6.4) and KP (57.7 versus 28.0). The Average Shortest Path (ASP) length is an indicator of a protein node's efficiency in transporting information in a PPI network The ASP length of essential proteins is significantly shorter than the ASP length of viable proteins (Fig 8A, Table 20). The betweenness centrality is an indicator of the centrality of a protein node in the PPI network. The betweenness centrality of essential proteins in each of the interaction networks is significantly higher than that of viable proteins (Fig 8B, Table 20). We also found significantly higher clustering coefficient values for essential proteins in K and KP networks compared to viable proteins. Essential proteins tend to have significantly high closeness centrality than viable proteins (Fig 8C). This difference was statistically significant for both networks (Table 21).
We wanted to identify protein nodes with large number of interactions (hubs) in the PPI network. We used the Hub object Analyser (Hubba) [41] to explore four additional network properties including: BottleNeck (BN), Edge Percolation Component (EPC), Maximum Neighbourhood Component (MNC) and Density of Maximum Neighbourhood Component (DMNC). These properties define probable hubs in the PPI network. Our investigation demonstrated that essential proteins tend to have high BN values in both K and KP networks (Fig  9). We further found that EPC and MNC of essential proteins are significantly higher than that of viable genes (Table 21). Although essential proteins exhibited high DMNC in the K network, the same trend was not observed for the KP network.
Housekeeping and enriched genes. Housekeeping genes are expressed at similar levels under all conditions, as they are required to maintain basic cellular functions [42]. In contrast, many genes are expressed only in certain conditions or environments, such as in individual tissues. We used the Pattern Gene Database [43] to test whether lethal and viable genes are also likely to be housekeeping or tissue-specific genes. Table 22 shows that viable genes are significantly more likely to be housekeeping or tissue enriched genes.

Discussion
Understanding what makes a gene essential shows which cellular, developmental and tissuespecific processes are crucial for mammalian development. Our non-culled dataset contained a total of 1,301 essential and 3,451 viable mouse genes, which were obtained from the MGI database. We included only targeted deletion null mouse phenotypes in our analysis. A lengthy literature search corrected numerous inaccurate annotations within this database, giving accurate gene lists to analyse. The presence of multiple copies of similar proteins could bias the analysis; we thereby removed redundant proteins from our dataset to generate non-redundant or culled datasets. Comparing culled and non-culled datasets showed whether redundancy affects particular gene properties, though in general all sets follow similar trends. We studied a wide range of gene and protein properties of Mus musculus genes, representative of different aspects of mouse biology, so that we could quantify their abilities to differentiate essential genes from viable genes. Our investigation focused on features that are attainable from existing databases and web-based tools. These properties fall into three categories: (1) genomic properties, which are based on gene sequence data. This group also included features such as evolutionary age and gene expression; (2) protein sequence properties, determined from protein sequence, including amino acid composition, enzyme class, post-translational modifications, signal peptides and transmembrane domains; (3) functional properties, which facilitate biological interpretations of gene functionality. These include GO annotations and PPIs. In total, we identified 75 features that show significant differences between essential and viable genes. These features, expressing different traits of mouse biology, are interrelated. Many (e.g. gene length, protein length, evolutionary age, gene expression, nuclear localization, PPI) are in broad agreement with those of previous studies on yeast [16][17][18] and bacteria [14,15], but have not been verified in mammals. In addition to previously evidenced features, we found a number of important novel features that are strongly associated with essential genes, summarised in Table 23.
Mouse essential genes are more likely to be longer in length, and have more transcripts than viable genes. Essential genes also tend to exhibit more exons and have a longer exon length. These results are in agreement with a prior study which showed that longer genes with a large number of exons tend to exhibit a higher degree of alternative transcripts compared to smaller genes with fewer exons [44]. Essential genes thus tend to encode complex proteins, having multiple domains and diverse cellular or tissue specialisations [45]. Essential genes also tend to have a significantly longer length of introns and a lower GC content. Intron and exon length is known to vary inversely with GC content [46,47]. GC content is also correlated with gene length [48] and recombination [49] in mammalian genomes.
Essential genes are expressed in greater proportions at the earlier stages of mouse development (pre-organogenesis stages) as compared to viable genes. Mouse genes that are expressed at early stages of development are more likely to be essential, as their disruption could affect downstream developmental events. We found that essential genes show a higher level of gene expression during development. This result is supported by previous studies that showed that highly expressed genes are likely to be essential and to evolve slowly [50,51].
Essential genes tend to have older evolutionary origins than viable genes and thus are evolutionarily more conserved [5,25], consistent with a previously reported observation that essential and highly expressed genes evolve more slowly than viable genes [52]. Overall, this analysis indicates that older mouse genes are more likely to be indispensible for fundamental cellular processes. Essential genes might be undergoing positive selection to retain their functionality, giving a lower mutation rate.  Numerous features derived from protein sequences differ between essential and viable gene products. We find that proteins encoded by essential genes tend to be longer in length and have greater molecular weight than proteins encoded by viable genes. This result is consistent with a prior study which stated that functionally essential proteins are more evolutionarily conserved and conserved proteins are, in general, longer in length [53]. Longer proteins contain more possible domains to mediate diverse cellular functionalities [45] and multiple protein-protein interactions, and our analysis with GO terms and protein domains supports this observation.
Essential proteins were found to have Ala, Asp, Glu, Lys, Gln and Ser residues in greater proportions. In contrast, viable proteins have higher proportions of Cys, Phe, Ile, Leu, Val and Trp. Cys can also form disulphide bonds, which are more common in extracellular proteins. A high Leu content is known to correlate negatively with the likelihood of being essential [20]. The enrichment of Lys in essential proteins agrees with our findings that they are likely to be more acetylated, as proteins are often acetylated on lysine residues [54]. The enrichment of acetylated proteins in essential datasets implies that they are involved in regulating proteinprotein interactions, gene expression and metabolic processes [55].
Essential proteins are likely to have more polar, charged, basic and acidic amino acids, whereas viable proteins have more aliphatic, aromatic and non-polar residues, because viable proteins are more likely to be membrane proteins. Viable proteins also tend to have Src Homology 2 (SH2), Src Homology 3 (SH3) and ion transport domains. This further establishes the propensity of viable proteins being membrane-bound, as these evolutionary conserved protein domains are common constituents of membrane proteins. Signal peptide motifs are found to be more frequent in viable proteins, as they are more likely to be secreted.
Proteins encoded by essential genes are more likely to function as ligases or transferases, consistent with results from a recent study [14]. Our analysis of GO biological process annotations reveals that proteins encoded by essential genes are often involved in regulating DNA replication, DNA repair and transferase activity, consistent with these proteins functioning as ligases and transferases. The enrichment for hydrolases in viable datasets is logical, as hydrolases are less critical to cellular function. Ligases also perform more complex chemistry than hydrolases, which may indicate their essentiality. Essential proteins are more likely to be phosphorylated, as phosphoproteins have critical roles in almost all cellular processes, including cell differentiation, gene transcription and cell division [56], confirmed by biological process GO terms. A greater number of N-glycosylated proteins in the viable datasets indicates the propensity of viable proteins to be membranebound or extracellular [57] and have a longer in vivo lifetime [58]. N-glycosylation can also aid folding of complex proteins [59]. Essential proteins are more likely to be engaged in control of gene transcription in the nucleus, as shown by GO and UniProt annotations and also by the presence of the zinc finger, C4 type (zf-C4) domain. Most occurrences of the zf-C4 domain are found within the DNA-binding regions of many nuclear receptors that function as transcription factors, which are enriched in the essential gene dataset.
Transmembrane proteins are common in the viable datasets, shown by various data, such as GO annotations, amino acid preferences and predicted number of transmembrane helices, consistent with their roles in cell communication, transport and signal transduction. These functions therefore appear to be less critical during mouse development.
Subcellular localisation has been recognised as a significant attribute for identifying essential genes in bacteria and yeast [15,16,18,35]. We found that proteins encoded by mouse essential genes are more likely to be intracellular. The majority of these proteins are located in the nucleus. This is expected, because almost one third of eukaryotic nuclear proteins are encoded by essential genes and are responsible for carrying out vital cellular processes like DNA replication, DNA repair and transcription [60,61]. Our analysis of GO biological process annotations further confirms this result. Proteins encoded by viable genes tend to be secreted (extracellular), as shown by the enrichment of signal peptide cleavage sites, fibronectin type III (fn3) domain and signal transducer activity. The fn3 protein domain is an evolutionarily conserved domain that is generally found in animal proteins, especially in extracellular proteins. Its main function is to mediate cell-cell signaling or interactions. These results agree with a recently reported study [14], which also showed the tendency of essential genes to encode nuclear proteins and viable genes to encode extracellular or membrane-bound proteins.
Unsurprisingly, GO data analysis shows that essential genes are more likely to be involved in developmental processes. These include the development of embryo, tissue, heart, nervous system, brain, lung, respiratory tube and blood vessel. Essential genes are enriched in Tbox domains, which are vital for heart development [62]. We observed a significant enrichment of essential genes in cell morphogenesis, cell division, cell proliferation, DNA replication, cell differentiation, DNA repair and transcription, all crucial processes for growth. The presence of homeobox domains further confirms their vital role in morphogenesis. In contrast, viable genes were associated with inflammatory response, apoptosis, behaviour and immune response, which are processes unlikely to be required in utero. Unlike viable genes, essential genes tended to participate in activities such as protein binding, DNA binding, transcription factor binding, transcription, and ATP binding. Viable genes were linked to transport, ion channels, signal transduction (with SH2 protein domains), enzyme binding, receptor binding, and lipid binding, consistent with their location in membranes and involvement with communication with other cells. A recently published study by Dickinson et al [14], analysed GO terms for 410 essential and 1143 viable genes. Of these, only 89 essential and 297 viable genes are also present in our datasets. Dickinson et al. also found the tendency of essential genes to participate in cell differentiation, cell proliferation, transcription and DNA binding. Confirmation of a consistent set of GO annotations with a different gene dataset adds strong support to our conclusions regarding the functions of essential and viable genes. Viable genes are more likely to be housekeeping or tissue enriched.
The correlation between PPI networks and gene essentiality has already been established in bacteria [15,16,39], yeast [15,35,39,[63][64][65], fly [66] and human [21,67]. Proteins that are highly connected in the PPI (hubs) tend to be essential and evolve slowly [68,69], and their absence disrupts cell viability [70]. Shorter ASP length and high values of closeness centrality and clustering coefficient show that essential proteins can quickly transfer information to other reachable protein nodes in the PPI network.
We conclude that mammalian essential genes are significantly different from non-essential genes based upon a number of features. Our manually curated datasets allowed for a large number of statistically significantly differing features to be identified. The interdependency of various features implies that multiple aspects of biology unite to determine whether a gene is essential or non-essential in mammals. The features we have identified as strongly associated with mouse essential genes can be compared to characteristics of essential genes in other organisms to gain insights into the evolution of essential functions. The wide variety of features we have identified associated with essential and non-essential genes will allow for improvements in predicting whether a gene is essential. Due to the genome similarities between mice and humans, future analyses may facilitate the identification of human genetic disease candidates and potential therapeutic targets.

Datasets
To construct the datasets for the current research, the phenotype information of knockout mice was collected from the MGI database [22] (http://www.informatics.jax.org/phenotypes. shtml, accessed on 1 November 2013). We included only null alleles of mouse genes that have known phenotypes resulting from single gene targeted (knockout) deletions, because mutations generated by other methods may not have created complete loss of function alleles, and thus a viable hypomorphic allele could be associated with a gene that produces a lethal null allele. Genes were included in the essential dataset if they produced lethality in either the homozygous or heterozygous state on any strain background, and were not separated further within the dataset. The phenotype of a mouse gene was marked as essential or lethal if it is associated with any essentiality annotation in the MGI (including prenatal, perinatal and postnatal annotations). The term 'prenatal essentiality' is a valid Mammalian Phenotype Ontology term which is defined in MGI as death of the mice anytime between fertilization and birth, whereas, 'perinatal essentiality' is defined as death any time between embryonic day E18.5 and postnatal day 1. We used 18 phenotypic annotations to classify a single-gene knockout phenotype as nonessential or viable (S10 Data). Since the majority of these terms refer to processes or tissues present only after birth, homozygous knockouts of these genes are evidence of a viable phenotype. We manually checked the literature for phenotypes of knockouts of genes for viability that were linked to the "adipose tissue", "abnormal skin morphology" and "abnormal skin physiology" terms, as these could be applied to embryos.
Our knockout datasets contained some ambiguous entries that have been annotated as both essential and viable in the MGI database. We manually checked phenotypes of these overlapped entries against the published literature and labelled them either as essential or viable. Each MGI gene symbol and identifier was further mapped to its corresponding Ensembl gene ID (http://www.ensembl.org) [71], UniGene expression clusters ID (http://www.ncbi.nlm.nih. gov/unigene) [24] and UniProt protein ID (http://www.uniprot.org/uniprot/) [28]. For some instances there were multiple UniProt protein IDs that correspond to one gene. For some of these cases, only one protein had the longest length and we included that in our dataset. For others, two or more protein IDs were found to have longest length. In these cases, to avoid bias due to annotation quality we included the longest length protein ID in our dataset that was marked as 'reviewed' in the UniProt annotations. Mouse protein sequences in FASTA format were downloaded from UniProt.
Non-redundant datasets. Redundancy was removed from our original essential and viable datasets by submitting the essential and viable protein sequences in FASTA format to Leaf [23] (http://leaf-protein-culling.appspot.com/). This gave four sets of non-redundant essential and viable proteins with protein pairs showing the maximum sequence similarities of 20%, 40%, 60% and 80% respectively. Protein sequences with <20% identity are structurally very different implying functional differences [72,73] so we therefore did not generate non-redundant datasets by removing proteins with <20% sequence identities. juvenile and adult. Since the total number of ESTs for a particular gene varies greatly between different developmental stages, we corrected the raw data to get gene expression in the form of transcripts per million (TPM). Eq 1 was used to estimate a TPM for the i th gene at j th developmental stage.
Evolutionary age. Evolutionary ages of mouse protein coding genes were determined by analysing the Ensembl (release 75) gene trees [75]. These gene trees represent the evolutionary processes by which genes diverged from their common ancestors. Ensembl runs a orthology and paralogy gene prediction pipeline that uses the TreeBeST method from the TreeFam methodology [78] to generate rooted phylogenetic trees. This pipeline merges tree topologies with the corresponding species trees inferred from the NCBI taxonomy and generates Ensembl genes trees with the tree internal nodes being annotated for duplication or speciation events.
Gene evolutionary ages were extracted from these Ensembl genes trees. We assigned two evolutionary ages to a mouse gene of our datasets: the age of the MRD event and the age of the evolutionarily most distantly related species, i.e., the age of the DCA that has an identified homolog to that gene.
Protein sequence properties. We retrieved the length of our essential and viable protein sequences by querying the UniProtKB database with their UniProt IDs. A script in Python was developed to compute the percentage frequencies of each of the 20 amino acid residues within protein sequences.
Pepstats (http://emboss.bioinformatics.nl/cgibin/emboss/pepstats) is a EMBOSS suite program [27] which outputs a report comprising statistics of a number of properties about a FASTA formatted protein sequence. These attributes include: molecular weight, number of residues, charge, isoelectric point, and amino acid composition. This program groups amino acids into nine categories: Tiny  and Acidic (B, D, E, Z). We used Pepstats to evaluate these sequence properties for our essential and viable protein sequences. The program was run with the default parameters setting. A Python script was written to extract features values from the output file generated by Pepstats.
Enzyme class. Primary EC numbers of mouse proteins were obtained from the definition lines (DE) of UniProtKB annotations by submitting UniProt IDs.
Post-translational modifications. Three post-translational modification (PTM) keywords 'Glycoprotein', 'Phosphoprotein' and 'Acetylation' were collected from the UniProtKB database for each protein of our datasets. The UniProt annotation 'Glycoprotein' is used for N-glycosylation sites.
We also collected information about the keyword 'Transcription'. It is a keyword in the biological process category representing proteins involved in regulating the process of transcription.
Signal peptides. Protein signal peptides were predicted using the SignalP program v4.1 (http://www.cbs.dtu.dk/services/SignaP/) [33]. This program uses artificial neural network (ANN) and hidden Markov model (HMM) algorithms to predict the amino acid composition and the cleavage site position of the signal peptide. A script in Python was written to extract the HMM probabilities generated by SignalP which is considered as the measure for signal peptide prediction.
Transmembrane domains. We extracted the total number of transmembrane domains in each mouse protein by querying the UniProtKB database. Transmembrane helices are annotated in the UniProt feature table line (FT) as TRANSMEM. UniProt also outputs the information about the transmembrane domain locations in a protein sequence. Subcellular location. Protein subcellular localizations were predicted from sequence data using the WoLF PSORT program (http://wolfpsort.org/) [37], chosen as it can make prediction on any protein sequence. WoLF PSORT predicts subcellular locations on the basis of known sorting signals, functional motifs and sequence features, such as amino acid composition. It outputs a report covering predicted locations with different confidence levels. We found prediction scores for six subcellular locations: nucleus, cytosol, plasma membrane, mitochondria, Golgi apparatus, peroxisome, and extracellular. We assigned a score of zero to a subcellular location if no prediction is made. We further collected information about all these six subcellular locations from the UniProtKB database. This feature is annotated as SUBCEL-LULAR LOCATION in the UniProt data file and is found in the comment lines (CC). The value of a subcellular location was set to 1 if found; otherwise, it was set to 0.
Gene ontology terms GO terms were obtained by using the 'Functional Annotation' tool of the web based application DAVID v6.8 (https://david.ncifcrf.gov/home.jsp) [36]. It integrates gene functional annotations with intuitive graphical displays to facilitate biological interpretations of any list of genes encoded by human, rat, mouse, or fly genomes. This program systematically associates a query gene list to their corresponding GO terms and highlights only the most pertinent terms among all along with their statistics. We extracted all possible GO terms for which the statistical test supported in DAVID has a p-value 0.05.

Protein-protein interactions
Mouse PPI data was downloaded from the Interologous Interaction Database (I2D) v2.3 [40] which is an integrated repository of known, experimental and predicted PPIs for human, mouse, rat, fly, yeast and worm genomes. To obtain high quality PPI data, we analysed all known and predicted mouse PPIs. The data obtained from the I2D database were imported into Cytoscape (v3.1.1) [76] to visualize and analyze PPI network as a graph. In this case, we removed all self-loops and duplicate edges. The 'network analyser' plugin of Cytoscape was further used to determine network properties including degree, the length of average shortest path, betweenness centrality, clustering coefficient, and closeness centrality. We further determined four other network properties including BN, EPC, MNC and DMNC by using a webbased service called Hub Object Analyzer (Hubba) (http://hub.iis.sinica.edu.tw/Hubba/) [41]. This system deciphers and visualizes hubs from the user-provided PPI networks. Query proteins are ranked in Hubba based on their topological features. Hubba also generates a subgraph for the top n ranked (n 100) hub along with their identifier.
PPI networks are usually characterized as undirected graphs. As an example, let G = (V, E) be an undirected graph representing a PPI network. In the graph G, nodes V represent proteins and edges E = {(a,b) | a,b 2 V} correspond to observed interactions between protein a and protein b. Graph topological feature are defined as: Degree. The most elementary property of a protein a is its degree or connectivity, which is the number of interactions a has to the other proteins in the network.
Average shortest path length (ASP). The shortest path measures the path with the minimum number of edges between proteins a and b. The ASP length therefore refers to the average over all shortest path length between all protein pairs. Betweenness centrality (BC). The betweenness centrality (BC) of a protein node a corresponds to the ratio of shortest paths passing through a [79,80] and is computed as follows: Here, σ bc denotes the number of shortest paths between proteins b and c; and σ bc (a) denotes the number of shortest paths between b and c to that go through protein node a.
Clustering coefficient (CCo). The clustering coefficient (CCo) of protein a (Eq 3) measures the ratio of the number of connections between all nodes within the neighborhood of a to the maximum number connections that could possibly present between them [81] Here, e bc denotes the number of connections between all neighbors b and c of a, and k a denotes the degree of a.
Closeness centrality (CC). The closeness centrality (CC) of the protein a corresponds to the reciprocal of the sum of average shortest path length between a and all the other nodes within the network [82] (Eq 4). It measures how close a protein node is to all the other nodes in the PPI network.
Here, d(a, b) is the length of the average shortest path between proteins a and b. BC, CCo and CC of each protein node are represented by a value between 0 and 1 where an isolated protein node has a value of 0 for these properties.
Bottleneck (BN). Let, T r be the shortest path tree derived from G considering protein node r 2 V as the root node. Protein b 2 V is a bottleneck node if at least n/4 nodes have their shortest path to r through b in T r . The BN score of the protein node b is defined to be the number of nodes r for which b is a bottleneck node in T r [83].
Edge percolation component (EPC). G' is a graph which is constructed n times from G by randomly removing a subset of edges. It is possible that proteins a and b are connected in G but not in G'. The EPC score [84] of the protein node a is computed using the following equation.
EPCðaÞ ¼ Density of maximum neighborhood component (DMNC). The DMNC for the protein a is calculated using the following equation: Here, E M denotes the number of edges and N denotes the number of protein nodes of DMNC(a); e is a constant which is equal to 1.7.

Statistical analysis
Statistical tests were carried out throughout using the statistics software package SPSS v20. First, the normality of different features was assessed, as many significant statistical tests including parametric tests depend upon normal data. All sequence features were tested for normality using a one sample Kolmogorov-Smirnov Test (K-S test). If a sequence property showed a normal distribution, a two-sample t-test with unequal variance analysis was used. The t-test was applied to test the null hypothesis that two samples of independent observations come from identical normal distributions with equal means. Statistical significance was determined at the 0.05 level. The statistical significance of each property was determined using the two-tailed nonparametric Mann-Whitney U test when the samples did not show normal distribution. The Chi-squared (χ 2 ) test was also carried out to check whether the frequencies of a particular feature in essential and viable gene differ from each other. We then applied the Bonferroni correction to calculate corrected p-values.