BACKGROUND
gene expression variation has been studied on three different levels: single cells across a common environment  <cit> , within one species across a variety of different environments  <cit> , and across different species/strains, which is often referred to as evolutionary variation  <cit> .
in this paper, we study genetic factors affecting gene expression variation within one species across many different environmental conditions.
broadly, the genetic factors affecting gene expression primarily include the binding of regulatory proteins to cis-elements in the upstream of the gene, as well as physical and genetic interactions with other genes.
with the availability of many gene expression profiles, protein interaction networks, and gene regulatory networks, it is now possible to study how gene expression variation is associated with both network features and genomic factors.
in the case of protein interaction networks, interaction degree, i.e., the number of interacting partners of a given protein, is one of many factors.
the presence or absence of tata box and the number of transcription factors  regulating a gene provide examples of genomic factors influencing gene expression variation.
many studies have focused on individual factors affecting gene expression variation.
for instance, newman et al.
<cit>  developed an experimental technique to study protein expression noise in single cells and showed that chromosomal distance to other genes and mrna-half life are associated with expression noise.
however, they did not find a relationship between protein expression noise and protein-protein interactions.
recently, using a more complete interaction dataset, batada et al.
<cit>  found that protein expression variation is negatively correlated with interaction degree when protein abundance was controlled using the data in newman et al.
<cit> .
this relationship continues to hold within the viable genes.
several groups investigated mrna expression variation within species.
for example, nelson et al.
<cit>  studied the relationship between the number of tissues or body parts , where the gene is expressed, and gene spacing in c. elegans and d. melanogaster.
they found that gene expression variation increases in relation to the intergenic distances between genes.
walther et al.
<cit>  found a positive correlation between the frequency of a gene's differential expression and the number of cis-regulatory elements of that gene in a. thaliana.
furthermore, several groups have studied gene expression variation, also known as evolutionary variation, across different species/strains.
using gene expression data from several yeast species  <cit> , as well as from different strains derived from mutation-accumulation experiments  <cit> , it was found, for instance, that the interspecies/interstrain variation of gene expression is significantly correlated with the presence/absence of the tata box in the promoter region.
lemos et al.
<cit>  studied the effect of protein-protein interactions and protein length on evolutionary variation .
they found that evolutionary variation is negatively correlated with protein-protein interactions in saccharomyces cerevisiae or drosophila melanogaster  <cit>  and negatively correlated with protein length in drosophila melanogaster  <cit> .
these studies highlighted the importance of protein interactions and gene regulatory regions on gene expression variation.
only a few studies, however, have integrated such different data sources in a way that collectively identifies and interprets the key factors affecting gene expression variation.
therefore, we conducted studies of proteomic and genomic factors marginally and collectively influencing gene expression variation across different perturbation conditions within one species: yeast.
protein interactions play an important role in gene expression variation.
protein-protein interactions are key biological events in a living cell, and proteins in a cell interact with each other to perform certain functions.
high throughput technologies, including yeast two-hybrid systems and mass spectrometry, have generated a large amount of protein interactions in yeast.
computational methods have also been developed to study the reliability of the observed interactions  <cit>  and to build reliable protein interaction networks.
these efforts have resulted in the development of several protein interaction databases, albeit with differing degrees of reliability, including mips  <cit> , dip  <cit>  and biogrid  <cit> .
from an evolutionary point of view, the expression profiles of neighboring genes of a target gene in a protein interaction network may put some constraints on the target gene's expression.
thus, in a protein interaction network, the interacting partners of a specific protein can affect the corresponding gene's expression.
therefore, protein physical interaction degree, i.e., the number of interacting partners of a given protein, can significantly affect gene expression variation.
in the present study, we show that gene expression variation decreases with protein interaction degree and that protein interaction degree accounts for 1–2% of the expression variation in model organism yeast, a result consistent with previous studies  <cit> .
another key factor affecting gene expression variation is gene essentiality.
genes can be classified into essential and non-essential genes based on the fitness phenotype of the yeast cell when the gene is deleted under normal growth conditions  <cit> .
essential genes are those that, when deleted, will render the yeast cell non-viable.
non-essential genes can be further classified into no-phenotype and toxicity-modulating genes based on the fitness phenotype of yeast cell when the gene is deleted under the conditions of four dna-damaging treatments  <cit> .
specifically, we define a gene's toxicity modulation degree as the number of dna-damaging treatments significantly affecting the deletion strain's fitness ,  <dig>   <dig>   <dig>  and 4).
the higher the toxicity modulation degree, the more important the gene is in relation to cell survival.
therefore, toxicity degree gives a quantitative measurement of a gene's importance to yeast cell survival.
we measure a gene's functional importance in relation to cell survival by the essentiality of the essential genes and the toxicity modulation degree of non-essential genes.
since the expression of genes important for cell survival are generally stable under many different stimuli and cannot fluctuate extensively, we hypothesize and show that expression variation of essential genes is lower than that of non-essential genes and decreases with toxicity degree within non-essential genes.
the number of cis-elements has been shown to be positively associated with gene expression variation  <cit> .
the number of cis-elements is usually approximated using computational approaches and many contain false positive and negative predictions.
theoretically, a given gene's expression pattern can become increasingly complex with the increasing number of transcription factors that regulate this gene, either directly or indirectly.
in this study, we hypothesize that the number of tfs is a significant predictor of expression variation and show that the number of tfs regulating a gene  accounts for 8–14% of its expression variation, much higher than that can be explained by the number of cis-elements .
this implies the importance of indirect trans-effect on expression variation.
the tata box is a conserved element in the eukaryotic promoter region and is usually bound by tata-binding proteins.
the presence of tata box has been shown to be one of the most important factors contributing to gene expression variation  <cit> .
further analysis of the individual genomic and proteomic factors affecting gene expression indicates that there might be two distinct mechanisms that specifically influence gene expression variation of tata-containing and non-tata-containing genes.
most importantly, we show that significant negative correlation between expression variation and toxicity degree is only present for tata-containing genes and that toxicity degree accounts for  <dig> – <dig> % of the expression variation.
in contrast, the relationship between expression variation and toxicity degree is absent for non-tata-containing genes.
the fact that tata-containing genes are enriched in stress-related genes  <cit>  may explain this difference.
although the number of tfs is significantly positively correlated with expression variation for both tata- and non-tata-containing genes, the association strength is higher for non-tata containing genes than for tata-containing genes.
these results imply that the mechanism influencing tata-containing gene expression variation is much more complicated than that in non-tata-containing genes.
for example, tata-containing genes were found more likely to be epigenetic regulated  <cit> .
thus, this study gives a more complete analysis of factors and their interaction affecting gene expression variation than may be found in previous studies.
RESULTS
we present our results based on the mips protein physical interaction data  <cit>  and the yeast gene expression profiles under  <dig> ca and na exposure conditions  <cit> .
the results based on three interaction datasets  and four other gene expression datasets   <cit> , environmental stress  <cit> , oxidative stress  <cit> , and a combined gene expression dataset over more than  <dig>  conditions  <cit> ) are given in the additional files  <dig> and  <dig>  we study the expression data individually in order to minimize the variation among different laboratories.
by doing so, we can also confirm whether the results based on different gene expression data are consistent.
consistency of results using a variety of different datasets adds confidence to the conclusions.
in this manuscript, we use genes and proteins interchangeably.
we declare statistical significance if a p-value is less than  <dig>  without adjusting for multiple comparisons.
in this study, we conducted an exploratory study of factors affecting gene expression variation.
as in many epidemiological studies, we did not adjust p-values for multiple comparisons.
therefore, some of our findings need to be further tested in other datasets.
gene expression variation versus protein interaction degree
we measured gene expression variation by the logarithm of the variance of the expression levels of each gene across the  <dig> ca and na exposure conditions  <cit> .
we then studied the relationship between the expression variation and protein interaction degree using the lowess function in r  <cit>  to fit the data.
on these bases, it was obvious that gene expression variation decreases with the degree of protein physical interaction .
the decreasing trend is especially significant when the interaction degree is relatively low .
in contrast, when the interaction degree is greater than  <dig>  the decreasing trend is not as obvious.
biologically, this may be explained by the fact that gene expression variation stabilizes when the interaction degree is above a given threshold.
another potential explanation is that interactions between proteins of high degrees are simply less reliable  <cit> .
the large number of less reliable interactions for proteins with high degrees can skew the true relationship between gene expression variation and interaction degrees.
since the true underlying mechanism of this phenomenon is not clear, we limited our further analysis to proteins with an interaction degree of no more than  <dig>  protein physical interaction degree has been found to be negatively associated with gene expression variation for single cells  <cit>  and evolutionary expression variation  <cit>  between different strains/species.
our result on the relationship between gene expression variation and interacting degree within s. cerevisiae is consistent with their findings.
accordingly, we then used linear regression to fit the expression variation for proteins with a maximal physical interaction degree of 20:
v = α + βd
where v is the gene expression variation and d is the interaction degree.
α and β are parameters.
the fitted line and the corresponding bar-plot for the expression variation are shown in figure 1b.
the gene expression variation is significantly negatively correlated with the protein interaction degree  .
the negative correlation between expression variation and interaction degree implies that protein with high interaction degrees do not tolerate extensive expression variation and such protein need more precise control on gene expression for an organism to function normally.
gene expression variation versus essentiality, toxicity modulation, and interaction degrees
as noted above, we divided genes into two classes: essential and non-essential genes.
essential genes are less likely to be perturbed than non-essential genes, as significant perturbations of essential genes will, for example, render the yeast cell non-viable.
we further classified the non-essential genes into five groups  according to the cell's fitness phenotype changes under four dna-damaging agents , the bulky alkylating agent 4-nitroquinoline-n-oxide , the oxidizing agent tert-butyl hydroperoxide , and 254-nm uv radiation) when the non-essential genes are knocked out  <cit> .
since toxicity degree reflects the functional importance of genes in relation to cell survival under several dna-damaging perturbations, we expected that gene expression variation would decrease as the toxicity degree increases.
since the number of genes with toxicity degree  <dig> was small , we combined them with the group having toxicity degree  <dig>  we referred to the essential genes as the group with toxicity degree  <dig>  figure 2a shows the bar-plot and the linear regression fit of the gene expression variation with respect to toxicity degree.
indeed, a significant negative association between gene expression variation and toxicity degree was observed .
in the study of the relationship between gene essentiality and evolutionary expression variation conducted by tirosh and barkai  <cit>  and choi et al.
<cit> , they found that essential genes tend to have lower variation than non-essential genes.
this is consistent with our result which demonstrates that the variation in essential genes  is lower than that of non-essential genes  .
however, our study differs from these two studies in that we further classified the genes according to their toxicity degree and found negative association between gene expression and toxicity degree.
we observed a positive correlation between interaction degree and toxicity degree  and, therefore, asked whether the observed negative correlation between gene expression variation and interaction degree is, conversely, caused by the positive correlation between interaction degree and toxicity degree.
we consequently studied the relationship between gene expression variation and protein interaction degree within gene groups stratified according to their toxicity degrees .
using ca and na exposure gene expression data  <cit> , we found a significant decreasing trend of gene expression variation with respect to interaction degree in all the strata except for the one with toxicity degree  <dig>  the corresponding  are , , , , and  for toxicity degrees  <dig>   <dig>   <dig>   <dig> and  <dig>  respectively.
the fraction of expression variation explained by the protein interaction degree seems to increase as the toxicity degree increases.
in our analyses, both toxicity and protein interaction degrees are negatively associated with gene expression variation.
hence, the more important a gene is to the survival of the yeast cell, the less variation there is in its expression levels across many different conditions.
similarly, the higher the interaction degree of a gene, the more stability is observed in its expression levels.
biologically, a gene is important to the cell's survival since it participates in many important biological processes.
any perturbation of this gene's expression will likely cause deleterious effect to the corresponding biological process and thus renders the cell non-viable.
an evolutionary consequence of this hypothesis is that genes important to cell survival appear to have robust expression levels.
expression variation versus gene regulatory regions: tata box, number of tfs, and toxicity degree
previous studies have established the relationship between gene expression variation and the regulatory regions, including the presence/absence of tata box  <cit> , the length of intergenic regions  <cit> , and the number of cis-regulatory elements  <cit> .
we therefore asked if the observed relationship between gene expression variation and toxicity degree are the same for tata-containing genes and non-tata-containing genes.
to accomplish this goal, we first stratified the yeast genes based on the presence/absence of tata boxes and reanalyzed the relationship between gene expression variation and toxicity degree.
consistent with previous findings  <cit> , it is clear that the gene expression variation of tata-containing genes is much higher than that of non-tata-containing genes for each fixed toxicity degree .
significant negative association between gene expression variation and toxicity degree was observed for the tata-containing group  .
however, the association between expression variation and toxicity degree for the non-tata-containing group is only marginally significant .
consistent with the other two gene expression datasets, chemostat  <cit>  and environmental stress  <cit> , as well as the combined gene expression data from  <cit>  , such a highly significant negative relationship between gene expression variation and toxicity degree could be observed in tata-containing genes, but not for the non-tata-containing genes.
thus, the effect of toxicity degree on gene expression variation is different for tata-containing genes versus non-tata-containing genes.
the relative small r <dig> value between gene expression variation and toxicity degree alone may be explained by the absence of association between them within the non-tata-containing genes.
previous studies identified the number of cis-element motifs contributing to gene expression variation in a. thaliana  <cit> , but landry et al.
<cit>  found that the number of cis-elements only marginally affects gene expression variation in yeast.
here we studied these two factors  in relation to gene expression variation.
we first analyzed the relationship between gene expression variation and number of cis-elements  <cit>  using the linear model.
similar to macisaac et al.
<cit> , we defined cis-elements according to binding probability  and classified them according to their conservation in two other yeast species.
the results are presented in table  <dig> for the different gene expression data.
for the ca and na exposure data  <cit> , the highest r <dig> is  <dig> %  when cis-elements that are conserved in at least one species are used.
we then studied the relationship between gene expression variation and the number of tfs that influence target gene expression based on the gene regulatory network developed in hu et al.
<cit> .
a highly significant positive correlation between gene expression variation and the number of tfs was observed  .
the fraction of variation explained by the number of tfs  is much higher than that by the number of cis-elements indicating that the number of tfs is a better predictor of gene expression variation than the number of cis-elements.
this result implies the importance of trans-effect for gene expression variation.
cis-elements are identified with three different criteria according to their conservation in two other species.
β is the linear coefficient in the linear model, the p-value is related to the null hypothesis that β ≠  <dig> versus β =  <dig>  and the r <dig> is the fraction of variation explained by the number of cis-elements.
overall analysis of factors affecting gene expression variation
as enumerated above, we have identified several factors influencing gene expression variation in s. cerevisiae.
in addition to the presence/absence of tata box identified in previous studies  <cit> , we found that gene expression variation decreases as both the protein interaction and toxicity modulation degrees increase.
these findings are consistent with other studies for expression variation of single cells  <cit>  and evolutionary expression variation across different strains/species  <cit> .
we also found that gene expression variation increases as the number of tfs or the number of cis-elements increases and that the number of tfs regulating a gene is a much better predictor of expression variation than the number of cis-elements.
we therefore studied the contribution of each factor and their interactions on gene expression variation by taking the other factors into consideration using the akaike information criterion   <cit> .
we retained the model with the smallest aic.
table  <dig> gives the results of factors included in the final linear model using the mips interaction data  <cit>  and the four expression profiles by stepwise selection with aic.
we then studied the effect of the selected factors on expression variation using linear regression, and the corresponding p-values and r <dig> values are given in table  <dig>  the results showed that protein interaction degree only explained less than 1% of variation when adjusted for the three other main factors and two interaction terms.
except for the oxidative stress gene expression dataset  <cit> , the consistently selected model contained four main factors including protein interaction degree, toxicity degree, the number of tfs, and tata box, and two interaction terms, i.e., interaction between tata and toxicity degree and tata and number of tfs.
the two interaction terms are found to be statistically significant across three expression datasets: ca and na exposure  <cit> , chemostat  <cit> , and environmental stress  <cit>   and explained  <dig> % –  <dig> % of the variation.
the four main factors include protein interaction degree , toxicity degree , number of tfs , and the presence of tata box .
the protein interaction data used in this analysis is based on the mips dataset.
the column marked with "√" indicates inclusion in the final linear model.
the multiple linear regression is based on the final linear model, respectively.
the p-value is related to the null hypothesis that β ≠  <dig> versus β =  <dig>  r <dig> is the variation explained by the model and each independent variable, respectively.
because we found that tata box interact with toxicity degree and the number of tfs influencing gene expression variation, we reanalyzed the contributions of toxicity degree and number of tfs on gene expression variation stratified by the presence/absence of tata box using the linear model.
the protein physical interaction degree was not included in this analysis because the interaction between tata box and physical interaction degree was not detected.
table  <dig> shows the different results for tata and non-tata gene sets.
first, the number of tfs is the most significant factor for gene expression variation in both tata- and non-tata-containing genes with p-values less than 5e- <dig> for all the datasets.
however, the association strength between the number of tfs and gene expression variation measured by the coefficient β for the non-tata-containing genes is about  <dig> -fold higher than the corresponding values for the tata-containing genes.
accordingly, the p-values for the non-tata-containing genes are about three orders of magnitude smaller than the corresponding p-values for the tata-containing gene set.
we noted that tata-containing genes tend to have higher number of tfs than non-tata-containing genes , but this cannot explain our observation that the association strength between expression variation and the number of tfs within the tata-containing genes is lower than that within the non-tata-containing genes.
one possible explanation is that the presence of tata-box weakens the effect of the number of tfs on gene expression during evolution.
second, the toxicity degree is only a significant contributor for gene expression variation within the tata-containing genes in the ca-na exposure  <cit>  and chemostat  <cit>  gene expression dataset, while it is not a significant contributor for gene expression variation in the non-tata-containing gene set.
within the environmental stress gene expression dataset  <cit> , the effect of toxicity degree on gene expression variation was not statistically significant.
however, we did observe a decreasing trend of gene expression variation with respect to toxicity degree within the tata-containing genes .
the linear model that includes the toxicity degree  and the number of tfs  is built for each of the four gene expression datasets, respectively.
r <dig> is the variation explained by each independent factor and the model, respectively.
β is the linear coefficient in the linear model, and the p-value is related to the null hypothesis that β ≠  <dig> versus β =  <dig>
the toxicity degree of a gene measures the tolerance of the yeast cell to different external stress conditions when the gene is knocked out.
therefore, it might be expected that the different relationship between gene expression variation and toxicity degree for tata-containing genes and non-tata-containing genes is due to the enrichment of stress-related genes in tata-containing genes, as found in  <cit> .
to test this hypothesis, we used a set of genes related to environmental stress response   <cit> .
if the hypothesis is true, we would expect a higher association between expression variation and toxicity degree within the esr genes than that within the non-esr genes.
however, our data shows that the β values are similar for the two groups of genes .
thus we cannot explain the interaction between tata and toxicity degree by the enrichment of esr genes in tata-containing genes.
the biological mechanisms underlying the observed interaction are not clear and need to be further studied.
the linear model is built for each of the four gene expression datasets, respectively.
β is the linear coefficient in the linear model and the p-value is related to the null hypothesis that β ≠  <dig> versus β =  <dig>
we also did the same analysis for the average gene expression variation across the four expression datasets  and the combined gene expression data of landry et al.
<cit> , and the results are presented as additional file  <dig>  the same conclusions can be obtained indicating the robustness of our results.
previous studies showed that tata- and non-tata-containing genes might recruit different coactivator complexes for gene expression  <cit> .
tata-containing genes were also found to be subject to greater nucleosomal regulation than non-tata-containing genes  <cit> .
basehoar et al.
<cit>  suggested that two distinct regulatory mechanisms may be present at tata- and tata-less promoters.
the results in table  <dig> support their findings.
the results based on the oxidative stress gene expression dataset  <cit>  are not consistent with the results based on the other three gene expression datasets.
this observation may be due to the relatively small gene expression variation in this data.
for example, the range of the variance of the expression levels within the oxidative stress dataset, , is much smaller than the corresponding ranges, , , and , for the ca and na exposure  <cit> , chemostat  <cit> , and environmental stress conditions  <cit> , respectively.
we also studied the contributing factors for gene expression variation using the dip  <cit>  and biogrid interactions  <cit> , and the results are given in additional file  <dig> and file  <dig>  similar conclusions as those based on the mips interaction data  <cit>  were obtained.
the consistency of the results using different combinations of protein interaction data sets and gene expression profiles showed the robustness of our conclusions.
however, the fraction of gene expression variation explained by all factors is less than 25%.
one possible explanation is that the measurement of gene expression changes and other factors, including the toxicity degree and interaction degree, are still very noisy.
we expect that the true r <dig> would be higher than that observed in this study.
CONCLUSIONS
we implemented a system-wide analysis of proteomic and genomic factors affecting gene expression variation.
among four different factors , tata-box and the number of tfs are the most important factors influencing gene expression variation.
the influence of tata-box on evolutionary gene expression variation has been extensively studied both computationally and experimentally  <cit> , and our results are consistent with their findings.
although it is intuitive that the number of tfs regulating a gene should have a significant effect on the gene's expression variation, the magnitude of its influence has not been studied in large scale expression datasets to the best of our knowledge.
our findings demonstrated that the gene regulation is a main factor affecting gene expression variation.
protein interaction degree and toxicity degree do not account for as much variation when compared to the influence of the number of tfs and the tata-box.
in our overall analysis, we also found the interactions between tata-box and toxicity degree as well as the number of tfs influence expression variation.
the further study stratified by tata-box indicated that tata-containing genes and non-tata containing genes behave differently in relation to the toxicity degree and the number of tfs.
the effect of the number of tfs on expression variation within the tata-containing genes is lower than that for the non-tata-containing genes.
on the other hand, toxicity degree is associated with expression variation within the tata-containing genes only.
these findings suggest that the regulatory mechanism might be more complicated for tata-containing genes than non-tata containing genes.

 gene expression variation has been studied on three different levels: single cells across a common environment, within one species across a variety of different environments, and across different species/strains, which is often referred to as evolutionary variation. in this paper, we study genetic factors affecting gene expression variation within one species across many different environmental conditions. broadly, the genetic factors affecting gene expression primarily include the binding of regulatory proteins to cis-elements in the upstream of the gene, as well as physical and genetic interactions with other genes. with the availability of many gene expression profiles, protein interaction networks, and gene regulatory networks, it is now possible to study how gene expression variation is associated with both network features and genomic factors. in the case of protein interaction networks, interaction degree, i.e., the number of interacting partners of a given protein, is one of many factors. the presence or absence of tata box and the number of transcription factors regulating a gene provide examples of genomic factors influencing gene expression variation. many studies have focused on individual factors affecting gene expression variation. for instance, newman et al. developed an experimental technique to study protein expression noise in single cells and showed that chromosomal distance to other genes and mrna-half life are associated with expression noise. however, they did not find a relationship between protein expression noise and protein-protein interactions. recently, using a more complete interaction dataset, batada et al. found that protein expression variation is negatively correlated with interaction degree when protein abundance was controlled using the data in newman et al.. this relationship continues to hold within the viable genes. several groups investigated mrna expression variation within species. for example, nelson et al. studied the relationship between the number of tissues or body parts, where the gene is expressed, and gene spacing in c. elegans and d. melanogaster. they found that gene expression variation increases in relation to the intergenic distances between genes. walther et al. found a positive correlation between the frequency of a gene's differential expression and the number of cis-regulatory elements of that gene in a. thaliana. furthermore, several groups have studied gene expression variation, also known as evolutionary variation, across different species/strains. using gene expression data from several yeast species, as well as from different strains derived from mutation-accumulation experiments, it was found, for instance, that the interspecies/interstrain variation of gene expression is significantly correlated with the presence/absence of the tata box in the promoter region. lemos et al. studied the effect of protein-protein interactions and protein length on evolutionary variation. they found that evolutionary variation is negatively correlated with protein-protein interactions in saccharomyces cerevisiae or drosophila melanogaster and negatively correlated with protein length in drosophila melanogaster. these studies highlighted the importance of protein interactions and gene regulatory regions on gene expression variation. only a few studies, however, have integrated such different data sources in a way that collectively identifies and interprets the key factors affecting gene expression variation. therefore, we conducted studies of proteomic and genomic factors marginally and collectively influencing gene expression variation across different perturbation conditions within one species: yeast. protein interactions play an important role in gene expression variation. protein-protein interactions are key biological events in a living cell, and proteins in a cell interact with each other to perform certain functions. high throughput technologies, including yeast two-hybrid systems and mass spectrometry, have generated a large amount of protein interactions in yeast. computational methods have also been developed to study the reliability of the observed interactions and to build reliable protein interaction networks. these efforts have resulted in the development of several protein interaction databases, albeit with differing degrees of reliability, including mips, dip and biogrid. from an evolutionary point of view, the expression profiles of neighboring genes of a target gene in a protein interaction network may put some constraints on the target gene's expression. thus, in a protein interaction network, the interacting partners of a specific protein can affect the corresponding gene's expression. therefore, protein physical interaction degree, i.e., the number of interacting partners of a given protein, can significantly affect gene expression variation. in the present study, we show that gene expression variation decreases with protein interaction degree and that protein interaction degree accounts for 1–2% of the expression variation in model organism yeast, a result consistent with previous studies. another key factor affecting gene expression variation is gene essentiality. genes can be classified into essential and non-essential genes based on the fitness phenotype of the yeast cell when the gene is deleted under normal growth conditions. essential genes are those that, when deleted, will render the yeast cell non-viable. non-essential genes can be further classified into no-phenotype and toxicity-modulating genes based on the fitness phenotype of yeast cell when the gene is deleted under the conditions of four dna-damaging treatments. specifically, we define a gene's toxicity modulation degree as the number of dna-damaging treatments significantly affecting the deletion strain's fitness, and 4). the higher the toxicity modulation degree, the more important the gene is in relation to cell survival. therefore, toxicity degree gives a quantitative measurement of a gene's importance to yeast cell survival. we measure a gene's functional importance in relation to cell survival by the essentiality of the essential genes and the toxicity modulation degree of non-essential genes. since the expression of genes important for cell survival are generally stable under many different stimuli and cannot fluctuate extensively, we hypothesize and show that expression variation of essential genes is lower than that of non-essential genes and decreases with toxicity degree within non-essential genes. the number of cis-elements has been shown to be positively associated with gene expression variation. the number of cis-elements is usually approximated using computational approaches and many contain false positive and negative predictions. theoretically, a given gene's expression pattern can become increasingly complex with the increasing number of transcription factors that regulate this gene, either directly or indirectly. in this study, we hypothesize that the number of tfs is a significant predictor of expression variation and show that the number of tfs regulating a gene accounts for 8–14% of its expression variation, much higher than that can be explained by the number of cis-elements. this implies the importance of indirect trans-effect on expression variation. the tata box is a conserved element in the eukaryotic promoter region and is usually bound by tata-binding proteins. the presence of tata box has been shown to be one of the most important factors contributing to gene expression variation. further analysis of the individual genomic and proteomic factors affecting gene expression indicates that there might be two distinct mechanisms that specifically influence gene expression variation of tata-containing and non-tata-containing genes. most importantly, we show that significant negative correlation between expression variation and toxicity degree is only present for tata-containing genes and that toxicity degree accounts for – % of the expression variation. in contrast, the relationship between expression variation and toxicity degree is absent for non-tata-containing genes. the fact that tata-containing genes are enriched in stress-related genes may explain this difference. although the number of tfs is significantly positively correlated with expression variation for both tata- and non-tata-containing genes, the association strength is higher for non-tata containing genes than for tata-containing genes. these results imply that the mechanism influencing tata-containing gene expression variation is much more complicated than that in non-tata-containing genes. for example, tata-containing genes were found more likely to be epigenetic regulated. thus, this study gives a more complete analysis of factors and their interaction affecting gene expression variation than may be found in previous studies.  we present our results based on the mips protein physical interaction data and the yeast gene expression profiles under ca and na exposure conditions. the results based on three interaction datasets and four other gene expression datasets, environmental stress, oxidative stress, and a combined gene expression dataset over more than conditions ) are given in the additional files and we study the expression data individually in order to minimize the variation among different laboratories. by doing so, we can also confirm whether the results based on different gene expression data are consistent. consistency of results using a variety of different datasets adds confidence to the conclusions. in this manuscript, we use genes and proteins interchangeably. we declare statistical significance if a p-value is less than without adjusting for multiple comparisons. in this study, we conducted an exploratory study of factors affecting gene expression variation. as in many epidemiological studies, we did not adjust p-values for multiple comparisons. therefore, some of our findings need to be further tested in other datasets. gene expression variation versus protein interaction degree we measured gene expression variation by the logarithm of the variance of the expression levels of each gene across the ca and na exposure conditions. we then studied the relationship between the expression variation and protein interaction degree using the lowess function in r to fit the data. on these bases, it was obvious that gene expression variation decreases with the degree of protein physical interaction. the decreasing trend is especially significant when the interaction degree is relatively low. in contrast, when the interaction degree is greater than the decreasing trend is not as obvious. biologically, this may be explained by the fact that gene expression variation stabilizes when the interaction degree is above a given threshold. another potential explanation is that interactions between proteins of high degrees are simply less reliable. the large number of less reliable interactions for proteins with high degrees can skew the true relationship between gene expression variation and interaction degrees. since the true underlying mechanism of this phenomenon is not clear, we limited our further analysis to proteins with an interaction degree of no more than protein physical interaction degree has been found to be negatively associated with gene expression variation for single cells and evolutionary expression variation between different strains/species. our result on the relationship between gene expression variation and interacting degree within s. cerevisiae is consistent with their findings. accordingly, we then used linear regression to fit the expression variation for proteins with a maximal physical interaction degree of 20: v = α + βd where v is the gene expression variation and d is the interaction degree. α and β are parameters. the fitted line and the corresponding bar-plot for the expression variation are shown in figure 1b. the gene expression variation is significantly negatively correlated with the protein interaction degree. the negative correlation between expression variation and interaction degree implies that protein with high interaction degrees do not tolerate extensive expression variation and such protein need more precise control on gene expression for an organism to function normally. gene expression variation versus essentiality, toxicity modulation, and interaction degrees as noted above, we divided genes into two classes: essential and non-essential genes. essential genes are less likely to be perturbed than non-essential genes, as significant perturbations of essential genes will, for example, render the yeast cell non-viable. we further classified the non-essential genes into five groups according to the cell's fitness phenotype changes under four dna-damaging agents, the bulky alkylating agent 4-nitroquinoline-n-oxide, the oxidizing agent tert-butyl hydroperoxide, and 254-nm uv radiation) when the non-essential genes are knocked out. since toxicity degree reflects the functional importance of genes in relation to cell survival under several dna-damaging perturbations, we expected that gene expression variation would decrease as the toxicity degree increases. since the number of genes with toxicity degree was small, we combined them with the group having toxicity degree we referred to the essential genes as the group with toxicity degree figure 2a shows the bar-plot and the linear regression fit of the gene expression variation with respect to toxicity degree. indeed, a significant negative association between gene expression variation and toxicity degree was observed. in the study of the relationship between gene essentiality and evolutionary expression variation conducted by tirosh and barkai and choi et al., they found that essential genes tend to have lower variation than non-essential genes. this is consistent with our result which demonstrates that the variation in essential genes is lower than that of non-essential genes. however, our study differs from these two studies in that we further classified the genes according to their toxicity degree and found negative association between gene expression and toxicity degree. we observed a positive correlation between interaction degree and toxicity degree and, therefore, asked whether the observed negative correlation between gene expression variation and interaction degree is, conversely, caused by the positive correlation between interaction degree and toxicity degree. we consequently studied the relationship between gene expression variation and protein interaction degree within gene groups stratified according to their toxicity degrees. using ca and na exposure gene expression data, we found a significant decreasing trend of gene expression variation with respect to interaction degree in all the strata except for the one with toxicity degree the corresponding are,,,, and for toxicity degrees and respectively. the fraction of expression variation explained by the protein interaction degree seems to increase as the toxicity degree increases. in our analyses, both toxicity and protein interaction degrees are negatively associated with gene expression variation. hence, the more important a gene is to the survival of the yeast cell, the less variation there is in its expression levels across many different conditions. similarly, the higher the interaction degree of a gene, the more stability is observed in its expression levels. biologically, a gene is important to the cell's survival since it participates in many important biological processes. any perturbation of this gene's expression will likely cause deleterious effect to the corresponding biological process and thus renders the cell non-viable. an evolutionary consequence of this hypothesis is that genes important to cell survival appear to have robust expression levels. expression variation versus gene regulatory regions: tata box, number of tfs, and toxicity degree previous studies have established the relationship between gene expression variation and the regulatory regions, including the presence/absence of tata box, the length of intergenic regions, and the number of cis-regulatory elements. we therefore asked if the observed relationship between gene expression variation and toxicity degree are the same for tata-containing genes and non-tata-containing genes. to accomplish this goal, we first stratified the yeast genes based on the presence/absence of tata boxes and reanalyzed the relationship between gene expression variation and toxicity degree. consistent with previous findings, it is clear that the gene expression variation of tata-containing genes is much higher than that of non-tata-containing genes for each fixed toxicity degree. significant negative association between gene expression variation and toxicity degree was observed for the tata-containing group. however, the association between expression variation and toxicity degree for the non-tata-containing group is only marginally significant. consistent with the other two gene expression datasets, chemostat and environmental stress, as well as the combined gene expression data from, such a highly significant negative relationship between gene expression variation and toxicity degree could be observed in tata-containing genes, but not for the non-tata-containing genes. thus, the effect of toxicity degree on gene expression variation is different for tata-containing genes versus non-tata-containing genes. the relative small r value between gene expression variation and toxicity degree alone may be explained by the absence of association between them within the non-tata-containing genes. previous studies identified the number of cis-element motifs contributing to gene expression variation in a. thaliana, but landry et al. found that the number of cis-elements only marginally affects gene expression variation in yeast. here we studied these two factors in relation to gene expression variation. we first analyzed the relationship between gene expression variation and number of cis-elements using the linear model. similar to macisaac et al., we defined cis-elements according to binding probability and classified them according to their conservation in two other yeast species. the results are presented in table for the different gene expression data. for the ca and na exposure data, the highest r is % when cis-elements that are conserved in at least one species are used. we then studied the relationship between gene expression variation and the number of tfs that influence target gene expression based on the gene regulatory network developed in hu et al.. a highly significant positive correlation between gene expression variation and the number of tfs was observed. the fraction of variation explained by the number of tfs is much higher than that by the number of cis-elements indicating that the number of tfs is a better predictor of gene expression variation than the number of cis-elements. this result implies the importance of trans-effect for gene expression variation. cis-elements are identified with three different criteria according to their conservation in two other species. β is the linear coefficient in the linear model, the p-value is related to the null hypothesis that β ≠ versus β = and the r is the fraction of variation explained by the number of cis-elements. overall analysis of factors affecting gene expression variation as enumerated above, we have identified several factors influencing gene expression variation in s. cerevisiae. in addition to the presence/absence of tata box identified in previous studies, we found that gene expression variation decreases as both the protein interaction and toxicity modulation degrees increase. these findings are consistent with other studies for expression variation of single cells and evolutionary expression variation across different strains/species. we also found that gene expression variation increases as the number of tfs or the number of cis-elements increases and that the number of tfs regulating a gene is a much better predictor of expression variation than the number of cis-elements. we therefore studied the contribution of each factor and their interactions on gene expression variation by taking the other factors into consideration using the akaike information criterion. we retained the model with the smallest aic. table gives the results of factors included in the final linear model using the mips interaction data and the four expression profiles by stepwise selection with aic. we then studied the effect of the selected factors on expression variation using linear regression, and the corresponding p-values and r values are given in table the results showed that protein interaction degree only explained less than 1% of variation when adjusted for the three other main factors and two interaction terms. except for the oxidative stress gene expression dataset, the consistently selected model contained four main factors including protein interaction degree, toxicity degree, the number of tfs, and tata box, and two interaction terms, i.e., interaction between tata and toxicity degree and tata and number of tfs. the two interaction terms are found to be statistically significant across three expression datasets: ca and na exposure, chemostat, and environmental stress and explained % – % of the variation. the four main factors include protein interaction degree, toxicity degree, number of tfs, and the presence of tata box. the protein interaction data used in this analysis is based on the mips dataset. the column marked with "√" indicates inclusion in the final linear model. the multiple linear regression is based on the final linear model, respectively. the p-value is related to the null hypothesis that β ≠ versus β = r is the variation explained by the model and each independent variable, respectively. because we found that tata box interact with toxicity degree and the number of tfs influencing gene expression variation, we reanalyzed the contributions of toxicity degree and number of tfs on gene expression variation stratified by the presence/absence of tata box using the linear model. the protein physical interaction degree was not included in this analysis because the interaction between tata box and physical interaction degree was not detected. table shows the different results for tata and non-tata gene sets. first, the number of tfs is the most significant factor for gene expression variation in both tata- and non-tata-containing genes with p-values less than 5e- for all the datasets. however, the association strength between the number of tfs and gene expression variation measured by the coefficient β for the non-tata-containing genes is about -fold higher than the corresponding values for the tata-containing genes. accordingly, the p-values for the non-tata-containing genes are about three orders of magnitude smaller than the corresponding p-values for the tata-containing gene set. we noted that tata-containing genes tend to have higher number of tfs than non-tata-containing genes, but this cannot explain our observation that the association strength between expression variation and the number of tfs within the tata-containing genes is lower than that within the non-tata-containing genes. one possible explanation is that the presence of tata-box weakens the effect of the number of tfs on gene expression during evolution. second, the toxicity degree is only a significant contributor for gene expression variation within the tata-containing genes in the ca-na exposure and chemostat gene expression dataset, while it is not a significant contributor for gene expression variation in the non-tata-containing gene set. within the environmental stress gene expression dataset, the effect of toxicity degree on gene expression variation was not statistically significant. however, we did observe a decreasing trend of gene expression variation with respect to toxicity degree within the tata-containing genes. the linear model that includes the toxicity degree and the number of tfs is built for each of the four gene expression datasets, respectively. r is the variation explained by each independent factor and the model, respectively. β is the linear coefficient in the linear model, and the p-value is related to the null hypothesis that β ≠ versus β = the toxicity degree of a gene measures the tolerance of the yeast cell to different external stress conditions when the gene is knocked out. therefore, it might be expected that the different relationship between gene expression variation and toxicity degree for tata-containing genes and non-tata-containing genes is due to the enrichment of stress-related genes in tata-containing genes, as found in. to test this hypothesis, we used a set of genes related to environmental stress response. if the hypothesis is true, we would expect a higher association between expression variation and toxicity degree within the esr genes than that within the non-esr genes. however, our data shows that the β values are similar for the two groups of genes. thus we cannot explain the interaction between tata and toxicity degree by the enrichment of esr genes in tata-containing genes. the biological mechanisms underlying the observed interaction are not clear and need to be further studied. the linear model is built for each of the four gene expression datasets, respectively. β is the linear coefficient in the linear model and the p-value is related to the null hypothesis that β ≠ versus β = we also did the same analysis for the average gene expression variation across the four expression datasets and the combined gene expression data of landry et al., and the results are presented as additional file the same conclusions can be obtained indicating the robustness of our results. previous studies showed that tata- and non-tata-containing genes might recruit different coactivator complexes for gene expression. tata-containing genes were also found to be subject to greater nucleosomal regulation than non-tata-containing genes. basehoar et al. suggested that two distinct regulatory mechanisms may be present at tata- and tata-less promoters. the results in table support their findings. the results based on the oxidative stress gene expression dataset are not consistent with the results based on the other three gene expression datasets. this observation may be due to the relatively small gene expression variation in this data. for example, the range of the variance of the expression levels within the oxidative stress dataset,, is much smaller than the corresponding ranges,,, and, for the ca and na exposure, chemostat, and environmental stress conditions, respectively. we also studied the contributing factors for gene expression variation using the dip and biogrid interactions, and the results are given in additional file and file similar conclusions as those based on the mips interaction data were obtained. the consistency of the results using different combinations of protein interaction data sets and gene expression profiles showed the robustness of our conclusions. however, the fraction of gene expression variation explained by all factors is less than 25%. one possible explanation is that the measurement of gene expression changes and other factors, including the toxicity degree and interaction degree, are still very noisy. we expect that the true r would be higher than that observed in this study.  we implemented a system-wide analysis of proteomic and genomic factors affecting gene expression variation. among four different factors, tata-box and the number of tfs are the most important factors influencing gene expression variation. the influence of tata-box on evolutionary gene expression variation has been extensively studied both computationally and experimentally, and our results are consistent with their findings. although it is intuitive that the number of tfs regulating a gene should have a significant effect on the gene's expression variation, the magnitude of its influence has not been studied in large scale expression datasets to the best of our knowledge. our findings demonstrated that the gene regulation is a main factor affecting gene expression variation. protein interaction degree and toxicity degree do not account for as much variation when compared to the influence of the number of tfs and the tata-box. in our overall analysis, we also found the interactions between tata-box and toxicity degree as well as the number of tfs influence expression variation. the further study stratified by tata-box indicated that tata-containing genes and non-tata containing genes behave differently in relation to the toxicity degree and the number of tfs. the effect of the number of tfs on expression variation within the tata-containing genes is lower than that for the non-tata-containing genes. on the other hand, toxicity degree is associated with expression variation within the tata-containing genes only. these findings suggest that the regulatory mechanism might be more complicated for tata-containing genes than non-tata containing genes. 
previous studies have shown that the presence of tata box, the number of cis-regulatory elements, gene essentiality, and protein interactions significantly affect gene expression variation. mutations in both cis-regulatory elements and transcription factors regulating a gene affect the gene's expression and thus we study the relationship between gene expression variation and the number of tfs regulating a gene. results showed that 1) gene expression variation is negatively correlated with the protein interaction degree in the protein interaction network, 2) essential genes tend to have less expression variation than non-essential genes and gene expression variation decreases with toxicity degree, and 3) the number of tfs regulating a gene is the most important factor influencing gene expression variation. identifying factors affecting gene expression variation is a challenging problem in genetics. nonetheless, the need to obtain a more complete understanding of such factors and how their interactions influence gene expression variation remains a challenge. since toxicity degree reflects a gene's importance to cell survival under dna-damaging conditions, we expect that it is negatively associated with gene expression variation. most importantly we study how these factors interact with each other influencing gene expression variation. using yeast as a model system, we evaluated the effects of four separate factors and their interactions on gene expression variation: protein interaction degree, toxicity degree, number of tfs, and the presence of tata box. in addition, the number of tfs regulating a gene was found to be an important factor influencing gene expression variation for both tata-containing and non-tata-containing genes, but with different association strength. moreover, gene expression variation was significantly negatively correlated with toxicity degree only for tata-containing genes. the finding that distinct mechanisms may influence gene expression variation in tata-containing and non-tata-containing genes, provides new insights into the mechanisms that underlie the evolution of gene expression. the growth rates of yeast cells under several dna-damaging conditions have been studied and a gene's toxicity degree is defined as the number of such conditions that the growth rate of the yeast deletion strain is significantly affected.
 identifying factors affecting gene expression variation is a challenging problem in genetics. previous studies have shown that the presence of tata box, the number of cis-regulatory elements, gene essentiality, and protein interactions significantly affect gene expression variation. nonetheless, the need to obtain a more complete understanding of such factors and how their interactions influence gene expression variation remains a challenge. the growth rates of yeast cells under several dna-damaging conditions have been studied and a gene's toxicity degree is defined as the number of such conditions that the growth rate of the yeast deletion strain is significantly affected. since toxicity degree reflects a gene's importance to cell survival under dna-damaging conditions, we expect that it is negatively associated with gene expression variation. mutations in both cis-regulatory elements and transcription factors regulating a gene affect the gene's expression and thus we study the relationship between gene expression variation and the number of tfs regulating a gene. most importantly we study how these factors interact with each other influencing gene expression variation.  using yeast as a model system, we evaluated the effects of four separate factors and their interactions on gene expression variation: protein interaction degree, toxicity degree, number of tfs, and the presence of tata box. results showed that 1) gene expression variation is negatively correlated with the protein interaction degree in the protein interaction network, 2) essential genes tend to have less expression variation than non-essential genes and gene expression variation decreases with toxicity degree, and 3) the number of tfs regulating a gene is the most important factor influencing gene expression variation. in addition, the number of tfs regulating a gene was found to be an important factor influencing gene expression variation for both tata-containing and non-tata-containing genes, but with different association strength. moreover, gene expression variation was significantly negatively correlated with toxicity degree only for tata-containing genes.  the finding that distinct mechanisms may influence gene expression variation in tata-containing and non-tata-containing genes, provides new insights into the mechanisms that underlie the evolution of gene expression. 