Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments

James H Bullard, Elizabeth Purdom, Kasper D Hansen, Sandrine Dudoit, James H Bullard, Elizabeth Purdom, Kasper D Hansen, Sandrine Dudoit

Abstract

Background: High-throughput sequencing technologies, such as the Illumina Genome Analyzer, are powerful new tools for investigating a wide range of biological and medical questions. Statistical and computational methods are key for drawing meaningful and accurate conclusions from the massive and complex datasets generated by the sequencers. We provide a detailed evaluation of statistical methods for normalization and differential expression (DE) analysis of Illumina transcriptome sequencing (mRNA-Seq) data.

Results: We compare statistical methods for detecting genes that are significantly DE between two types of biological samples and find that there are substantial differences in how the test statistics handle low-count genes. We evaluate how DE results are affected by features of the sequencing platform, such as, varying gene lengths, base-calling calibration method (with and without phi X control lane), and flow-cell/library preparation effects. We investigate the impact of the read count normalization method on DE results and show that the standard approach of scaling by total lane counts (e.g., RPKM) can bias estimates of DE. We propose more general quantile-based normalization procedures and demonstrate an improvement in DE detection.

Conclusions: Our results have significant practical and methodological implications for the design and analysis of mRNA-Seq experiments. They highlight the importance of appropriate statistical methods for normalization and DE inference, to account for features of the sequencing platform that could impact the accuracy of results. They also reveal the need for further research in the development of statistical and computational methods for mRNA-Seq.

Figures

Figure 1
Figure 1
Comparison of differential expression statistics: ROC curves. (a) All DE statistics, no gene filtering. (b) GLM-based likelihood ratio statistics and t-statistics, before and after removing genes with fewer than 20 reads in either Brain or UHR. In both plots, a gene was declared non-DE if its qRT-PCR absolute log-ratio was less than 0.2 and DE if its absolute log-ratio was greater than 2.0. Note that we require a true positive to be differentially expressed in the same direction according to both mRNA-Seq and qRT-PCR (see Table 1 and Methods).
Figure 2
Figure 2
Differential expression statistics, by length. Boxplots of the ranks of DE statistics vs. gene lengths for UI genes at least 250 bp-long and with non-zero counts in both Brain and UHR. (a) Delta method t-statistics. (b) Delta method t-statistics weighted by the inverse of the square root of gene length.
Figure 3
Figure 3
Impact of base-calling calibration method on read-mapping. Barplots of average read counts per lane with and without phi X calibration, for each of the four biological sample (Brain, UHR) and flow-cell (F2, F3) combinations. Reads are classified into three nested categories: purity-filtered perfectly matching reads (FPM); purity-filtered reads with either 0, 1, or 2 mismatches (FMM); unfiltered reads with either 0, 1, or 2 mismatches (MM).
Figure 4
Figure 4
Comparison of biological, library preparation, and flow-cell effects. Boxplots of estimated log-fold-changes for UHR vs. Brain biological effects (GLM 2 in [Additional file 2: Supplemental Table S4]), flow-cell effects adjusting for biology (GLM 4), library preparation effects within flow-cell (GLM 7). Estimates are presented for total-count (black) and upper-quartile (blue) normalization.
Figure 5
Figure 5
Impact of highly-expressed genes. (a) Cumulative percentage of total read count for Brain (green) and UHR (purple) samples, starting with the gene with the highest read count (across the seven Brain or UHR lanes). Cumulative read counts are marked for the 5, 10, 20, and 30 percent most highly expressed genes. (b) Running value of the UHR/Brain expression fold-change for unnormalized counts, starting with the gene with the lowest total count across all 14 lanes. Horizontal lines correspond to: the ratio of the counts for all genes (black), the ratio of the counts for the POLR2A gene (red), and the ratio of the per-lane upper-quartile of counts for genes with reads in at least one lane (blue).
Figure 6
Figure 6
Comparison of mRNA-Seq and microarray differential expression calls to qRT-PCR: ROC curves. Genes common to all three platforms and present for both qRT-PCR and sequencing (see [Additional file 2: Supplemental Section S6]). evaluated and declared DE if their qRT-PCR absolute log-ratio was (a) greater than 2 or (b) greater than 0.5; genes were declared non-DE if their absolute log-ratio was less than 0.2. The GLM-based likelihood ratio test was used for the sequencing data. Two normalization procedures are presented for mRNA-Seq: total-count (black) and upper-quartile (blue) normalization. Microarray data were normalized using RMA (gray). Note that we require a true positive to be differentially expressed in the same direction according to both mRNA-Seq and qRT-PCR (see Table 1 and Methods).
Figure 7
Figure 7
Comparison of mRNA-Seq and microarray differential expression measures to qRT-PCR. Difference scatterplots comparing the estimates of UHR/Brain expression log-ratio from qRT-PCR to those from (a) mRNA-Seq, using the standard total-count normalization, and (b) microarrays, using the standard RMA normalization. Shown are the genes shared between all three platforms, present in both Brain and UHR according to both mRNA-Seq and qRT-PCR (see [Additional file 2: Supplemental Section S1]), and having absolute qRT-PCR expression log-ratio less than 4. Horizontal lines in (a) represent the median UHR/Brain log-ratio for the sequencing data after the standard total-count normalization (black), POLR2A normalization (red), quantile normalization (yellow), upper-quartile normalization (blue); horizontal lines in (b) show the median UHR/Brain log-ratio for the microarray data after the standard RMA normalization (black) and POLR2A normalization (red).
Figure 8
Figure 8
Comparison of normalization procedures: Goodness-of-fit of Poisson model. The multiplicative Poisson model (GLM 1 in [Additional file 2: Supplemental Table S4]) is fit to the seven Brain lanes in the MAQC-2 experiment after (a) total-count, (b) POLR2A, (c) upper-quartile, and (d) quantile normalization. Goodness-of-fit statistics are computed and displayed in χ2 quantile-quantile plots. Genes with goodness-of-fit statistics in the top quantiles of the χ2-distribution are displayed using colored plotting symbols: red (1, 5]%, purple (.1, 1]%, gold [0, .1]%. Similar plots for UHR show the same patterns.

References

    1. Chiang DY, Getz G, Jaffe DB, O'Kelly MJT, Zhao X, Carter SL, Russ C, Nusbaum C, Meyerson M, Lander ES. High-resolution mapping of copy-number alterations with massively parallel sequencing. Nature Methods. 2009;6:99–103. doi: 10.1038/nmeth.1276.
    1. Dohm JC, Lottaz C, Borodina T, Himmelbauer H. Substantial biases in ultra-short read data sets from high-throughput DNA sequencing. Nucleic Acids Research. 2008;36(16):e105. doi: 10.1093/nar/gkn425.
    1. Hoen PAC, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RHAM, de Menezes RX, Boer JM, van Ommen GJB, den Dunnen JT. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Research. 2008;36(21):e141. doi: 10.1093/nar/gkn705.
    1. Lee A, Hansen KD, Bullard J, Dudoit S, Sherlock G. Novel low abundance and transient RNAs in yeast revealed by tiling microarrays and ultra high-throughput sequencing are not conserved across closely related yeast species. PLoS Genetics. 2008;4(12):e1000299. doi: 10.1371/journal.pgen.1000299.
    1. 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. PNAS. 2008;105(51):20179–20184. doi: 10.1073/pnas.0807121105.
    1. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays. Genome Research. 2008;18(9):1509–1517. doi: 10.1101/gr.079558.108.
    1. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods. 2008;5(7):621–628. doi: 10.1038/nmeth.1226.
    1. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, Snyder M. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008;320(5881):1344–1349. doi: 10.1126/science.1158441.
    1. 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(7221):470–476. doi: 10.1038/nature07509.
    1. MAQC Consortium. The MicroArray Quality Control (MAQC) project shows inter-andintraplatform reproducibility of gene expression measurements. Nature Biotechnology. 2006;24(9):1151–1161. doi: 10.1038/nbt1239.
    1. Oshlack A, Wakeffeld MJ. Transcript length bias in RNA-seq data confounds systems biology. Biology Direct. 2009;4(14)
    1. Illumina. Sequencing Analysis Software User Guide For Pipeline Version 1.3 and CASAVA Version 1.0 T. Illumina, Inc.; 2008. [Part # 1005359 Rev. A]
    1. Canales RD, Luo Y, Willey JC, Austermiller B, Barbacioru CC, Boysen C, Hunkapiller K, Jensen RV, Knight CR, Lee KY, Ma Y, Maqsodi B, Papallo A, Peters EH, Poulter K, Ruppel PL, Samaha RR, Shi L, Yang W, Zhang L, Goodsaid FM. Evaluation of DNA microarray results with quantitative gene expression platforms. Nature Biotechnology. 2006;24(9):1115–1122. doi: 10.1038/nbt1236.
    1. Illumina. Preparing Samples for Sequencing mRNA. Ilumina, Inc.; 2009. [Part # 1004898 Rev. A]
    1. Bentley DR, Balasubramanian S, Swerdlow HP. Accurate whole human genome sequencing using reversible terminator chemistry. Nature. 2008;456(7218):53–59. doi: 10.1038/nature07517.
    1. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology. 2009;10(3):R25. doi: 10.1186/gb-2009-10-3-r25.
    1. Ewing B, Green P. Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Research. 1998;8(3):186–194.
    1. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data. Biostatistics. 2003;4(2):249–264. doi: 10.1093/biostatistics/4.2.249.
    1. Taub MA. PhD thesis. Department of Statistics, UC Berkeley; 2009. Analysis of high-throughput biological data: some statistical problems in RNA-seq and mouse genotyping.
    1. Durinck S, Bullard J, Spellman PT, Dudoit S. GenomeGraphs: integrated genomic data visualization with R. BMC Bioinformatics. 2009;10:Article 2. doi: 10.1186/1471-2105-10-2.
    1. Lu J, Tomfohr JK, Kepler TB. Identifying differential expression in multiple SAGE libraries: an overdispersed log-linear model approach. BMC Bioinformatics. 2005;6:165. doi: 10.1186/1471-2105-6-165.
    1. Robinson MD, Smyth GK. Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007;23(21):2881–2887. doi: 10.1093/bioinformatics/btm453.
    1. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4(1465-4644 (Print)):249–64. doi: 10.1093/biostatistics/4.2.249.

Source: PubMed

Szponzorok és közreműködők

Egészségi állapot

Kábítószer-beavatkozások

3
Iratkozz fel