Model and inference

We denote by k_ijl the number of reads overlapping counting bin l of gene i in sample j. We interpret k_ijl as a realization of a random variable K_ijl. The number of samples is denoted by m, i.e., j = 1, …, m.

We write μ_ijl for the expected value of the concentration of cDNA fragments contributing to counting bin l of gene i, and relate the expected read count E(K_ijl) to μ_ijl via the size factor s_j, which accounts for the depth that sample j was sequenced: E(k_ijl) = s_jμ_ijl. Note that s_j depends only on j, i.e., the differences in sequencing depth are assumed to cause a linear scaling of the read counts. We estimate the size factors with the same method as in DESeq (Anders and Huber 2010; for details, please see Supplemental Note S.1).

A generalized linear model

We use generalized linear models (GLMs) (McCullagh and Nelder 1989) to model read counts. Specifically, we assume K_ijl to follow a negative binomial (NB) distribution:

where α_il is the dispersion parameter (a measure of the distribution’s spread; see below) for counting bin (i, l), and the mean is predicted via a log-linear model as

The negative binomial distribution in Equation 1 has been useful in many applications of count data regression (Cameron and Trivedi 1998). It can be seen as a generalization of the Poisson distribution: For a Poisson distribution, the variance v is equal to the mean μ, while for the negative binomial, the variance is v = μ + αμ², with the dispersion α describing the squared coefficient of variation in excess of the Poisson case. Lu et al. (2005) and Robinson and Smyth (2007) motivated the use of the NB distribution for SAGE and RNA-seq data; we briefly summarize their argument in Supplemental Note S.2.

We fit one model for each gene i, i.e., the index i in Equation 2 is fixed. The linear predictor μ_ijl is decomposed into four factors as follows: represents the baseline expression strength of gene i. is (up to an additive constant) the logarithm of the expected fraction of the reads mapped to gene i that overlap with counting bin l. is the logarithm of the fold change in overall expression of gene i under condition ρ_j (the experimental condition of sample j). Finally, is the effect that condition ρ_j has on the fraction of reads falling into bin l.

To make the model identifiable, constraints on the coefficients are needed; see Supplemental Note S.3.

Of interest in this model are the effects and . If one of the is different from zero, this indicates that the counting bin it refers to is differentially used. A value of different from zero indicates an overall differential abundance that equally affects all counting bins, i.e., overall differential expression of the gene. Before we describe the analysis-of-deviance (ANODEV) procedure to test for these effects, we need to discuss the aspect of dispersion.

Parameter fitting

For a fixed choice of the dispersion parameter, the NB distribution is a member of the exponential family with respect to the mean. Hence, the iteratively reweighted least square (IRLS) algorithm, which is commonly used to fit GLMs (McCullagh and Nelder 1989), allows fitting of the model (1, 2) if the dispersion α_il is given.

Ordinary maximum likelihood estimation of the dispersion is not suitable, because it has a strong negative bias when the number of samples is small. The bias is caused by not accounting for the loss of degrees of freedom that arises when estimating the coefficients. Robinson and Smyth (2008) reviewed alternatives and derived an estimator based on the work of Cox and Reid (1987) and Smyth and Verbyla (1996). Cox and Reid suggested modifying the profile likelihood for the parameter of interest (here, the dispersion) by dividing out a term containing the Fisher information for the other parameters as an approximation to conditioning on the profiled-out parameters. This works if the parameter of interest is approximately independent from the other parameters with respect to Fisher information, which is the case for the NB likelihood with respect to its parameters mean and dispersion. However, calculating the Cox-Reid correction term for dispersion estimation in GLMs is not straightforward. The (to our knowledge) best method has been proposed by McCarthy et al. (2012). The authors have been using it in their edgeR package (Robinson et al. 2010a) since September 2010 (version 1.7.18). We make use of this approach to estimate the dispersion for each counting bin; details are provided in Supplemental Note S.4.

Two noise components

It is helpful to decompose the extra-Poisson variation of K_ijl into two components: variability in gene expression and variability in exon usage. If the expression of a gene i (i.e., the total number of transcripts) in sample j differs from the expected value for experimental condition ρ_j, the values μ_ijl for all of the counting bins l of gene i will deviate from the values expected for condition ρ_j by the same factor. We denote this the variability in gene expression. By variability in exon usage, we refer to variability in the usage of particular exons or counting bins. The dispersion parameter α_il in Equation 1 with respect to the model of Equation 2 contains both of these parts. However, if we replace Equation 2 with

i.e., instead of fitting one parameter for the effect of each condition ρ on the expression, we fit one parameter for each sample j, the gene expression variability is absorbed by the model parameters and we are only left with the exon usage variability. Hence, we use model (3) to increase power in our test for differential exon usage. This is possible because we test for an interaction effect. If the aim were to test for a main effect such as differential expression, dispersion estimation would need to be based on model (2).

We fit the model (3) for each gene i separately and use the Cox-Reid dispersion estimator of McCarthy et al. (2012), as described above, to obtain a dispersion value for each counting bin l in the gene.

Information sharing across genes

If only a few replicates are available, as is often the case in high-throughput sequencing experiments, we need to be able to deal with the fact that the dispersion estimator for a single counting bin has a large sampling variance. A commonly used solution is to share information across estimators (Tusher et al. 2001; Lönnstedt and Speed 2002). We noted that there is a systematic trend of dispersions as a function of the mean, and consider the relationship

This relation appears to fit many data sets we have encountered in practice. (See also Di et al. 2011 for a comparison of approaches to model mean-variance relations in RNA-seq data.) To obtain the coefficients a₀ and a₁, we regress the dispersion estimates for all counting bins from all genes on their average normalized count values with a gamma-family GLM. To ensure robustness of the fit, we iteratively leave out bins with large residuals until convergence is achieved (Huber 1981).

Figure 2 shows a scatterplot of dispersion estimates against average normalized count values , together with the fit α(μ). For many counting bins, the difference between the sample estimate and the fitted value is compatible with a χ² sampling distribution (indicated by the dashed lines). Nevertheless, there are sufficiently many bins with a sample estimate so much larger than the fitted value that it would not be justified to only rely on the fitted values. Hence, for the ANODEV (see below), we use as dispersion value α_il the maximum of the per-bin estimate and the fitted value . On average, this overestimates the true dispersion and costs power, but we consider this preferable to using either only the fitted values or the sample estimates, both of which carry the risk of producing many undesirable false positives. More sophisticated alternatives for this step, which usefully interpolate between the two extremes, and perhaps incorporate further covariates besides μ, might become available in the future.

Open in a separate window

Figure 2.

Dependence of dispersion on the mean. Each dot corresponds to one counting bin in the data of Brooks et al. (2010) (discussed in detail in the Results section); (x-axis) normalized count, averaged over all samples; (y-axis) estimate of the dispersion. The bars at the bottom denote dispersion values outside the plotting range (in particular, those cases in which the sample dispersion is close to zero). (Solid red line) The regression line; (dashed lines) the 1-, 5-, 95-, and 99-percentiles of the χ² distribution with 4 degrees of freedom scaled such that it has the fitted mean.

Analysis of deviance

We test for each counting bin whether it is differentially used between conditions. More precisely, we test against the null hypothesis that the fraction of reads overlapping with a counting bin l, of all the reads overlapping with the gene, does not change between conditions. To this end, we fit for each gene i a reduced model with no counting-bin–condition interaction:

and, separately for each bin l′ of gene i, a model with an interaction coefficient for only this bin, but as in Equation 5, main effects for all bins l,

Here, δ_ll′ is the Kronecker delta symbol, which is 1 if l = l′ and 0 otherwise. We compute the likelihood of these models using the dispersion values α_il as estimated from model (3), with the information-sharing scheme presented earlier. Comparing the fit (6) for counting bin l′ of gene i with the fit (5) for gene i, we get an analysis-of-deviance P-value p_il′ for each counting bin by means of a χ² likelihood-ratio test. Note that we test against the null hypothesis that none of the conditions influences exon usage, and hence, if there are more than two different conditions ρ, we aim to reject the null hypothesis already if any one of the conditions causes differential exon usage.

Differential exon usage, as treated here, cannot be distinguished from overall differential expression of a gene if the gene only consists of a single counting bin or if all but one of its counting bins have zero counts. Hence, we mark all counting bins with zero counts in all samples, and all bins in genes with less than two nonzero bins, as not testable. Furthermore, we skip counting bins with a count sum across all samples below a threshold chosen low enough that a significant result would be unlikely, to speed up computation. Such filtering can also improve power (see Bourgon et al. 2010).

Note that we perform one test for each counting bin, always fitting an interaction coefficient only for the single bin l′ under test. Therefore, it is valid that a read that overlaps with several exons is counted for each of these exons: In each test, for the purpose of estimating and testing the interaction coefficient, any given read is only considered at most once.

Additional covariates

The flexibility of GLMs makes it easy to account for further covariates. For example, if in addition to the experimental condition ρ_j we wish to account for a further covariate τ_j, we extend model (3) as follows:

When testing for differential exon usage, the extra term is added to both the reduced model (5) and the full model (6).

An example is provided in the next section with Equation 9.

Visualization

The DEXSeq package offers facilities to visualize data and fits. An example is shown in Figure 3, using the data discussed in the next section. Data and results for a gene are presented in three panels. The top panel depicts the fitted values from the GLM fit. For this plot, the data are fitted according to model (2), with the y coordinates showing the exponentiated sums:

The tildes indicate that a decomposition of the linear predictors has been used that separates the effects of expression and isoform regulation, as described in Supplemental Note S.3.

Open in a separate window

Figure 3.

The treatment of knocking down the splicing factor pasilla affects the fourth exon (counting bin E004) of the gene Ten-m (CG5723). (Top panel) Fitted values according to the linear model; (middle panel) normalized counts for each sample; (bottom panel) flattened gene model. (Red) Data for knockdown samples; (blue) control.

For genes with differential overall expression, it can be difficult to see the evidence for differential exon usage in a plot based on Equation 8. For these cases, the software offers the option to average over the expression effects. Supplemental Figure S1 shows this for the pasilla gene.

Analysis of the chimpanzee data of Brawand et al.

While the preceding application was on a controlled experiment with a cell culture under sharp treatment, in this section, we analyze data from an observational study with complex subject-to-subject variation (Brawand et al. 2011). This data set includes RNA-seq from prefrontal cortex samples from six chimpanzees and cerebellum samples from two further chimpanzees. We used DEXSeq to test for exon usage differences between these two brain tissue types.

We aligned the RNA-seq reads (GEO accessions GSM752664–GSM752671) from these samples to the chimpanzee genome (CHIMP2.1.4 from Ensembl 64) using GSNAP 2012-01-11 (Wu and Nacu 2010). Prior to alignment, we trimmed all reads to a common length of 76 nt, single-ended. The trimming was necessary to make the data comparable across samples; DEXSeq itself has no length limitation and can deal with any read length.

At 10% FDR, DEXSeq found significant differential exon usage for 866 counting bins in 650 genes. The result table, with plots for all genes with significant differential exon usage, can be found at http://www-huber.embl.de/pub/DEXSeq/chimp/testForDEU.html.

Exploration of this hit list reveals interesting differences between the tissues. For example, one of the top hits, the gene PRKCZ (protein kinase C zeta; ENSPTRG00000000042) expresses its first four exons only in cerebellum but not in the prefrontal cortex (Supplemental Fig. S4). Inspecting the Pfam (Finn et al. 2010) and SMART (Letunic et al. 2012) databases of protein domains reveals that these four exons encode the heterodimerization domain PB1. This suggests the hypothesis that the gene product loses its ability to bind to its partner protein in the prefrontal cortex. Indeed, these two isoforms are well studied (for review, see Hirai and Chida 2003). The long isoform protein product PKZζ is widely expressed and is activated by a second messenger, PARD6A, which removes the protein’s auto-inhibition by binding to the PB1 domain. The truncated protein, denoted PKMζ, is specific to the brain and, due to the lack of the PB1 domain, constitutively active. It plays a major role in long-term potentiation and memory formation. In this context, it is noteworthy that, as our analysis shows, its expression is confined to certain brain regions.

Another example is provided by the gene PLCH2 (phospholipase C eta 2; ENSPTRG00000000051), for which DEXSeq indicated differential usage of counting bin E011 (fourth exon). According to SMART and Pfam, this exon contains an EF hand, a calcium binding helix–loop–helix motif. Here, we are not aware of prior work on the isoform(s) lacking this exon. We can speculate that the shorter isoform’s activity might no longer depend on calcium concentration, on which PLCH2’s enzymatic activity normally depends strongly (Nakahara et al. 2005). Furthermore, Zhou et al. (2008) studied the activation of PLCH2 by Gβγ complexes and found that the EF hand domain of PLCH2 is required for this interaction. Another hypothesis might hence be that the observed tissue-specific usage of the fourth exon serves to modulate the regulation of PLCH2 by G-proteins.

For the gene ENSPTRG00000000130, DEXSeq reports increased usage of the second exon in the cerebellum and of the second-to-last exon in the prefrontal cortex. This gene codes for precortistatin, a protein that gets cleaved to give rise to the neuropeptide cortistatin, which (in human) comprises the last 17 amino acids of the full protein’s C terminus (de Lecea et al. 1997), which are contained in the last exon. While the overall expression differences seen in the data agree with the known main location of cortistatin—the cortex, the observed differential exon usage is intriguing and more difficult to interpret: The affected parts of the protein are not part of the final product. Nevertheless, the presence or absence of these parts could affect the efficiency of the cleavage process or the stability of the mRNA, to coregulate the tissue-specific expression.

These three examples illustrate how a DEXSeq analysis can serve as a starting point for hypothesis formation. We picked these three genes by inspecting the first 10 genes with significant differential exon usage, as sorted by numerical Ensembl gene ID (not by P-value); that is, in essence we inspected a mere 10 randomly chosen hits. The richness of the biology seen indicates that many novel insights into gene function and regulation may be expected from the analysis of tissue-specific exon usage patterns.

Analyses based on Fisher’s test

To test for differential isoform regulation, Wang et al. (2008) and Brooks et al. (2010) used 2 × 2 contingency tables and a Fisher’s exact test. In this approach, the contingency table’s rows corresponded to control and treatment, the cells in one column contained the numbers of reads supporting inclusion of an exon (i.e., reads overlapping the exon), and the cells in the other column gave the numbers of reads supporting exclusion (e.g., in the case of cassette exons, reads straddling the exon). In the study of Wang et al. (2008), each row corresponded to a single sample, while Brooks et al. (2010) summed up the number of reads from their replicates. The MISO method (Katz et al. 2010) proposed a different way of setting up the contingency table. In all cases, the contingency tables did not contain information on sample-to-sample variability (Baggerly et al. 2003), and, therefore, one should expect the results to contain an inflated number of false positives.

As an example, Supplemental Figure S3 shows the gene Lk6, for which Brooks et al. (2010) reported differential use of its alternative first exons, while our analysis did not call a significant differential use. Clearly, the average expression strength of exon E002 is different between the conditions. However, examining the counts from the individual biological replicates reveals that the variance within the treatment groups is large compared with this difference, and hence, the data do not support the claim of a significant effect of the treatment.

Heterogeneity of dispersions

In our model, we allow the counting bins of a gene to have different dispersion values. The gene RpS14b (Fig. 5) exhibits very different variability for its three exons and thus illustrates the need for this modeling choice.

Open in a separate window

Figure 5.

Ribosomal protein gene RpS14b (from the Brooks and coworkers data) is shown here as an example for a gene with heterogeneous dispersion. The first exon has zero count in the paired-end samples untreated 2, in the single-end sample treated 2, and in the paired-end sample treated 3, and large nonzero counts in the four other samples. Colors are as in Figure 3.

The first exon also illustrates the value of replicates and the importance of making use of their information. This exon had between 252 and 416 (normalized) counts in four of the samples and no counts in three. However, this difference cannot be attributed to the treatment because both the control and the treatment group contained samples with zero counts as well as samples with several hundreds of counts. Hence, the reason for the difference in read counts for this exon cannot be the knockdown of pasilla and is likely some other difference between the samples’ treatment that was not under the experimenters’ control.

If one just adds up or averages the samples in a treatment group, as done in the contingency table method, one would only see a sizeable difference, as in the upper panel of the figure, and might call a significant effect. It is also crucial that the test for differential exon usage does not rely on the fitted dispersion (solid line in Fig. 2) only, because the effect size would seem significant if one did not take note that the actual observed within-group variance is so much larger that the fitted value is implausible. The maximum rule discussed in the section on information sharing ensures this.

Comparison with Cuffdiff

Cufflinks (Trapnell et al. 2010) is a tool to infer gene models from RNA-seq data and to quantify the abundance of transcript isoforms in an RNA-seq sample. In addition to this, the Cuffdiff module allows testing for differences in isoform abundance. Cuffdiff, as described in Trapnell et al. (2010), compares a single sample with another one and does not attempt to account for sample-to-sample variability. The latter is also true for the version described by Roberts et al. (2011), which allows processing of replicate samples, but uses this for the assessment only of bias, not of variability. Hence, the same drawbacks may be expected as discussed earlier for the Fisher-test-based methods. More recently, starting with version 1.0.0, Cufflinks attempts to assess overdispersion and account for it.

We compared the three knockdown samples of the Brooks and coworkers data set against the four control samples with version 1.3.0 of Cuffdiff. With nominal FDR control at 10%, Cuffdiff reported differential splicing for only 50 genes, and thus showed less power than our approach.

To test the control of false-positive rates, we made use of the fact that there were four replicates for the untreated condition. We formed one group from samples 1 and 3 and another group from samples 2 and 4. We tasked both DEXSeq and Cuffdiff with comparing between the two groups at a nominal FDR of 10%. Because this is a comparison between replicates, ideally no significant calls should be made. Note that each group contained one single-end and one paired-end sample, i.e., the blocking caused by the library type was balanced between the groups. In this mock comparison, DEXSeq found eight genes significant, compared with 159 in the comparison of treatment versus control. Surprisingly, Cufflinks found 639 genes in the mock comparison, many more than the 37 genes that it found in the proper between-groups comparison. Supplemental Note S.6 describes further tests, which confirmed Cufflinks’ difficulty with providing type I error control in this data set.

We also performed the same type of comparison on a data set with quite different characteristics and experimental design, the chimpanzee data of Brawand et al. (2011). In a comparison of the six chimpanzee prefrontal cortex (PFC) samples with the two cerebellum samples, Cuffdiff 1.3.0 reported 114 genes at 10% FDR, again showing less power then DEXSeq (650 genes; see above).

We then used the five PFC samples from male chimpanzees to assess type I error rates. Both tools were tasked to compare any combination of two samples versus two other samples. DEXSeq in each case found substantially fewer genes in these mock comparisons than in the proper comparison (with one exception, always less than 1/65). Cuffdiff, however, always found more than twice as many genes in the mock comparisons than in the true one. For details, see Supplemental Note S.6. Also see Supplemental Note II, which contains the exact commands used for all computations performed for this article.

Comparing exon or isoform usage

The interpretation of the results of our method is straightforward when a single exon of a gene with many exons is called differentially used. However, if many exons within a gene are affected, the interpretation is more complex. For instance, consider a gene with two isoforms, a long one with n exons and a short one consisting of only the first n/2 exons. If an experimental condition increased the number of long transcripts at the expense of the short ones, without changing the total number, one might expect an analysis to indicate differential usage for the last n/2 exons. However, our method cannot distinguish this situation from one in which the gene is overall down-regulated, while the first n/2 exons are more strongly used.

Hence, if differential exon usage is detected within a gene, we can safely conclude that this gene is affected by alternative isoform regulation. However, the test’s output with regard to which of the counting bins are affected can be unreliable if the isoform regulation affects a large fraction of the exons. In practice, the assignment to counting bins is reliable as long as only a small fraction of counting bins in the gene is called significant.

Methods that attempt to estimate not just the abundance of exons but of isoforms, such as the method of Jiang and Wong (2009), Cufflinks (Trapnell et al. 2010) and MMSeq (Turro et al. 2011), may be able to circumvent this issue. Of these, only Cufflinks/Cuffdiff offers the functionality of comparing between samples. We commented on Cuffdiff in the preceding section. (Note added after revision: Recently, Glaus et al. 2012 published BitSeq, another method to identify differential expression of isoforms.)

Apart from the lack of tools for inferring differential expression at the transcript level, there can be concrete advantages in per-exon analysis. If, for example, several transcripts have most exons in common and differ by only a few exons, their abundance estimates will contain substantial correlated uncertainties that reduce the power for inference of differential expression. The remedy would be to disregard the reads that inform about the shared parts of the transcripts and to focus on those reads in which they differ. Hence, an exon-centric analysis might be a crucial component even of a transcript-level method.

In addition, it is not clear that inference about transcripts is always more useful for biological interpretation than inference at the per-exon level. After all, we have knowledge about the functional differences of multiple translated isoforms of a gene for only a small number of proteins. If currently a researcher finds that a gene of interest expresses different transcripts in different conditions, her further analysis will typically start with assessing the difference between the two transcripts. She might find, for example, that they differ in the presence of certain exons and ask which regulatory signals or functional domains these exons may contain. Therefore, we expect that a method such as ours that pinpoints the location of the differences by focusing on specific exons will be valuable for biological interpretation, and sometimes perhaps more valuable than a transcript-centric approach. This expectation is supported by the three examples discussed in the analysis of the chimpanzee data. The next step will be to leverage, in a systematic and automated way, databases with annotation for parts of gene products, e.g., information on protein domains provided by resources such as Pfam (Finn et al. 2010), SMART (Letunic et al. 2012), and PROSITE (Sigrist et al. 2010), or predicted miRNA target sites.

Junction reads

Junction reads are reads whose genomic alignment contains a gap because they start in one exon, end in another exon, and “jump” over the intron in between and possibly over skipped exons. In DEXSeq, such reads are counted for each counting bin with which they overlap, i.e., they appear multiple times in the count table. However, because we test for each exon separately, this does not affect the validity of the test.

Junction reads contain additional information that is especially valuable when inferring gene models and the positions of splice junctions. Unless one works with a very well annotated model system, this information should be used when defining the counting bins, by parsing the spliced alignments with appropriate tools.

Furthermore, junction reads give evidence for connections between counting bins and thus are crucial for isoform deconvolution tools such as Cufflinks and MMSeq. For our exon-by-exon test, however, leveraging this information is not essential, and also not straightforward. In the method presented, we essentially consider for each sample the ratio of the number of reads overlapping with an exon to the number of reads falling onto the whole gene. Alternatively, one could consider the ratio of the number of reads skipping over the exon under consideration to the total count. We anticipate that the latter would offer a moderate increase in power in cases in which the counting bin is much shorter than the typical read length. It may be an interesting future extension to the DEXSeq method to switch to this scheme for bins that are short compared with the read length.

Implementation

We implemented DEXSeq as a package for the statistical programming language R (R Development Core Team 2009) and have made it available as open source software via the Bioconductor project (Gentleman et al. 2004). See the Bioconductor web page for downloading instructions. DEXSeq can be used on MacOS, Linux, and Windows.

For the preparation steps, namely, the “flattening” of the transcriptome annotation to counting bins and the counting of the reads overlapping each counting bin, two Python scripts are provided, which are built on the HTSeq framework (Anders 2011). The first script takes a GTF file with gene models and transforms it into a GFF file listing counting bins, and the second takes such a GFF file and an alignment file in the SAM format and produces a list of counts. The R package is used to read these counts, estimate the size factors and dispersions, fit the dispersion-mean relation, and test for differential exon usage. After the analysis has been performed, all the results are available, together with the input data, in an object derived from the ExpressionSet class, Bioconductor’s standard container type for data from high-throughput assays. The results provided include for each counting bin the following data: the conditional-maximum-likelihood estimate for the dispersion, the dispersion value actually used in the test (which may be different, due to the information sharing across genes), the P-value from the test for differential exon usage, the Benjamini-Hochberg-adjusted P-value, and the fit coefficients describing the fitted log₂ fold change between treatment controls (or, if there are more than two conditions, for pairs of conditions as chosen by the user). Other R or Bioconductor functionality can be used for downstream analyses of these results. If required, the other coefficients as described in Supplemental Note S.3 are also available.

Furthermore, DEXSeq can create a set of HTML pages that contains the results of the tests, and, for each gene, plots like Figures 3 and ​and55 and Supplemental Figures S1 and S3. The HTML output allows interactive browsing of the results and facilitates sharing of the results with colleagues by uploading the files to a web server.

The DEXSeq package provides functions on different levels. In the simplest case, a single function is called that runs all the steps of a standard analysis. To give experienced users the possibility to interfere with the workflow, functions are also provided to run each step separately, to run some steps only for single genes, and to inspect intermediate and final results.

The use of the package is explained in the vignette (a manual with a worked example) and documentation pages for all functions.

Because the DEXSeq method relies on fitting GLMs of the NB family, a reliable IRLS fitting function is required. We use the function nbglm.fit (McCarthy et al. 2012) from the statmod package, which offers better performance and convergence than older implementations.

Fitting GLMs for many genes and counting bins is a computationally expensive process. When running on a single core of a current desktop computer, the analysis of the Brooks and coworkers data presented here takes several hours. However, the method lends itself easily to parallelization: We use the multicore package (Urbanek 2011) to distribute the computation on several CPU cores.

The complete workflow used to perform all calculations for this article is documented in Supplement II.