RNA seq data: normal models and missing data Rhonda Bacher Spring 2012 Honors Undergraduate Thesis Bachelor of Science in Statistics
RNA seq data: normal models and missing data Rhonda L. Bacher 1,2,3 1 Department of Sta tistics, University of Florida, Gainesville, Florida, USA 2 Department of Mathematics, University of Florida, Gainesville, Florida, USA 3 University of Florida Genetics Institute, Cancer/Genetics Research Complex, Gainesville, Florida, USA. Abstract RNA Seq is a tool used for assessing gene expression based on read counts from high throughput sequencing. Many analysis methods to detect differential expression thus far have focused on using the Poisson distribution or the negative binomial distribution fo r analysis of differential expression. Looking at both the underlying data and the measurement on which the analysis is performed we propose that RNA Seq data can be considered continuous. The underlying libraries are made from a solution of mRNA quantifi ed by concentration. This solution is sampled and sequencing technology is used to estimate the number of molecules in the sample for a particular gene. Normalization techniques, which result in non integer values, are often applied to RNA Seq data in orde r to account for systematic effects on the total number of counts, for example the length of the exon/transcript and the total number of reads mapped to the reference per sample. Examination of four different experiments and using a general linear model wi th a normal distribution reveal the residuals do conform to underlying assumptions when data are complete. When some observations are missing the residual assumptions are often violated. The effect from these missing observations on the model and significa nce tests can be best avoided through an imputation technique. Introduction The application of next generation sequencing technologies to sequence cDNA fragments fro m samples of RNA is called RNA s eq (Mortazavi, Williams, Mccue, Schaeffer, & Wold, 2008; Z. Wang, Gerstein, & Snyder, 2010) RNA seq is poised to become the dominant technology in th e analysis of gene expression. Its ability to quantify transcription levels has led to the use of RNA Seq to measure tissue specific expression (Mortazavi et al., 2008) allele specific expression (McManus et al., 2010; Graze et al., 2012) and alternative transcript isoforms (E. T. Wang et al., 2008) As with all new technologies, analytical and s tatistical approaches need to be developed. The statistical issues are varied and complex and as yet no consensus ha s been reached. RNA Seq data have unique properties in which issues with missing data (Degner et al., 2009) technical error (Hillier et al., 2008; McIntyre et al., 2011; Quail et al., 2008) normalization (Bullard, Purdom, Hansen, & Dudoit, 2010) and a ppropriate models for differential expression must be considered. Sta tistical solutions are needed for these problems in order to realize the full potential of RNA Seq to make accurate biological inferences. To date, much discussion on the analysis of RN A Seq data has focused on modeling read counts as a discrete variable. In the presence of counting data, a Poisson model is a natural choice (Jiang & Wong, 2009; Marioni, Mason, Mane, Stephens, & Gilad, 2008) The Poisson employs the same parameter for its mean and variance, a high
mean represen ts a high variance. Yet, the Poisson distribution may not adequately describe the variation in RNA Seq data (Robinson & Smyth, 2007) As a result, other approaches such as the negative binomial (Anders & Huber, 2010; Hardcastle & Kelly, 2010; Robinson, McCarthy, & Smyth, 2010) beta binomial (Zhou, Xia, & Wright, 2011) the generalized Poisson model (Srivastava & Chen, 2010) a two stage Poisson (Auer & Doerge, 2011) Poisson linear model (Li, Jiang, & Wong, 2010) and an over dispersed Poisson (Auer & Doerge, 2010 ) have been used to describe variation in the data. Gene by gene tests of differential expression using likelihood based approaches have been developed for each of thes e models (Auer & Doerge, 2011; Bullard Purdom, Hansen, & Dudoit, 2010; Langmead, Hansen, & Leek, 2010; Marioni et al., 2008; Srivastava & Chen, 2010; Zhou et al., 2011) Sharing information across genes has been ibution (Anders & Huber, 2010; Hardcastle & Kelly, 2010; Robinson et al., 2010) Selecting a model for RNA Seq data analysis is based on how well the model describes the distribution of rea d counts. However, biological differenc es affect the number of counts, such as length of the transcript (Bullard et al., 2010; Mortazavi et al., 2008; Oshlack & Wakefield, 2009) In addition, technical variation during the sequencing process such as the total number o f reads mapped to the reference per sample also affect the total number of reads. Normalization is used to correct these differences, making it a critical step in RNA Seq data analysis (Anders & Huber, 2010; Bullard et al., 2010; Langmead et al., 2010; Robinson & Oshlack, 2010) Normalizing as a function of r eads per kilobase of exon model per million mapped reads (RPKM) (Mortazavi et al., 2008) is widely used. Scaling methods (Bullard et al., 2010; Robinson & Oshlack, 2010) have also been proposed. In addition to these methods, a logarithmic transformation has been shown to stabilize the variance seen in genes having a high number of reads (Robinson & Oshlack, 2010) The primary reason for viewing the data as discrete are the initial read cou nts. Yet, when the data generating process and measurement are examined the assumption that RNA Seq data are discrete is called into question. The underlying data from RNA Seq are the molar concentration of mRNA (Mortazavi et al., 2008) Molar concentration is continuous. The read c ounts themselves span from counts of 0, 1 to 6 and 7 digit numbers. The normalization methods employed produce non integer values which are often used as the dependent variable in the analysis. Additional transformations may further smooth the distributio n (Bullard et al., 2010; Marioni et al., 2008; Robinson & Oshlack, 2010) The combination of large range, non integer values, and an underlying continuous process suggests that a continuous distribution is worth considering. A continuous approach to RNA Seq data analysis has been suggested (Langmead et al., 2010; Marioni et al., 2008 ) It has also been found that using generalized linear model derived likelihood ratio statistics for detecting differential expression perform as (Bullard et al., 2010) The likelihood ratio statistics are a flexible tool and have been reported to detect differential expression among low count genes (Bullard et al., 2010) Langmea d et al., introduces a RNA Seq data analysis pipeline for identifying differential expression named Myrna. The Myrna pipeline utilizes
a generalized linear model which implements both a Gaussian and Poisson distribution. They apply both of these functions to data from the HapMap project and compare p values from the tests of differential expression; they demonstrate that the Gaussian is more flexible than the Poisson, provides a better fit of the data, and reduces bias in detecting differential expression. A likelihood ratio statistic combined with a raw data permutation procedure is implemented. A n additional complexity in the data becomes apparent when counts are observed in one sample from a particular condition but not all samples (Bullard et al., 2010) Are these transcripts truly not expressed? Unlike microarray experiments, in RNA Seq it is possible to fail to sample low abundant transcripts due to the low sampling fraction and the competition among genes (McIntyre et al., 2011) Even in technical replicates there are differences in detection, especially among low counts (McIntyre et al., 2011) As in proteomics, missingness is not random but rather is related to abundance (Liu, Sadygov, & Yates, 2004; Oberg et al., 2008; P. Wang, Tang, Zhang, Whiteaker, & Paulovich, 2006) I hypothesize that varying degrees of missingness affects the model and many of the issues around finding an appropriate model fit are due to missing data. Though many techniques have b een developed to handle normalization and tests for differential expression, there has not yet been a procedure detailing how to handle missing data in RNA Seq. There are a variety of statistical methods for analyses with missing data. Oberg et al. conclud e that fitting their ANOVA model with a censoring procedure works well with missing data in proteomics. Censoring of data is a procedure well known in survival analysis, where the missing data does not have to be excluded from the analysis (Kleinbaum & Kle in, 2005). The caveat here is the data are not missing completely at random or even missing at random, but missingness as a function of abundance. In which case, precautions must be taken to insure that inferences are not heavily affected due to the underl ying processes (Egleston & Wong, 2009) For RNA Seq this implicates that tests for differential expression may ultimately be affected by missing data. This thesis seeks to address the following questions: 1) Can Normal theory apply to RNA seq data? 2) Does missing data affect model fit? These questions are addressed by sy stematically examining four different RNA Seq datasets. A general linear model with Normal assumptions is fit to each and model fit is evaluated for the following categories: complete data, mostly complete data, large amounts of missing and extremely miss ing. Models are evaluated where missing data are dropped from the model, set to a constant, and imputed using a nave procedure. Materials and Methods I begin by looking at the distribution within a gene, since differential expression is tested gene by g ene. The Poisson model is a logical choice when considering the data for all genes. However, the distribution within a gene is not obviously Poisson (figure 1), and even appears approximately normal. I fit a general linear model to the RNA Seq data and eva luate the model fit. R esidual analysis is a common method for determining the
model fit. Residual analysis is o ften performed visually, including residual plots across variables ( Kutner et al., 2005 ) However, with tens of thousands of genes being tested fo r each set of data, individual gene by gene visual analyses are not feasible. Instead, a test for normality for each set of residuals can be used to determ ine how many genes are deviating from model assumptions. I compare the number of genes with residuals inconsistent with normality across our various response variables and consider the cause of these inconsistencies. Datasets Two way ANOVA Hi Seq of Sorghum Wilfred Vermerris, Saballos A., Murray S., Rooney W., and Kresovich S. performed Hi Seq on Sorghu m tissue samples. There are two types of leaf tissue and two types of stem tissue which have either low or high expression of sugar. They use two biological replicates of each treatment resulting in 16 different samples. Four samples are bar coded and run in one lane; this is done in four lanes. RNA Seq of Cegs Sergey Nuzhdin at University of Southern California performed RNA Seq on two genotypes of Drosophila r208 and r303, with three biological replicates of female heads, flies were either mated or virg in. All replicates were multiplexed and sequenced on ten lanes. One Way ANOVA RNA Seq of Fru Network Justin Dalton and Michelle N. Arbeitman performed RNA Seq on Berlin a lab strain of Drosophila melanogaster The samples were extracted from adult fly hea ds for three biological replicates for each male and female. RNA Seq of Dros. Mel. Sergey Nuzhdin at University of Southern California sequenced male and female flies, each with three biological replicates. There are two technical replicates for each biolo gical replicate. Models RPKM is widely used in RNA Seq analyses (Mortazavi et al., 2008) ; adjusts for the number of total mapped reads per sample and transcript length. The mean variance re lationship seen in RNA Seq renders the log transformation appropriate ( Kutner et al., 2005 ). I use the natural log of RPKM as the response variable in our analysis. A general linear model applied by gene for each of the datasets. Two way ANOVA model: Y i j 1 X 1 2 X 2 + B 3 X 1 X 2 Sorghum Model:
Where Y ij = ln (RPKM), X 1 = sugar, X 2 = tissue, i=1,2,..,8 and j= 1,2. Cegs model: Where Y ij = ln (RPKM), X 1 =treatment, X 2 One way ANOVA model: Y ij 1 X 1 Dmel model: Y ij = ln (RPKM), X 1 =sex, i=1,2 and j=1,2,3. Fru model: Y ij = ln (RPKM), X 1 =genotype, i=1,2 and j=1,2,3. For datasets containing technical replicates I averaged across the replicates. Analysis The model is fit separately for each gene. Genes which contain zero coverage across all samples are removed from the analysis. In addition, genes which have zero coverage in all samples for at least one treatment are also removed as no F test can be calculated in this case. The assumption th at the residuals are normally distributed is tested using the Shapiro Wilk test. Small values of the p value(less than .05) suggest that the data are not consistent with normality. If the proportion of genes with residuals not consistent with normality is large this may suggest that the model is not appropriate. Influential points in the data are identified by removing each point and calculating the effect removing the point has on the overall model. Two useful statistics in determining influence are DFFI between the fitted value and the predicted value once the i th point is removed from the data when fitting the model. A general cutoff for DFFITS is 2, but 2* is a recommended cutoff where p are the number of predictors in the model and n are the th point being removed from the data on all of the fitted values. The point F distribution (Kutner et al., 2005) T he residuals themselves are also useful in determining influence. Studentized residuals are able to account for the possibility of the residuals having different variances. The studentized residuals are useful in looking for candidate influential points which have extreme Y values or are outliers. The externally studentized residuals or the studentized deleted resid uals, which is the i th residual when the i th point is removed from the data, are an additional resource. An observation is influential if the absolute value of the residual is large (Kutner et al., 2005)
Some statistics on the residuals such as skewness and kurtosis reveal the shape of the distribution for the data. The measures of skewness and kurtosis of the residuals is compared to the values of skewness and kurtosis for a normal distribution. Skewness measures the symmetry surrounding the peak of the distribution. A normal distribution is symmetric thus it has zero skew. Kurtosis measures the peak of the distribution. Some distributions have very flat peaks with thin tails (low kurtosis) while other have narrow peaks with heavy tails (high kurtosis). A normal distribution has a kurtosis of about three (Lovric, 2011) Fitting a linear model on ln (RPKM) causes the zero observations to drop out. When missing observations are present, the model is only analyzing the data that is nonzero. In the case that all samples in a treatment are zero, there is no way of analyzing this data appropriately. It is also not appropriate to analyze the group of missing more than 50% of data. However, for smaller amounts of missing data, the deletion of zero observations may not be optimal, particularly if the missingness is a re sult of a low abundance level. I can resolve these missing data in two ways. The addition of a c onstant is a common solution, it corrects the presence of zero read counts observed in RNA Seq data, for each analysis, let C = 100. I also apply a nave imputation method. The imputation calculates a variance component using the values that are present and the mean for each treatment combination (Lovric, 2011). Since our within gene distributions are approx imately normally distributed, we draw imputed values from a normal distribution. Can the zero observations be the cause of non normality? Genes are divided into five categories describing how much missing data each gene contains. The categories of missin g are no missing, missing less than 20%, missing between 20% 35%, missing between 35% 50%, and missing more than 50%. Examination of residuals is performed for each category. An F test is carried out for each gene for all methods to determine significance Since comparisons are across methods and no biological inferences made, a nominal p value of .05 (Fisher, 1925) is used. Results Choice of model Since we are testing a Normal model and differential expression is tested gene by gene we begin by showing the distribution within a gene is approximately normally distributed in Figure 1. We look at the distribution within a gene that has no missing data since this comprises the majority of genes and notice that there does not appear to be overdispersion. Mul tiplexing, where many samples are sequenced in one lane, may lead to low counts and as a result produce more missing values. Summing or averaging over technical replicates is encouraged. We found no difference in our model fits using the average versus usi ng the sum. To sum over technical replicates however, the design must be balanced, meaning an equal number of samples for each treatment combination.
To determine whether missing data are a source of non normality, we compare analysis for differing categ ories of missing data. Table 1 describes the percentage of genes that are not able to be analyzed due to missing treatment means. After removing these un analyzable data, the remaining genes are categorized into five categories of missingness. For all the datasets the majority of the data falls into the no missing category. Although the percentage overall of genes in the missing between 20% and 35% category is small, this is not a trivial number of genes. To explain, in the Dmel dataset we start with 60,277 genes, of which 7,751 are not able to be included in the analysis. Of the 52,526 genes remaining, 88.2% or 46,329 genes contain no missing data and 1,643 genes contain between 20% and 35% missing data. Residual Analysis A visual inspection of the residua ls cannot be performed for such a large dataset. We rely on finding the number of genes inconsistent with normality for each category of missing. Table 2 shows the percent of residuals inconsistent with normality. This table can be compared to the performa nce of the model when ln (RPKM), ln (RPKM+100) or imputed ln (RPKM) values are used as the response variable in Tables 3 and 4. Although we perform the tests using the Shapiro Wilk test for normality, the results are not qualitatively different when using any of the other normality tests such as Kolmogorov Smirnov, Kramer von Mises, or Anderson Darling (data not shown). Agreement between methods Kappa statistics represent a measurement of agreement (Fleiss, 1981). Kappa statistics are listed for each metho d and for each category of missing (Table 6 ). The agreement for every method in every da taset gets worse as more data are missing. The imputation method for the most pa rt has better agreement with ln (RPKM) than does ln ( RPKM+100). T he imputation method cal ls many more genes significant versus the other two methods (Table 5) Conclusions & Discussion As RNA Seq technology becomes the dominant source of measuring transcriptome expression levels, analysis techniques will need to be resolved. Here we conside r a general linear model and find that the overall model performs well for each of the datasets. However, the analysis is limited by missing data. We define data as missing due to the ambiguity as to the origin of the zero coverage. Either the gene is not expressed or the transcript is expressed but not sequenced. Since missingness is related to gene abundance, it is difficult to distinguish between these two cases. Consider the situation where we have zero expression in one treatment but coverage in anothe r treatment. This is interesting biologically, but is not analyzable. Figure 2 shows that when zero coverage is present in some samples, the average of the remaining measurements of coverage is also low. We are not seeing the extreme case of completely abs ent in one treatment but highly expressed in another. The best way to examine these genes with high amounts of missing data is to increase the amount of sequencing per sample, that is to increase the sampling fraction.
The group of data that has more than 50% missing but not complete zero coverage in any one treatment should also be avoided in the analysis. Even though the model is able to run on this data and we can look at residuals; the tests for gene significance often fail and very few are statistical ly significant. Imputing data with such extreme missingness is not recommended. Similarly, adding a constant builds a point mass at that value causing the model to have clear under dispersion and perform poorly. Fortunately, this group does not constitute a large portion of the data. There are implications of using ln (RPKM) as a response variable where the zeros are dropped from the analysis. This response variable is appropriate when there is no missing data and the model performs well. However, as more missing data are present and zeros drop out and the model is analyzing less data. This results in a drop in the percent of genes found to be significantly different at a nominal significance level of .05 a s more missing data are present. It is unclear if d ropping the zero values is biologically sensible, again due to the propensity for low abundance transcripts to have zero values. I attempt to resolve this issue by adding a constant to RPKM prior to the natural log transformation. This method in itself cr eates problems. In fact, we see a decrease in the performance of t he model as more missing data are present. More genes are inconsistent with nor mality when more missing data are present. Adding the constant also results in calling fewer genes significant especially in the categories with more than 35% missing data These issue s are best explained by describing the addition of a constant as creating a point mass at the value ln (RPKM+100). For the no missing group of data this method has considerable agreem ent with the ln (RPKM) method. However, t he agreement between the methods when missing data is present is moderate and the agreement reduces as more missing data are present; the kappa statistics get smaller. I also perform a nave imputation method which uses the overall variance and treatment means to predict values for each treatment group. The imputation method results in calling more genes significant versus the other methods for every category of missing. The agreement between the imputation methods and ln (RPKM) tends to be slightly improved; however the agreement get s weaker as more missing data are present. These findings are consistent across all four datasets, representing different species and various designs and it is unlikely to be an artifa ct of any one dataset. I conclude that missingness is an issue in analyzing RNA Seq data. Although I have only examined a very simple imputation I can already say that it is a better to impute than to add a constant. With the imputation we are not affectin g the analysis of the no missing data, meaning we are not introducing any bias, which is important since no missing tends to be the majority of the data included in the analysis. In adding a constant to RPKM, there is an implicit assumption that these miss ing values have no variance. Biologically this is not reasonable, making the imputation a more appropriate approximation to what is expected. The agreement between the imputation and ln (RPKM) is substantial when 20% or less of the data are missing. The 35 % missing group has only moderate
agreement, due to ln (RPKM) method identifying fewer genes significant since less data are analyzed. The missing data should not be ignored and dropped from an analysis, clearly they affect statistical significance. In an analysis where missing data are ignored, the genes with more than 35% of data are likely to be excluded since the model using ln (RPKM) will fail to be significant in many cases. Imputing the missing values for this category may improve inferences. Fina lly, there is no evidence that normal models fail to fit for complete data. On the contrary, for complete data normal models fit very well indicating that continuous methods can be used for RNA Seq analysis. The breadth of normal theory means that accomm odating larger studies with more complex designs will be more straightforward than if non care fully considered. Imputation is preferable to ignoring the zero counts or adding a constant. More sophisticated multiple imputation techniques will undoubtedly improve model fits further, but perhaps the best solution is to increase the sampling fraction. FIGURES AND TABLES Figure1: The distribution of ln (RPKM+100) across all genes (top) and within one gene with no zero coverage (bottom) for all four datasets. The distribution of ln (RPKM +100) within a gene approximates a normal distribution with no obvious signs of over dispersion.
Figure 2: Removing the zeros from all treatments and ta king the mean of RPKM. The mean among the nonzero R PKM is smaller when more data are missing. Plotted on a log scale. Percent of genes in each category of missing Cegs Fru Dmel Sorghum No Missing 30.97% 84.55% 88.20% 91.40% Missing less than 20% 51. 11% 6.33% 4.00% 7.67% Missing between 20% 35% 13.85% 3.99% 3.13% .83% Missing between 35% 50% 3.33% 3.09% 2.68% .10% Missing greater than 50% .74% 2.05% 1.99% -Zero coverage for at least one treatment 37.42% 7.46% 6.74% 21.68% Zero covera ge for all treatments 2.24% 11.01% 6.56% 8.59% Table 1: Genes are classified based on the percentage of missing observations. Genes that have zero coverage for at least one treatment are removed prior to the analysis; these percentages are in the bottom two panels.
Percent of residuals failing normality in each category of missing using ln(RPKM) Cegs Fru Dmel Sorghum Overall 2.55% 2.45% 2.54% 5.41% No Missing 4.17% 2.77% 2.73% 5.21% Missing less than 20% 1.93% .13% .05% 5.69% Missing between 20% 35% 1.38% 1.82% 2.80% 18.75% Missing between 35% 50% .77% 0.00% 0.00% 57.14% Missing greater than 50% 7.45% 0.00% 0.00% -Table 2: Using ln(RPKM) as the response variable in the linear model on the non zero observations the percent of residuals incons istent with normality separated by percent of missing data. Percent of residuals failing normality in each category of missing using ln(RPKM + 100) Cegs Fru Dmel Sorghum Overall 11.92% 4.05% 3.38% 12.59% No Missing 5.32% 2.74% 2.41% 11.27% Missing l ess than 20% 12.00% 2.71% 1.90% 27.17% Missing between 20% 35% 19.07% 8.74% 6.03% 23.30% Missing between 35% 50% 34.02% 4.17% 3.48% 14.29% Missing greater than 50% 48.84% 53.05% 44.92% -Table 3: Percent of residuals failing normality for each categor y of missing for the general linear model with ln(RPKM+100) as the response. Percent of residuals failing normality in each category of missing using imputed ln(RPKM ) values Cegs Fru Dmel Sorghum Missing less than 20% 2.86% 3.37% 2.67% 4.39% Missing between 20% 35% 2.27% 4.09% 3.35% 3.98% Missing between 35% 2.83% 3.59% 3.19% 5.00%
50% Missing greater than 50% 10.63% 0.00% 0.00% -Table 4: Percent of residuals failing normality for each category of missing for the general linear model with imputed ln(RPKM) as the response. Percentage of genes significant in each category of missing for all three methods. Cegs Fru Dmel Sorghum Method A B C A B C A B C A B C No Missing 7.56 8.06 -13.81 14.00 -45.18 44.77 -88.31 86.83 -20% Missing 6.37 11.75 14.07 6.75 5.38 11.91 9.14 9.66 15.66 65.72 65.16 74.63 35% Missing 5.77 1.52 24.92 8.34 8.29 23.90 8.22 11.38 26.11 32.39 42.61 64.20 50% Missing 4.71 .09 38.90 13.43 0.00 47.59 4.83 0.00 39.82 20.00 9.52 57.14 Table 5: Percentage of genes statis tically significant (pvalue < .05) for each dataset separated by categories of missing. Method A is ln (RPKM), method B is ln (RPKM +100), and method C is the imputation (no values are imputed in the no missing group) Kappa statistics for each dataset Cegs Fru Dmel Sorghum Compari son A B C A B C A B C A B C 20% Missing .5395 .3653 .4520 .6751 .5003 .4411 .6724 .6228 .5944 7745 3733 6348 35% Missing .2951 .1200 .0426 .4491 .1613 .2821 .4002 .2722 .4019 4215 5924 5201 50% Missing .1442 .0017 .0027 .2918 --.1424 --.2 857 -.1137 Tab le 6: Kappa statistics for a comparison of statistical significance. Comparison of statistical significance for each dataset comparing methods and separated by category of missing. Comparison A is between the imputation method and ln(RPKM). Comparison B i s between the ln(RPKM+100) method and ln(RPKM). Comparison C is b etween the two methods imputation of ln(RPKM) and ln(RPKM+100). Acknowledgements I would like to thank Rita M. Graze 2, 3 Wilfred Vermerris 3, 4 Ann L. Oberg 5 Sergey V. Nuzhdin 6 Justin D alton 7 Michelle N. Arbeitman 7 (in no particular order) for their support and use of their data. In addition, I thank Lauren M. McIntyre 1 2, 3 for guidance and support throughout this research.
1 Department of Statistics, University of Florida, Gainesvil le, Florida, USA 2 Department of Molecular Genetics and Microbiology, University of Florida, Gainesville, Florida, USA 3 University of Florida Genetics Institute, Cancer/Genetics Research Complex, Gainesville, Florida, USA. 4 Department of Agronomy, Univer sity of Florida, Gainesville, Florida, USA 5 Department of Health Sciences Research, Division of Biomedical Statistics and Informatics, Mayo Clinic, Rochester, Minnesota, USA 6 Molecular and Computational Biology, Department of Biological Sciences, Univer sity of Southern California, Los Angeles, California, USA 7 Department of Biomedical Sciences, College of Medicine, Florida State University, Tallahassee, Florida, USA References Anders, S., & Huber, W. (2010 ). Differential expression analysis for sequence count data. Genome biology 11 (10), R106. BioMed Central Ltd. doi:10.1186/gb 2010 11 10 r106 Auer, P. L., & Doerge, R. W. (2011). A Two Stage Poisson Model for Testing RNA Seq Data. Statistical Applications in Genetics and Molecular Biology 10 (1). doi:10.2202/1544 6115.1627 Bullard, J. H., Purdom, E., Hansen, K. D., & Dudoit, S. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA Seq experiments. BMC bioinformatics 11 94. doi:10.1186/1471 2105 11 94 Degner, J. F., Marioni, J. C., Pai, A. a, Pickrell, J. K., Nkadori, E., Gilad, Y., & Pritchard, J. K. (2009). Effect of read mapping biases on detecting allele specific expression from RNA sequencing data. Bioinformati cs (Oxford, England) 25 (24), 3207 12. doi:10.1093/bioinformatics/btp579 Egleston, B. L., & Wong, Y. N. (2009). Sensitivity analysis to investigate the impact of a missing covariate on survival analyses using cancer registry data. Statistics in medicine 2 8 (10), 1498 1511. Wiley Online Library. doi:10.1002/sim.3557.Sensitivity Fisher, R. A. Statistical methods for research workers. Edinburgh:Oliver & Boyd, 1925. F leiss J.L. (1981). Statistical Methods for Rates and Proportions. 2nd ed. New York: John Wile y & Sons. Hardcastle, T. J., & Kelly, K. a. (2010). baySeq: empirical Bayesian methods for identifying differential expression in sequence count data. BMC bioinformatics 11 422. doi:10.1186/1471 2105 11 422 Hillier, L. W., Marth, G. T., Quinlan, A. R., D ooling, D., Fewell, G., Barnett, D., Fox, P., et al. (2008). Whole genome sequencing and variant discovery in C elegans. Nature Methods 5 (2), 183 188. doi:10.1038/NMETH.1179
Jiang, H., & Wong, W. H. (2009). Statistical inferences for isoform expression in RNA Seq. Bioinformatics (Oxford, England) 25 (8), 1026 32. doi:10.1093/bioinformatics/btp113 Kleinbaum, David G., and Mitchel Klein. Introduction to Survival Analysis Springer New York, 2005. eBook. Kutner, Michael H., Christopher J. Nachtsheim, John N eter and William Li. 2005. Applied linear statistical models. 5th ed. Boston: McGraw Hill. Langmead, B., Hansen, K. D., & Leek, J. T. (2010). Cloud scale RNA sequencing differential expression analysis with Myrna. Genome biology 11 (8), R83. doi:10.1186/gb 2010 11 8 r83 Li, J., Jiang, H., & Wong, W. H. (2010). Modeling non uniformity in short read rates in RNA Seq data. Genome biology 11 (5), R50. doi:10.1186/gb 2010 11 5 r50 Liu, H., Sadygov, R. G., & Yates, J. R. (2004). A Model for Random Sampling and Es timation of Relative Protein Abundance in Shotgun Proteomics proteolytic digestion and liquid chromatography in com 76 (14), 4193 4201. Lovric, Miodrag. International Encyclopedia of Statistical Science.Berlin, Heidelberg: Springer, 2011. Marioni, J. C., Mason, C. E., Mane, S. M., Stephens, M., & Gilad, Y. (2008). RNA An assessment of technical reproducibility and comparison with gene expression arrays. Genome Research 1509 1517. doi:10.1101/gr.079558.108. McIntyre, L. M., Lopiano, K. K., Morse, A. M., Amin, V., Oberg, A. L., Young, L. J., & Nuzhdin, S. V. (2011). RNA BMC genomics 12 (1), 293. BioMed Central Ltd. doi:10.1186/1471 2164 12 293 McManus, C. J., Coolon, J. D., Duff, M. O., Eipper Mains, J., Grave ley, B. R., & Wittkopp, P. J. (2010). Regulatory divergence in Drosophila revealed by mRNA seq. Genome research 20 (6), 816 25. doi:10.1101/gr.102491.109 Mortazavi, A., Williams, B. A., Mccue, K., Schaeffer, L., & Wold, B. (2008). Mapping and quantifying m ammalian transcriptomes by RNA Seq. Nature Methods 5 (7), 1 8. doi:10.1038/NMETH.1226 Oberg, A. L., Mahoney, D. W., Eckel Passow, J. E., Malone, C. J., Wolfinger, R. D., Hill, E. G., Cooper, L. T., et al. (2008). Statistical analysis of relative labeled ma ss spectrometry data from complex samples using ANOVA. Journal of proteome research 7 (1), 225 33. doi:10.1021/pr700734f
Oshlack, A., & Wakefield, M. J. (2009). Transcript length bias in RNA seq data confounds systems biology. Biology direct 4 14. doi:10 .1186/1745 6150 4 14 Quail, M. A., Kozarewa, I., Smith, F., Scally, A., Stephens, P. J., Durbin, R., Swerdlow, sequencing system. Nature Methods 5 (12), 1005 1010. doi:10.1038/NmETH. 1270 Robinson, M. D., & Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA seq data. Genome biology 11 (3), R25. doi:10.1186/gb 2010 11 3 r25 Robinson, M. D., & Smyth, G. K. (2007). Moderated statistical tests fo r assessing differences in tag abundance. Bioinformatics (Oxford, England) 23 (21), 2881 7. doi:10.1093/bioinformatics/btm453 Robinson, M. D., McCarthy, D. J., & Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (Oxford, England) 26 (1), 139 40. doi:10.1093/bioinformatics/btp616 Srivastava, S., & Chen, L. (2010). A two parameter generalized Poisson model to improve the analysis of RNA seq data. Nucleic acids research 38 (17), e170. doi:10.1093/nar/gkq670 Wang, E. T., Sandberg, R., Luo, S., Khrebtukova, I., Zhang, L., Mayr, C., Kingsmore, S. F., et al. (2008). Alternative isoform regulation in human tissue tr anscriptomes. Nature 456 (7221), 470 6. doi:10.1038/nature07509 Wang, P., Tang, H., Zhang, H., Whiteaker, J., & Paulovich, A. G. (2006). Spectrometry Data Martin Mcintosh. Cancer Research 326 315 326. Wang, Z., Gerstein, M., & Snyder, M. (2010). RNA : a revolutionary tool for transcriptomics. October 10 (1), 57 63. doi:10.1038/nrg2484.RNA Seq Zhou, Y. H., Xia, K., & Wright, F. a. (2011). A powerful and flexible approach to the analysis of RNA sequence count data. Bioinformatics (Oxford, England) 27 (1 9), 2672 8. doi:10.1093/bioinformatics/btr449