A fundamental assumption, common to the vast majority of high-throughput transcriptome analyses, is that the expression of most genes is unchanged among samples and that total cellular RNA remains constant. As the number of analyzed experimental systems increases however, different independent studies demonstrate that this assumption is often violated. We present a calibration method using RNA spike-ins that allows for the measurement of absolute cellular abundance of RNA molecules. We apply the method to pooled RNA from cell populations of known sizes. For each transcript, we compute a nominal abundance that can be converted to absolute by dividing by a scale factor determined in separate experiments: the yield coefficient of the transcript relative to that of a reference spike-in measured with the same protocol. The method is derived by maximum likelihood theory in the context of a complete statistical model for sequencing counts contributed by cellular RNA and spike-ins. The counts are based on a sample from a fixed number of cells to which a fixed population of spike-in molecules has been added. We illustrate and evaluate the method with applications to two global expression data sets, one from the model eukaryote Saccharomyces cerevisiae, proliferating at different growth rates, and differentiating cardiopharyngeal cell lineages in the chordate Ciona robusta. We tested the method in a technical replicate dilution study, and in a k-fold validation study.
Accurate transcriptome measurements are central to understanding the fundamental mechanisms of gene expression. A main challenge presented by the RNA-seq method for digitizing information about cellular RNA content—both its composition and abundance—is correcting noise, errors, and biases introduced in the process of making the measurement. An important step in typical library preparation for sequencing is random fragmentation of the molecules to be sequenced. The actual unit that is being digitized during RNA-seq is therefore not the RNA molecules directly, but their fragments, whose number depends on transcript length and RNA abundance. However, it is the molecular abundance of a transcript that is of interest, rather than the distribution of fragments over its gene model. The introduction of Reads per Kilobase of exon model per Million mapped reads (RPKM) as the unit measurement by Mortazavi et al. [1] addressed this problem. The work of Mortazavi et al. [1] and of Tarazona et al. [2] both addressed the problem of varyng sample size (sequencing depth). Other researchers looked more carefully at fragment biases and developed a maximum likelihood algorithm to estimate the true differences [3]. The adoption of “read-counts” overlapping each transcript [4] instead of FPKMs (Fragments per Kilobase of exon model per Million mapped reads) [5] has allowed a more intuitive interpretation of the data. Among the models that have been proposed for RNA-seq count data [6–9], the intuitively appealing negative binomial model [10] has become the most popular. The negative binomial probability mass function can be thought of as a mixture distribution arising from concatenating biological noise in transcript abundance (described by a gamma probability density function) and sampling noise (described by a Poisson distribution) in the compilation of corresponding sequencing counts. Applications of the negative binomial distribution and methods of hypothesis testing have been reviewed recently in [11]. Widely used normalization, statistical modeling, and hypothesis-testing methods are implemented in widely used R packages, edgeR [12], EDASeq [13], and DESeq2 [14]. Methods that focus particularly on the removal of unwanted variation from unspecified extraneous, nuisance sources are implemented in the R package RUVSeq [15]. A common assumption among all these methods is that, while the proportions of some transcripts vary across conditions/treatments, most transcripts do not vary between experimental conditions, and the total abundance of cellular RNA remain more-or-less fixed. Lovén et al. [16] were the first to demonstrate an experimental system in which the central assumption of transcriptome equivalence among conditions is not satisfied. The researchers discovered that in cells overexpressing the oncogene cMyc, 90% of all transcripts are also overexpressed. To overcome the problem and allow comparison of expression levels between normal and cMyc overexpressing cells, the group incorporated external spike-ins in their samples, which were then used as a de facto invariant pool of RNAs. The spike-in approach had been previously used successfully in microarray experiments and its use in RNA-seq was facilitated by the development of external RNA spike-in mixes by the External RNA Controls Consortium (ERCC) [17,18]. The need to normalize high-throughput RNA and DNA counts, in general, by the use of spike-in standards was recently explained and validated in the wide-sweeping paper of [19]. Even more recently, external RNA spike-ins were used in single cell RNA-seq experiments [20]. The ERCC external RNA spike-in mix 1 (Ambion) that was used in this study, consists of 92 different synthetic RNAs at 22 different concentrations spanning six orders of magnitude (30,000–0. 01 amol/μL). Instead of relying on normalization methods that aim to match the cellular and spike-in RNA read-count distributions, we take advantage of the digital nature of the RNA-seq output and we use the spike-ins as calibrators of known absolute abundance in the samples. We demonstrate that our calibration/normalization model is applicable in two different model organisms (the unicellular eukaryote Saccharomyces cerevisiae and the multicellular chordate Ciona robusta, three experimental setups, a growth rate regulation and a dilution study in yeast, as well as an embryonic differentiation and cell lineage specification study in Ciona), and two library preparation protocols. We perform hypothesis testing and detect global amplification of gene expression in both organisms.
We used three distinct experimental setups, two representing different cases in which the assumption of constant transcriptome sizes is violated. We added a fixed, known amount of external RNA spike-ins to the sample of cells. In the rest of the paper, we refer to spike-in abundance per cell to mean the ratio of spike-in amount (molecules) added to the sample divided by the number of cells in the sample. In all libraries we avoided the use of poly-dT for reverse transcription as it has been previously shown to be incompatible with quantitation through external spike-ins for RNA-seq [28]. In all cases the filtered aligned reads were converted to counts using the function featureCounts from the package Rsubread (R, Bioconductor) and strand information. GR samples and samples for the dilution study were prepared essentially as described in the Borodina et al. [29] directional RNA-seq protocol. We modified the protocol by using UMI adaptors [30,31] that were used to eliminate PCR duplicates from the results. RNA was extracted from ten million cells after the addition of 2 μl of 1: 20 dilution of ERCC spike-in Mix 1 stock [32] in the lysis buffer. All samples were distributed in two lanes and sequenced on an Illumina HiSeq 2000, with 100 nt-long, single end reads. The RNA-seq data (fastq files) were first filtered for residual rRNA reads. The data were aligned to the latest version of the yeast genome (sacCer3) and spike-in sequences, and filtered for mapping quality using Bowtie with optimized parameters. The aligned reads were then processed with a custom script that removes PCR duplicates based on the combination of mapping coordinates and UMI adaptor barcodes. For the samples of three different in vivo Ciona cell lineages, 800 cells were directly sorted into lysis buffer from RNAqueous-Micro Total RNA Isolation Kit (Ambion) containing a fixed amount of ERCC spike-in mix 1. Total RNA extraction was performed according to the manufacturer’s instructions, followed by depletion of rRNA in the samples. The quality and quantity of total RNA in all stages were measured using Agilent RNA 6000 Pico Kit (Agilent) on Agilent 2100 Bioanalyzer. cDNA were synthesized using the SMART-Seq v4 Ultra Low Input RNA Kit (Clontech). RNA-Seq Libraries were prepared and barcoded using Ovation Ultralow System V2 1-16 (NuGen). The samples were pooled in one lane and sequenced on an Illumina HiSeq 2500, with 50 nt-long, single end reads. The RNA-seq reads were mapped to the Ciona genome (v. 2008, ghost database) using Tophat2 with default parameters. The mapped reads were assigned to Ciona KH gene models (v. 2013). All computer programming, including data analysis, Monte Carlo simulations, hypothesis testing, and figure preparation were done in the R programming language and environment. For pairwise testing for differential gene expression, we used DESeq function in the DESeq2 package [14]. In some cases hypothesis testing followed the application of an RUV normalization method in a suite of 3 R functions, RUVr, RUVs and RUVg, in the RUVseq package [15]. For maximum likelihood estimation of parameters we used the R function nlm, a general R function that minimizes a supplied objective function over its parameters. We used it to minimize minus log likelihood of the observed data in the context of a model. The first argument of the nlm function is the name of the name of the R function that computes minus log likelihood for the problem at hand. We wrote such a function for each maximum likelihood problem we considered.
We use external RNA spike-ins as a calibration tool to normalize RNA-seq counts by introducing the variable relative yield (Fig 1). We parametrize a multinomial model for sampling noise, conditioned on native RNA abundances, with library size, relative yield coefficients, and known absolute abundances of spike-in molecules. In the context of this model, we derive a maximum likelihood estimation of nominal RNA abundance, proportional to absolute abundance, for each endogenous transcript in our sample. In the remainder of this paper, we often omit the qualifier “nominal” when we mean nominal abundance, and use the phrase “absolute abundance” when we mean molecules or attomoles (per cell or per sample). In S1 Appendix, we show that the maximum likelihood estimator of RNA abundance is a library dependent scaling of counts by a factor that is proportional to total spike-in counts; it also depends on the known abundance of a reference spike-in, and fraction of overall spike-in counts contributed by this reference spike-in. The expected proportion of counts for a given molecule (native RNA or spike-in) i represented in library j depends on the product of its abundance (attomoles or molecules per cell) in the original sample ni, j and its relative yield coefficient αi; (Fig 1B (ii) ). For spike-in molecules, ni, j are known. We define the relative yield coefficient of a molecule (spike-in or RNA) to be the ratio of its yield coefficient to that of a reference spike-in. By yield coefficient of a molecule we mean the expected number of fragments per molecule contributed by that molecule to a total RNA-seq library of fixed size and prepared according to a fixed protocol. The relative yield coefficient captures specific properties of an RNA molecule such as transcript length and GC content. By convention, we assign the index 1 to the reference spike-in; consequently α1 = 1, by definition (Fig 1B (i) ). For the sake of generality, we do not assume that relative yield coefficient is proportional to transcript length as in [33]. The relative yield coefficient of a spike-in is related to its FPKM within an RNA-seq spike-in library as follows. The relative yield coefficient multiplied by abundance in the original sample, and divided by length is proportional to FPKM. In S1 Appendix, we express the expected proportion of counts for each molecule, indexed by i, in library j in terms of all zi, j = αini, j, in a multinomial joint distribution of counts, and then solve for the maximum likelihood estimator of zi, j. Because we do not estimate the relative yield coefficients, αi for cellular RNA molecules in the present paper, we cannot disentangle here their relative yield coefficients αi and absolute molecular abundances, ni, j (see Discussion). Consequently, we refer to zi, j as a nominal abundance. The corresponding terms for spike-in molecules do not depend on library index j when the same amounts of spike-ins are used in each sample; so, for spike-ins, we can write more simply zi = αini. For sake of clarity, we mention that, by our convention, the indices for the s spike-in molecules are i = 1,2, …, s, and the molecule indices for the q detected native RNA molecules are i = s + 1, s + 2, …, s + q. As shown in S1 Appendix, the derived maximum likelihood values of abundances, zi, j for the RNA molecule i in library j are given by z i, j = y i, j ν j for i = s + 1, s + 2, …, s + q, and j = 1,2, …, r (1) where yi, j is the count (sequencing reads for RNA molecule i in library j and νj is the maximum likelihood calibration constant for library j. The calibration constant νj is given by νj=deff1LjSIn1, (2) where f1 is the proportion of spike-in counts across all libraries contributed by the reference spike-in, n1 is the attomoles or molecules per cell, depending on one’s choice of units, for the reference spike-in, and L j SI is the size (total counts) of spike-in library j. In S1 Appendix we extend the estimation of zi, j to a full statistical model, including biological variation, linking cellular RNA abundance to RNA-seq counts. The νj calibration constant in Eq (1) is qualitatively like the dimensionless “technical” (library) size factor sj of [34], but with an explicit relationship to absolute abundance, because 1/νj is on the scale of attomoles or molecules per cell. The relationship between the two factors is discussed thoroughly in S8 Appendix. The numerator on the right-hand side of Eq (2), according to our statistical model, is the expected number of counts from the reference spike-in, in replicate j, given the spike-in library size L j SI. As shown in S2 Fig, the “expected” counts given by the model for the reference spike-in, closely approximate the actual number. Therefore, Eq (1) says that the inferred abundance of RNA transcript i in replicate j is given by the counts for this transcript multiplied by a scale factor that is the attomoles or molecules per cell, per count of the reference spike-in. If it were known that, in fact, RNA transcript i on average yields twice as many aligned counts per amol as the reference spike-in, the abundance of RNA transcript i in the 107 cells (from which our sample came) would be given by zi, j/2. The normalization in Eq (2) is reminiscent of RPM (reads per million mapped reads) normalization, but the denominator involves the total number of reads in the spike-in library only. It makes intuitive sense, because read depth scales spike-in counts and endogenous RNA counts the same way [17]. Therefore, dividing an endogenous RNA count by spike-in counts derived from a fixed number of molecules in the original biological sample simultaneously normalizes for read depth and provides a measure proportional to the molecular abundance of the endogenous RNA in question. [17] applied this sort of normalization to ERCC spike-in counts, but also divided by spike-in length to obtain FPKM (fragments per kilobase per million mapped reads). Under the assumption that count scales with molecular length, FPKM would be proportional to spike-in abundance, in the absence of any molecular biases, e. g. , GC content (see Fig 2 in [17]). As shown in S1 Appendix, the derived maximum likelihood values of abundances for the s spike-in molecules are given by z i = n 1 f 1 f i for i = 1, … s, (3) where fi is the empirical fraction of total spike-in counts, across all libraries, that is accounted for by spike-in molecule i, and f1 is that of the reference spike-in. Because the absolute abundances of the spike-ins, ni are known, Eq (3) and the definition zi = αini, allows us to estimate the spike-in relative yield coefficients as α i = z i n i = n 1 n i f i f 1 (4) We chose as the reference spike-in the one that contributes the largest fraction of overall spike-in counts. Any of the top few spike-ins would do just as well. We compute the spike-in relative yield coefficients, αi, describe their statistical properties and model the spike-in relative yield coefficients in terms of their biophysical properties in S6 Appendix. However, we do not yet have a large enough repertoire of spike-in molecules nor accurate enough biophysical models to use the model relative yield coefficients for spike-ins to estimate those of native RNA molecules. Thus, the computed spike-in αi terms in the present paper are not participating in the estimation of native RNA abundances. In the Discussion we analyze how αi could be used in the estimation approach. Following the recommendation of [35], we prepared diagnostic relative-log-expression (RLE) plots for all three of our data sets (dilution study, yeast GR study, Ciona embryonic differentiation study) to help in evaluating our maximum likelihood (νj) calibration method. We found unwanted variation in the total inferred RNA abundance within conditions for all three data sets, and the variations are similar across data sets. We ascribe this variation to technical errors in one or more steps in the preparation of samples to be sequenced: variation in RNA extraction efficiency, error in cell count, dilution and/or volume errors in preparation of the spike-ins added to the cellular RNA. We refer to these errors collectively as library preparation errors. Details of results and analyses are presented in S1 Fig and S2 Appendix. Accordingly, we derived a library-specific scale factor, δj (S2 Appendix), much like the total RNA correction factor, ξj, in the single-cell RNA-seq study of [36] (S8 Appendix). The corrected nominal abundance values are computed as zi, j/δj, and we performed statistical analyses and hypothesis testing on these corrected values. In S8 Appendix we discuss, and in S5 Fig. we illustrate, the similarities and differences of this noise reduction to a removal of unwanted variation (RUV) method RUVr in the RUVSeq R package [15] that is based on residuals in a generalized linear model. For this dilution study we prepared 3 replicate libraries with a “high” spike-in aliquot dilution, and 3 with a“low” dilution as described in Materials and Methods. In the absence of any library preparation noise, the molar amount of spike-ins are 4. 44 times larger in samples that were added to low dilution spike-in aliquots compared to those with the high-dilution aliquots. The generalization of the νj normalization in this case is to simply add the subscript j (indexing library) to nref to give ν j = f 1 L j SI / n 1, j. Perfect performance of our method, would give identical mean abundances for libraries prepared with the high- and low-volume aliquots. In Fig 3 the MA plot shows less than ideal performance in that the ordinates in the scatter plot are offset a bit (by 0. 28) from zero. This corresponds to a mean fold difference of 1. 2 rather than 1. This discrepancy might be due, at least in part, to sample preparation handling errors. Modeling RNA abundances within a condition as gamma-distributed random variables results in counts with the familiar negative binomial distribution, as shown in S1 Appendix. Previous applications of the negative binomial distribution for modeling RNA sequencing counts and various methods of hypothesis testing have recently been reviewed in [11]. Our model for sequencing counts is formally equivalent to that in [34], but with an important distinction. The library-specific size factor of [34] (based on genes), written according to our notation, is s j = median i y i, j (∏ k = 1 r y i, k) 1 / r. (5) The size factors sj is a dimensionless constant that carries with it an implicit assumption of fixed total amount of cellular RNA. Furthermore, in practice, the median of sj across libraries is of order 1. In contrast, although our calibration constants νj can be thought of as “size factors, ” they are proportional to the total spike-in library size L j SI in Eq (2), and 1/νj has dimensions of attomoles, or molecules per cell, depending on the units one chooses to use for the spike-in abundances ni (i ∈ {1,2, …s}). Our νj calibration factor is closely related to the extension of [34], by [36] to estimate cellular RNA abundance with the use of a “technical size factor” resembling that in Eq (5), but computed using spike-in counts only. Hypothesis testing methods currently in the literature that are based on the negative binomial distribution of counts with library size factors sj, such as DESeq [34], can be used following our νj normalization in the manner described in S4 Appendix.
Our focus in this paper is on deriving a nominal cellular molecular abundance that can be converted to absolute abundance by the transcript’s relative yield coefficient, which could be measured in separate experiments. In this study however, we do not attempt to measure the relative yield coefficient values, or estimate the absolute number of molecules per cell for each transcript within a condition. The current work allows us to say, that, for example, RNA transcript A has x times more molecules per cell, on average, in condition 1 compared to condition 2, even if the corresponding RNA-seq libraries were prepared in different batteries of experiments, different studies, or even prepared in different laboratories. Such a conclusion about what might be called, an absolute ratio of abundances, can be drawn without knowing the relative yield coefficient of transcript A. In the section that follows, we discuss the links between our work and methods by which these relative yield coefficients might be measured. In this manuscript we offer RNA abundance estimates that are proportional to absolute transcript abundance. For this we assign a (relative) yield coefficient value of 1 to a reference spike-in, arbitrarily chosen from among those that contribute a sizable fraction of total spike-in counts. Our nominal abundance of an RNA molecule is based on the temporary assumption that this molecule has the same yield coefficient as the reference spike-in. If our calibration method is supplemented with additional data on the effect that a broad range of transcript physicochemical characteristics has on library preparation and sequencing, a more realistic relative yield coefficient could be assigned to each RNA molecule of interest. A technical statement of the outstanding problem is that our inferred nominal abundances zi, j do not disentangle true absolute molecular abundance, ni, j, and the corresponding relative yield coefficient, αi; because, by definition, zi, j = αi ni, j. However, once one measures absolute cellular abundance of transcript i in a preparation of cells from which library j was derived (ni, j), the relative yield coefficient becomes known, at, least in the idealized situation ignoring various sorts of noise, because αi = zi, j/ni, j. For example, ni, j might be measured by single-cell Fluorescence In Situ Hybridization (FISH) methods, performed on a large population of cells from which library j was derived. Statistical methods taking into account biological noise and technical noise could be used to compute a confidence interval for αi, provided ni, j could be estimated. Likelihood methods could be used to integrate data across several libraries in the estimation of αi. In principle, once αi is estimated from one or more libraries and a population of cells from which those libraries were derived, this estimate could be used for other libraries (prepared using the same protocol), past, present, and future, to allow the determination of absolute cellular molecular abundances of transcript i. Modeling, like that presented in S6 Appendix and S2 Fig, and like that of [17] could also play a vital role in estimating relative yield coefficients, especially if a wider array of synthetic spike-ins covering a large gamut of physical properties were designed and utilized. Our methods have the potential of facilitating statistical modeling of RNA counts because of the explicit relationship between our nominal abundances and absolute molecular, cellular abundances of RNA. In principle, variation in counts as a consequence of true biological variation in random attomoles, N, and variation in counts due to variation in relative yield coefficient across transcripts with nearly identical mean abundances, μN, could be disentangled. Our approach lays the groundwork for investigating, testing, and modeling how the physical properties—e. g. , length, GC content, folding energy—determine the relative yield coefficient of spike-ins and native RNA transcripts alike. Empirical measurements of relative yield coefficients, as we have defined them, and biophysical modeling could facilitate progress in making the connection between sequencing counts and the underling molecular cellular abundances of the corresponding transcripts. Our work follows up on and extends the work of [15,16,36,43,44]. Our inference method is linear and global for each library, like that of [19], [36] and [45]. We showed that our global (library specific) νj calibration constants are closely related to the Anders and Huber-like “technical” size factors of [36], which are based on spike-in counts. We called their normalization constants s j SI, and we showed that they are proportional to our νj normalization constants in the cases of 2 of our data sets with large library sizes, as predicted by theory (S8 Appendix). An important difference is that the s j SI calibration constants are on a dimensionless scale, on the order of 1, and do not allow one to infer absolute abundances of transcripts once their relative yield coefficients become known. [16] applied loess normalization to ERCC spike-in counts to derive a normalization function that they then applied to the counts corresponding to native RNA. Our analysis and rigorous testing of our theory and methods suggest that a local nonlinear transformation, such as loess normalization of the count data is not needed for our RNA-seq data. It seems likely that any local nonlinear fitting of counts to make replicate spike-in libraries as similar as possible would involve overfitting the data. Our work has some important features in common with the HTN method of [46], particularly, the assumptions underlying their Eq (1) and our Eqs S1 Appendix (2) and (3). These equations explicitly allow for differences in total RNA abundance across conditions. In addition, both normalization methods are global and linear. However the HTN method of [46]: relies on having de facto housekeeping genes rather than experimentally-added spike-ins; does not include a model for biological noise; assumes that relative yield is simply proportional to transcript length; is focused primarily on testing for differential gene expression; and does not provide estimates of absolute RNA abundance. Their global scale factor for a given library is determined by minimizing the sum over spike-ins of the square differences between the spike-in counts in that library and those of a library chosen to be the reference library. That scale factor is then used for the native RNA counts within the same non-reference library. It can be shown that this library-by-library normalization procedure, in the limit as library size (native RNA and spike-ins) approaches infinity, will give an abundance measure that is proportional to our z abundance values based on νj normalization. A quite different suite of normalization methods, called RUV (removal of unwanted variation), was introduced by [15,36,43,44] and applied with great effect to many different data sets. The methods involve singular value decomposition (SVD) variant of factor analysis to compute a factor matrix W, which is used to model nuisance sources of variation that are unrelated to the experimental design. The factor matrix W is included, in addition to a design matrix, in a generalized linear model for normalized counts. One qualitative way of thinking about the W matrix is that is adds columns to the original design matrix for explanatory variables that one didn’t originally know about. Although this method is widely effective at reducing unwanted variation in RNA-seq data, it does not allow one to infer absolute cellular molecular RNA abundances, even if the factor matrix is computed based on spike-ins or an invariant gene set (S8 Appendix), as the authors are well aware. The simple reason is that proportion of spike-in count is tightly correlated with the biological phenomenon of interest the change of total RNA abundance with condition. However, we showed that results of our maximum likelihood normalization method can be improved, with respect to clustering and detection of differential gene expression, by applying an an RUV method based on residual, RUVr (RUVSeq package [15]) after νj normalization. We obtained closely similar results by a simpler method involving a correction factor δj for each library that was based on our discovery in a dilution study with technical replicates that we seem to have some noise in the actual overall amount of spike-ins added to the cellular RNA. We tentatively ascribed these to dilution/volume errors in handling the stock spike-in mixture. This finding highlights the importance of replacing pipetting methods for handling the spike-ins with more accurate robotic methods. The continuing discovery of examples in which there are gross transcriptome differences between cellular states, has established a need for spike-in controls in RNA-seq experiments [19]. Despite some criticisms [15], external RNA spike-ins have been adopted in several recent studies alongside methods developed to use them for RNA-seq quantitation [16,19,36,46,47]. The model presented in this work lends itself for both absolute and relative RNA quantitation, dependent on the experimental ability to accurately isolate a fixed number cells for library preparation. In both cases, we offer evidence that our approach provides reproducible results in a wide variety of conditions and has a strong predictive power. In conclusion, the presented model allows for improved unbiased RNA-seq quantitation in any experimental setup using external RNA spike-ins.
We present a complete statistical model for the analysis of RNA-seq data from a population of cells using external RNA spike-ins and a maximum-likelihood method for genome-wide estimation of transcripts per cell. The model includes biological variability of cellular transcript number and sampling noise. We derive an unbiased estimator of transcripts per cell for every transcript, given by simply multiplying the count by a library-dependent, but transcript-independent, scale factor. This is a nominal estimate that can be converted to an absolute estimate by dividing by the transcript’s relative yield coefficient, measured in a separate experiment. A negative binomial probability mass function with novel normalization (size) factors allows for parametric testing of hypotheses concerning dependence of the absolute abundance of each transcript on experimental condition. Our method integrates information from every RNA-seq experiment across all replicates and experimental conditions to determine the calibration constants. We test the method with a dilution study and a k-fold cross-validation study. We illustrate our method with applications to two independent data sets from yeast and the sea squirt that were derived by different library preparation protocols. We show that our methods detect genome-wide amplification of expression, and we compare our method to others.
cdna libraries
sequencing techniques
statistics
rna extraction
fungi
mathematics
forms of dna
dna libraries
molecular biology techniques
dna
rna sequencing
cell enumeration techniques
cellular structures and organelles
extraction techniques
research and analysis methods
total cell counting
gene expression
complementary dna
molecular biology
ribosomes
yeast
biochemistry
rna
statistical models
eukaryota
ribosomal rna
cell biology
nucleic acids
genetics
biology and life sciences
physical sciences
non-coding rna
organisms
A complete statistical model for calibration of RNA-seq counts using external spike-ins and maximum likelihood theory