Identifying stably expressed genes from multiple RNA-Seq data sets

RNA-Seq data collection and processing

Overview of the RNA-Seq data sets

We examined RNA-Seq data from 49 Arabidopsis experiments stored on the NCBI GEO repository (see more details below). After screening, we retained data from 211 biological samples in 24 experiments. To illustrate our methods for finding stably expressed genes, we divided the experiments into three groups: the seedling group contains 60 Arabidopsis seedling samples from 9 experiments; the leaf group contains 60 Arabidopsis leaf samples from 5 experiments; the multi-tissue group contains 91 samples from 10 experiments on multiple tissue types (shoot apical, root tip, primary root, inflorescences and siliques, hypocotyl, flower, carpels, aerial tissue, epidermis, seed). Table 1 summarizes the basic information about the three groups (see Table S1 for more details).

Table 1:

Summary statistics for the three groups of Arabidopsis samples.

Group	# Experiments	# Treatments	# Samples	# Genes
Seedling	9	27	60	22,207
Leaf	5	28	60	20,967
Multi-tissue	10	39	91	23,611

DOI: 10.7717/peerj.2791/table-1

To find stably expressed genes in each group, we processed the raw sequencing data and summarized the results as count matrices of mapped RNA-Seq short reads (see details below). We removed genes with low mean numbers (less than 3) of mapped read counts for all experiments. Such genes tend to be more prone to sequencing noise, less interesting to biologists, and also cause convergence issues when fitting statistical models. Many other researchers (such as Anders et al., 2013) recommend removing such genes before analyzing RNA-Seq data. The number of remaining genes in each group is also summarized in Table 1. Figure 1 shows the numbers and overlap of the genes after this step.

Count normalization

As explained in the introduction, count normalization is needed when analyzing RNA-Seq data to (1) adjust for differences in sequencing depths or library sizes; (2) to adjust for the apparent changes in relative read frequencies of non-DE genes that occur as a consequence of changes in relative read frequencies of truly DE genes.

For the second type of adjustment, we follow Anders and Huber’s method (Anders & Huber, 2010) for estimating normalization factors. Let y_ij denote the read count for ith gene of the jth sample (m genes and n samples in total). We first create a pseudo-reference sample where each gene’s expression value is the geometric mean expression over all real samples for that gene, (1)${\displaystyle y_{i,0}={\left(\prod _{j=1}^n{y}_{i,j}\right)}^{1∕n},i=1,\dots ,m.}$ Next we calculate the median fold-change in relative frequency between each sample j and the pseudo-reference sample, (2)${\displaystyle R_j^{\prime }=\mathrm{median}\left(\frac{y_{1,j}∕N_j}{y_{1,0}∕N_0},\dots ,\frac{y_{m,j}∕N_j}{y_{m,0}∕N_0}\right),}$ where N_j is the library size for sample j (the sum of RNA-Seq counts mapped to all genes retained in each sample). Finally, the normalization factor R_j for sample j is calculated as (3)${\displaystyle R_j=\frac{R_j^{\prime }}{{\left(\prod _{j=1}^n{R}_j^{\prime }\right)}^{1∕n}}.}$ Using the estimated normalization factors, the relative frequencies will be computed as y_ij∕N_jR_j, which we will call the normalized relative frequency for gene i in sample j. The assumption made here is that the median fold change between normalized relative frequencies in two samples should be 1. In other words, this normalization method assumes that the majority of genes are not DE. The NBPSeq package (Di, Schafer & Di, 2013) has an inbuilt function for this procedure and it will be used for count normalization in this paper. With the estimates from Eq. (3), we see that the median fold change in normalized relative frequencies between each sample and the pseudo-reference sample will be set to 1: (4)${\displaystyle \mathrm{median}\left(\frac{y_{1,j}∕N_j{R}_j}{y_{1,0}∕N_0{R}_0},\dots ,\frac{y_{m,j}∕N_j{R}_j}{y_{m,0}∕N_0{R}_0}\right)=1,}$ where $R_0={\left({\prod }_{j=1}^n{R}_j^{\prime }\right)}^{-1∕n}$.

We can apply Eq. (2) to a subset of reference genes to estimate normalization factors. In doing so, effectively, the median fold change in Eq. (4) among the reference genes will be set to 1 in each sample j. Other normalization methods may make different assumptions than Anders and Huber’s, but some assumptions of a similar nature seem unavoidable. For example, the TMM method of Robinson & Oshlack (2010) is based on a similar principle: assuming the majority of the genes are not DE. The TMM method can be applied to a subset of genes selected based on an initial screening of mean expression level and fold changes. In TMM method, one can also specify certain quantile (instead of the median) of the fold changes to be 1.

In this paper, we will identify stably expressed genes from multiple data sets based on numerical measure and use them as reference for estimating normalization factors (from Eqs. (2) and (3)). However, to identify the stably expressed genes, we first need a set of initially estimated normalization factors. To tackle this circular dependence, we use a one-step iteration method to estimate the normalization factors:

• 1.

  First, we use all the genes to calculate the initial normalization factors;

• 2.

  Then, we fit a GLMM to each gene and estimate the total variance measure, incorporating the initial normalization factors as an offset term (see ‘Poisson log-linear mixed-effects regression model and the total variance measure of expression stability’);

• 3.

  Next, we select the top 1,000 stably expressed genes based on the total variance measure estimated from step 2 above, and use them as reference genes to recalculate the normalization factors.

In practice, this one-step method seems to be adequate and further iterations will only slightly change the set of 1,000 stably expressed genes. For example, for the multi-tissue group of experiments, if we were to run one more iteration of steps 2 and 3, there would be 946 overlapping genes between the top 1,000 genes from the first iteration and those from the second iteration.

Poisson log-linear mixed-effects regression model and the total variance measure of expression stability

We fit a Poisson log-linear mixed-effects regression model to the RNA-Seq counts mapped to each gene and measure gene expression stability using a total variance measure. Let Y_ijkl be the number of RNA-Seq reads mapped to gene i in sample j from treatment group k in experiment l. We will fit regression models to each gene separately and suppress subscript i from the model equations. For each gene, we fit a Poisson log-linear mixed-effects regression model (5)${\displaystyle Y_{jkl}\sim \mathrm{Poisson}\left({\mu }_{jkl}\right),}$ (6)${\displaystyle log\left({\mu }_{jkl}\right)=log\left(R_{jkl}{N}_{jkl}\right)+\xi +{\alpha }_l+{\beta }_{k\left(l\right)}+{ϵ}_{jkl},}$ which is a specific type of generalized linear mixed model (GLMM, McCulloch & Neuhaus, (2001)). In Eq. (6), N_jkl and R_jkl are the library size and normalization factor discussed in ‘Count normalization’. We will call R_jklN_jkl the normalized library size. The parameter ξ is a fixed-effect term for the baseline log mean of the relative counts (counts divided by the normalized library sizes). The values α, β, and ϵ represent the experiment effect, the treatment effect (nested within each experiment), and the sample effect respectively. We view α, β and ϵ as random effects and assume that they are independent and follow normal distributions: (7)${\displaystyle {\alpha }_l\sim N\left(0,{\sigma }_{\text{experiment}}^2\right),{\beta }_{k\left(l\right)}\sim N\left(0,{\sigma }_{\text{treatment}}^2\right),{ϵ}_{jkl}\sim N\left(0,{\sigma }_{\text{sample}}^2\right),}$ where ${\sigma }_{\text{experiment}}^2$, ${\sigma }_{\text{treatment}}^2$ and ${\sigma }_{\text{sample}}^2$ are called variance-components—they quantify the overall variances of the corresponding random effect terms.

The sample effect ϵ represents the extra-Poisson variation in read counts among samples in the same treatment group and ${\sigma }_{\text{sample}}^2$ plays a similar role as the over-dispersion parameter in a negative binomial model (Anders & Huber, 2010; Di et al., 2011). The experiment effect, α, accounts for all sources of variation at the experiment level, including differences in lab personnel and conditions, day light hours, age of the plants, temperature, sequencing platform, and other unidentified sources. The contributions from these different experiment-level sources are often difficult to separate statistically. We treat the experiment effect α as a random effect because we view the collected experiments as a random sample from the pool of all Arabidopsis RNA-Seq experiments. We also treat the treatment effect β as a random effect. In a DE test, β is usually considered as a fixed-effect term. Here for evaluation of expression stability, we are not interested in the specific levels of the individual β’s and focus more on the overall variation of β under a range of treatment types.

We define the stability measure as the estimated total variance, (8)${\displaystyle {\hat{\sigma }}^2={\hat{\sigma }}_{\text{sample}}^2+{\hat{\sigma }}_{\text{treatment}}^2+{\hat{\sigma }}_{\text{experiment}}^2.}$ The parameters $\left(\xi ,{\sigma }_{\text{experiment}}^2,{\sigma }_{\text{treatment}}^2,{\sigma }_{\text{sample}}^2\right)$ are estimated using the glmer() function of the R package lme4 (Bates, Mächler & Bolker (2015), version 1.1.7), which uses a Gaussian-Hermite quadrature to approximate the likelihood function. We rank all the genes according to their values of ${\hat{\sigma }}^2$ in increasing order (smallest first), and consider highly ranked (e.g., top 1,000) genes to be stably expressed.

Normal models (Eq. (7)) are commonly assumed for the random effects in the GLMM settings. The normality assumption is likely a simplification of reality, yet it is a good starting point and should be adequate for finding genes with low total variation—the stably expressed ones.

Other stability measures

The assessment of gene expression stability depends on the specific stability measure used. Czechowski et al. (2005) and Dekkers et al. (2012) used the coefficient of variation (CV) measure, computed as standard deviation divided by mean, to find stably expressed genes from microarray data.

The M-value in geNorm (Vandesompele et al., 2002) is a well-cited measure. For a set of m₀ genes, the M-value measure works as follows: First, the pairwise variation between gene i₁ and gene i₂ is calculated as the standard deviation of the log fold changes between their expression levels across all the n samples: ${\displaystyle V_{i_1,i_2}=\mathit{st.dev}\left\{log\left(\frac{y_{1,i_1}}{y_{1,i_2}}\right),\dots ,log\left(\frac{y_{n,i_1}}{y_{n,i_2}}\right)\right\}.}$ Next, the M-value for gene i is defined as the average pairwise variation between gene i and all other genes ${\displaystyle M_i=\frac{\sum _{k\ne i}{V}_{i,k}}{m_0-1}.}$

In the Results section, we compare the M-value to the total variance measure on RNA-Seq data from the multi-tissue group experiments, and compare the stably expressed genes identified from these two measures to those identified from microarray data using the CV measure.

Stably Expressed Genes

Using the total variance, ${\hat{\sigma }}^2$, from the GLMM (see Eq. (6) in ‘Poisson log-linear mixed-effects regression model and the total variance measure of expression stability’) as a stability measure, we identified stably expressed genes from the three groups of experiments described in ‘RNA-Seq data collection and processing’: the group of seedling experiments, the group of leaf experiments, and the group of experiments on different tissue types (see Table 1 for a summary). As we mentioned in the Introduction, absolutely stably expressed genes may not exist. Choosing different sample sets as reference allows us to identify stably expressed genes for different biological contexts.

In Tables S2–S4, we summarize the top 1,000 most stably expressed genes in each group. In Fig. 2, we provide the histograms of the mean Count Per Million (CPM) for the 1,000 most stably expressed genes identified in each group. For each gene, the CPM is computed as (9)${\displaystyle \frac{\text{count}\times 10^6}{\text{normalized library size}}}$ in each sample and the mean is computed over all samples.

The lists of the top 1,000 genes in the three groups share 104 genes in common (see Table S5 for more details). These genes are stably expressed under a wide range of experimental conditions and in different tissue types, and thus may be worth further study. This list of 104 genes has significant overlap with the top 100 stably expressed genes identified by Czechowski et al. (2005) from a developmental series of microarray samples: 9 out of these 104 genes (see Table S6 for details),

Comparison to house-keeping genes and stably expressed genes identified from microarray data

Czechowski et al. (2005) discussed the expression stability of house-keeping genes and showed that the house-keeping genes are not stably expressed according to their numerical measure. In particular, they compared the expression profiles of five traditional house-keeping genes (AT1G13440, AT3G18780, AT4G05320, AT5G12250, AT5G60390) and five genes (AT1G13320, AT5G59830, AT2G28390, AT4G33380 and AT4G34270) that they identified as stably expressed according to the CV measure from a developmental series of microarray experiments (see Fig.1 of that paper). In Fig. 3, we compare the expression profiles of these 10 genes from Czechowski et al. (2005) to the expression profiles of five genes (AT1G64840, AT1G75420, AT2G32910, AT3G51310, AT5G48340) that we randomly selected from the top 100 most stably expressed genes identified from the multi-tissue group RNA-Seq data according the total variance ${\hat{\sigma }}^2$. For each of the 15 genes, Fig. 3 shows the expression levels measured in CPM over 91 samples in the eight experiments in the multi-tissue group, and Table 2 summarizes the variance components estimated from the GLMM in ‘Poisson log-linear mixed-effects regression model and the total variance measure of expression stability’.

The five house-keeping genes show large total variation with all three variance-components relatively large as compared to the other 10 genes. This is consistent with Czechowski’s observation that house-keeping genes are not necessarily stably expressed according to a numerical measure. Three of the five stably-expressed genes identified by Czechowski et al. (2005) are among the top 1,000 stably-expressed genes according to our stability measure, the total variance ${\hat{\sigma }}^2$. Czechowski et al. (2005) identified those five genes from microarray data and different experiments. It is not too surprising those genes might not be the most stable in RNA-Seq experiments: the two technologies differ in many aspects including coverage and sensitivity.

The house-keeping genes identified in Czechowski et al. (2005) tend to have higher CPM. This is partly due to a selection preference: the authors there intentionally found genes with higher CPM for use as references so that they can be observed in most of the experiments. As we will explain later, we suggest using a collection of 100–1,000 genes as reference gene set for normalization, we did not specifically target for genes with high CPM.

Factors affecting stability ranking

The previous two subsections demonstrate that when using a numerical measure to quantify gene expression stability, the outcome is dependent on (1) the biological context reflected in the reference sample set used and (2) the technology used for measuring gene expression. It should also be intuitive, and we will further clarify in the second half of this subsection, that the stability ranking is also dependent on (3) the specific numerical measure used. In this section, we will first compare the lists of stably-expressed genes identified under different scenarios where one or more of the above three factors differ. We then further discuss the subtle roles played by the specific stability measure and the reference gene set by comparing the total variance ${\hat{\sigma }}^2$ measure from the GLMM (see Eq. (6)) to the M-value measure used in the geNorm method (Vandesompele et al., 2002).

We look at an additional five lists of stably expressed genes identified under different scenarios and examine how each of these five lists overlaps with the the top stably-expressed genes identified from the multi-tissue group of RNA-Seq experiments according to the total variance measure ${\hat{\sigma }}^2$ (see ‘Poisson log-linear mixed-effects regression model and the total variance measure of expression stability’). The five lists are:

• L₁:

  100 top stably expressed genes from the multi-tissue group according to the M-value in geNorm (applied to (count + 1)) of Vandesompele et al. ;

• L₂:

  100 top stably expressed genes from the seedling group according to the total variance ${\hat{\sigma }}^2$ from the GLMM;

• L₃:

  100 top stably expressed genes from the leaf group according to the total variance ${\hat{\sigma }}^2$ from the GLMM;

• L₄:

  100 stably expressed genes identified from a developmental series of microarray experiments by Czechowski et al. (2005) using the CV measure (see ‘Other stability measures’);

• L₅:

  50 stably expressed genes identified by Dekkers et al. (2012) from microarray seed experiments using the CV measure.

For each list L_i above, we measure how it overlaps with the top stably expressed genes (the reference set) from the multi-tissue group using the recall percentage (10)${\displaystyle \frac{\#\left\{L_i\cap \text{reference set}\right\}}{\#\left\{L_i\right\}}\times 100,}$ where #{} denotes the number of elements in the list. In Fig. 4, we plot the recall percentage versus the number of top stably-expressed genes we selected as reference from the multi-tissue group.

We have the following observations:

• 1.

  The list L₁ is identified from the same set of RNA-Seq experiments as the reference sets, but using a different stability measure (M-value in geNorm). This list has significant overlap with the top stably-expressed genes identified using the total variance measure: 29 and 98 out of the 100 genes from the list L₁ are among the top 100 and 1,000 most stably-expressed genes, respectively, from the multi-tissue group identified using the total variance measure.

• 2.

  The lists L₂ and L₃ are identified from different sets of RNA-Seq experiments (leaf and seedling experiments) using the same stability measure as used for the reference sets. The lists L₄ and L₅ are identified from microarray experiments (a developmental series and a seed group) and using the CV measure. The overlapping (recall) percentages are still statistically significant, but much less than in the case of L₁. This shows that differences in tissue type and in measuring technology both influence the expression stability ranking, and to comparable degrees. The lists L₃ and L₅ have the least overlapping percentages with the reference sets. These lists are identified from a leaf group and a seed group respectively. Our understanding is that the leaf group and the seed group are more biologically homogeneous than the multi-tissue group and thus provide very different biological contexts for evaluating expression stability.

When applied to the same set of samples, the M-value and total variance measure ${\hat{\sigma }}^2$ give similar expression stability ranking: the rank correlation is 0.97 (see also, observation 1 above). We point out that the reason is because the M-value and normalization step needed for computing our total variance measure have similar fundamental assumptions. The basic principle behind the M-value is that the expression ratio of two stably-expressed genes should be identical in all samples. In formula, it means that the expression values of two stably-expressed genes i₁, i₂ in any two samples j₁, j₂ should satisfy (11)${\displaystyle \frac{y_{i_1,j_1}}{y_{i_2,j_1}}=\frac{y_{i_1,j_2}}{y_{i_2,j_2}}.}$ Our total variance measure ${\hat{\sigma }}^2$ is estimated from normalized data. The basic assumption in the normalization step is that majority of genes are not DE. In formula, it means for any stably-expressed gene i₁, its expression level as measured by the relative frequency should be stable across all samples, (12)${\displaystyle \frac{y_{i_1,j_1}}{S_{j_1}}=\frac{y_{i_1,j_2}}{S_{j_2}},}$ where S_j1 to S_j2 are the normalized library sizes (i.e., R_jN_j in Eq. (6)). This implies for any two stably-expressed genes i₁ and i₂ (13)${\displaystyle \frac{y_{i_1,j_1}}{y_{i_1,j_2}}=\frac{y_{i_2,j_1}}{y_{i_2,j_2}}=\frac{S_{j_1}}{S_{j_2}}.}$ The first equation in Eq. (13) is equivalent to Eq. (11). (In practical application of both methods, the stability of any single gene is evaluated by comparing its expression to a set of reference genes. See the Method ‘Count normalization’ for more details.)

In practice, the geNorm program (Vandesompele et al., 2002) is frequently used to rank a set of reference genes identified from other methods. An iterative elimination procedure is used along with the M-value to determine the final ranks of the expression stability: after each iteration, the gene receiving the largest M-value will be removed and a new set of M-values will be computed for the remaining genes, and the iteration will go on until there are only two genes left. We did not use such an iterative procedure in the comparisons above (i.e., we only computed one set of M-values for all genes). We provided some comments about the iterative elimination procedure in the Appendix.

Sources of variation

For each gene, the GLMM (Eq. (6) of ‘Poisson log-linear mixed-effects regression model and the total variance measure of expression stability’) allows us to decompose total count variance into between-sample, between-treatment and between-experiment variance components. The estimated variance components tell us how much each component contributes to the overall count variation. Table 3 summarizes the percentages—averaged over all genes—of the total variance attributable to each of the three components for three groups of RNA-Seq samples (seedling, leaf and multi-tissue groups in ‘RNA-Seq data collection and processing’). Fig. 5 shows the histograms of the percentages. Figure 6 shows the stacked bar plot of variance components estimated from the multi-tissue group for 20 genes randomly selected from the top 1,000 stably expressed genes and 20 genes randomly selected from 23,611 genes. As expected, the between-experiment variance component, on average, explains the largest proportion of the total variation. The between-experiment variation is relatively smaller among the leaf samples, indicating that the leaf samples are more homogeneous. There is more variation in the relative percentages of total variance explained by the between-sample and between-treatment variance components. In principle, the between-treatment variation will be greater when there is a higher proportion of DE genes or when the samples are more homogeneous. In practice, the between-sample variance depends greatly on what samples are used as biological replicates.

Reference gene set for normalization

Once we have ranked the genes according to our numerical stability measure (i.e., the total variance measure, ${\hat{\sigma }}^2$), one application is to use an explicit set of most stably expressed genes as reference genes for count normalization. This new approach allows investigators to prescribe a specific biological context for evaluating gene stability by choosing the most relevant reference samples and experiments when computing the stability measure. For example, the most stably expressed genes identified from the multi-tissue group and those identified from the seedling group will provide different biological contexts. In contrast, existing normalization approaches are often applied to the single data set under study, and thus provide a single, narrow context.

Even under a specific biological context, it is almost impossible to know whether the genes in any reference set are absolutely stably expressed, even though commonly used normalization methods often enforce some assumptions on the reference gene set: for example, when we use Anders and Huber’s method to estimate the normalization factors based on a subset of reference genes, roughly speaking, the median fold change among the reference genes will be set to 1 (see ‘Count normalization’ for more details). A subtle point we want to make is that since it is impossible to know how well such or similar assumptions on DE hold for a reference gene set, we can improve the interpretability of the DE test results by making the reference gene set explicit: we can slightly change our perspective and interpret all DE results as relative to the reference gene set. For example, a fold change of 2 inferred from the GLMM model can be interpreted as the fold change of a gene is 2 times the true (but often unknowable) median fold change of the reference genes. When one estimates the normalization factors based on all genes, one is effectively specifying an implicit set of genes as a reference set. Our proposal is to make the reference set explicit and interpret DE results as relative to the reference gene set.

Interpreting the DE results as relative to an explicit reference set is especially beneficial when one wants to compare DE results from an experiment to previously published results. When the interest is in comparing different experiments, we recommend using a common reference set. For example, when two RNA-Seq data sets are separately normalized with different reference sets, a fold change of two observed in one experiment may not be directly comparable to a fold change of two observed in the other. This concern can be alleviated by using a common set of reference genes. We use a toy example to illustrate this point in Table 4 where we examine the mean counts for 5 genes in two two-group comparison experiments. If we use different reference gene sets for count normalization in the two experiments, for example, we use genes 1–3 as reference in experiment 1, but use genes 3–5 as reference in experiment 2, we may conclude that gene 3 is not DE in either experiment. If we use a common reference gene set—either genes 1–3 or genes 3–5—for normalization, however, we will be able to discover, in either case, that the DE behavior of gene 3 is different in the two experiments. Note that the DE conclusion in both experiments will depend on the reference genes used: if genes 1–3 are used as reference, gene 3 is not DE in experiment 1, but will be DE in experiment 2; if genes 3–5 are used as reference, gene 3 will be considered DE in experiment 1, but not DE in experiment 2. The point is, in either case, we will notice that the DE behavior of gene 3 is different between the two experiments. This information will be lost if one uses different reference sets to assess DE in the two experiments.

In practice, we recommend using the top 1,000 most stably expressed genes for estimating normalization factors. The key is to avoid using too few (e.g., less than 10) or too many (e.g., using all genes) reference genes: intuitively, using too few, the estimates will be unstable; using too many, the results may be subject to influence from highly unstable genes. Our simple simulations suggest that using between 100 to 10,000 genes seems to give stable results. In the first set of three examples, we used Anders and Huber’s method (see Eq. (2)) to estimate normalization factors for samples in each of the seedling, leaf and multi-tissue groups of experiments (see ‘RNA-Seq data collection and processing’). We used the top 10, 100, 1,000, and 10,000 stably expressed genes identified earlier (see ‘Stably expressed genes’ for details) as reference gene sets. Figure 7 shows the pairwise scatter plots and correlation coefficient between the normalization factors when different numbers of top stable genes are used as reference. A stronger correlation indicates the normalization factors estimated from the two settings are highly consistent. The plots and correlation coefficients suggest using between 100 and 1,000 genes tend to give similar normalization factor estimates. We also used the top 10, 100, 1,000, and 10,000 stably expressed genes identified from the multi-tissue group as reference set for estimating normalization factors for a set of 48 root samples from a new experiment (GSE64410, Vragović et al. (2015)). The largest Pearson correlation 0.993 is between the normalization factors estimated using the top 100 and top 1,000 stably expressed genes as reference. Based on the above observations, using 1,000 most stably expressed genes as reference seems to be a reasonable heuristic rule.

An example

In this part, we illustrate the effect of using different reference gene sets for computing normalization factors on a real data set and explain the implication on DE analysis.

Wang et al. (2012) performed RNA-Seq experiments using 10-day old seedlings to investigate the role of Arabidopsis SNW/Ski-interacting (SKIP) protein on transcriptome-wide changes in alternative splicing. Two biological replicates each from wild type (Col-0) and skip-2 mutant were compared. We retrieved and processed the raw RNA-Seq data from this experiment using our pipeline (see ‘RNA-Seq data collection and processing’, accession number GSE32216). For this data set, the normalization factors for the four samples (two wild types followed by two mutants) estimated using all genes, (0.84, 0.62, 1.38, 1.39), differ markedly from the normalization factors, (0.71, 0.54, 1.59, 1.63), estimated using the 1,000 reference genes that we identified using the total variance measure from the seedling group (see ‘Stably expressed genes’).

The implication on DE analysis is that if we use the 1,000 stably expressed genes for normalization, we will expect to see more under-expressed genes and less over-expressed gene in the mutant group relative to the wild type group. The two sets of estimated normalization factors reflect different assumptions: roughly speaking, when using all genes to compute the normalization factors, the assumption is that median fold change among all genes is 1; when using the 1,000 reference genes to compute the normalization factors, the assumption is that the median fold change among the set of 1,000 genes is 1. It is difficult to know which assumption is more reasonable without additional biological insights. However, the benefit of using an explicit set of 1,000 genes as reference is the improved interpretability by making the reference gene set, and thus the implied assumption, more explicit. Furthermore, if one wants to compare the DE results from this experiment to other DE results from the collection of seedling experiments, then one should use a common reference set of genes for count normalization.