A longitudinal feature selection method identifies relevant genes to distinguish complicated injury and uncomplicated injury over time

Experimental data

The raw data for the injury microarray experiment by [14] were downloaded from the Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/). The accession number is GSE36809. This experiment was hybridized on Affymetrix HGU133 plus2 chips. Because we focus on identifying genes that present longitudinal expression change pattern between the trauma patients with complication and those without complication, only patients with uncomplicated recoveries and patients with complicated recoveries were considered here.

Based on the duration of recovery, uncomplicated recovery represents recovery within 5 days while complicated recovery includes recovery after 14 days, no recovery by 28 days, or death [14]. Thus, the potential time points for an uncomplicated recovery are days 1/2, 1, 4, 7 and 14, while those for a complicated recovery are days 1/2, 1, 4, 7,14, 21, and 28. Here, we collected 18 complicated patients with 7 full measures and 25 uncomplicated patients with 5 full measures and used the resulting dataset to train final models. Next, the data for patients with complications were truncated at 14 days in order to make comparisons with uncomplicated patients at each single time point. Lastly, we used the rest of complicated and uncomplicated patients as a test set to validate the results given by the proposed method. There were 50 uncomplicated patients and 23 complicated patients in the test set, and the time points considered in the test data were days 1/2, 1, 4, and 7 since early discharge of patients occurred before day 14. The characteristics of patients in the training set and the test set may be different since the patients in the test set were those had been discharged early from the hospitals. Of note, since different pre-processing procedures may impact the data analysis, the summary expression values of the experimental data provided by the GEO were downloaded and used directly. No pre-processing procedures were carried out.

Statistical methods

In this subsection, we first gave an introduction to the SAMGSR algorithm, and then provided a detailed description on the proposed procedure, which is referred to as the longitudinal SAMGSR method herein.

SAMGSR

The SAMGSR method extends the SAMGS method [15] by providing an additional reduction of significant gene sets into respective core subsets. This reduction step may approximately result in a 90% of reduction in the size of selected genes, in an effort to improve predictive performance and allow biological patterns to become more obvious. The SAMGSR method consists of two major steps [1]. The first step is the SAMGS process in which an SAMGS statistic is calculated. This statistic is the squared sum of SAM statistics over all genes within the specific gene set. Using a permutation test by perturbing phenotype labels to calculate a p-value for the SAMGS statistic, the statistical significance of a gene set is determined. In the second step, only those genes within significant gene sets identified by the first step are considered. The SAMGSR method reorders the genes within gene set j based on the magnitude of its SAM statistic and gradually partitions the entire gene set into two subsets: the reduced subset R_k which includes the first k genes with largest SAM statistic, and the residual subset \( {\overline{R}}_k \) being the complement of R_k for k = 1,…, |S-1|. Here S is the size of gene set j. Let c_k be the SAMGS p-value of the residual subset \( {\overline{R}}_k \). The optimal size of reduced set R_k is the smallest k such that c_k is larger than a pre-specified cut-off, e.g., 0.05. Conceptually, the significance level of a gene within a gene set is determined by the magnitude of its SAM statistic. It implies that if in a gene set |SAM_i| > |SAM_j| for genes i and j, gene j cannot enter the reduced subset unless gene i is inside the reduced subset.

Modification to SAMGSR for longitudinal data

here, d_t is the SAM statistic [16] of gene g (g = 1,…,G) at time point t (t = 1,…, T). In the SAM statistic, \( {\overline{x}}_d(t) \) and \( {\overline{x}}_c(t) \) are the sample averages of gene g at time point t for the diseased and control group, respectively. Parameter s(t) is a pooled standard deviation that is estimated by pooling the expression values of all samples at the specific time point t, while s_0t is a small positive constant used to offset the small variability in microarray expression measurements and thus to avoid the denominator of the SAM statistic being zero. Both s(t) and s_0t are specific for individual time points because the variability of expression measurements may differ over time.

From the above equation, it is observed that each gene’s expression profiles over time are combined into a gene set. Namely, a gene set represents one specific gene over different time points. Our rational is that a gene’s expression values for the same individual over time are correlated, mimicking a gene set/pathway. This method first calculates the SAMGS statistics for all genes to select the relevant genes and then determines exact time point(s) where the expression values of the specific gene differ between two phenotypes.

In the method, c_k is regarded as a tuning parameter. Using the sequence of 0.05, 0.1, …, 0.5, the optimal value of c_k corresponds to the one associated with the minimum 5-fold cross-validation (CV) error. Figure 1 elucidates graphically the longitudinal SAMGSR algorithm. Of note, the essential difference between the SAMGSR method and the longitudinal SAMGSR algorithm is what a gene set refers — while one corresponds to a set of genes and a single time point, the other contains only a single gene but many time points.

Performance statistics

In this study, we use four metrics - Belief Confusion Metric (BCM), Area Under the Precision-Recall Curve (AUPR), Generalized Brier Score (GBS) and the misclassified error – to evaluate the performance of a classifier. Our previous study [17] described those metrics in detail. Briefly, they all range from 0 to 1. For the first two metrics the closer to 1 the better a classifier is, whereas a value of 0 is optimal for the GBS and the misclassified error.

Since an evaluation on individual time points using the selected statistical metrics might be unfair for the SAMGSR extension in that its tendency to identify those genes that are insignificant at isolated time points but significant jointly over time, we used the following steps to calculate overall performance statistics. Specifically, for those methods incapable of providing the final model simultaneously with the selection of relevant genes, e.g., the longitudinal SAMGSR method, a linear support vector machine (SVM) model using the selected genes was fit to estimate the β coefficients at individual time points. Then, the posterior probabilities of the samples belonging to the diseased group were calculated at each time point on the basis of the SVM models, and the averages of those probabilities over time were taken and used to obtain these performance statistics. Figure 2 shows graphically how to calculate BCM, AUPR, GBS and the error rate. For those methods that are able to provide final models, e.g., LASSO, no extra SVM fitting was needed.

Statistical language and packages

Statistical analysis was conducted in the R language version 3.2.1 (http://www.r-project.org), and R codes for the longitudinal SAMGSR algorithm are available in an additional file (see Additional file 1).

Injury application

In order to evaluate the predictive performance of EDGE and make a more fair comparison between EDGE and our proposed method, we conducted the EDGE analysis on the training set and evaluated its performance on the test set. The results of EDGE analysis are presented in Table 1. Compared with the performance statistics from the longitudinal SAMGSR extension, EDGE does not show any superiority. Furthermore, there were 1083 genes identified by the EDGE method. In contrast, the overall unique number of selected genes is less than one hundred with our SAMGSR extension. The superiority of the SAMGSR extension over the EDGE method in terms of model parsimony is obvious. In terms of computing time, a single run of the simple SAMGSR algorithms takes 4.03 min on a Mac Pro equipped with a 2.2 GHZ dual-core processor and 16GB RAM. A single run of EDGE takes 4.2 min if the bootstrapping method is chosen to estimate the null distribution, which is more suitable than the asymptotic normality method for this specific application because of its moderate sample size.

The other feature selection algorithms evaluated were the original SAMGSR method [1], penalized SVM [18] and LASSO [19]. Of note, these three methods were applied separately on individual time points since they are incapable of carrying out longitudinal feature selection. The results are presented in Table 1 as well. Firstly, the comparison between the SAMGSR extension and the SAMGSR separately at each time point was made. While the longitudinal SAMGSR extension accounts for the correlations among the expression values of one specific gene over time, the application of SAMGSR at each individual time point considers the membership of genes (i.e., the specific gene belongs to which gene sets, thus might interplay with other genes inside those gene sets). The results of this comparison indicate that the SAMGSR extension is superior to the separate SAMGSR procedure. Incorporation of the pathway information inside gene sets, the clusters of genes that might be potentially co-expressed/co-regulated together, did not result in the separate SAMGSR method having substantially superior performance. One possible explanation relates to the quality of the pathway database itself. The canonically curated databases on pathways/gene sets are biased toward well-studied diseases such as cancers, and with few work has investigated on traumatic injury using gene expression profiles, the gene-to-gene interaction information contained inside those curated pathways may have no or limited meaning for injury. Therefore, we conclude that the consideration of coordinated effects existing among one gene’s expression profiles over time results in a better predictive performance.

Secondly, we compared the longitudinal SAMGSR algorithm with well-known conventional feature selection algorithms, namely, LASSO and penalized SVM separately at each time point. In these comparisons, we observed that the SAMGSR extension has at least comparable predictive performance but being superior in terms of overall model parsimony. For example, the longitudinal SAMGSR method identifies 97 unique genes while LASSO selects 147 genes. Moreover, we observed that in the 97- gene signature identified by longitudinal SAMGSR, there is a substantial proportion of overlaps for all 5 time points together (25.77%), while the number of genes being significant only at one specific time point is one half of this number (Fig. 3). Again, this highlights the ability of our SAMGSR extension to identify genes that present mild but concordant change over time points between two different phenotypes. In contrast, though at each individual time point LASSO selects the smallest number of genes, there is almost no overlap among those genes, resulting in a 147-gene list.

Regarding computing time, LASSO is the most efficient with one single run only taking less than 1 s since the implementation of LASSO in the R glmnet package [20] calls upon Fortune language and utilizes a fast-updating strategy. These two factors make substantial contributions to the time efficiency of LASSO. On the other hand, penalized SVM is the least time efficient one, which is unsurprising since the cross-validation process is automatically a part of model construction inside the penalized SVM modeling.

Simulations

Here, p_i = exp.(logit_i)/(1 + exp.(logit_i)). Then a random variable Y_i that follows a Bernoulli distribution with the parameter of p_i was simulated, it has two possible values with 1 indicating the sample belongs to the complicated injury group and 0 for the uncomplicated injury group. Notably, we considered one gene (i.e., F13A1) whose significance arises from its joint association accumulated from time point 1 to time point 4 and the other (i.e., GSTM1) whose association with the outcome is only at the third point time in this simulation.

The aim of this simulation was to illustrate the advantage possessed by the SAMGSR extension, namely, by incorporating the accumulated effect of genes over time, it can recognize genes with a coordinated change accumulating over time but only mild or moderate change at each time point.

Although in the second simulation the number of relevant time points was less than that in the first one, the number of selected genes by the longitudinal SAMGSR algorithm was dramatically larger in the second simulation. This might be because the relevant genes in the second simulation were highly correlated with other genes compared to the first simulation. The highly correlated structure between relevant features and irrelevant ones produced a large number of redundant features that the SAMGSR extension cannot exclude. To our best knowledge, however, many feature selection algorithms, especially the filter methods [11], suffer from this drawback. As illustrated in our previous work [21, 22], an additional reduction step using a statistical method such as bagging [23] may alleviate this problem.