Robust background normalization method for one-channel microarrays

Multiple-probe robust gene expression index: multi-RGX

As described in the previous section, the major aim of our study is to suggest a new gene expression index. Moreover, we aim that the new index can estimate the true signals more accurately, since we model not only gene specific intensities, but also both gene and probe intensities. Finally, multi-RGX covers a wider range of distributions for the signals, hence, does not ignore any information in the data. Accordingly, it suggests a robust method to model the true signals.

Hereby, in this novel model, called as multi-RGX (multiple-probe robust gene expression index), we assume that the intensities of the PM and MM probes on the logarithmic scale have the density from the LTS, similar to the idea of RGX. On the other hand, different from this, it does not follow an iterative process in the optimization of the likelihood function, rather, it works with the explicit estimates of the desired parameters.

In a standard Affymetrix chip, accepting that there are m probes, and in each probe, there exist n genes coming from the PM and MM probes, it is assumed that both PM and MM have true signals in such a way that PM covers the complete true signals, whereas, MM includes only a fraction of them. To define the amount of the underlying true signals in MM, a fraction term p is defined. Hereby, the intensities coming from the PM probes are distributed as LTS with a mean containing a constant term and another term, which stands for the true signal changing per gene. Additionally, its variance is constant over genes and probes. Indeed, this assumption may be considered strong in modelling. However, from the application of this assumption in FGX and RGX indices, it has been found that it is competitive in terms of the accuracies of the signals with respect to the gene specific variance as BGX implements. Moreover, from alternative modellings via the gene specific variance, it has been observed that none of the estimators can have explicit solutions under a gene specific variance [14], and the estimates can be obtained via optimization steps which cause an increase in the computational time. On the other hand, different from PM, MM has a mean containing the same constant term and the part of the true signal. Hence, in modelling of the PM and MM intensities for each probe j (j=1, …, m) and each gene i (i=1, …, n), we consider the following distributional assumption on the log-scale:

where S_i shows the true signal for the gene i and μ_H denotes the constant background signal. Furthermore, σ² indicates the constant variance and p stands for the fraction value within the range of zero to one that implies the ratio of the true signal in the MM probes. Thereby, the corresponding likelihood of the model in Equation (1) is found via

L(Si, μH, p|a, b)=∏i = 1n∏j = 1mf(PMij)f(MMij),L=∏i = 1n∏j = 1m12πσ2e− (PMij − Si − μH)22σ212πσ2e− (MMij − pSi − μH)22σ2,

which is proportional to

where v≥4 and k=2v−3.

Then, the first derivatives of the function lnL are taken with respect to each parameter as stated below:

But it is observed that the underlying partial derivatives cannot produce unique solutions due to the repeated nonlinear functions in their formulas. Therefore, we suggest to use the modified maximum likelihood approach (MMLE) [15] that is based on approximating the nonlinear functions via the first order Taylor expansions around the ith population quantile t_(i) of the student-t distribution with the (2v−1) degrees of freedom. In the calculation, we apply the order statistics and define the partial derivations of the objective log-likelihood function in Equation (2) in terms of small linear sub-pieces. Moreover, it is found that the estimators based on MMLE are asymptotically efficient and equivalent to MLE [15, 16].

In this study, we observe that the nonlinearity in the partial derivatives is caused by the following expressions such that

z PMij=(PMij−Si−μH)σ and zMMij=(MMij−pSi−μH)σ.

Then, by approximating the common nonlinear functions g(z)=z1+z2k for PM as

g(z PMij)=(PMij−Si−μH)σ1+(PMij−Si−μH)2kσ2

and for MM as

g(z MMij)=(MMij−pSi−μH)σ1+(MMij−pSi−μH)2kσ2;

their first-order Taylor expansions can be derived by

g(z PMi(j))=αj+βjzPMi(j) and g(zMMi(j))=αj+βjzMMi(j),

in which

αj=2tj3k(1+tj2k) and βj=(1−tj2k)(1+tj2k)2.

In these expressions, zPMi(j) and zMMi(j) show the ordered probes in increasing magnitude for each gene i in the PM and MM standardized intensities, respectively. Accordingly, g(zPMi(j)) and g(zMMi(j)) are their associated linearized functions via the Taylor series. Thereby, the partial derivatives of MMLE are as follows:

Finally, by setting Equations (7)–(10) to zero, we get the optimal MML estimates of parameters.

As a result, μ^H is found as:

On the other side, from Equation (9), we obtain the estimate of the true signal for each gene i as:

However, by substituting Equation (12) into Equation (11), we can define μ^H in an alternative way as follows:

where PM ¯.j=∑i = 1nPMi(j)n and MM ¯.j=∑i = 1nMMi(j)n.

The estimate of the common fraction term p is computed by four roots and the optimal one which maximizes the log-likelihoods is found as the following equation:

p^=(SSMM −SSPM )+(SSPM −SSMM )2+4SSPM,MM 22SSPM,MM ,

where

SSPM =∑i = 1n∑j = 1mβj[∑j = 1mβjPMi(j)−∑j = 1mβjPM ¯.j]2,SSMM =∑i = 1n∑j = 1mβj[∑j = 1mβjMMi(j)−∑j = 1mβjMM ¯.j]2,SSPM,MM =∑i = 1n∑j = 1mβj[(∑j = 1mβjPMi(j)−∑j = 1mβjPM ¯.j)    ×(∑j = 1mβjMMi(j)−∑j = 1mβjMM ¯.j)]

Finally, from Equation (10), we can get MMLE of σ by:

σ^=B+B2+4nmCnm,

Adjusting the degree of freedom as 2nm−(n+2):

σ^=B+B2+4nmCnm−n−2,

in which

B=vk[∑i = 1n∑j = 1mαj(PMi(j)−MMi(j))],C=vk[∑i = 1n∑j = 1mβj(PMi(j)−S^i−μ^H)2 + ∑i = 1n∑j = 1mβj(MMi(j)−p^S^i−μ^H)2].

Observed Fisher information matrix and estimators for variances and covariances

The variance of the asymptotic maximum likelihood estimator can be found from the Fisher information matrix [17] and as stated previously that MMLE is asymptotically equivalent to it, resulting in the same feature of estimators.

Accordingly, the Fisher information matrix I is generated from the second partial derivatives of the loglikelihood function with respect to each parameter as stated below:

I11=2vkσ2∑i = 1n∑j = 1m[1−(PMi(j)−Si−μH)2kσ2(1+(PMi(j)−Si−μH)2kσ2)2  +1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2],I22=2vkσ2∑i = 1n∑j = 1mSi21−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2,Iii=2vkσ2∑j = 1m1−(PMi(j)−Si−μH)2kσ2(1+(PMi(j)−Si−μH)2kσ2)2  +2vp2kσ2∑j = 1m1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2 ,I12=I21=−∂2l∂μH∂p=2vkσ2∑i = 1n∑j = 1mSi1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2 ,I1i=Ii1=2vkσ2∑j = 1m1−(PMi(j)−Si−μH)2kσ2(1+(PMi(j)−Si−μH)2kσ2)2+2vpkσ2∑j = 1m1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2,I2i=Ii2=−2vkσ∑j = 1m(MMi(j)−pSi−μH)σ(1+(PMi(j)−Si−μH)2kσ2)2+2vpkσ2∑j = 1mSi1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2,Iik=Iki=−∂2l∂Si∂Sk=0,

where i, k=1, …, n, present the gene and j=1, …, m, shows the probe indicator, respectively. Moreover, l stands for the loglikelihood function lnL. Then, the I matrix is derived via

I=[∂2l∂μH2∂2l∂μH∂p∂2l∂μH∂S1…∂2l∂μH∂Sn∂2l∂p∂μH∂2l∂p2∂2l∂p ∂S1…∂2l∂p∂Sn∂2l∂S1∂μH∂2l∂S1∂p∂2l∂S12…0……………∂2l∂Sn∂μH∂2l∂Sn∂p0…∂2l∂Sn2].

In order to find the variance-covariance matrix, we partition the above matrix into the four sub-matrices. The details of this partition can be found in [18] and its application can be seen in [2, 3]. Once I is reorganized with respect to its sub matrices, the variances of the model parameters can be derived from simultaneous solutions of partial derivatives from the reorganized I. For an illustration, we show the variance of μ^H as follows:

V(μ^H)=1C0[2vkσ2∑i = 1n∑j = 1mSi21−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2− ∑i = 1n−2vkσ∑j = 1m(MMi(j)−pSi−μH)σ(1+(PMi(j)−Si−μH)2kσ2)2+2vpkσ2∑j = 1mSi1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)22vkσ2∑j = 1m1−(PMi(j)−Si−μH)2kσ2(1+(PMi(j)−Si−μH)2kσ2)2+2vp2kσ2−∑i = 1n1−(MMi(j)−pSi−μH)2kσ2(1+(MMi(j)−pSi−μH)2kσ2)2],

in which

C0=(−∂2l∂μH2−∑i = 1n(∂2l/∂Si∂μH)2− ∂2l/∂Si2)× (−∂2l∂p2−∑i = 1n(∂2l/∂Si∂p)2− ∂2l/∂Si2)− (−∂2l∂p∂μH−∑i = 1n(− ∂2l/∂Si∂μH)(− ∂2l/∂Si∂p)− ∂2l/∂Si2).

Application via a real dataset

In the analysis, we use a benchmark Affymetrix spike-in dataset. This dataset consists of 11 spike-in genes coming from four bacterial ancestor genes and 59 arrays with 10,864 probe pairs.

For the analysis, the 16 spike-in probe pairs are taken. Here, the gene expression values of the individual cRNA fragments, which are hybridized to U95A GeneChip arrays at the same concentration are analyzed under 14 different concentration levels [0.0, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, and 1024.0 pM (picoMolar)]. Accordingly, every gene is described by 14 probes in each array, resulting in 224 observations per array.

Finally, for the analysis, we use the R programme language with our original codes for all calculations.

Hereby, in the estimation of the model parameters for each array, we check the possible shape parameters for the LTS distribution from 2 to 40 with the jump size 0.05. In the computation, we find that most of the array can be optimally defined under the shape parameter 40 or close to 40, which practically indicates the normal distribution.

As presented in Figure 1A, while the level of concentrations increases, a nonlinear pattern is seen in the structure of the estimated signals, like other methods. Figure 1B shows the associated graph via all other methods.

Figure 1:

Average estimated spike-in signal versus concentrations on the nominal logarithmic scale (log₂) for Affymetrix data via (A) multi-RGX and (B) MAS5.0, MBEI, RMA, mgMOS, GC-RMA, multi-mgMOS, FGX and RGX methods.

Furthermore, the average estimated signals via multi-RGX versus nominal log-concentrations are plotted for each gene and the graph is reported in Figure 2.

Figure 2:

Signal versus concentrations on log-scale for Affymetrix data.

Then to evaluate the performance in accuracy for each level of concentrations, we regress the spike-in genes which have the same concentration according to the distinct log-concentrations and apply the simple linear regression method. The computed R² is called as the signal detect R² and its associated slope term is named as the signal detect slope. Then, we compute the R² values and slope terms for each gene separately and later calculate the mean of both values individually. We call this measure as simply R². Hereby as seen in Table 1, the signal detect R², signal detect slope and R² are found as 0.79, 1.04 and 0.92, respectively. From the results, it is seen that R² is comparable with respect to its alternatives. Yet, under all concentrations simultaneously, the array indicates relatively more linear relation between estimated signals and concentrations on the nominal log-scale on the average, similar to multi-mgMOS, which has one of the perfection results for this type of modelling. But different from multi-mgMOS, our novel index is fully frequestist, i.e. based on the maximum likelihood approach. On the other hand, it indicates a stronger linearity between the relation of signals and concentrations with respect to the FGX and RGX outputs. This can be seen another advantage of multi-RGX regarding other choices in the analyses.

Besides, in order to assess the computational demand for alternative approaches, we compare the time for calculation via multi-RGX with other records in Table 2. It is seen that multi-RGX is significantly faster than BGX and multi-mgMOS. But with respect to the time for FGX, it is relatively slower. The reason is that we need to compute the optimal concomitant of probe sets, which converge to the true order in two or three iterations and this calculation causes extra time that is observed in tabulated values. Moreover, the calculation of FGX is based on the mean probe level, on the contrary, the computation of multi-RGX depends on both probe and gene specific values, which complicate the calculation.

Application via simulated dataset

To assess the quality of our novel gene expression index, when the data are far from normal or have outliers, we simulate two datasets from 10,000 Monte Carlo runs. Each Monte Carlo run is conducted for 10 genes under 20 probe pairs, resulting in 200 observations for every generated array. In the first data, we generate the PM and MM values under the normal distribution and in the second data, we produce the values under the mean shifted normal density. Then, we mix these two sets in order to get two location-mixture datasets with different ratios of outliers. Hereby, the mathematical description of the first set is presented as below:

and the second location mixture can be described as the following structure:

where N(·,·) indicates the univariate normal with the given parameters.

Here, the first data possess a large number of outliers. On the other hand, the second dataset owns a relatively moderate number of extreme observations. Moreover, in all simulations, we assume that every gene is measured under specific concentrations. Thereby, the data are simulated under S_i setting to S_i =2, …, 11, for the gene i=1, …, 10, respectively, on the original scale for each presumed concentration level. Furthermore, μ_H , p, σ are equated as 1, 0.7 and 1, in order. Finally, the shift in the location, δ, is set to δ=10 for both datasets.

In the assessment of accuracies for the estimated results, we select three criteria: the average error (AE), mean absolute error (MAE) and the root mean square error (RMSE) [19].

In Tables 3 and 4, we list the estimated values, i.e. S_i , μ_H , p and σ, of FGX, RGX and multi-RGX for each gene and model parameters, in order, for the expression in Equations (15)–(16), respectively, by computing the associated absolute errors.

From Tables 3 and 4, we observe that there is a large amount of bias in the estimates whose cause can be explained as follows: The first dataset is generated for 10 genes on the original scale, representing very low concentrations. Then, we take the nominal log-scale of all these simulated measurements in order to calculate the estimated signals. By this way, indeed, we evaluate the performance of all three models under very low concentrations since the estimates of all indices are problematic under these levels. Furthermore, as observed in Figure 1B, FGX, RGX and multi-RGX are able to measure the amount of the noisy signals on the array when there is no effect on the gene. Accordingly, from the outcomes we can see that this amount is approximately 21 units as measured under μ^H, and if we subtract this shift from all estimated signals, we can find closer estimates of the signals with respect to their true values.

Additionally, in Table 5, we present the mean absolute error (MAE) and the root mean square error (RMSE) for the estimated signals via three selected models and in Table 6, we list the real and CPU time of each index in the calculation of the estimated model parameters for these two mixture data.

From all the tabulated values, it is seen that the estimates of multi-RGX are more accurate than the ones computed via FGX and RGX. But the computational demand of multi-RGX is more than other two alternatives in the real time. However, this difference is indeed can be tolerable if we compare the results of the CPU time. This output shows that the real time of multi-RGX can be even improved with the help of an effective programming since the calculation of the estimates is almost performed under the same computer time in the end.

Additionally, to evaluate the performance of each alternative in a large dataset, we extent the calculation via the simulated data by 10,000 and 20,000 genes. Finally, we compute AE, MAE and RMSE for estimated S_i in which the first data have 10,000 genes (i=1, …, 10,000) and the second data have 20,000 genes (i=1, …, 20,000). In contrast, for the remaining model parameters p, μ_H and σ, we compute directly AE, MAE and RMSE, which are unique per study.

Hereby, in the simulation we set the true value of p, μ_H and σ to 0.7, 1 and 1, in order. For the true signals S_i , we initially generate values according to the number of total genes ranging from 2 to 11 from uniform distributions, assuming that the genes are measured under 10 distinct concentrations similar to the previous dataset. Also, the intensities are measured on the original scale due to the fact that we aim to evaluate the performance of all three methods under the most problematic range of intensities, i.e. the worst scenario for the comparison. Then, alike other estimates, we calculate the average of all model selection criteria, i.e. AAE, MAE and RMSE, for all signals. Thereby, the average of these criteria for signals is found by using 10,000 and 20,000 values for the first and second dataset, respectively. Whereas, they are directly computed for estimates of p, μ_H and σ. The results are presented in Tables 7–9 for the first large dataset (with 10,000 genes) whose mixture proportions are arranged according to Equations (15)–(16), in order. Moreover, the outcomes of the second large dataset (with 20,000 genes) are listed in Tables 9 and 10 whose mixture model is generated with respect to Equations (15)–(16), respectively.

From all findings with large datasets, we observe that the performance of all three indices become very close to each other. In contrast, under small numbers of genes, indicating no unique best model, we find that FGX, RGX and multi-RGX estimates work well in the bench-mark dataset and for the simulated data with less numbers of genes, multi-RGX performs better than its strong alternatives. Hence, our suggested method can be seen as a promising approach to estimate the true signals for one-channel microarray datasets.