Gene ontology analysis for RNA-seq: accounting for selection bias
© Young et al.; licensee BioMed Central Ltd. 2010
Received: 1 December 2009
Accepted: 4 February 2010
Published: 4 February 2010
We present GOseq, an application for performing Gene Ontology (GO) analysis on RNA-seq data. GO analysis is widely used to reduce complexity and highlight biological processes in genome-wide expression studies, but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts. Application of GOseq to a prostate cancer data set shows that GOseq dramatically changes the results, highlighting categories more consistent with the known biology.
Next generation sequencing of RNA (RNA-seq) gives unprecedented detail about the transcriptional landscape of an organism. Not only is it possible to accurately measure expression levels of transcripts in a sample , but this new technology promises to deliver a range of additional benefits, such as the investigation of alternative splicing , allele specific expression  and RNA editing . However, in order to accurately make use of the data, it is vital that analysis techniques are developed to take into account the technical features of RNA-seq output. As many of the specific technical properties of RNA-seq data are not present in previous technologies such as microarrays, naive application of the same analysis methodologies, developed for these older technologies, may lead to bias in the results.
In RNA-seq experiments the expression level of a transcript is estimated from the number of reads that map to that transcript. In many applications, the expected read count for a transcript is proportional to the gene's expression level multiplied by its transcript length. Therefore, even when two transcripts are expressed at the same level, differences in length will yield differing numbers of total reads. One consequence of this is that longer transcripts give more statistical power for detecting differential expression between samples . Similarly, more highly expressed transcripts have a greater number of reads and greater power to detect differential expression. Hence, long or highly expressed transcripts are more likely to be detected as differentially expressed compared with their short and/or lowly expressed counterparts. The fact that statistical power increases with the number of reads is an unavoidable property of count data, which cannot be removed by normalization or re-scaling. Consequently, it is unsurprising that this selection bias has been shown to exist in a range of different experiments performed using different analysis methods, experimental designs and sequencing platforms . When performing systems biology analyses, failure to account for this effect will lead to biased results.
One simple, but extremely widely used, systems biology technique for highlighting biological processes is gene category over-representation analysis. In order to perform this analysis, genes are grouped into categories by some common biological property and then tested to find categories that are over represented amongst the differentially expressed genes. Gene Ontology (GO) categories are commonly used in this technique and there are many tools available for performing GO analysis - for example, EasyGO , GOminer , GOstat , GOToolBox , topGO , GSEA , DAVID  (see supplementary data 1 in Huang da et al.  for a more complete list). Although these tools have some differences in methodology , they all rely on similar underlying assumptions about the distribution of differentially expressed (DE) genes. Specifically, it is assumed that all genes are independent and equally likely to be selected as DE, under the null hypothesis. It is this assumption that makes the standard GO analysis methods unsuitable for RNA-seq data due to the bias in detecting differential expression for genes with a high number of reads. It follows then that when using a standard analysis, any category containing a preponderance of long genes will be more likely to show up as over-represented than a category with genes of average length. Similarly, categories with highly expressed genes are also more likely to be found as over-represented than categories of lowly expressed genes. Having identified this issue, it is possible to compensate for the effect of selection bias in the statistical test of a category's enrichment among differentially expressed genes.
This paper will be concerned with developing a statistical methodology that enables the application of GO analysis to RNA-seq data by properly incorporating the effect of selection bias. Using published RNA-seq data, we will show that accounting for this effect leads to significantly different results, which agree much better with previous microarray studies and the known biology than the results of an uncorrected analysis.
Because it is so widely used, we choose to focus our category testing on GO analysis. The techniques developed here apply more generally, however, to any analysis that looks for over-representation of some category of interest amongst differentially expressed genes. For example, alternative analyses might look for over-representation of KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, gene sets in the Molecular Signatures Database , or gene lists derived from earlier microarray experiments.
To illustrate the methodology, the GOseq technique was applied to an experiment examining the effects of androgen stimulation on a human prostate cancer cell line, LNCaP . This published data set includes more than 17 million short cDNA reads obtained for both the treated and untreated cell line and sequenced on Illumina's 1G genome analyzer. These reads were re-mapped back to the reference genome and the number of reads per gene was recorded (see Additional file 1 for details). Examples of the methodology throughout the paper are made using this data set and we show that the GOseq method makes a substantial difference to the results of the GO analysis, which are more consistent with prior knowledge of the system.
Standard methods for testing over-representation of a GO category assume that, under the null hypothesis, each gene has equal probability of being detected as DE. Under this assumption, the number of genes associated with a category that overlap with the set of DE genes follows a hypergeometric distribution. Hence the GO test can be conducted using Fisher's exact test, which uses the hypergeometric distribution, or Pearson's chi-square test, which is a computationally convenient approximation . Because of the existence of selection bias, genes with more reads are more likely to be detected as DE, violating the assumption behind the hypergeometric distribution. Therefore, these standard test distributions should not be used for GO analysis with RNA-seq data. Instead, we need a new test that corrects for selection bias.
The GOseq method
In order to correct for selection bias in category testing, we propose the following three-step methodology. First, the genes that are significantly DE between conditions are identified. The GOseq method works with any procedure for identifying DE genes. Second, the likelihood of DE as a function of transcript length is quantified. This is obtained by fitting a monotonic function to DE versus transcript length data. Finally, the DE versus length function is incorporated into the statistical test of each category's significance. This final step takes into account the lengths of the genes that make up each category.
In the second step, a probability weighting function (PWF) is estimated from the data. The PWF quantifies how the probability of a gene selected as DE changes as a function of its transcript length. To estimate this trend function, each gene is assigned a binary value (zero or one), according to whether or not the gene is found to be DE. A monotonic spline with 6 knots is then fitted to this binary data series using the transcript length of each gene as predictor (see Materials and methods). Monotonicity is imposed as the power to detect DE using statistical tests increases monotonically with an increasing number of reads. Figure 2 also shows the resulting probability weighting function for the prostate cancer data set. The PWF then forms the null hypothesis for our enrichment test.
Unlike GO analysis for microarray data, the null probability distribution does not conform to a standard distribution, precluding an analytical solution for determining the probability of a category being over-represented among DE genes. However, the P-value for each category can be computed using a resampling strategy. For the final step of the GOseq method, resampling was performed by randomly selecting a set of genes, the same size as the set of DE genes, and counting the number of genes associated with the GO category of interest. This random selection weights the chance of choosing a gene by its length or read count, from the previously fitted probability weighting function. The resampling is repeated many times and the resulting distribution of GO category membership is taken to approximate the shape of the true probability distribution. This sampling distribution allows calculation of a P-value for each GO category being over-represented in the set of DE genes while taking selection bias into account.
The Wallenius approximation
Although accurate, random sampling is very computationally intensive, particularly when a fine granularity in P-value is needed to distinguish between large numbers of categories. In order to make the calculation in the final step of the GOseq procedure computationally manageable, we implemented an alternative approximation technique, based on an extension of the hypergeometric distribution known as the Wallenius non-central hypergeometric distribution . This distribution extends the hypergeometric distribution to the case where the probability of success and failure differ. The GOseq implementation of the approximation assumes that all genes within a category have the same probability of being chosen, but this probability is different from the probability of choosing genes outside this category. The mean of the probability weightings for each gene within/outside the category is defined as the common probability of choosing a gene within/outside that category. While the Wallenius approximation is obviously a simplification, it is significantly closer to the true distribution than the standard hypergeometric distribution.
Although the accuracy lost by using the Wallenius approximation is not negligible, the gain in computational efficiency is dramatic. Furthermore, the ability to differentiate the most highly over-represented categories from one another (Additional file 1) makes the Wallenius approximation an attractive alternative, particularly when the range of the probability weighting function is moderate.
Comparisons of GOseq with the standard GO analysis
Gene Ontology categories ranked in the top 25 in the standard method but not in length bias adjusted GOseq
Average gene length in category
Zinc ion binding
Metal ion binding
Regulation of transcription, DNA-dependent
Gene Ontology categories ranked in the top 25 in length bias adjusted GOseq but not in the standard method
Average gene length in category
Post-translational protein modification
Small conjugating protein ligase activity
Regulation of protein metabolic process
Measuring GOseq's accuracy by comparison with microarray data
Transcript length bias versus read count bias
As transcript length bias is a technical effect, it is always necessary to correct for it when performing category testing on RNA-seq data. However, the bias in power to detect DE in longer genes arises from an increase in the total number of reads for each gene, where the number of reads is given by transcript length multiplied by expression level. Therefore, there may be circumstances where it is desirable to correct for the effect of expression level on power to detect DE, in addition to the contribution from transcript length, that is, total read count bias (for further discussion see Additional file 1). The GOseq method is capable of handling both types of bias. We have mainly focused on transcript length bias as it will always need to be accounted for and because the decision to account for expression level or not ultimately depends on the questions the user wishes to answer.
Top ten Gene Ontology categories using the standard method
Median read count of genes in category
Top ten Gene Ontology categories using GOseq adjusted for total read count bias
Median read count of genes in category
Integral to plasma membrane
Glucuronosyl transferase activity
Integral to membrane
Protein amino acid dephosphorylation
Positive regulation of transcription from RNA polymerase II promoter
Biological relevance of selected categories
To determine the effect of GOseq on the ability to draw biologically meaningful conclusions, the top ten GO categories using the standard hypergeometric method and the read count adjusted GOseq method were compared in the prostate cancer data set (Tables 3 and 4). We found that categories identified by GOseq are more consistent with previous studies looking at the relationship of androgens with prostate cancer. The role of androgens in prostate cancer is well supported, with androgen required in rodent prostate cancer induction models, and castration prior to puberty being protective against prostate cancer. Androgen is thought to be responsible for promotion of prostate cancer progression through enhancing the androgen regulated processes of growth and cellular activity. In normal prostate, androgen supports the secretary epithelial, which turns over at a rate of 1 to 2% of cells per day, and most prostate cancers are derived from these cells . Based on this biological knowledge the prior expectation is that there will be an increase in cellular activity, proliferation and secretion in LNCaP cells in response to androgen. Previous microarray experiments have shown that LNCaP cells retain androgen responsiveness and that most genes upregulated are involved in the production of seminal fluid .
In the standard analysis, the top ten GO categories indicate a change in intracellular genes, including nuclear and DNA binding genes (Table 3). The top ten categories using GOseq with total read count adjustment indicate significant changes at membranes and extracellular space, transcriptional upregulation, and in cell cycle genes (Table 4). The categories identified as most significant by GOseq therefore better match the known biology of androgen response in prostate cancer. The genes that are significant only with the length bias adjustment (Table 2) include four categories consistent with increased translation and protein production, vesical mediated transport consistent with secretion, and two categories related to nucleosomes consistent with increased replication. The category of small conjugating protein ligase activity is supported by the previously reported up-regulation of ubiquitin ligases UBE2C and HSPC150 .
Possibility of true biological trends in gene length
A key assumption of GOseq is that longer genes are not of biologically greater interest than shorter genes, per se. This assumption is supported by microarray data, where no systematic trend between gene length and differential expression has been observed . The authors find it hard to imagine that any genuine biological process could induce a trend in differential expression versus gene length comparable in magnitude to the technical trend that is removed by GOseq. Nevertheless, users should be aware that any biological trend in DE versus gene length will be adjusted for by GOseq.
Using the Wallenius approximation
The reduction in computation afforded by the Wallenius approximation is predicated on the assumption that the variance in the probability weighting within categories is low. Hence, the approximation will perform better for a probability weighting function with a low range of probabilities. When the range in probabilities as a function of gene length or read count is large, random sampling performs better than the Wallenius approximation, even at a relatively small number of replicates. The probability weighting function for read count bias has a larger range in probabilities compared to the PWF for length bias in the prostate cancer data set (Figure 2) and so random sampling may be more appropriate when accounting for total read count bias.
GOseq and other technologies
GOseq with total read count adjustment is relevant for other tag-based next generation expression profiling technologies, such as SAGE or CAGE, which have no transcript length bias . Even without length bias, statistical power to detect differential expression still depends on the expression level of each transcript and correcting for this bias will generally be desirable in this context.
GO analysis of microarray expression data has so far ignored the possibility of selection bias, but such bias clearly does exist. It is well known that fold changes are more precisely estimated for microarray probes at higher intensity levels, so selection bias is likely to exist as a function of intensity. Furthermore, most microarray platforms have multiple probes for some genes. Genes with more probes will have a great chance of being selected as DE, assuming the analysis is conducted probe-wise. The methodology developed here for RNA-seq could easily be adapted to GO analysis of microarray data, and would likely yield benefits in terms of biological relevance.
In order to implement the GOseq method, we developed a set of freely available R functions, which includes functions for calculating the significance of over-representation of each GO category amongst DE genes. These functions give the user the option of selecting which type of bias they wish to compensate for (transcript length bias or total read count bias). The option of using random sampling or the Wallenius approximation is also available. The ability of the user to supply their own categories for unbiased testing is also included. The software is freely available and can be downloaded from our website .
Here we have developed a statistical framework for GO analysis for use with RNA-seq data. It is mathematically indisputable that all commonly used criteria for judging DE interact with gene length and read count. This provides a well understood causal model of why length bias exists and why it needs to be accounted for. The GOseq method is able to account for such biases when performing GO analysis. We find that the new method makes a substantial difference to the categories identified as the most significant. We show that the GOseq method is able to recover well established microarray results more readily than existing methods of GO analysis of RNA-seq data. Furthermore, using an androgen treated prostate cancer data set we find that the most significant categories identified using GOseq match the known biology better than existing methods.
Materials and methods
The prostate cancer data set
The LNCap cell line was treated with androgen. Mock treated and treated cell lines were sequenced using the Illumina GA 1 . Raw 35-bp RNA-seq reads were provided and mapped to the human genome using Bowtie. Each mapped read was associated with an ENSEMBL gene. A Poisson exact test [15, 16] was used to determine differential expression between treated and mock-treated LNCap cells.
The liver versus kidney data set
Genome-wide expression was measured in liver and kidney using RNA-seq on the Illumina GA I and hybridization of the same samples to Affymetrix HG-U133 Plus 2.0 arrays. The sample preparation and data analysis was designed to maximize the similarity between the microarray and RNA-seq experiments (see Marioni et al. ). Differential expression between kidney and liver was determined using an empirical Bayes modified t-statistic on the microarray platform and P-values for DE were downloaded from their website. For the RNA-seq experiment, the data were normalized using TMM normalization  and a negative binomial exact test was used to determine DE . To test the GOseq method, we used the genes called DE from the microarray experiment to calculate the significance of over-representation of each GO category using the standard GO analysis methods. We also calculated P-values for each GO category being over-represented among genes that were DE in the RNA-seq data, using both the GOseq and hypergeometric methods. GOseq's ability to outperform the hypergeometric method, as measured by its ability to reproduce the results of the microarray GO analysis, was quantified by calculating a P-value for the difference in the two methods being due to chance. To do this, a NULL was chosen under which both methods were equally likely to correctly recover each microarray GO category, with this likelihood given by a binomial distribution.
The probability weighting function
To calculate the PWF, a cubic spline with a montonicity constraint is fitted to the binary data series where a value of 1 refers to a DE gene and 0 refers to a non-DE gene. This fit can be calculated against either the length of the gene or the read counts of a gene. We used the R function pcls in the mgcv package to generate the fit.
Calculating significance of categories
Random samples of genes are created by selecting a subset of genes from the experiment, with each gene weighted by the probability derived from the PWF. Each random sample contains the same number of genes as the set of DE genes. For each sample the number of genes with a given GO category is calculated. Many samples are generated in order to produce a null distribution from which the P-value for the significance of a category can be estimated.
To implement the Wallenius approximation, we used the BiasedUrn package in R. The 'odds' parameter is defined as the relative probability of genes within a category to the genes outside the category. The 'odds' ratio is calculated by taking the mean of the values from the PWF for each gene in the set and dividing by the mean of the values from the PWF for genes outside the set.
probability weighting function.
We thank Mark Robinson for providing insightful discussion and some R code and to Yoav Gilad for critical reading of the manuscript. We thank Gene Yeo and Michael Lovci for providing raw RNA-seq data.
- Fu X, Fu N, Guo S, Yan Z, Xu Y, Hu H, Menzel C, Chen W, Li Y, Zeng R, Khaitovich P: Estimating accuracy of RNA-Seq and microarrays with proteomics. BMC Genomics. 2009, 10: 161-PubMedPubMed CentralView ArticleGoogle Scholar
- Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB: Alternative isoform regulation in human tissue transcriptomes. Nature. 2008, 456: 470-476.PubMedPubMed CentralView ArticleGoogle Scholar
- Wang X, Sun Q, McGrath SD, Mardis ER, Soloway PD, Clark AG: Transcriptome-wide identification of novel imprinted genes in neonatal mouse brain. PloS One. 2008, 3: e3839-PubMedPubMed CentralView ArticleGoogle Scholar
- Wahlstedt H, Daniel C, Enstero M, Ohman M: Large-scale mRNA sequencing determines global regulation of RNA editing during brain development. Genome Res. 2009, 19: 978-986.PubMedPubMed CentralView ArticleGoogle Scholar
- Oshlack A, Wakefield MJ: Transcript length bias in RNA-seq data confounds systems biology. Biol Direct. 2009, 4: 14-PubMedPubMed CentralView ArticleGoogle Scholar
- Zhou X, Su Z: EasyGO: Gene Ontology-based annotation and functional enrichment analysis tool for agronomical species. BMC Genomics. 2007, 8: 246-PubMedPubMed CentralView ArticleGoogle Scholar
- Zeeberg BR, Feng W, Wang G, Wang MD, Fojo AT, Sunshine M, Narasimhan S, Kane DW, Reinhold WC, Lababidi S, Bussey KJ, Riss J, Barrett JC, Weinstein JN: GoMiner: a resource for biological interpretation of genomic and proteomic data. Genome biology. 2003, 4: R28-PubMedPubMed CentralView ArticleGoogle Scholar
- Beissbarth T, Speed TP: GOstat: find statistically overrepresented Gene Ontologies within a group of genes. Bioinformatics. 2004, 20: 1464-1465.PubMedView ArticleGoogle Scholar
- Martin D, Brun C, Remy E, Mouren P, Thieffry D, Jacq B: GOToolBox: functional analysis of gene datasets based on Gene Ontology. Genome Biol. 2004, 5: R101-PubMedPubMed CentralView ArticleGoogle Scholar
- Alexa A, Rahnenfuhrer J, Lengauer T: Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006, 22: 1600-1607.PubMedView ArticleGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005, 102: 15545-15550.PubMedPubMed CentralView ArticleGoogle Scholar
- Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009, 4: 44-57.PubMedView ArticleGoogle Scholar
- Li H, Lovci MT, Kwon YS, Rosenfeld MG, Fu XD, Yeo GW: Determination of tag density required for digital transcriptome analysis: application to an androgen-sensitive prostate cancer model. Proc Natl Acad Sci USA. 2008, 105: 20179-20184.PubMedPubMed CentralView ArticleGoogle Scholar
- Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25: 25-29.PubMedPubMed CentralView ArticleGoogle Scholar
- Robinson MD, Smyth GK: Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007, 23: 2881-2887.PubMedView ArticleGoogle Scholar
- Robinson MD, Smyth GK: Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics. 2008, 9: 321-332.PubMedView ArticleGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009, 26: 139-140.PubMedPubMed CentralView ArticleGoogle Scholar
- Benjamini Y, Hochberg Y: Controlling the false discovery rate: a pratical and powerful approach to multiple testing. J R Stat Soc Series B (Methodological). 1995, 57: 289-300.Google Scholar
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5: 621-628.PubMedView ArticleGoogle Scholar
- Wallenius KT: Biased sampling: the non-central hypegeometric probability distribution. 1963, PhD: Stanford UniversityGoogle Scholar
- Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008, 18: 1509-1517.PubMedPubMed CentralView ArticleGoogle Scholar
- Feldman BJ, Feldman D: The development of androgen-independent prostate cancer. Nat Rev Cancer. 2001, 1: 34-45.PubMedView ArticleGoogle Scholar
- DePrimo SE, Diehn M, Nelson JB, Reiter RE, Matese J, Fero M, Tibshirani R, Brown PO, Brooks JD: Transcriptional programs activated by exposure of human prostate cancer cells to androgen. Genome Biol. 2002, 3: RESEARCH0032-PubMedPubMed CentralView ArticleGoogle Scholar
- t Hoen PA, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, de Menezes RX, Boer JM, van Ommen GJ, den Dunnen JT: Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 2008, 36: e141-View ArticleGoogle Scholar
- GOseq software. [http://bioinf.wehi.edu.au/software/goseq/]
- The R Project for Statistical Computing. [http://www.r-project.org/]
- Robinson MD, Oshlack A: A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology. 2010.Google Scholar
- Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R25-PubMedPubMed CentralView ArticleGoogle Scholar
- Ensembl BioMart. [http://www.biomart.org]
- Barnard GA: Contribution to the discussion of Professor Bartlet's paper. J Roy Stat Soc Series B (Methodological). 1963, 25: 294-Google Scholar
- Efraimidis P, Spirakis PG: Weighted random sampling with a reservoir. Information Processing Lett. 2005, 97: 181-185.View ArticleGoogle Scholar
- Smyth GK: Limma: linear models for microarray data. Bioinformatics and Computational Biology Solutions using R and Bioconductor. Edited by: Gentleman R, Care V, Dudoit S, Irizarry R, Huber W. 2005, New York: Springer, 397-420.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.