No counts, no variance: allowing for loss of degrees of freedom when assessing biological variability from RNA-seq data

2.1 Definition of the reduced residual d.f.

Let μ^gi denote the fitted value for library i in the GLM for this gene. Assume that there is a set of libraries Z_g for which μ^_gi = 0 in each library i ∈ Z_g. For each of these libraries, the count y_gi must also be zero as all counts are non-negative and any positive count would result in μ^_gi > 0. Knowing that μ^_gi = 0 means that y_gi = 0, which can be determined without the need to know any other count. Thus, libraries in Z_g will not contribute any d.f. to the fit. The canonical definition for the residual d.f. fails to consider this subtlety. Consequently, the computed value for d may be larger than the true residual d.f. The latter is denoted as d_g and will be referred to as the “reduced residual d.f.” in the following text.

The value of d_g can be determined by identifying and ignoring the libraries in Z_g. In practice, this is done by removing all libraries with fitted values below some arbitrarily small value like 10⁻⁴. This ensures that libraries in Z_g are not overlooked due to numerical imprecision of the model fit. The effective number of libraries for gene g becomes n−|Zg|. Similarly, a refined design matrix X_g can be constructed by removing all rows in X corresponding to i ∈ Z_g. When X_g is fitted to the counts for the remaining libraries, the reduced residual d.f. is defined as

where the column rank of X_g is simply the number of independent coefficients in X_g. This can be obtained by performing a QR decomposition on X_g and counting the number of non-zero diagonal elements in the resulting R matrix. The above expression for d_g is appropriate as it explicitly ignores the libraries that do not contribute any d.f. to the model fit. This avoids any overestimation of the residual d.f. that might occur with the canonical definition. Of course, d_g = d if Z_g is empty as the rank of X_g is equal to p.

2.2 Overview of the canonical QL framework

Lund et al. (2012) assume that σg2∼σ02χd0−2 across all genes g, where d₀ is the prior d.f. and σ02 is a constant scaling factor. This means that σ^g2∼σ02F(d,d0). Both d₀ and σ02 can be simultaneously estimated from the distribution of σ^g2 across all g (Smyth, 2004; Phipson et al., 2016). The shrunken dispersion in the EB framework is defined as

which effectively squeezes σ^g2 towards the estimated σ02 for each gene.

The QL F-test is used to test a pre-specified null hypothesis by comparing nested designs. One or more columns are chosen and removed from X to obtain the null design matrix X₀. The moderated F-statistic uses the shrunken dispersion and is defined as

where t is the difference in the number of coefficients between designs and LR_g is the likelihood ratio, i.e. the difference in the total residual deviances between GLMs fitted with X₀ and X. Under the null hypothesis, F_g will be F-distributed on t and d + d₀ d.f.’s. This can be used to compute a p-value to reject or accept the null for each gene. In practice, a lower bound for the p-value is defined by using the likelihood ratio test after fitting a Poisson GLM to the counts for each gene. This reflects the fact that RNA-seq data should exhibit at least Poisson-level variance due to sequencing noise (Marioni et al., 2008).

2.3 Redefined statistics with the true residual d.f.

The canonical definitions are only correct for gene g when d = d_g. If these two values are not equal, all occurrences of d should be replaced with d_g. For simplicity, assume that all libraries i ∈ Z_g have true means close to zero. This is reasonable in most replicated designs, where a fitted value of zero from a large true mean would require zero counts for all replicate observations (which is unlikely). The assumption means that the libraries in Z_g do not contribute to any sampled instance of D_g, i.e. the unit deviance from each library will always be zero. Thus, Dg|σg2∼σg2χdg2 as only those libraries contributing d.f. are considered. The correct estimator of the QL dispersion for gene g is defined as

This expression ensures that E(σ^g(2)2|σg2) is still equal to σg2 when Z_g is not empty. EB shrinkage should be performed using the distribution of σ^g(2)2 across all genes, instead of σ^g2. This ensures that the correct scaling factor σ0(2)2 and prior d.f. d0(2) are estimated. The shrunken dispersion for each gene is similarly redefined as

while the moderated F-statistic is redefined as

This is F-distributed on t and d_g + d₀₍₂₎ d.f.’s, and can be used to compute a p-value as previously described.

3.1 The canonical method underestimates the dispersion

Consider a one-way layout with two replicate libraries in each of four conditions A₁, A₂, B₁ and B₂. NB-distributed counts for 10,000 genes were generated using a dispersion of 0.05 and a mean of μ_c for each library in condition c. For the first 5000 genes, μA1=μA2=0 and μB1=μB2=200, while for the last 5000 genes, μA1=μA2=200 and μB1=μB2=0. This simulation design ensures that each library has a non-zero total number of reads for downstream normalization. (A constant value of 200 is only chosen for simplicity, and can be replaced with other non-zero values specific to each gene without affecting the results.) Methods in the edgeR package v3.16.5 were then used to estimate a common NB dispersion across all genes and to fit a NB GLM to the counts for each gene. QL dispersion estimates σ^g2 and σ^g(2)2 were computed for each gene as previously described, and the mean and standard error of each value was computed over all genes.

The expected value of the QL dispersion in this simulation is σg2=1 for each gene. This is because counts are exactly NB-distributed such that no QL-based modification of the mean-variance relationship is required. Estimation of the mean σ^g2 across all genes yields a value of 0.508 whereas the mean σ^g(2)2 is 1.016, with negligible standard errors in both cases. The latter is closer to the expected value for σg2, indicating that the reduced definition is correct. This difference in behaviour is due to the presence of zero counts such that d_g = 2. In contrast, d = 4 as one d.f. is provided from each condition in the experimental design. This drives the two-fold underestimation of the QL dispersion across all genes when d is used.

3.2 Assessing type I error control

3.3 Care is required when the reduced residual d.f. is zero

Some additional care is required when dealing with genes where d_g = 0. These genes provide no information on the variability of the counts and have undefined σ^g(2)2. This means that they cannot be used to estimate σ0(2)2 or d0(2). However, the null hypothesis can still be tested for these genes. As d_g = 0, the expression for the shrunken dispersion σ~g(2)2 simply collapses to σ0(2)2. This means that Fg(2) and a p-value can be computed.

To demonstrate, consider a one-way layout containing three conditions A, B₁ and B₂ for 10,000 genes. Condition A contains two replicates whereas B₁ and B₂ contain one replicate each. For the first 50% of genes, μA=0 and μB1=μB2=200. This simulates non-DE genes between B₁ and B₂ for which there are no residual d.f.’s. For the remaining genes, μA=200 and μB1=μB2=0. This setup ensures that all library sizes are non-zero and that σ0(2)2 and d0(2) can be estimated from genes with defined QL dispersions. Counts for each library were sampled from a NB distribution with the specified mean. The QL framework was applied using the reduced residual d.f., and the null hypothesis μB1=μB2 was tested for the first 50% of genes. The observed type I error rate was computed at a nominal threshold of 1%. This was repeated for 20 simulation iterations, yielding a mean error rate of 0.91% with a standard error of 0.09. Similar results were obtained at thresholds of 0.1% (0.14 ± 0.03%) and 10% (9.72 ± 0.19%). Error rates close to the nominal threshold indicate that the QL framework remains valid when d_g = 0.

3.4 Effect on log-transformed counts in linear models

3.4.1 Overview

An alternative approach to analyzing count data involves fitting a linear model to log-CPM values (Law et al., 2014). Consider an example data set where all libraries have 10⁶ reads. Any zeroes will be converted into a lower bound for the log-CPM, e.g. log2(0+0.5)=−1 where the 0.5 represents a continuity correction. If any fitted value of the linear model is equal to the lower bound, the corresponding log-CPM must also be equal to the lower bound, as no smaller log-CPM can exist. Thus, no d.f. will be provided by those libraries, which means that the canonical definition of the residual d.f. may be incorrect. This can be remedied for each gene by removing all libraries with fitted values equal to the lower bound. Linear modelling and variance estimation will then be performed using the reduced d.f. as described above for GLMs.

To demonstrate, NB-distributed counts for a four-condition design were simulated as described in Section 3.1. The voom function was applied to log-transform the counts and to compute precision weights from a fitted mean-variance trend. A linear model was fitted to the log-CPM values with the precision weights using limma v3.30.9, and an EB strategy was applied to robustly shrink the variances (Phipson et al., 2016). The mean variance estimate across all genes was recorded. The F-test was applied to compute p-values against the null hypothesis, i.e. μB1=μB2 for the first 5000 genes, and μA1=μA2 for the last 5000 genes (again, this distinction avoids a trivial result when means and counts are fixed at zero for some libraries). As the null is always true, the distribution of p-values across all genes should be uniform. The analysis was repeated after removing libraries with μc=0 from the dataset and design matrix prior to fitting, i.e. all libraries in A₁ and A₂ for the first 5000 genes, or in B₁ and B₂ for the last 5000 genes.

The mean sample variance across all genes was estimated as 0.065 and 0.132 using the canonical and reduced definitions, respectively. This fold difference is consistent with the differences in the residual d.f.’s, i.e. d = 4 and d_g = 2. Type I error control was subsequently lost with the canonical method (Figure 1), consistent with the liberalness observed in simulations for NB-based GLMs. Uniformity of the p-values was restored by removing the libraries in the offending conditions. This ensured that the reduced residual d.f. was used during modelling. In practice, removal of offending libraries is complicated by the presence of variable library sizes, resulting in a variable lower bound that is difficult to define for any given library. For simplicity, such cases will not be considered here.

Figure 1

Histograms of p-values generated by voom and limma on simulated count data, before and after removal of libraries with true means of zero.

4.1 Overview

To determine the relevance of the simulation results, analyses with the canonical and reduced residual d.f.’s were compared on a real RNA-seq dataset. Data was generated from a Pax5 knock-out experiment in pro-B cells by Drs. Rhys Allan and Steve Nutt in the Molecular Immunology division of the Walter and Eliza Hall Institute of Medical Research. This dataset contains two replicate libraries for the knock-out (KO) condition and one replicate for the wild-type (WT) condition.

4.2 Processing steps for real RNA-seq data

Paired-end sequencing was performed by the Beijing Genomics Institute on an Illumina HiSeq 2000. Each library consisted of approximately 9 million pairs of 90 bp reads. Reads were aligned to the mm10 build of the mouse genome using subread v1.4.6 (Liao et al., 2013) in paired-end mode with unique mapping and tie splitting by Hamming distance. Read pairs were summarized into gene counts using the featureCounts function (Liao et al., 2014) in the Rsubread package v1.18.1. The number of fragments mapped to exons was counted for each gene in the NCBI mouse build 38 annotation. Note that each read pair corresponds to a single cDNA fragment and is counted no more than once. Reads with MAPQ scores below 10 were ignored to avoid non-uniquely/poorly mapped reads. Approximately 64% of read pairs in each library were counted into genes.

Genes were filtered to remove those with a count sum across all libraries below 10. This removes lowly expressed genes that are unlikely to be DE and is roughly independent of the p-values when library sizes are similar. A differential expression analysis was performed using edgeR to compare the WT and KO conditions. Briefly, a trended NB dispersion was estimated and used for GLM fitting. Offsets were defined from the log-transformed effective library sizes, after using TMM normalization (Robinson and Oshlack, 2010) to remove composition biases. Raw QL dispersions were estimated with the canonical residual d.f., as previously described. A second mean-dependent trend was fitted to the raw QL estimates, and EB shrinkage was performed towards this trend (Lund et al., 2012). Trend fitting is necessary here to empirically model non-NB mean-variance relationships, but was not required for the simulations where counts were exactly NB-distributed for simplicity. Finally, the QL F-test was used to compute a p-value for each gene.

Genes were considered to be significantly DE if the false discovery rate (FDR) was <5% after applying the Benjamini and Hochberg method to all p-values. This was repeated using the reduced residual d.f. and its associated statistics.

4.3 Differences between definitions are present in real data

The QL framework was applied with the reduced and canonical residual d.f.’s to identify DE genes between WT and KO cells. For the canonical analysis, a cluster of low dispersion estimates was observed (Figure 2A). These correspond to 139 genes that have zero counts in the two KO samples. Here, there are no residual d.f. so the deviance is always zero. As d = 1, σ^g2 is incorrectly defined as zero for these genes. In comparison, σ^g(2)2 is undefined when d_g = 0 and should not participate in EB shrinkage. Removal of the corresponding genes increased the estimated prior d.f. from 2.71 to 5.17 and increased the estimated QL dispersions (Figure 2B). This is consistent with inflation of heteroskedasticity and underestimation of the QL dispersion when d is used instead of d_g in the simulations.

Figure 2

QL dispersion estimates for the Pax5 dataset. (A) Histogram of log-transformed gene-wise estimates using the canonical d.f. A small prior was added to avoid undefined values. The black bar marks genes with a reduced residual d.f. of zero. (B) Loess-fitted trend in the estimates for all genes, using the average abundance as the covariate. Results are shown before (full) and after (dashed) removing genes with no residual d.f.

Hypothesis testing was also performed to identify DE genes using the reduced and canonical methods. At a FDR of 5%, the number of DE genes increased from 1926 to 2809 when d_g was used instead of d. Reduced detection with the canonical d.f. is consistent with the conservativeness observed in Table 1, where the drop in d₀ outweighs the liberalness caused by the underestimation of σg2 and σ02. Note that the 139 genes with undefined σ^g2 were identified as DE in both the reduced and canonical methods. This is not unexpected, given that the counts are zero in the KO samples and equal to or greater than 10 in the WT sample. However, it does illustrate that incorrect estimation of the dispersion for some genes will affect inferences for the entire data set. Almost 900 DE genes were uniquely detected when the reduced method was used instead of the canonical method, despite the fact that all of them have correctly defined d = d_g = 1. This unique set includes genes such as Tnfaip3 (adjusted p-value of 0.035), Id2 (0.023) and Nfatc2 (0.021) that have been implicated in B-cell differentiation (Peng et al., 2001; Becker-Herman et al., 2002; Chu et al., 2011). Such genes would not have been identified as putative downstream targets of Pax5 if the canonical method was used.