Statistical calibration of qRT-PCR, microarray and RNA-Seq gene expression data with measurement error models

The accurate quantification of gene expression levels is crucial for transcriptome study. Microarray platforms are commonly used for simultaneously interrogating thousands of genes in the past decade, and recently RNA-Seq has emerged as a promising alternative. The gene expression measurements obtained by microarray and RNA-Seq are, however, subject to various measurement errors. A third platform called qRT-PCR is acknowledged to provide more accurate quantification of gene expression levels than microarray and RNA-Seq, but it has limited throughput capacity. In this article, we propose to use a system of functional measurement error models to model gene expression measurements and calibrate the microarray and RNA-Seq platforms with qRT-PCR. Based on the system, a two-step approach was developed to estimate the biases and error variance components of the three platforms and calculate calibrated estimates of gene expression levels. The estimated biases and variance components shed light on the relative strengths and weaknesses of the three platforms and the calibrated estimates provide a more accurate and consistent quantification of gene expression levels. Theoretical and simulation studies were conducted to establish the properties of those estimates. The system was applied to analyze two gene expression data sets from the Microarray Quality Control (MAQC) and Sequencing Quality Control (SEQC) projects.

1. Introduction. The transcriptome of a cell is the entire set of mRNA molecules or transcripts produced by DNA transcription under certain biological or environmental conditions. Systematic profiling of the transcriptome cannot only provide a dynamic characterization of the cell's molecular constitution, but also shed light on gene functional annotation, regulatory mechanisms and transcriptional networks underlying various biological processes of the cell.
Since the mid-1990s, DNA microarray has served as the leading experimental platform for transcriptome study [Schena et al. (1995), Lockhart et al. (1996)]. Despite its huge success, microarray is known to suffer from some limitations such as reliance on existing knowledge of transcript sequences, high level background noise and limited dynamic range of detection.
A second platform called the quantitative real-time Reverse Transcription Polymerase Chain Reaction (qRT-PCR) can also be used in gene expression studies. Although well acknowledged as the most reliable gene expression measurement technology, the throughput of qRT-PCR is low, that is, the number of genes or transcripts that can be measured in a single qRT-PCR experiment is limited. In addition, qRT-PCR is also subject to variations caused by various biological, technical and experimental factors involved in qRT-PCR experiments [Bustin (2002), Bustin and Nolan (2004)]. Currently, qRT-PCR is used either for detecting and quantifying a small number of specific transcript targets or as a gold standard for validating hits or findings from other high-throughput platforms such as microarray [Applied Biosystems (2006)].
Recently, RNA-Seq has emerged as a new experimental platform for transcriptome profiling, and it is believed to overcome the major limitations of microarray [Wang, Gerstein and Snyder (2009)]. Although promising, data generated from RNA-Seq experiments still demonstrate excessive variability [Schwartz, Oren and Ast (2011)]. Therefore, RNA-Seq data need to be normalized before used in downstream transcriptome analysis. A number of methods have been proposed in the literature for normalizing RNA-Seq data, including RPKM [Mortazavi et al. (2008)], quantile-based procedures [Bullard et al. (2010)], TMM [Robinson and Oshlack (2010)], mseq [Li, Jiang and Wong (2010)], GPseq [Srivastava and Chen (2010)], POME [Hu et al. (2012)] and DESeq [Anders and Huber (2010)]. Nonetheless, the problem of how to statistically characterize and normalize RNA-Seq data has not been fully settled with satisfaction and demands further investigation [Mak (2011)].
The strengths and weaknesses of qRT-PCR, RNA-Seq and microarray can be summarized as follows. In terms of throughput capacity, RNA-Seq is the highest, microarray is the second, and qRT-PCR is the lowest, whereas, in terms of accuracy, the order from the highest to the lowest is qRT-PCR, RNA-Seq and microarray. An ideal platform for transcriptome profiling should combine the accuracy of qRT-PCR with the high-throughput capacity of RNA-Seq. Unfortunately, such a platform is not currently available yet. A natural question is whether it is possible to use statistical methods STATISTICAL CALIBRATION OF GENE EXPRESSION DATA 3 to combine the strengths of the three platforms and generate gene expression measurements of a higher quality at the genome-wide scale. The answer turns out to be positive, and the methodology that can be applied is statistical calibration.
Statistical calibration is typically used for the scenario where p instruments are available for measuring the same quantity. Among the instruments, one is more accurate but at the same time more expensive than the others. In order to make a less accurate instrument generate more reliable measurement results for general use, it needs to be calibrated by the most accurate instrument through statistical analysis. There are two different types of calibration, which are absolute calibration and comparative calibration. Absolute calibration assumes the most accurate instrument gives the true value of the targeted quantity, whereas comparative calibrate assumes the most accurate instrument is also subject to measurement error. The literature on statistical calibration is primarily focused on absolute calibration, and a variety of statistical methods have been developed; see Osborne (1991) for a comprehensive review. On the other hand, comparative calibration is more challenging than absolute calibration and is mainly discussed in the literature on measurement error models; see Fuller (1987) and Cheng and Van Ness (1999) for more thorough discussions.
The platforms for measuring gene expression levels using qRT-PCR, microarray and RNA-Seq represent three different instruments, with the qRT-PCR platform being the most accurate and the most expensive. In this article, we propose to use a system of measurement error (ME) models to model gene expression measurements obtained by the qRT-PCR, microarray and RNA-Seq platforms, and to further use the system to calibrate RNA-Seq and microarray measurements with qRT-PCR measurements. A two-step approach is used to estimate the parameters of the system of ME models, and the statistical properties of the resulting estimates are discussed. Through both theoretical and simulation studies, we show that the calibrated gene expression measurements are more consistent and accurate than those by any of the individual platforms. Furthermore, we apply the system of ME models to calibrate gene expression data generated by qRT-PCR, RNA-Seq and microarray from both the Microarray Quality Control (MAQC) and Sequencing Quality Control (SEQC) projects (detailed information about the data is given in Section 4.1), and show that the resulting calibrated measurements provide more accurate quantification of gene expression levels and lead to more discoveries in gene differential expression analysis.
In Section 2 we define the system of measurement error (ME) models and describe the two-step approach for estimating the parameters in the system. In Section 3 we present some simulation results. In Section 4 we describe the gene expression data from the MAQC and SEQC projects and discuss the results from applying the system of ME models to analyze the data. 2. Measurement error models and calibration. Let T denote the collection of all the genes or transcripts in the transcriptome under study. As discussed in the Introduction, RNA-Seq has the capacity to measure the expression levels of all genes in T , microarray can target thousands of genes, and qRT-PCR generally will measure at most hundreds of genes in a single experiment. Let A, B and C be the collections of genes measured by qRT-PCR, microarray and RNA-Seq, respectively. We assume that A ⊂ B ⊂ C ⊂ T , and the cardinalities of A, B and C are n, m and l, respectively. By the assumption, n genes (i.e., those in A) are measured by all three platforms, m − n genes (i.e., those in B − A) are measured by both microarray and RNA-Seq, and l − m genes (i.e., those in C − B) are measured only by RNA-Seq. We further assume that commonly used platform-specific methods are used to process and normalize raw data generated by each platform to produce gene expression measurement data. Some examples are the PLIER method [Affymetrix Inc. (2005)] for microarray, the RPKM method for RNA-Seq and the delta-delta C t method for qRT-PCR. Following the convention in transcriptome studies, the log-2 transformation is further applied to the normalized gene expression measurements of the three platforms, and we refer to the resulting normalized transformed values as the qRT-PCR, microarray and RNA-Seq gene expression measurements, respectively, in the rest of this article.
Each type of gene expression measurement data can be generated by one lab or multiple labs. Based on the number of labs involved in generating one type of gene expression measurement data, we distinguish two scenarios, which are the single-lab scenario and the multi-lab scenario. In the singlelab scenario, each type of gene expression measurement data (e.g., RNA-Seq data) is generated by a single lab. However, the expression measurement data of different types are not necessarily generated from the same lab. In the multi-lab scenario at least one type of gene expression data is generated by more than one lab. Due to limited space, we treat only the single-lab scenario in this section, and the multi-lab scenario will be discussed in Section S.7 of the supplementary material [Sun, Kuczek and Zhu (2014)].
2.1. Single-lab scenario. Suppose that the genes in A, B − A and C − B are labeled by integers from 1 to n, n + 1 to m, and m + 1 to l, respectively. We assume each platform generates one expression measurement for each gene. When technical replicates are available, their average value is used as the expression measurement of a gene. For 1 ≤ j ≤ n, let X j denote the expression measurement of gene j by qRT-PCR. Similarly, let Y j denote the expression measurement of gene j by microarray for 1 ≤ j ≤ m, and Z j the STATISTICAL CALIBRATION OF GENE EXPRESSION DATA 5 expression measurement of gene j by RNA-Seq for 1 ≤ j ≤ l. Note that X j , Y j and Z j are all in log-2 scales.
As discussed in the Introduction, gene expression measurements produced by the three platforms are all subject to measurement errors. We propose to use the following measurement error models to characterize X j , Y j and Z j , respectively, The above models are not simple linear models because none of the terms on the right-hand sides of the equations is observed, and terms β 2 µ j and β 3 µ j are cross terms of unknown quantities. Note that n < m < l, that is, the three platforms measure different subsets of genes. Next, we will discuss the terms in these three models in detail and then propose the estimation method.
In model (2.1a), µ j is the true expression level of gene j and ε 1j is the random error due to the qRT-PCR platform. Depending on the normalization method, qRT-PCR can lead to either absolute quantitation or relative quantitation of a gene's expression level [Pfaffl (2004), Chapter 3]. Therefore, the interpretation of µ j depends on whether absolute or relative quantitation is used in an experiment. For simplicity, in this article, we do not further distinguish these two quantitation methods and simply refer to µ j as the true expression value of gene j. We assume that for 1 ≤ j ≤ n, ε 1j 's are i.i.d. N (0, σ 2 1 ). The variance of X j is σ 2 1 , representing the reproducibility of the qRT-PCR platform. According to the model, the qRT-PCR measurement X j is unbiased with respect to µ j .
In model (2.1b), ε 2j is the random error due to the microarray platform. For 1 ≤ j ≤ m, we assume ε 2j 's are i.i.d. N (0, σ 2 2 ). The other two terms gives the mean measurement, that is, E(Y j ) = α 2 + β 2 µ j . In other words, the expectation of the microarray measurement of gene j is assumed to be a linear function of µ j , where α 2 and β 2 are the intercept and slope, respectively. The variance of Y j is σ 2 2 . Note that due to the differences in technologies and platform specific normalization methods, gene expression measurements from different platforms are in different scales and may have shifts. Therefore, α 2 and β 2 represent the shift and scale of the microarray measurements relative to the qRT-PCR measurements. When comparing the reproducibilities of the qRT-PCR and microarray platforms, Y j needs to be transformed to Y j = (Y j − α 2 )/β 2 , which is of the same scale as X j . The variance of Y j , which is σ 2 2 /β 2 2 , is referred to as the reproducibility of the microarray platform.
In model (2.1c), the mean of the RNA-Seq measurement is also assumed to be a linear function of µ j , that is, E(Z j ) = α 3 + β 3 µ j , and ε 3j is the random error due to the RNA-Seq platform. We assume that for 1 ≤ j ≤ l, ε 3j 's are i.i.d. N (0, σ 2 3 ). The intercept α 3 and the slope β 3 represent the shift and scale of RNA-Seq measurements relative to the qRT-PCR measurements. The variance of Z j is σ 2 3 . Similar to the microarray platform, we refer to σ 2 3 /β 2 3 as the reproducibility of the RNA-Seq platform. Different measurement platforms generally have different dynamic ranges of detection. This is also the case for the qRT-PCR, microarray and RNA-Seq platforms. The models proposed above may hold for only a range that fits all three platforms. Therefore, when applying the models in practice, a proper range of expression levels needs to be used, and genes with extremely low or high expression levels need to be excluded. Random measurement errors are typically heteroscedastic, that is, the variance of a random measurement error depends on the magnitude of the targeted quantity. Notice that the errors in the models above are assumed to homoscedastic. The justification for assuming homoscedastic errors is twofold. First, the gene expression measurements X j , Y j and Z j are the log-2 values of the original measurements produced by the three platforms, respectively, and the log-2 transformation is known to mitigate heteroscedasticity. Second, as previously discussed, the proposed models will be applied to genes with expression levels within a certain proper range, which further alleviates the concern of heteroscedasticity. Nonetheless, when applying the models, diagnostic analysis always needs to be performed.
When µ j 's in models (2.1a)-(2.1c) are assumed to be unknown fixed quantities, the measurement error models are said to be functional; and when µ j 's are assumed to be i.i.d. random variables, the models are said to be structural. In this article, µ j 's are gene expression values and considered fixed. Therefore, the models defined above are functional. Following the convention in the literature, we call the µ j 's incidental parameters k and the other parameters (α i , β i for 2 ≤ i ≤ 3, and σ 2 i for 1 ≤ i ≤ 3) structural parameters. In general, when p measurement platforms or instruments are used to measure the same quantity, it requires p measurement error models to characterize the measurement results. For ease of discussion, the p measurement error models are said to form a system of ME models of order p. When the measurement error models are structural, the system is said to be structural; and when the models are functional, the system is said to be functional. Hence, models (2.1a)-(2.1c) defined above form a functional system of order 3, and any two of them form a functional system of order 2.
2.2. Parameter estimation. The statistical literature on measurement error models is focused primarily on systems of order 2. Reiersøl (1950) showed that structural systems of order 2 are not identifiable under the nor-STATISTICAL CALIBRATION OF GENE EXPRESSION DATA 7 mality assumption, unless additional information such as the ratio between the variances of the two measurement instruments is available. For structural systems of order higher than 2 (e.g., p = 3), Barnett (1969) showed they become identifiable with no further information needed, and the maximum likelihood estimates (MLEs) of the structural parameters can be obtained.
Parameter estimation for functional systems of ME models is more difficult than that for structural systems due to the presence of the incidental parameters, µ j 's. First, under the normality assumption on the involved measurement errors, the likelihood function of a functional system becomes unbounded [Kendall and Stuart (1973)]. Solari (1969) showed that the likelihood function of a functional system of order 2 does not have a local maximum. We have found that this is also generally true for functional systems of order 3. Therefore, the maximum likelihood method fails to produce proper estimates for both structural and incidental parameters of a functional system of ME models.
The relationship between structural and functional systems was discussed in the early literature on measurement error models. Much attention has been given to the connection between the identifiability of structural systems and the existence of proper estimates of the structural parameters in functional systems. Gleser (1983) showed that when a structural system is identifiable, consistent estimates of the structural parameters in the corresponding functional model exist. Further, he suggested that in most cases, a good way to find a consistent estimator for the functional model is to directly verify whether the consistent estimator in the structural model is also consistent in the corresponding functional model.
Based on the discussion above, a two-step approach can be used to estimate the parameters of the functional system of order 3 defined above for the qRT-PCR, microarray and RNA-Seq measurements. In the first step, the genes with measurements from all three platforms (i.e., genes in A) are used to obtain the MLEs of the structural parameters under the assumption that these µ j 's are i.i.d. N (µ, σ 2 ). In the second step, the estimates obtained in the first step are used in place of their corresponding structural parameters in the functional ME system, and then the generalized least squares method is used to obtain estimates of the incidental parameters (or gene expression levels).
Structural parameters. We use the expression measurements of the genes in A to estimate the structural parameters. Let X . , Y . and Z . be the averages of X j , Y j and Z j over all genes in A. The sample variances and covariances for {X j } 1≤j≤n , {Y j } 1≤j≤n and {Z j } 1≤j≤n are denoted as S xx , S yy , S zz , S xy , S xz and S yz , respectively. Following the first step of the two-step approach discussed previously, the estimates of the structural parameters are given as Note that when deriving the above estimates, the nonnegativity constraints on the variance estimates were not enforced due to two considerations. First, we can enforce the constraints by using more complicated estimates such as the ones given in Carter and Fuller (1980). However, in practice, the constraints are usually satisfied automatically and, thus, the benefit of using the more complicated estimates is minimal. Second, gene expression data are often noisy, and the true values of the variance components are unlikely to be close to 0. The violation of the constraints by the variance components estimates indicates either the estimates are not reliable due to insufficient sample size or the model assumptions are invalid and need to be reconsidered. The enforcement of the constraints may miss the opportunity of identifying these potential pitfalls. For convenience, let θ = (α 2 , α 3 , β 2 , β 3 , σ 2 1 , σ 2 2 , σ 2 3 ) ⊤ be the vector of the structural parameters andθ the estimate of θ. The following proposition establishes the asymptotic distribution ofθ under the functional system of ME models as n goes to ∞. Proposition 1. Assume the functional system of ME models is true and the following limits exist, −→ means converge in distribution and Γ θ is the variance-covariance matrix with its explicit expression given in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.1.2.
The proof of Proposition 1 is outlined in Section S.1.1 of the supplementary material [Sun, Kuczek and Zhu (2014)]. The assumptions on the mean and variance of µ j 's in Proposition 1 appear to be reasonable, because µ j 's are the true expression levels of genes in a given cell line or tissue sample.

9
Note that ∆ appears in the denominator of the variances ofα 2 ,α 3 ,β 2 and β 3 . Therefore, the estimates of the structural parameters will attain higher accuracy when µ j 's in A are more spread out. The asymptotic property of θ is stated as n goes to ∞. In practice, this property will hold approximately when n is sufficiently large. In our simulation study, we found that when n ≥ 150, the Mean Square Errors (MSE) ofθ are close to zero (see Figure S1 and Figure S2 in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.4) and the standard errors ofθ are within 10% of the true parameter values.
Incidental parameters. How to best estimate the incidental parameter µ j depends on whether gene j is in A, B − A or C − B. We consider these three cases separately. For each case, we follow the second step of the two-step approach discussed previously. First, we replace the structural parameters by their estimates, and then we apply the generalized least squares method to obtain the estimate of µ j .
Depending on whether gene j is in A, B − A or C − B, given the structural parameters, the measurement error models involving µ j can be rewritten as a linear model of µ j (see the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.2). Replacing the structural parameters θ by their estimatesθ and applying the generalized least squares method, the estimates of µ j for j in A, B − A and C − B are given as follows: We callμ xyz j ,μ yz j andμ z j the calibrated estimates of µ j for j in A, B − A and C − B, respectively. The quality of the calibrated estimates of µ j depends on how accurate the estimates of the structural parameters are, which further depends on n, the number of genes in A. The properties ofμ xyz j ,μ yz j andμ z j as n goes to ∞ are given in the following proposition.
Proposition 2. Assume the functional system of the ME models is true and the limits in (2.3) exist.
From Proposition 2, as n goes to ∞, the calibrated estimatesμ xyz j ,μ yz j andμ z j are asymptotically unbiased estimates of µ j for j in A, B − A and C − B, respectively. The variances of the calibrated estimates do not converge to zero asymptotically, instead, they converge to γ A , γ B−A and γ C−B , respectively. When n is sufficiently large, the second terms in (2.5a), (2.5b) and (2.5c) are negligible, and the variances of the calibrated estimates of µ j for j in A, B − A and C − B can be approximated by γ A , γ B−A and γ C−B , respectively. When n is small or moderate, the second terms in the above expansions may not be negligible.
It is clear that γ B−A is larger than γ A , and γ C−B is larger than γ B−A and γ A . The order implies that when n is sufficiently large, the calibrated estimate of µ j based on only the measurements of RNA-Seq is less accurate than the calibrated estimate of µ j based on the RNA-Seq and microarray measurements, which, further, is less accurate than those calibrated estimates based on all three platforms.
We further compare the three types of calibrated estimates with the qRT-PCR measurement X j , microarray measurement Y j and the RNA-Seq measurement Z j . Note that X j , which is presumably the most accurate measurement among X j , Y j and Z j , is an unbiased estimate of µ j , and its reproducibility is σ 2 1 . Because γ A is less than σ 2 1 ,μ xyz j is more accurate than X j when n is sufficiently large. Therefore, by combining the measurements from all three platforms, we can obtain a more accurate estimate of the expression level of gene j. Because γ B−A is less than either of σ 2 2 /β 2 2 and σ 2 3 /β 2 3 , which are the reproducibilities of Y j and Z j , respectively,μ yz j is more accurate than Y j and Z j when n is large. γ C−B is equal to σ 2 3 /β 2 3 , which is the reproducibility of Z j . Therefore,μ z j has the same reproducibility as the RNA-Seq measurement Z j . The advantage for usingμ z j instead of Z j is that µ z j is in the same scale asμ xyz j andμ yz j so that the calibrated estimates in C − B are scale compatible with those in A and B − A.
The calibrated estimates {μ xyz j : j ∈ A}, {μ yz j : j ∈ B − A} and {μ z j : j ∈ C − B} can also be considered gene expression measurements normalized by the qRT-PCR platform. As discussed above, the standard errors of these three types of calibrated estimates are different. Let ϕ j,A , ϕ j,B−A and ϕ j,C−B denote the first order expansions of the variances ofμ xyz j ,μ yz j andμ z j in (2.5a), (2.5b) and (2.5c), respectively. By replacing the structural parameters with their corresponding estimates, we obtain the estimated variances and standard errors of the three types of calibrated estimates, which areφ j,A and φ j,A ,φ j,B−A and φ j,B−A , andφ j,C−B and φ j,C−B , respectively. When used for gene expression analysis such as detecting differentially expressed genes, both the calibrated estimates and their standard errors need to be used.

3.
Simulation. An extensive simulation study was conducted to evaluate the performance of the proposed calibration method. In this section we present the simulation results under the single-lab scenario. We first report the simulation results on the accuracy of the calibrated estimates or calibrated expression levels, and then we present the performance of the calibrated estimates when used in gene differential expression (DE) analysis.
3.1. Accuracy of estimates. In this section we report the simulation results for three settings of the system of ME models under the single-lab scenario. The structural parameters of the three settings are listed in Table 1. Setting 1 is set to resemble the results from real data analysis in Section 4. In setting 2 the reproducibilities of microarray and RNA-Seq are mildly worse than that of qRT-PCR technology; and in setting 3 the reproducibilities of microarray and RNA-Seq are much worse that that of qRT-PCR. The incidental parameters (i.e., µ j 's) for each model setting were randomly drawn from N (0, 25). For each model setting, we independently generated a training data set of 300 genes and a testing data set of 1000 genes. Both data sets contain the measurements of the genes by all three platforms. The measurements of the first n genes in the training set were used to estimate the structural parameters (n was varied from 20 to 300), and then the three types of calibrated estimates of the expression levels of the genes in the testing set were obtained. For each combination of model  setting and n, this procedure was repeated 200 times to calculate the MSEs of the estimates for the purpose of performance evaluation and comparison.
We first examine and compare the variances of the calibrated expression levelsμ xyz j ,μ yz j andμ z j . Recall that Var(μ xyz j ) = γ A + n −1 ω A (θ, µ j ). When n is large, γ A is the dominant term and the other higher order terms are negligible; but when n is moderate (e.g., n ∈ [20, 50]), the second term n −1 ω A (θ, µ j ) is not negligible. Further notice that γ A depends on the structural parameters only, but n −1 ω A (θ, µ j ) also depends on the incidental parameters, particularly (µ j −μ) 2 as shown in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.3. This implies that the variance ofμ xyz j increases as µ j deviates away from the meanμ when n is moderate. The variances ofμ yz j andμ z j behave in the same way as the variance ofμ xyz j . Figure 1 from bottom to top shows the scatter plots of the sample variances ofμ xyz j ,μ yz j andμ z j versus the incidental parameters µ j for the first model setting and three different sizes of training data subsets (black for n = 20, blue for n = 50, and red for 300). For convenience, we refer to the curve generated by plotting the variance of a type of estimate against the incidental parameter as the variance curve of the type of estimate. The bottom three curves are the variance curves ofμ xyz j for n = 20, 50 and 300, respectively, the middle three curves are those ofμ yz j , and the top three curves are those ofμ z j . For each type of estimate, the variance curve for n = 20 is above that for n = 50, which is above that for n = 300, indicating that as n increases, the variance of the respective estimate decreases. Furthermore, all three variance curves for n = 20 demonstrate a stronger quadratic pattern than those for n = 50 and n = 300, indicating their dependence on (µ j −μ) 2 . The variance curves for n = 50 show a much mitigated quadratic pattern, and the variance curves for n = 300 become flat. Comparing the variance curves, it is clear that under the same n, the variance curve ofμ xyz j is lower than that ofμ yz j , which is lower than that ofμ z j , indicating that, in terms of accuracy, the order of the three types of estimates, from the best to the worst, isμ xyz j ,μ yz j andμ z j . These results confirm the properties of the variances of the three types of estimates discussed in Section 2.

STATISTICAL CALIBRATION OF GENE EXPRESSION
We further examine the average MSEs of the three types of calibrated estimates over all of the genes in the testing data set. For convenience in discussion, we refer to the scatter plot of the average MSE of a type of estimate over all of the genes in the testing data set versus the training data size n as the aMSE curve of the type of estimate. Note that X j is also an unbiased estimate of µ j , and we also compare the average MSEs of X j with the three types of estimates. Figure 2 demonstrates the aMSE curves of the three types of estimates {μ xyz j }, {μ yz j } and {μ z j } as well as that of the original qRT-PCR measurements {X j }. Plots (a), (b) and (c) in Figure 2 from the left panel to the right panel correspond to the three model settings given in Table 1, respectively.
The aMSE curves of {X j } are flat because they do not depend on n. Overall, the aMSE curves of {μ xyz j }, {μ yz j } and {μ z j } decrease as n increases, and they become flat after n is greater than 100, indicating that the structural parameters are accurately estimated and the benefit of further increasing n in the training data set becomes negligible.
In all of the plots, the aMSE curves of {μ xyz j } are always the lowest, indicating that {μ xyz j } give the most accurate quantification of the gene expression levels. The aMSE curves of {μ yz j } are always below the aMSE curves of {μ z j }, indicating that {μ yz j } is more accurate than {μ z j }. However, the aMSE curve of {μ yz j } is not always below that of {X j }. Plots (a), (b) and (c) demonstrate three different cases where when n is sufficiently large (e.g., ≥100), the aMSE curves of {μ yz j } are below, close to and above that of {X j }. In general, the performance of {μ yz j } relative to {X j } depends on the model settings. In all of the plots, the aMSE curves of {μ z j } are always above that of {X j }, indicating that the calibrated estimate based on RNA-Seq measurement alone cannot outperform the qRT-PCR measurement.
As discussed in Section 2, the structural parameters are fundamentally different from the incidental parameters. As n increases, the estimates of the structural parameters will converge to their respective targets as shown in Proposition 1. The simulation results about the structural parameters can be found in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.4.
3.2. Differential expression. In this section we report the simulation study on the performance of calibrated estimates when used in DE analysis. In order to make the simulation study convincing, we simulated the RNA-Seq data from a popularly used simulator called the Flux Simulator [Griebel et al. (2012)]. The Flux Simulator mimics the pipeline of RNA-Seq experiments and generates RNA-Seq short reads. It allows researchers to investigate the properties of RNA-Seq data and the analysis tools in silico. In this study, we considered two biological conditions, which are referred to as conditions 1 and 2. Totally, 5000 genes were considered, and among them 500 genes were set to be differentially expressed. Let µ 1j and µ 2j denote the log-2 values of the true expression levels of gene j under the two conditions, respectively, and d j = µ j1 − µ j2 the difference between µ 1j and µ 2j . Under condition 1, µ 1j 's were randomly assigned by the Flux Simulator using its default setting. When gene j is not differentially expressed between the two conditions, we set d j to be 0; otherwise, d j is randomly generated from the uniform distribution on [−2, −0.5] ∪ [0.5, 2]. Then under condition 2, µ 2j 's were set to be µ 1j + d j . For each biological sample, 5 million reads were generated and mapped back to the reference genome. Finally, the log-2 RPKM values were calculated and used as the RNA-Seq expression measurements. qRT-PCR and microarray data were generated from the system of ME models, with the true expression levels on the two conditions set the same as those used for the RNA-Seq data. The structural parameters related to qRT-PCR and microarray (i.e., α 2 , β 2 , σ 1 and σ 2 ) were set the same as setting 1 in Table 1. We generated qRT-PCR data for 500 genes and microarray data for 3000 genes. Therefore, in this simulation study, there are 500 genes in set A, 3000 genes in set B, and 5000 gene in set C.
After generating the expression measurements, we applied the system of ME models for each biological condition. The calibrated expression levels and their variances areμ xyz j andφ j,A ,μ yz j andφ j,B−A , andμ z j andφ j,C−B , for j ∈ A, B − A and C − B, respectively. In what follows, for convenience, we do not further distinguish the notation for these three types of calibrated estimates. Instead, we useμ j1 andμ j2 to denote the calibrated expression levels of gene j under condition 1 and condition 2, respectively, and useφ j1 andφ j2 to denote the estimated variances ofμ j1 andμ j2 , respectively.
We performed DE analysis for the two conditions using the z-test based on the RNA-Seq measurements (Z j 's) and the calibrated measurements (μ j 's) separately, and compared the results. When performing the z-test based on the calibrated measurements for gene j, the p-value for testing H 0 : µ j1 = µ j2 was where Z follows the standard normal distribution. Based on the calculated p-values, the standard Benjamini-Hochberg procedure [Benjamini and Hochberg (1995)] was used to identify differentially expressed genes at a given false discovery rate (FDR). The z-test based on the RNA-Seq measurements was similarly conducted. Figure 3 depicts the Receiver Operator Characteristic (ROC) curves comparing the performances of the calibrated expression measurements and the RNA-Seq expression measurements in DE analysis. The ROC curve of the calibrated expression levels (black line) is above that of the RNA-Seq expression measurements (red line), indicating the DE results based on the calibrated expression measurements have higher true positive rate than those based on the uncalibrated RNA-Seq expression measurements at any level of false positive rate. For example, at the false positive rate level of 0.05, the true positive rate of the uncalibrated RNA-Seq measurements was 0.565, while that of the calibrated expression measurements was 0.641. Therefore, the calibrated expression measurements outperform the uncalibrated RNA-Seq measurements in this example.

Application and results.
4.1. Data sets. We applied the system of ME models to analyze the qRT-PCR, RNA-Seq and microarray gene expression data of two RNA samples, namely, the human brain reference RNA sample (or, in short, the Brain sample) and the human universal reference RNA sample (or the UHR sample), generated by the Microarray Quality Control (MAQC) and Sequencing Quality Control (SEQC) projects.
The qRT-PCR data were generated by one single lab using TaqMan R Gene Expression Assays and can be downloaded from Gene Expression Omnibus (GEO) with the series number GSE5350 (ftp://ftp.ncbi.nih.gov/ pub/geo/DATA/supplementary/series/GSE5350/). The qRT-PCR data were normalized by the delta-delta Ct method and originally contain the qRT-PCR measurements of 1001 genes; see the supplementary material Section S1.2 of Bullard et al. (2010) for detailed information. Among the 1001 genes, 7 genes have multiple entries with distinct expression values but under the same RefSeq ID. To avoid ambiguity, these genes were removed. Each gene has 4 technical replicates for each of the UHR and Brain samples.
The RNA-Seq data include measurements generated from two RNA-Seq experiments conducted in two different labs using the Illumina Genome Analyzer. The data from the first lab [Bullard et al. (2010)] can be downloaded from the NCBI Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih. gov/sra/) under the accession number SRA010153, and the data from the second lab can be downloaded from the same website under the accession number SRA008403. For convenience, we will use the accession numbers to refer to these two data in the rest of the article. The RPKM value is calculated lane-by-lane for each gene. For each of the Brain and UHR samples, those genes that received no reads in at least one lane of the two RNA-Seq experiments were excluded. See the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.7, for more information about the two RNA-Seq data sets.
The microarray data were generated by five different microarray experiments conducted in five different labs using Affymetrix U133 Plus2.0. For convenience, we label the labs as MA1 to MA5 in the rest of the article. The original probe-level data can be downloaded from Gene Expression Omnibus (GEO) with the series number GSE5350 (ftp://ftp.ncbi.nih.gov/pub/ geo/DATA/supplementary/series/GSE5350/). Data from each lab have 5 replicates for both UHR and Brain samples, and was normalized using the PLIER [Affymetrix Inc. (2005)] method. Detailed information about these data sets can be found in the website of the MAQC projects (http://www. fda.gov/). For each replicate in each lab, the average probe-level measurement of a gene is considered the gene's expression level intensity.
We integrated the three different types of gene expression data as follows. First, the log-2 transformation was applied to all three types of gene expression data. Then for each gene and each lab, the mean measurement across technical replicates is used as the gene's expression value. Furthermore, we exclude genes with extremely low or extremely high expression levels (i.e., those with qRT-PCR expression measurements below −6 or above 4 in log-2 scale) in A; the remaining genes in A are used to estimate the structural parameters. For the Brain sample, there are 409 genes with expression data from all three platforms (A), 6419 genes with only RNA-Seq and microarray measurements (B − A), and 5949 genes with only RNA-Seq measurements (C − B); and for the UHR sample, there are 477 genes with measurements by all three platforms (A), 7187 genes with measurements only by RNA-Seq and microarray (B − A), and 6892 genes with only measurements by RNA-Seq (C − B).
Two schemes were used to analyze the integrated data. First, we considered the single-lab scenario, in which each platform has data from one lab. In total, there are ten possible combinations, and we applied the system of ME models (2.1a)-(2.1c) to each combination. The analysis results of individual combinations were similar, and we report only the results for the combination that includes the RNA-Seq data from SRA010153 and the microarray data from MA1. Second, we applied the system of ME models for the multi-lab scenario to analyze the entire data set for each RNA sample. Due to limited space, the results from the multi-lab scenario are presented and briefly discussed in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.8.

4.2.
Diagnostics of model assumptions. The system of ME models (2.1a)-(2.1c) imposes the normality and homoscedasticity assumptions on the measurement errors. We checked these assumptions using genes in A.
After the system of ME models was fitted, we calculated the residuals by e 1j = X j −μ j , e 2j = Y j −α 2 −β 2μj and e 3j = Z j −α 3 −β 3μj corresponding to the measurement errors due to the qRT-PCR, microarray and RNA-Seq platforms, respectively. To check the normality assumption, we generated the QQ plots for the residuals and did not detect significant violation of the assumption. To check the homoscedasticity assumption, we generated residual plots and constructed approximate 95% confidence intervals of the Box-Cox transformation. Because the diagnostic results are similar for the two RNA samples, we present only those from the UHR sample as an example. The residual plots corresponding to the qRT-PCR, microarray and RNA-Seq platforms are presented in Figure S3 in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.8. The plots do not demonstrate strong heteroscedastic patterns. The 95% confidence intervals of Box-Cox transformation for the residuals from the three platforms are presented in Figure S4 in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.8. The confidence intervals corresponding to the qRT-PCR and RNA-Seq platforms both contain 1, and 1 is on the boundary of the confidence interval corresponding to the microarray platform. These results together indicate that there does not exist significant violation of the homoscedasticity assumption imposed on the platform measurement errors. 4.3. Structural parameters. The estimates of the structural parameters for the Brain and UHR samples are given in Table 2. The standard errors of the estimates are also reported in parentheses in the table. The estimates of the structural parameters can be used to compare the three platforms in terms of the quality of measurements they provide.
From models (2.1a)-(2.1c), the qRT-PCR measurements are unbiased with respect to µ j 's, and the intercepts α 2 and α 3 represent the shifts of In both samples,σ 2 1 <σ 2 3 /β 2 3 <σ 2 2 /β 2 2 , indicating that qRT-PCR has the best reproducibility, microarray has the worst reproducibility, and RNA-Seq is slightly worse than qRT-PCR but much better than microarray. 4.4. Incidental parameters. After the estimates of the structural parameters were obtained, we used the formulas obtained in Section 2 to estimate the gene expression levels, and the standard errors of the calibrated gene expression values were also calculated. These calibrated gene expression values together with their standard errors are expected to have a higher quality than the original measurements by the three platforms and lead to better results in downstream gene expression analysis.
We compare the calibrated gene expression levels with their original qRT-PCR, microarray and RNA-Seq measurements, and use the genes in the Brain and UHR samples that have measurements by all three platforms as an illustrative example.  From Figure 4, in the Brain sample, in terms of the linear correlation coefficient with the best measurements {μ xyz j }, microarray is the worst (0.8112) and qRT-PCR is the best (0.9403) among the three platforms. This is consistent with the findings reported in the literature that in terms of the quality of the gene expression data produced by the three platforms, qRT-PCR is most accurate and microarray is least accurate. The calibrated expression levels {μ yz j } have a higher correlation coefficient with {μ xyz j } (0.9551) than the qRT-PCR, microarray and RNA-Seq measurements. In the UHR sample, in terms of the linear correlation coefficient with the best measurements {μ xyz j }, again microarray is the worst (0.8553) and qRT-PCR is the best (0.9583) among the three platforms. The calibrated expression levels {μ yz j } have a higher correlation coefficient with {μ xyz j } (0.9645) than qRT-PCR, microarray and RNA-Seq measurements. 4.5. Gene differential expression. To demonstrate that calibrated expression levels can lead to better gene differential expression results, we used the RNA-Seq measurements (Z j 's) and the calibrated measurements (μ j 's), separately, to perform gene DE analysis for the Brain and UHR samples.
We carried out the DE analysis using the same z-test procedure as described in Section 3.2. Controlling the FDR at 0.01, the numbers of identified genes by the two types of measurements are reported in Table 3. In total, 331 genes were detected to be differentially expressed in the UHR and Brain samples using the calibrated expression measurements, whereas only 158 genes were detected using the original RNA-Seq measurements. These two groups of genes share 153 genes. Therefore, the calibrated expression measurements led to the detection of almost all the genes (153 out of 158 genes) detected by the original RNA-Seq measurements. In addition, the former detected 173 more genes than the latter. The total number of genes detected by each type of measurements was further broken down according to whether a gene belongs to A, B − A or C − B, and the results are also reported in Table 3. From the table, 27, 132 and 14 more genes were detected by the calibrated expression measurements in A, B − A and C − B, respectively, than by the original RNA-Seq measurements. In this case, there does not exist a gold standard to further verify the selected genes by the calibrated expression measurements and the RNA-Seq expression measurements. However, based on the simulation study in Section 3.2, we believe that the DE analysis based on the calibrated measurements has higher true positive rate and can lead to more discoveries than the uncalibrated RNA-Seq measurements. A list of the identified genes by the calibrated expression measurements is given in the supplementary material [Sun, Kuczek and Zhu (2014)], Section S.9.
5. Discussion. A system of functional measurement error (ME) models was proposed to calibrate the microarray and RNA-Seq measurements of gene expression levels by qRT-PCR. Due to limited space, the design issue of the proposed approach was not discussed in this article. The success of the proposed approach hinges on the genes that are measured by all three platforms, and the major bottleneck is the relative low throughput of the qRT-PCR platform. Therefore, the design issue is centered on the qRT-PCR platform with respect to two questions. The first question is how many genes should to be measured by qRT-PCR, and the second question is which genes should be measured. For the first question, based on the theoretical, simulation and real data application results in this article, it seems that at least 150 genes are needed to ensure that the bias and variances (i.e., the structural parameters) can be accurately estimated and the calibrated gene expression levels (i.e., the incidental parameters) can reach their best possible accuracy. Vendors of the qRT-PCR platform such as Life Technologies now offer assays and services to measure a sufficiently large number of genes simultaneously. This makes our proposed approach feasible in practice. For the second question, based on our theoretical results (e.g., Proposition 1), the true expression levels of the genes selected to be measured by qRT-PCR should be as spread out as possible. These two questions and the design issue in general will be addressed in a future publication.
As discussed in the Introduction, RNA-Seq data are also found to be subject to excessive variability and various methods have been proposed to normalize RNA-Seq data. In this article, only the RPKM method and the resulting RPKM measurements were used as the RNA-Seq measurements. Clearly, other normalization methods and their resulting measurements can also be considered. The variance component due to the RNA-Seq platform estimated under the multi-lab scenario (i.e., ς 2 2 ) in this article corresponds to the RPKM method. If another normalization method (e.g., the "mseq" method) is used, then ς 2 2 will correspond to that normalization method. Therefore, the system of ME models can be used to compare different normalization methods.
The functional system of ME models of order 3 can be extended to that of order p > 3, and the two-step approach to parameter estimation can also be extended in a straightforward manner. In this article, the bias of the measurement of a platform is assumed be a linear function of the true expression level. In general, nonlinear models or even nonparametric models can be considered. Furthermore, covariates can also be incorporated into the system of functional ME models. These possible extensions of the proposed approach will be investigated in the near future.