A framework for significance analysis of gene expression data using dimension reduction methods

Smoker transcriptome

We first look at a data set that investigates the effects of smoking on human airway transcripts [15]. Data are available from three groups of subjects according to their smoking history. The original article reports significant differences both between a group of smokers compared to subjects that have never been smoking, and between former smokers and the group of never smokers. A PCA of the expression data X reveals that there is a large source of variation in the data that does not correspond to the available group information. In figure 1, a PCA score plot of the two first components is given, where the objects are coloured according to the classes "Never", "Former", and "Current". Unreported experimental factors, for instance a change in laboratory procedures, may be responsible for the separation of objects into two clusters along the first component. In the corresponding loading plot in figure 2, the 725 most significant genes from the PCA analysis are given in green. These genes span mainly the unreported X-variation and are not important for the smoking-induced differences of expression. The Venn diagram in figure 3 confirms that the genes detected by PCA are to a large degree irrelevant for classification of these data.

The PCA-analysis indicates that the group information should be allowed to guide the decomposition if any genes related to smoking are to be found. We know that a minimum of K - 1 independent variables are needed for a linear separation of K classes. For instance, one gene measured for several subjects may in theory be used to differentiate between two classes, while a minimum of two genes are needed for assignments into three classes. Due to co-regulation and inter-dependencies between genes, such single, descriptive features are usually both uninformative and hard to find. However, dimension reduction offers independent linear combinations of genes that may be used instead. If Y is a binary matrix holding the class information for each object, a good linear separation based on the expression data may be obtained with components that span all the Y-related variance in X. In Bridge-PLSR, this variance may be completely described in a minimum number of components. The score plot from a two-component Bridge-PLSR model is given in figure 4. Here, three groups corresponding to the predefined classes are revealed. This model explains 54% of Y, whereas an ordinary PLSR for comparison explains 45% in two components. Jack-knife is performed based on a full leave-one-out cross-validation, and the significant features are given in the loading plot in figure 5. The green spots correspond to features called significant by both SAM and Bridge-PLSR, whereas the red spots indicate significant genes not found by SAM. Neighbouring genes in the loading plot have similar profiles, whereas genes in opposite quadrants are negatively correlated. Also, the arrays in a specific area of the score plot have a high expression value for the genes in the corresponding area of the loading plot. As the significant genes span mainly the first component, they discriminate well between smokers and never-smokers, according to the score plot. Along the second component, which separates former smokers from the rest, fewer significant genes are found. This indicates that this correlation is weaker and more susceptible to random variations in the data. It is seen, however, that Bridge-PLSR finds more genes than SAM along this direction. Any significant genes along the second component describes smoking-induced transcriptomic changes among the former smokers.

A Venn diagram which compares the significant genes found with Bridge-PLSR, SAM and Limma is given in figure 6. All methods agree on the significance calls for the majority of the genes, however, Bridge-PLSR finds an additional 185 genes which are undetected by the other methods. To verify that these genes are not primarily false positives, the genes are annotated and their biological significance with regard to smoking is evaluated. A comprehensive list of our findings is given as supplementary information (see additional file 2). Of the genes identified as differentially expressed by Bridge-PLSR, only 66 have unknown or poorly known biological function. An additional 14 have well understood biological function often related to the natural system, but with no readily apparent link to smoking or lung damage. Of the remaining genes, 49 are known to be involved in regulation of cell growth, cell cycle, or apoptosis, processes known to be affected by smoking [16, 17]. Seven of these genes as well as another 14 genes on the list not related to cell growth, have been shown to be involved in various forms of lung cancer. As the majority of lung cancer sufferers are smokers or previous smokers, these genes can be linked to smoking directly. Another interpretation is that these genes change early in the cancer development process. Seven genes are identified that have been reported to take part in lung development, this is not surprising as smoking can cause dramatic changes in the airways. In addition to cancer, we identify 17 genes that are related to other diseases in the lungs and bronchia, such as asthma, inflammation, fibrosis, or response to metallic toxins. Another large group of genes belongs to the protein life cycle, particularly the ribosome, or the ubiquitin cycle. This could point to an overall change in the rate of protein turnover caused by smoking. This is consistent with the results of Tomfohr et al. [18], who analysed the same data set with gene set enrichment analysis. Five genes are also found that are known to be expressed in lungs, but for which we find no link to cancer development or lung damage in the literature. Overall it seems that of the genes for which biological background information is available, the majority relate to biological processes that are influenced by smoking. The additional genes found by Bridge-PLSR are therefore unlikely to be false discoveries.

For prediction purposes, ordinary PLSR seems to be superior to Bridge-PLSR. The best predictive PLSR model is obtained for two components, with an explained Y-variance of 32%. The best predictive Bridge-PLSR model is a one-component solution with a cross-validated explained variance in Y of 20%. The residual mean squared error of prediction even increases slightly when two components are used. Caution must therefore be used when interpreting the second component. They are both included here because the main goal of the analysis is significance testing. Hotelling's T²-test using only one component corresponds to an ordinary one-sample t-test, and is, in this case, only able to separate between the groups "Never" and "Current". In order to find genes that might be relevant for former smokers, the second component must be investigated. The jack-knife is believed to reject most of the genes that cause the residual error to increase, as those genes have a high variance across cross-validation segments. For prediction purposes, the model with the lowest residual error might be used instead. It should be noted, however, that even if the Bridge-PLSR or PLSR-model with one or two components, respectively, is optimal for prediction, neither will be able to predict former smokers.

A ridge parameter γ controls the influence of Y in the decomposition. A small value of this parameter ensures that the first Bridge-PLSR components span the directions of maximum covariance between X and Y. This is ideal from the point of view of significance analysis, as the major components will span all the Y-relevant variance. However, if γ = 0, the maximum number of components is limited by the rank of Y, which for these particular data is two. For explorative reasons, γ is set to a default value of 0.01 throughout this work in order to enable interpretation also of later components. For the first two components, this has no visible effect compared to a γ-value of zero. Components three and four, however, are almost identical to the first two PCA components in figure 1 (not shown). This is a further verification that the Bridge-PLSR is able to separate the Y-relevant structures in X from technical and other irrelevant variation. It is possible to optimise both the ridge parameter and the number of components to minimise the residual error of the bilinear model. This would require an extra layer of validation, for instance by cross model validation [19], to obtain unbiased error estimates. As the primary goal in this work is significance analysis, such an optimisation is out of scope. However, a supplementary response surface is provided for reference, along with the ordinary PLSR results (see additional file 3). It can be seen that the cross-validated error variance is robust with regard to γ for small values of the parameter.

Analysis by ordinary PLSR is performed for the sake of comparison. This model uses two components to separate never and current smokers, and additional two or three components are needed to describe the former smokers. The optimal number of components according to the cross-validation is two, but the majority of the significant genes found by this model span only the irrelevant X-information (not shown). Because of the irrelevant, dominating structure in X, Bridge-PLSR is a much better choice for significance analysis in this case.

Spiked-in data

A further analysis is performed with Affymetrix spike-in data in a constant human background [20]. The spiked-in transcripts hybridise to 42 known probe sets, however, several authors report that also other probe sets show differential expression due to sequence or function similarity. Sheffler et al. [21] remove 56 such genes as outliers prior to analysis. Two independent groups [22, 23] suggest that additional 22 or 23 probe sets, respectively, should be added to the original list of known spike-ins. These lists differ in a single probe set, which is known to have identical sequence to one of the spiked-in genes for 5 of 11 probe pairs. We include the 22 extra genes which are suggested by both groups in our analysis, yielding a total of 64 spike-in probe sets. Bridge-PLSR with jack-knife finds all these genes among the 66 most significant ones. An estimated FDR of 0.01 results in 69 significant genes, corresponding to a true FDR of 7%. The five false positives are the probe sets "203173_s_at", "204890_s_at", "204891_s_at", "217019_at", and "200822_x_at". The maximum squared correlation coefficients between the expression values of these probe sets and the concentration profiles are, respectively, r² ∈ {0.92, 0.88, 0.86, 0.33, 0.29}. For comparison, r² for the list of true positives range from 0.39 to 0.74. In passing, it is noted that the three false positives with high correlation are included also in the outlier-list of Sheffer et al. [21].

Rat liver study

A data set investigating liver damage in rats after exposure to toxic doses of paracetamol is also analysed [24]. In addition to design information of dose and time of exposure, several pathological variables based on blood and tissue chemistry are available for the rats examined. There is an abundance of possible comparisons in these data, and analysis by Limma or SAM is therefore difficult. We perform an analysis by SAM in which two groups of "high" and "low" dosage are defined for each of the four time-points of the design. The drastic treatment some of these rats receive induces profound genetic responses, leading to a long list of differentially expressed features. Testing for all 28 contrasts defined by the 8 groups, SAM assigns significance to 50% of the genes at a significance level of 0.001. This result verifies something that might be expected – that poisoning of an organism can affect several biological processes that in turn involve many genes. Investigation of more complex relationships based on these methods is not straightforward.

PLSR is well suited for handling both multiple comparisons and continuous responses. The score plot from a Bridge-PLSR model with all available response information in Y is given in figure 7. The spots are coloured according to the administered drug-doses. The low doses are clearly separated from the toxic doses along the first component, but the classes within those two groups are mixed together. The dashed lines envelope three groups that are related to the time from drug exposure to necropsy. The lower group consists of the 6-hour design points, the upper group corresponds to the 48-hour points, while in the middle group the 18- and 24-hour points are mingled together. The second component is thus related to time since exposure. A leave-one-out cross-validation is performed also for these data. The optimal model for prediction is obtained after four components, with a residual variance of 47%. However, we are interested mainly in the genes which are differentially expressed across dose and time. We therefore use a two-component model in the significance analysis, with a prediction error of 58%. For these data with much structure in Y, the differences between the results from Bridge-PLSR and ordinary PLSR are small. All predictive components describe the same phenomena, and the residual errors are similar, both in a calibrated and validated sense. For reference, the full response surface including several values γ of is given as supplementary information (see additional file 3).

A helpful plot is the Y-loading plot in figure 8. Here, each spot corresponds to a response-variable, and the predictabilities are given by the associated colour-bar. The colour of a spot reflects the cross-validated, explained variance for that variable, given in percent of the true response's variation. The variables with red spots are very well predicted while the blue spots represent variables that cannot reliably be modelled based on the expression data. Two variables which have high predictive abilities and are correlated with the dose are ALT (alanine aminotransferase) and AST (aspartate aminotransferase). Increased levels of these enzymes in the blood stream are diagnostic markers for liver damage. The large degree of explained variance for these variables indicates that the first component actually describes liver injury directly, rather than just the concentration of paracetamol in the rats. Because the individual responses to the drug vary, these markers are better suited to model the genetic response to overdose than the design is. This information should be included also in the interpretation of the scores in figure 7. It is seen that the response to overdose peaks between 18 and 24 hours after exposure, as these groups span the first component better than the 6- and 48-hour groups. After 48 hours, metabolism may have disposed of some of the paracetamol, the liver-cells may have started to repair, and the levels of ALT and AST in the blood have started to decrease. Other correlations seen in figure 8 may also be of interest. Liver damage seems to be associated with increased concentrations of total bile acids and blood urea nitrogen, and there is also an indication that the level of cholesterol decreases. The reduced glutathione level in the liver tissue increases with time after poisoning, and the level of blood albumin seems to decrease accordingly.

An estimated FDR significance cutoff of 0.1% yields 25% significant features for the Bridge-PLSR model. The X-loading plot is given in figure 9, with the significant genes marked by green spots. The majority of these features span the component describing liver injury, but some are also found that vary with time. The significant genes found by SAM are scattered more uniformly over the same loadings, making any systematic trends hard to detect (not shown).

Data sets

The smoker data [15] contain expression levels for M = 7939 probe sets measured for N = 75 subjects on Affymetrix HG-U133A GeneChips (Affymetrix, Inc., USA). The subjects are divided into 34 current smokers, 18 former smokers and 23 individuals that have never been smoking. The expression data are arranged in a matrix X of size N rows by M columns, and the class information is contained in a binary array Y. Each class is represented by a column in Y, where a value of one assigns an object to that class. X is logarithmically transformed prior to the analysis. No other pre-processing or scaling is performed, that is, we use the data directly as provided by Spira et al. [15].

The spike-in data [20] are obtained from a designed experiment where 14 groups of three transcripts each were spiked in at 14 different concentration levels in three replicates. There is thus a total of N = 42 objects in the experiment, corresponding to specific concentration mixtures. Expression levels were measured on a modified Affymetrix HG-U133A array, and concentrations ranging from zero to 512 pM were arranged in a Latin square setup. The 42 spike-in genes were measured in a constant human background, 22300 probe sets in total. The data are pre-processed using the gcrma-package [31], which is part of the Bioconductor project [32]. As groups of three transcripts are spiked in at the same amounts, the response matrix Y contains 14 concentration profiles, one for each group. Several authors have remarked that some of the background genes show consistent differential expression across the experiments [21–23]. Some of those genes contain sequences which match the spike-ins exactly. We therefore adopt the extended spike-in gene-list given by McGee and Chen [22], yielding a total of 64 entries. The list is identical to the one suggested by Lo and Gottardo [23], except the probe set "205397_x_at" is removed. This probe set has identical sequence to the originally included "205398_at" for 5 of 11 probes. Due to memory limitations, there is a need to filter away some genes before the analysis. We therefore include only those normalised probe sets with a standard deviation above 0.4, which results in a total of M = 14.490 variables in X. Potential alternatives to filtering would have been to apply a leave-two-out cross-validation, or to use fewer Y-variables or components in the analysis.

In the liver toxicity study [24], male rats of the inbred strain Fisher 344 were administered non-toxic (50 or 150 mg/kg), moderately toxic (1500 mg/kg) or severely toxic (2000 mg/kg) doses of paracetamol in a controlled experiment. Necropsies were performed at 6, 18, 24 and 48 hours after exposure. RNA analysis was performed on an Agilent platform (Agilent Technologies, USA), and four biological replicates were obtained for each of the 16 design-points. All spots flagged as bad or not found by the image analysis program are removed. In addition, spots are removed if they have total intensity less than 200, foreground to background intensity less than 1.5, spot diameter less then 25 pixels, or more than 50% of pixels saturated. As the background estimates available in the Agilent image analysis report are local backgrounds measured in pixels near the spots, it is doubtful if this is a useful estimate of the background in the spot [33]. Therefore, no background correction is done. The filtered data are transformed to log₂-ratios and normalised using global loess to improve colour balance and remove global intensity effects [34]. Ratios for spots that have been removed in the filtering are imputed using k-nearest neighbours with k = 10 [35]. Two histological and ten clinical variables containing markers for liver injury are available for each object. The expression data are arranged in a matrix X of N = 64 objects and M = 7445 expression levels. For analysis by SAM, the two non-toxic levels are pooled in a group of low dose, and the two toxic levels are pooled in a group of high dose. The data are further divided into groups for each of the four time-points, giving a total of 8 groups to be tested for differential expression. In the bilinear analysis, a matrix Y of 14 columns is constructed instead that contains the design information of dose and time as well as the phenotypic responses available. A description of all the variables in Y are given in figure 8. The Y-variables SDH, BILE, AST and ALT are log-transformed in order to get responses which are more closely bell-shaped.

Bridge-PLS regression

PLSR finds linear combinations (scores) of the original features spanning directions of high covariance between X and Y. In the most common PLSR implementation, the data are deflated each time a component is found, by subtraction of that component from the remaining X-data. This ensures that each score vector is orthogonal to all previous components. It also enables better model fit, as several components are allowed in the modelling of each individual response. Cross-validation is usually applied to find the number of components which yields the best predictions [3]. However, the best prediction model does not necessarily contain the most informative components. If a single response is modelled (i.e. if Y has only one column), qualitatively useful information is expected from the vector of covariances between the variables in X and the response [36, 37]. The covariances are given by the so-called loading weights corresponding to the first component, and this vector is therefore well suited for significance analysis. Subsequent loading weights are less interpretable, as they, due to the deflation step, may be influenced by information in X unrelated to Y. In order to obtain informative weights for multiple responses, it is possible to make one PLSR-model for each individual response. Such an approach was taken e.g. in [9], where cell cycle regulated genes were found using the first loading weight from two single-response PLSR-models.

Bridge-PLSR has no deflation step, and the loading weights are used directly in the model rather than calculating new loadings [14]. If the ridge-parameter γ is set to zero, Bridge-PLSR is simply a singular value decomposition of the covariance matrix X ' Y, with a subsequent prediction relation between X and Y. At the other extreme, when γ = 1, Bridge-PLSR becomes an ordinary principal component regression [3]. Any intermediate value of γ corresponds to some weighted average between those two methods. For significance analysis, and for interpretation of components in relation to Y, the optimal value of γ is zero. The loadings will then collectively explain the same Y-information as the set of loading weights from all one-component, individual response PLSR-models. For multiple correlated responses, Kvalheim [38] recommends orthogonalisation of Y before calculating such single-response models. This ensures that the models will be close to orthogonal to each other as long as the predictive relation between X and Y is good. These steps are all circumvented with Bridge-PLSR, as the Bridge-PLSR loadings are orthogonal by default. In addition, a single, predictive Bridge-PLSR model is more useful and appealing to work with than multiple single-response PLSR-models.

When γ = 0, the maximum number of components equals the sum of independent responses in Y. If all these components are used, the residuals will be orthogonal to Y, and the loadings (= loading weights) will together account for all the Y-related information in X. Because only one component is available for each orthogonal response when γ = 0, the predictive ability of Bridge-PLSR is often less than that of an ordinary PLSR. However, a tuning of the γ-parameter may increase the predictive ability of the model. This is expected for instance in cases where Y is large or noisy compared to X. It is also seen from the response surfaces provided in supplementary information that a larger γ-value may stabilise the decomposition such that early local minima are avoided. However, the response surfaces provided are not validated by an external validation and should not be used for optimisation. Any such tuning of parameters should be accompanied for instance by cross model validation to avoid optimistically biased results [19, 39].

For the analysis in this work, γ is set to the small default value of 0.01. It is seen from the response surfaces in supplementary information that this yields the same predictive ability (down to two decimals) as if γ were zero. Also, the components obtained using the two values of γ are qualitatively equal, such that the influence on the significance analysis is negligible. The γ-parameter is given this value only to enable interpretation of more than two components for the smoker data, although these components are not used in the jack-knife. Bridge-PLSR is algorithmically simple, and the models are validated by cross-validation or test sets in a straight forward manner.

Significance analysis

A cross-validation using an appropriate dimension reduction method is performed, and loadings from each sub-model with left-out objects is obtained. The perturbed loadings are used to estimate the stability of each feature. The optimal number of components A for prediction, and the overall predictability of the model may be estimated as described in [3]. Only components that lead to a significant decrease in the cross-validated prediction error should then be included. In this work, we concentrate more on significance analysis than on optimisation of prediction models, and interpretation of score plots is therefore at least as important when A is found.

Many of the dimension reduction methods yield orthonormal loadings, in which cases the relative importance of the components are reflected by the length of the scores. For the significance analysis, it is important that the loadings themselves are weighted by the variance of their components. All perturbed loadings from the cross-validation also have to be aligned properly. A bilinear model can be written on the form $\widehat{X}=T_A{P}_A^T=T_{A}C{C}^{-1}{P}_A^T$, where C is any invertible matrix of rank A. This means that the components have full rotational freedom, and it is often observed that a component is flipped 180° in one sub-model relative to another. Components of similar eigenvalues may also be sorted differently in different sub-models. A Procrustes rotation of the perturbed models toward the calibration model will correct for these inconsistencies [3]. We perform un-centred and un-scaled rotations of all loadings prior to significance testing.

Each variable (gene) is tested under the null hypothesis that its corresponding loading-parameter is zero in the A-dimensional space spanned by the components. It follows that the choice of dimension reduction is very important for the outcome of the test. Figures 1 and 2 show that components spanning only irrelevant X-variation will yield irrelevant genes. On the other hand, when the components describe some Y-information of interest, we directly find the features important for the different responses (e.g. figures 4 and 5). Hotelling's T²-test is able to assess the significance of parameters in several dimensions simultaneously. We have modified the test to make it more suited for jack-knife of bilinear model parameters. Some of these modifications are described in [3, 12] for the regular one-sample t-test. The variation of the perturbed sub-models around an estimated true model describes the stability for each gene. Traditionally, the mean of the sub-models represents the centre points, but the median or even the full calibration model may also be used. Because the perturbed models are not independent of the full model, the deviations are summed instead of averaged. The latter point makes a difference if the actual value of the T²-score is interpreted, not if the T²-score is used simply as a ranking of the genes. Finally, shrinkage of the gene-specific variances toward a total variance across genes is included [25]. This is to avoid poorly estimated variances due to the small number of objects in microarray experiments. A shrinkage factor of zero means that every gene is evaluated individually, while a factor of one treats all genes as having equal variances. The actual value of the parameter is a choice that must be made on the basis of the number of replicates and the overall data quality. Poor quality data and few replicates may justify increased shrinkage.

The vector p _i is of length A and represents the true loadings, β ∈ [0, 1] is the shrinkage factor, $P_i^{CV}$ is a matrix of size K-by-A holding the loadings from the K sub-models and 1 is a vector of ones. S ^tot is the mean or median of the scaled, gene-specific variance-covariance matrices $S_i^{gene}$ across all genes. The scaling $g=\frac{K-1}{K}$ is taken from [3] and is close to one in most cases. In this work, p _i is the mean value over sub-models, the coefficient β = 0.1 and the median of covariances is used to find S ^tot .

FDR estimation

Every feature i is associated with a score $T_i^2$. The denominator holds the total number of genes at least as significant as this based on the calibration model. The number of false positives at the same significance level is not known, and must be estimated. The estimated count for gene i is found by reshuffling and denoted $\widehat{FP}$ _i . The rows of Y are shuffled R times to give several matrices $Y_r^{\ast }$, r ∈ {1,..., R}, with random object-information only. A jack-knife based on one of the randomised matrices is expected to give no true positives, especially if the number of significant features is small compared to the total number of features tested. Here, $\widehat{FP}$ is found as the median number of significant outcomes over R = 1000 shuffled responses, and all features with ${\widehat{q}}_i$ smaller than a pre-set significance level α are called significant. Estimation of $\widehat{FP}$ when no Y is available is difficult, because it is not possible to shuffle X for each feature separately without simultaneously disrupting all other systematic structure in the data. For the smoker-data, we have therefore used the same number of significant features for PCA as we found for Bridge-PLSR. The FDR significance level is set to α = 0.05 for the smoker-data and to α = 0.01 for the spike-in data. Because the relevant signal in the paracetamol data set is very strong, a significance level of α = 0.001 is enough to get a large number of significant features for these data.

SAM and Limma analysis

SAM and Limma are used to identify differentially expressed genes in the smoker data set. Both methods are part of the Bioconductor software (Version 1.9) [32], and we use the packages limma and siggenes for Limma and SAM, respectively. Limma is performed with the functions lmFit and eBayes with default settings. We control Benjamini and Hochberg's FDR at a level of α using the function decideTests. SAM is run with the function sam.dstat with 1000 permutations, the parameter delta = 0.8 and the remaining parameters at default setting.