Asthma-predictive genetic markers in gene expression profiling of peripheral blood mononuclear cells.

PURPOSE: We sought to identify asthma-related genes and to examine the potential of these genes to predict asthma, based on expression levels.
METHODS: The subjects were 42 asthmatics and 10 normal healthy controls. PBMC RNA was subjected to microarray analysis using a 35K array; t-tests were used to identify genes that were expressed differentially between the two groups. A multiple logistic regression analysis was applied to the differentially expressed genes, and area under the curve (AUC) values from receiver operating characteristic (ROC) curves were obtained.
RESULTS: In total, 170 genes were selected using the following criteria: P≤0.001 and ≥2-fold change. Among these genes, 57 were up-regulated and 113 were down-regulated in asthmatics versus normal controls. A multiple logistic regression analysis was done using more stringent criteria (P≤0.001 and ≥5-fold change), and eight genes were selected as candidate asthma biomarkers. Using these genes, 255 models (2(8)-1) were generated. Among them, only 85 showed P≤0.05 by multiple logistic regression analysis. Based on the AUCs from ROC curves for the 85 models, we found that the best model consisted of the genes MEPE, MLSTD1, and TRIM37. The model showed 0.9928 of the AUC with 98% sensitivity and 80% specificity.
CONCLUSIONS: MEPE, MLSTD1, and TRIM37 may be useful biomarkers for asthma.

INTRODUCTION

Asthma is a common and heterogeneous respiratory disease characterized by intermittent airway obstruction and respiratory symptoms that are related to chronic airway inflammation and remodeling.1 Pathological features of airway remodeling include goblet cell hyperplasia, subepithelial fibrosis, collagen deposition, mucosal gland hyperplasia, smooth muscle hypertrophy, and changes in the extracellular matrix.1

Inflammation and remodeling are the main causes of airway hyperresponsiveness and chronic airway obstruction, characteristic features of asthma. These changes are attributed to the altered expression of genes associated with transcriptional pathways, inflammatory processes, apoptosis, and cell proliferation.2

However, the complex nature of the asthma phenotype, together with genetic heterogeneity and environmental influences, has made it difficult to uncover genetic aspects of this common condition.3 Thus, the search for asthma-specific genes is an important issue because such genes may promote the early detection, positive prognosis, and treatment of the disease. To this end, large-scale, high-throughput, whole-genome studies are needed to understand the genomic contribution to asthma.

Although many approaches have been used to investigate relationships between diseases and genes, the high-throughput microarray is one of the most important. Microarray technology allows the monitoring of gene expression on a genomic scale. It also enables the classification of genes causing various diseases and for a diagnosis to be made at the genetic level.4,5

Several types of tissues have been used in human microarray studies of asthma, including airway epithelial cells and airway inflammatory cells from patients and healthy subjects.6-10 Nasal mucosal cells are a good choice for genomic studies.11

Although an analysis of gene expression in the airways can provide important information, obtaining airway tissues, including airway epithelial cells, using bronchoalveolar lavage (BAL) or brushing is difficult, especially in severe asthma12; thus, peripheral blood mononuclear cells (PBMCs) have been used as an alternative.13-16 PBMCs contain lymphocytes, monocytes, and dendritic cells, which are recognized sources of allergic response mediators.

Over the past decade, using gene expression profile analysis, the diagnosis of multifactorial complex diseases and their subphenotypes have been attempted. Although the diagnostic value of the expression of a single gene is regarded as minimal, multi-gene analyses can increase the predictability, and possibly lead to usable biomarkers.

In the present study, we identified candidate genes affecting asthma development using an mRNA expression chip. Additionally, we evaluated the diagnostic value of differential gene expression levels in the discrimination of asthmatics from normal controls through an additional analysis of the area under the curve (AUC) values from receiver operating characteristic (ROC) curves, sensitivity, and specificity.17 Finally, the selected genes were analyzed according to asthma severity.

MATERIALS AND METHODS

Study subjects

The subjects were recruited from the Genome Research Center for Allergy and Respiratory Diseases at Soonchunhyang University Hospital. All subjects were Korean. The protocol was approved by the local ethics committee of Soonchunhyang University Hospital. Written informed consent was obtained from all subjects.

All patients met the definition of asthma given in the Global Initiative for Asthma (GINA) guidelines.1 Each patient showed airway reversibility, as documented by a positive bronchodilator response of a greater than 15% increase in FEV1 and/or airway hyperreactivity to less than 10 mg/mL methacholine. Normal controls were recruited from hospital personnel; they gave negative answers to a screening questionnaire for respiratory symptoms and had an FEV1 ≥75% of the predicted value, a PC20 methacholine ≥10 mg/mL, and a normal chest X-ray.

Subjects were assessed using a skin-prick test for 24 common inhaled allergens, including house dust mites, Alternaria, Aspergillus, pollen, dogs, cats, and cockroaches. Atopy was defined as having a wheal reaction with the allergen extract equal to or greater than that with histamine (1 mg/mL) or 3 mm in diameter. Total IgE was measured using the UniCAP system (Pharmacia Diagnostics, Uppsala, Sweden). Peripheral venous blood was obtained when the subjects were in a stable state. No patient had an exacerbation or systemic steroid treatment within 6 weeks of the study. Clinical profiles of the asthma patients and control subjects are summarized in Table 1.

Table 1

Clinical profiles for asthma patients and normal controls, mild asthma and moderate to severe asthma subjects

Mild asthma and moderate to severe asthma are separated by FEV1 value of 80% predicted and P values were acquired by t-test or Fisher's exact test.

NS, non-smoker; SM, smoker; ES, ex-smoker.

cDNA synthesis and microarray hybridization

A human oligonucleotide microarrays (Genomictree Inc., Daejeon, Korea) containing 34,560 oligonucleotide probes, representing human genes and gene transcripts, were used for gene expression analysis. The 34,560 oligoncleotide probes were purchased from (Qiagen, Germantown, MD, USA). The synthesis of target cDNA probes and hybridization were performed according to previously described.8 Briefly, each 50 µg of total RNA was mixed with 2 µg of oligo-24N (dT) (GenoTech, Daejeon, Korea) in 15.4 µL of RNase free water and incubated at 65℃ for 10 minutes. After incubation, the single stranded cDNA was synthesized in the presence of Cy3-dUTP or Cy5-dUTP (1 mM each, NEN Life Science Products, Boston, MA, USA) at 42℃ for 2 hours. Total RNA from test samples and common reference RNA pooled from 10 normal samples were labeled with Cy5 and Cy3, respectively. Both of Cy3 and Cy5-labeled cDNA were purified using PCR purification kit (Qiagen) as recommended by manufacturer. The purified cDNA was resuspended in 100 µL of hybridization solution containing 5X SSC, 0.1% SDS, 30% formamide, 20 µg of Human Cot-1 DNA, 20 µg of poly A RNA and 20 µg of Yeast tRNA (Invitrogen, Grand Island, NY, USA). The hybridization mixtures were heated at 100℃ for 2-3 minutes and directly pipetted onto microarrays. The arrays hybridized at 42℃ for 12-16 hours in the humidified hybridization chamber (Genomictree Inc.). The hybridized miroarrays were washed with 2×SSC/0.1% SDS for 5 minutes, 0.1×SSC/0.1% SDS for 10 minutes, and 0.1×SSC for 2 minutes two times. The washed microarrays were immediately dried using the microarray centrifuge (Genomictree Inc.).

Statistical methods

For microarray analysis, background-corrected values for each probe on the BeadChip array were extracted using Bead-Studio (Illumina, San Diego, CA, USA); P values were computed from the background model characterized by the chance that the target sequence signal was distinguishable from the negative controls on the same chip. Normalization was performed for the extracted values. The statistical significance of the microarray data was calculated by t-tests using 'R' (ver. 2.8.1; http://www.r-project.org/). If the distribution of the mRNA expression levels showed normality and equality in the variances in each group (i.e., asthma and normal controls), a classic two-sample t-test was used. If the distribution of mRNA expression levels showed normality and inequality in the variances for each group, the Welch two-sample t-test was used. If the distribution of mRNA expression levels did not show normality in each group, then the Mann-Whitney U-test was applied. For the diagnostic gene expression values used to discriminate the asthmatics from the normal controls, we selected candidate genes with P<0.001 and a fold-change ≥5 between the two groups. Fold changes in gene expression were obtained by dividing the mean value for the asthmatics by that for the normal controls. After filtering by P value and fold-change, a power set was made using the candidate genes, and a multiple logistic regression analysis was performed with elements in the power set. ROC curves for all models were obtained and AUCs were calculated19,20 to select disease marker genes.

RESULTS

Comparison of gene expression between the asthmatics and normal controls

To identify genes that may be related to asthma, we applied a high-throughput gene expression microarray consisting of 15,054 featured genes on RNA samples, which were obtained from normal controls (n = 10) and asthma patients (n = 42). All values in the microarray cell were normalized and used for feature selection. The overall approach is depicted in Fig. 1.

Fig. 1

Gene expression profiling strategy and general workflow.

To evaluate overall differences in gene expression levels in PBMCs between asthmatics and normal controls, we calculated gene expression as shown on a volcano plot. To identify differentially expressed genes between the asthmatics and normal controls, we applied two types of t-tests and the Mann-Whitney U-test to the difference in mean expression level between the two groups. We selected genes having a P value<0.001 and a change of 2-fold or greater. Volcano plots of significance against the fold-change values for each gene in the PBMCs revealed that the expression levels were quite different between the asthmatics and normal controls (Fig. 2A). Using the criteria of P≤0.001 and a fold change ≥2, we identified 57 genes showing a significant increase in expression and 113 genes showing a significant decrease in expression in the asthmatics compared with the normal controls (see Table S1 in additional file 1). A heat map of these differentially expressed genes (n=170) revealed distinctive RNA expression profiles for the PBMCs (Fig. 2B). The up- and down-regulated genes were well-classified by this hierarchical clustering method.

Fig. 2

(A) Volcano graph. The x-axis represents the logarithm of the fold-change value, while the y-axis represents the negative logarithm of the P value. A and B denote areas satisfying the following criteria: ≥2-fold change and P≤0.001. (B) Hierarchical clustering graph and heat map.

Because 170 genes would be too many to handle as biomarkers, we used the more stringent criteria of P≤0.001 and a fold change ≥5 to the differences in gene expression between the asthmatics and normal controls. Consequently, we identified eight genes (Fig. 3).

Fig. 3

Volcano graph. The x-axis represents the logarithm of the fold-change value, while the y-axis represents the negative logarithm of the P value. A and B denote areas satisfying the following criteria: ≥5-fold change and P≤0.001.

Development of genetic biomarkers for the diagnosis of bronchial asthma

The genes that met these criteria are listed in Table 2. The expression of ZFP161 and NOX5 was increased, while that of LMAN1, MEPE, MLSTD1, TRIM37, KNS2, and CCT5 was decreased in the PBMCs of the asthmatics versus the normal controls.

Table 2

List of genes meeting the criteria of P≤0.001 and a fold change of ≥5-fold

*Fold change (asthma/normal control).

Selection of the best discriminative model using multiple logistic regression analysis and ROC curve analysis

Using the eight genes as candidate biomarkers to discriminate between asthmatics and normal controls, we performed a multiple logistic regression analysis for all elements of the power set and measured the AUCs from ROC curves. First, we made 255 (28-1) models for all elements of the power set and measured the P values of the variable for each model. We separated the 255 models into eight groups (Groups 1-8). Group n indicates models made of n genes for multiple logistic regression analysis. Among the 255 models, only 85 showed P≤0.05. All data for the AUC and P values are presented in Supplementary Table S2 (see additional file 2). As the number of genes increased, the P value decreased (Fig. 4). Only models comprising fewer than three genes showed significant P values (i.e., <0.05).

Fig. 4

(A) Distribution of the average AUCs for the top five in group n. (B) Distribution of the average log of the P values for the top five in group n. The dashed line indicates the cut-off value (P<0.05).

The AUC and ROC were calculated for Groups 1, 2, and 3. Among them, single genes or each combination of two or three genes having the five highest AUC values are presented in Fig. 5A. The best AUC value (0.9928) was observed for a combination of three genes, MEPE, MLSTD1, and TRIM37 showing asymptotic P values: 0.000001, asymptotic 95% confidence interval (lower bound: 0.977, upper bound: 1). The ROC curves and AUC values for the other genes and combinations are presented in Supplementary Figure S1 (see additional file 4).

Fig. 5

(A) Values of the AUCs for each of the eight genes, two-gene combinations, and three-gene combinations. (B) ROC curve of the best model, consisting of MEPE, MLSTD1, and TRIM37 (P value: 0.000001; asymptotic 95% confidence interval lower bound: 0.977, upper bound: 1; AUC: 0.9928).

Discriminating power of the combination of MEPE, MLSTD1, and TRIM37 between asthmatics and normal controls

To evaluate the discriminating power of the combination of MEPE, MLSTD1, and TRIM37 between asthmatics and normal controls, we calculated the sensitivity and specificity using a contingency table of 42 asthmatics and 10 normal controls. As shown in Table 3, the sensitivity and selectivity were 0.98 and 0.80, respectively, while the accuracy was 0.942. To evaluate the diagnostic accuracy of the three-gene combination in an independent data set, we applied three-fold cross-validation (CV). The average results of three-fold CV for sensitivity, specificity, and accuracy were 1, 1, and 1, respectively.

Table 3

Contingency table for the best model

*Fold change (asthma/normal control).

Analysis of the selected genes according to asthma severity

We divided the asthma group into mild asthma and moderate-to-severe asthma according to their FEV1% (threshold: 80), and we analyzed the genes for each paired group as follows: normal controls versus mild asthma, normal controls versus moderate-to-severe asthma, and mild asthma versus moderate-to-severe asthma. By this analysis, we could predict genes involved in asthma development. Our complete results are presented in Supplementary Table S3 (see additional file 3). Table 4 shows the results for the eight genes.

Table 4

The selected 8 genes analysis according to asthma severity

Based on our results, TRIM37 is predictive for the occurrence of asthma, while NOX5, LMAN1, KNS2, and MLSTD1 are predictive for mild asthma. CCT5 is predictive for moderate-to-severe asthma.

DISCUSSION

We identified genes related to asthma using a microarray analysis of PBMCs. Using these data, we found that a combination of three genes, MEPE, MLSTD1, and TRIM37, showed the best discriminating power between asthmatics and normal controls. To our knowledge, this is the first reported study to select disease markers using ROC analyses of PBMCs. In our microarray analysis, when a P value threshold of <0.001 and a fold change of ≥2 were applied, a total 170 genes were selected. Of these genes, 57 were up-regulated and 113 were down-regulated in PBMCs. The genes included inflammatory and immune response genes, such as NOX5, MALT1, TNFRSF10C, GRK5, CXCL3, RELA, CD40, ABR, RELB, REBB2, PGLYRP1, CD82, RPE, CFTR, KITLG, and IKBKG, which would be expected to be expressed differentially in asthmatics versus normal controls. These findings suggested that our gene chip study was performed adequately.

There have been several previous reports of gene expression analyses in asthma. Previous gene array studies have shown that CLCA121,22, SerpineB27,21,22, MUC5AC6,21, AGR221, CPA36,22, and tryptase6,22 are over-expressed in bronchial epithelial cells or the airway mucosa of asthmatics. Such studies have used tissues obtained from bronchial epithelial cells or the mucosa of healthy controls and subjects with allergic asthma. Genes such as NFκB p65 subunit (RELA) and scinderin (SCIN) were expressed at similar levels in our and these other studies, suggesting that the expression of certain genes in the airways is reflected in the peripheral blood cells of similar subjects. However, many genes were expressed differently in the present study than in previous reports. This may be attributable to the different tissue sources used, such as airways and peripheral blood cells. For mechanistic or pathogenesis studies, airway tissue, such as bronchial mucosa, would be the best target. However, obtaining airway tissues or airway epithelial cells using BAL or brushing is not convenient for asthmatics, especially those with severe asthma.12 Thus, PBMCs were used for the development of diagnostic gene markers in the present study.

First, we narrowed the number of candidate gene markers from the 170 genes that were initially selected, because 170 would be too many to handle. By applying more stringent criteria (P<0.001 and ≥5-fold change), eight candidate genes were selected for modeling. Using these genes, a multiple logistic regression analysis was performed to identify disease marker genes. Because we did not know the complex dependencies among the eight selected genes,23 we made all 255 (28-1) possible models using the genes and measured the P values to examine their validity. ROC curves and the AUC were measured to assess the predictability of the gene markers for asthma.

Of the 255 models, only 85 showed a P value <0.05. Furthermore, the number of genes analyzed in combination was inversely correlated with the P values. Those models made of more than three genes showed insignificant statistics when a P value <0.05 was applied, so we compared the AUC values only for those genes or combinations of genes having P values <0.05. After performing a multiple logistic regression analysis, the best model (MEPE, MLSTD1, TRIM37) and the second best model (CCT5, MEPE, TRIM37) were found to show almost perfect classification. The AUC values were 0.9928 and 0.9904, respectively. Considering that an AUC value of '1' indicates perfect classification, 0.9928 and 0.9904 indicate very high predictability for the markers. The sensitivity and specificity for the best model were 0.98 and 0.80, respectively, while the accuracy was 0.942. Thus, these models are sufficient to discriminate between asthmatics and normal controls. For three-fold CV, the results showed perfect classification. This means that the genes are very effective for independent data. Thus, the genes may be used as disease markers for asthma in PBMCs. The approach of selecting gene markers using ROC curves has been applied to select biomarkers in other diseases, including gastric cancer, knee osteoarthritis, and lung cancer.24-26

In previous studies, the most up-regulated gene and/or down-regulated gene have been selected and subjected to logistic regression analyses and measurements of the AUC from ROC curves. However, such models might not have the best discriminative performance as disease markers. In contrast, in the current study, we selected several candidate genes that satisfied our criteria. We applied a multiple logistic regression analysis for all possible (28-1) models and selected the one with the highest AUC. In our experiment, a model consisting of three genes showed the best performance.

We also analyzed the eight genes selected according to asthma severity. Our results show that TRIM37 was predictive of the occurrence of asthma, while NOX5, LMAN1, KNS2, and MLSTD1 were predictive of mild asthma, and CCT5 was predictive of moderate-to-severe asthma.

Additionally, we analyzed the effect of asthma mediation on the expression levels of the three candidate genes. We found no difference in terms of mRNA level between the treated and untreated subjects using t-tests and multiple logistic regression analyses, indicating that there was no effect of mediation on the expression of the three genes.

To date, the functions of the genes used as diagnostic markers of asthma have not been characterized. MEPE is an inhibitor of mineralization in situ and in cell cultures where altered expression is associated with oncogenic osteomalacia and hypophosphatemic rickets. MLSTD1 is related to fatty alcohol synthesis in mammals.27 CCT5 is associated with cell proliferation, the cell cycle, morphological changes, and apoptosis.28 TRIM37 encodes a member of the tripartite motif (TRIM) family, which is involved in diverse cellular functions, such as developmental patterning and oncogenesis. We searched various databases to determine any relationship between asthma and the selected genes; however, we were unable to find any connection. Thus, the functional relation of these genes with asthma remains to be solved; however, these genes do highlight the complex mechanism(s) underlying the pathogenesis of asthma. Furthermore, they comprise a set of potential biomarkers for use in discriminating between asthmatics and normal individuals. The three genes presented here as biomarkers were selected by computational methods alone. Additional functional studies are needed to determine the exact mechanism or indirect relationship among these genes with asthma.

In summary, a 35K whole-genome mRNA expression study of PBMCs from asthmatics and normal controls resulted in the selection of 170 genes based on criteria of P≤0.001 and a fold change ≥2. By applying more stringent criteria (P≤0.001 and a fold change ≥5), we selected eight genes as candidate asthma biomarkers. Using the AUCs from ROC curves, we identified a genetic biomarker consisting of three genes (MEPE, MLSTD1, TRIM37) having a diagnostic accuracy of 0.9928 AUC, with 98% sensitivity and 80% specificity. This marker may prove to be useful diagnostically after validation using a larger number of samples.

ADDITIONAL FILES

Additional file 1 - Table S1, Excel file, 170 genes (P≤0.001 and fold change ≥2).

Additional file 2 - Table S2, Excel file, 255 models (AUC, P-value).

Additional file 3 - Table S3, Excel file, gene list according to asthma severity.

Additional file 4 - Figure S1, PowerPoint file, ROC curves and AUCs for all models.