A genome-by-environment interaction classifier for precision medicine: personal transcriptome response to rhinovirus identifies children prone to asthma exacerbations.

OBJECTIVE: To introduce a disease prognosis framework enabled by a robust classification scheme derived from patient-specific transcriptomic response to stimulation.
MATERIALS AND METHODS: Within an illustrative case study to predict asthma exacerbation, we designed a stimulation assay that reveals individualized transcriptomic response to human rhinovirus. Gene expression from peripheral blood mononuclear cells was quantified from 23 pediatric asthmatic patients and stimulated in vitro with human rhinovirus. Responses were obtained via the single-subject gene set testing methodology "N-of-1-pathways." The classifier was trained on a related independent training dataset (n = 19). Novel visualizations of personal transcriptomic responses are provided.
RESULTS: Of the 23 pediatric asthmatic patients, 12 experienced recurrent exacerbations. Our classifier, using individualized responses and trained on an independent dataset, obtained 74% accuracy (area under the receiver operating curve of 71%; 2-sided P = .039). Conventional classifiers using messenger RNA (mRNA) expression within the viral-exposed samples were unsuccessful (all patients predicted to have recurrent exacerbations; accuracy of 52%). DISCUSSION: Prognosis based on single time point, static mRNA expression alone neglects the importance of dynamic genome-by-environment interplay in phenotypic presentation. Individualized transcriptomic response quantified at the pathway (gene sets) level reveals interpretable signals related to clinical outcomes.
CONCLUSION: The proposed framework provides an innovative approach to precision medicine. We show that quantifying personal pathway-level transcriptomic response to a disease-relevant environmental challenge predicts disease progression. This genome-by-environment interaction assay offers a noninvasive opportunity to translate omics data to clinical practice by improving the ability to predict disease exacerbation and increasing the potential to produce more effective treatment decisions.

BACKGROUND AND SIGNIFICANCE

The advent of increasingly rapid and cost-effective genome-wide measures (omics methods), is transforming our understanding of human pathophysiology and has stimulated numerous initiatives in precision medicine. In particular, our understanding of single-gene disorders is advancing rapidly. However, genes alone do not determine the onset or progression of complex diseases – environmental factors, including lifestyle and specific stimulus exposure, also contribute substantially. Already, a number of stimulation-based bioassays have been developed for diagnosis and treatment in clinical medicine. These include in vivo bronchial challenge tests (eg, inhaled methacholine challenge for asthma) and in vitro antibiotic-sensitivity bacteriograms, though available tests rely on available pathognomonic (diagnostic) biomarkers and have not yet leveraged the power of big-data genomics. Tissue-stimulation assays coupled with omics measurements provide a minimally invasive opportunity to understand an individual’s personal genome-by-environment response to a stimulus relevant to disease state or risk of progression. While cohort statistics can identify the most common stimulus-provoked molecular responses between clinical populations, there are no unbiased tools for assessing a whole transcriptome response that can be scaled down to a single patient without requiring inference from reference sets based on well-powered, cross-patient comparisons. Tailoring patient treatment according to the results of a given in vivo or in vitro stimulation thus remains an unmet challenge. Nonetheless, diseases with complex inheritance that are precipitated by acute environmental challenge are conspicuous candidates for assay development to provide clinical insight into using this strategy. Focusing on dynamic gene expression changes within a subject additionally gets to the heart of the matter more effectively than conventional static classifiers derived from population-averaged expression.

We have used an approach that assesses changes to transcript expression between 2 paired samples from a single subject in response to stimulus exposure. This is followed by analyses of biological process involvement through gene set enrichment. Expression changes are determined using a suite of single subject–oriented N-of-1-pathways analytical tools uniquely designed for paired data analysis without use of a reference set. In the accompanying experimental assay focused on asthma, live peripheral blood mononuclear cells (PBMCs) were drawn from patients, aliquoted, and exposed in vitro to either human rhinovirus (HRV) stimulus (stimulated) or vehicle control (unstimulated). We have called the provoked transcriptional response a viral regulome specific to a personal genome, or a “virogram.”

Epidemiological studies have identified viral infection as the most prominent environmental risk factor for severe asthma exacerbation episodes at all ages, with HRV infection (aka the common cold) specifically contributing to 48%–65% of pediatric cases., These severe asthmatic exacerbations are a major cause of morbidity, require acute hospitalizations, and in rare cases can lead to death. Though the role of viral infection is well established, not all children with respiratory infections become asthmatic or experience significant exacerbations later in life. This supports an additional significant role for non-Mendelian genetic predisposition interacting with viral exposure. Identifying individuals at significantly increased risk for exacerbation before a clinically severe episode occurs would provide benefit at the individual case management and larger population health levels.

We hypothesize that clinically useful pathway classifiers can be derived from dissecting the interaction of an individual patient’s aggregate genomic risk with relevant environmental or laboratory assay exposure (genome × environment classifier; G × E classifier). In this study, we designed a robust dynamic transcriptional signature that can identify asthmatic patients at risk of exacerbation with a single blood test.

MATERIALS AND METHODS

The framework of our study is shown in Figure 1 . We used a published microarray dataset measuring mRNA expression from an in vivo HRV exposure study in non-asthmatic volunteers (training set; see Table 1). These data were used to build classifiers based on differentiating responses associated with symptomatic and asymptomatic infection. The G × E single-subject classifier built from our analysis of individualized within-subject responses was compared to cohort-based classifiers constructed using more standard methods. In our validation set (see Table 1), PBMCs were isolated from asthmatic patients and stimulated in vitro with the same virus to provoke transcriptional changes. RNA was harvested and transcript expression quantified by microarray from paired unstimulated and stimulated cells. Single subject–derived and cohort-derived viral response classifiers were then explored for their ability to predict exacerbation episodes in the asthmatic validation set.

Figure 1.

Overview of the experimental and analytic design of the study. We hypothesize that a classifier based on the pathway-level transcriptional responses that differ between symptomatic and asymptomatic responses to HRV infection in healthy patients can predict which asthmatic patients will have exacerbations over a 1-year follow-up period, based on those patients’ transcriptomic responses to an in vitro HRV stimulation assay. Panel A illustrates the development of a classifier using innovative features. Both shared and unique genetic and nongenetic variables influence transcript expression and healthy vs disease state within an individual. Exposure to a relevant stimulus (here, HRV infection) reveals relevant pathway gene sets whose genome and environmentally informed responses (G × E response pathways) can be used to predict individual prognosis. Panel B describes the development of classifier from the training set using data from PBMCs of healthy volunteers exposed to HRV in vivo. For each patient, paired microarrays analyzing gene expression before and after HRV exposure were compared using N-of-1-pathways analysis to identify significant Gene Ontology biological process (GO-BP) features describing each response. Responses in asymptomatic patients were then compared to responses in symptomatic patients to develop the classifier. Panel C describes the laboratory stimulation assay of asthmatic patients' PBMCs in this study. As in the training set, paired microarrays for each patient were used to determine the response using N-of-1-pathways. The classifier developed in panel B was then applied to individual responses and used to predict recurrent exacerbation.

Table 1.

Descriptions of datasets

Dataset		Training set	Validation set
Purpose		Learn classifier	Predict outcome
Source	Authors	Zaas et al. Cell Host and Microbe 2009	Present study
	Data Source	GSE17156 (downloaded 9/17/2014)	GSE68479
Platform	Microarray	Affy. Human Gene U133A 2.0	Affy. Human Gene 1.0ST33297
	Probe	22277	33297
Protocol		Inhaled HRVPBMC samples drawn before and 48 h after HRV inoculation	PBMCs isolated, then incubated in vitro with HRV (stimulated) or vehicle (unstimulated)
Subjects	Total	19 healthy adult volunteers	23 pediatric asthmatic patients
		□ 10 symptomatic for common cold	□ 12 recurrent exacerbations of asthma (hospitalizations and/or emergency room visits)
		□ 9 asymptomatic	□ 11 no exacerbation of asthma
Samples	Samples:	38 RNA microarrays	46 RNA microarrays
		□ 19 PBMC drawn prior to infection	□ 23 PBMC unstimulated
		□ 19 PBMC drawn after HRV inhalation	□ 23 PBMC HRV-stimulated in vitro

Gene set analysis for classifier construction was performed using Gene Ontology (GO) Biological Process annotations., To provide additional clinical translational value, we also introduced a novel visualization tool using star plots to help physicians make treatment decisions from the mechanism-anchored metrics revealed for each patient.

Description of training set to derive responsive features for development of an asthmatic exacerbation classifier

For our training set, we downloaded published microarray data from the Gene Expression Omnibus (GEO; GSE17156) and Table 1. These data were derived from PBMC samples collected from 19 healthy, non-asthmatic volunteers exposed in vivo to HRV via intranasal inoculation. Each volunteer donated 2 blood samples, the first drawn at baseline prior to challenge, the second drawn 48 h after viral exposure. Of these 19 volunteers, 10 individuals developed symptoms of acute respiratory infection after exposure to HRV (symptomatic), while the other 9 had no symptoms (asymptomatic). Note: though the original study included 20 subjects, only 19 pairs of matched mRNA expression datasets were deposited in GEO.

Asthmatic subjects accrued for validation set and clinical phenotype determination

All human studies were carried out according to protocols approved by the Institutional Review Boards of the University of Arizona and the University of Wisconsin (IRB 0600000587 for the TREXA clinical trial and IRB06035604 for the BADGER clinical trial). Twenty-three pediatric asthma patients who previously participated in the Treating Children to Prevent Exacerbations of Asthma (TREXA) and Best Add-on Therapy Giving Effective Responses (BADGER) clinical trials were asked to participate in a yearlong extension study 2 years after the clinical trials ended. Based on exacerbation status, the pediatric asthma patients were classified in 2 groups during the clinical trial: (1) 11 patients who did not receive oral corticosteroids during the year of the extension study were coded as no exacerbation, and (2) 12 patients who received 2 or more courses of oral corticosteroids during the year of the study were coded as RE. (For more details on patient characteristics, see Supplementary Table S1.) Importantly, none of the 23 asthma patients received steroids during the month prior to the extension study when blood samples were collected. We statistically tested for obvious confounding factors between the 2 phenotypes, such as age, gender, family history, etc., and found no evidence of group-specific demographic or clinical differences (Supplementary Table S1).

PBMC sample extraction, HRV stimulation, and RNA expression quantification in the in vitro study to develop the validation set

PBMCs were isolated from blood samples collected from the 23 pediatric asthma patients enrolled in the extension study. Each individual’s PBMCs were divided into 2 cultures, nonstimulated and stimulated with human rhinovirus RV-16 provided by James E Gern, University of Wisconsin-Madison. Specifically, 100 µl of the working stock RV-16 was added to 2 mL of PBMC suspension, yielding a final concentration of 2.5 × 106 plaque-forming units/mL. These samples were then incubated at 37°C and 5% CO2 for 24 h (HRV-stimulated). In parallel, the second aliquot of patient PBMCs was incubated under the same culture conditions, but exposed only to vehicle (unstimulated). Following this, RNA was isolated from each sample, quality was assessed, and then RNA was amplified, tagged, and hybridized on Affymetrix Human Gene 1.0 ST microarrays according to standard operating procedures. Data were deposited in GEO (GSE68479).

Processing of microarray datasets

Within both the training and validation sets, each pair of matched samples (HRV-stimulated and unstimulated) was normalized using Robust Multiple-array Average (2 paired samples at a time to avoid bias in single-subject analyses) from Affymetrix Power Tools. For the training set, each pair of samples derived from each patient was obtained in a log2-normalized format and further corrected for batch effects using ComBat. In all, 12 157 genes were found in common across both datasets.

Pathway (gene set) annotations using gene ontology biological processes

To obtain mechanistic profiles of responsive pathways for each patient, virus-responsive genes were aggregated into gene sets according to GO annotations of biological processes (GO-BPs), using the org.Hs.eg.db package (Homo sapiens) of Bioconductor, available for R statistical software. Hierarchical GO terms were retrieved using the org.Hs.egGO2ALLEGS database (downloaded on May 15, 2013), which contains a list of GO-BP terms along with all of their child nodes. Gene sets comprising a minimum of 15 and a maximum of 500 genes were included, as documented in our previous studies.,,,,

Transcriptome and pathway (gene set) analysis of gene-by-environment interactions using N-of-1-pathways framework

To determine personally responsive gene set involvement for each patient in each dataset (training and validation sets, Table 1), the N-of-1-pathways Wilcoxon framework was applied using each patient’s gene expression microarrays derived from paired samples (HRV-stimulated and nonexposed control PBMCs). Briefly, N-of-1-pathways is designed to create personal profiles of responsive pathways (gene sets) and is based on 3 principles: (1) the individual patient is the sole source/unit of observation, and any statistical descriptions and inferences are intended to relate only to that patient; (2) significance and interpretation are derived from gene sets; and (3) gene set–level information is used to answer questions of clinical importance. Principle 1 allows for detection of individual signals that traditional cohort-level studies may overlook. Principle 2 anchors the results in mechanism, and this affords dimension reduction and interpretation. Principle 3 provides quantitative and qualitative measures to address questions related to patient care.

As described previously,, each patient’s paired gene expression was transformed via N-of-1-pathways Wilcoxon into pathway-level metrics (annotated using GO-BP). Specifically, for each pathway gene set, a Wilcoxon signed-rank test was conducted between the gene expression measurements. This test produces a P-value that quantifies whether the central gene expression is different between the pair of samples (higher or lower expressed). In effect, observations from a patient in both datasets (training and validation) were transformed from paired, whole-transcriptome measurements to a simple list of pathway-associated P-values for each patient.

The P-values obtained from the training set were further manipulated to serve as inputs for the classifier algorithms. The output of N-of-1-pathways was transformed into a matrix of P-values representing all patients and pathways (23 patients by 3055 pathways). This matrix was further transformed into a ternary matrix by considering a pathway gene set as responsive if its nominal P-value was <5% and signed to indicate the direction of response. Specifically, the ternary matrix contains the values −1 (significantly lower in HRV-stimulated), 0 (not significant responsive), and 1 (significantly higher).

Principal component analysis

The principal component analyses (PCAs) were computed using the FactoMineR package in R (with default parameters) based on the ternary matrix derived from the training set. A projection over the first component was used to assess the statistical significance of the separation shown by the PCA between 2 phenotypes (eg, symptomatic vs asymptomatic). In particular, a Mann-Whitney U test was used to assess the degree of separation of the 2 phenotypes. Projections over the first and third principal components were used to visualize the phenotypic separation.

Classifier modeling using gene-by-environment response pathways as features

The N-of-1-pathways ternary matrix (pathways × patients) that pertains to the analysis of the training set was used for training the classifier through machine learning procedures with feature selection. The training set was transformed into an .arff file using Weka v.3.6.11, and then fed into Weka algorithms for further analysis. Briefly, feature selection was performed on the training set (clinical assay ternary matrix) using a chi-squared test. Using the training data alone, chi-squared test of association was performed using the pathway-associated ternary values and the phenotype (symptomatic vs asymptomatic). The pathways were then ranked in ascending order of P-values. The pathway gene sets with the smallest P-values from the chi-squared test were selected for further modeling based on the best prediction performance by varying the number of features in the classifiers listed below (assessed via leave-one-out cross-validation; Supplementary Table S2). Five different classifier models were assessed for performance on the training set using the default parameters of Java software: random forest classifier, naïve Bayes, decision tree, support vector machine, and nearest neighbor. The random forest classifier model was further assessed for superior performance on the training set using receiver operating characteristic (ROC) curve based on prediction scores extracted from Weka. The area under the curve (AUC) and significance tests were computed with the default parameters of the verification package in R. The fully specified classifier derived from the training set was then applied once on the validation set, and the accuracy of predictions was measured from the ROC curve.,

Visualization of personal gene-by-environment interaction via pathways gene set responsive to stimulation using star plots

Star plots were computed in R using the stars function in the default graphics package. Each patient has an individual star plot, where each spike represents the response of a specific pathway (GO-BP; method: N-of-1-pathways framework). The spikes within each star plot represents the 20 gene set features of the classifier described in Table 4. Each individual star plot corresponds to a single validation set individual’s N-of-1-pathways scores (n = 23 patients; Table 1). An individual pathway is considered upregulated when the N-of-1-pathways test computes a concordant increase in expression of the genes within the pathway in the HRV-stimulated sample as compared to the control. The median radius of the star plot corresponds to no changes in the pathway, while upregulated pathways are assigned values in the white range and downregulated pathways are assigned values in the gray range. Direction-adjusted (signed) P-values were transformed into 4 states: (1) significantly upregulated (reaching the border of the star plot), (2) significantly downregulated (center of the star plot), (3) upregulated not reaching statistical significance (midrange radius of the white zone of the star plot), and (4) downregulated not reaching statistical significance (midrange radius of the gray zone of the star plot). For grouping of related pathways in the star plot, similarities between GO-BP terms were determined using Jiang’s information theoretic similarity, which ranges from 0 (no similarity) to 1 (perfect match).

Table 4.

Gene ontology (GO) gene sets of the G × E classifier obtained in the training set

Class of GO-BP gene sets	GO ID	GO term
I. Acquired Immunity	GO:0002768*	immune response–regulating cell surface receptor
	↪GO:0002429	immune response–activating cell surface receptor
	↪GO:0050851	antigen receptor–mediated signaling pathway
	↪GO:0050852	T cell receptor signaling pathway
II. Innate Immune Response	GO:0045087*	innate immune response
	↪GO:0060337	type I interferon-mediated signaling pathway
	(p)↪GO:0034340	response to type I interferon
	(i)↪GO:0071357	cellular response to type I interferon
III. Morphogenesis	GO:0050807	regulation of synapse organization
	GO:0001658	branching involved in ureteric bud morphogenesis
	GO:0060688	regulation of morphogenesis of a branching structure
	GO:2000027	regulation of organ morphogenesis
IV. Response to Stimulus	GO:0043279	response to alkaloid
	GO:0050795	regulation of behavior
	GO:0009581	detection of external stimulus
V. Chromatin Organization	GO:0006325*	chromatin organization
	↪GO:0016568	chromatin modification
	↪GO:0016569	covalent chromatin modification
	↪GO:0016570	histone modification
	GO:0006913	nucleocytoplasmic transport

Twenty GO biological processes (GO-BPs) responsive to HRV stimulation were selected by the best classifier after evaluation in the validation set (Materials and Methods section). GO-BPs were organized into 5 categories by an unbiased information theoretic similarity score (Materials and Methods section), and these classes were manually assigned representative names. Of note, some GO-BPs were ordered according to their GO hierarchy when available, with the parent term annotated with an asterisk and located above the child term. For example, GO:0002768 is the parent of GO:0002429. Legend: ↪ = “is a” (parent-child relationship); (p)↪ = “part of” (parent-child relationship); (i)↪inferred as “is a” (parent-child relationship).

RESULTS

Generation of the predictive classifier for rhinovirus stimulation response

As an exploratory analysis, we first computed the principal components resulting from the transformed ternary representation of all N-of-1-pathways scores from the training set (see Materials and Methods section; Principal Component Analysis). The idea was to investigate whether these scores carry information to discriminate asymptomatic from symptomatic individuals prior to fitting a complex classifier. Indeed, the first component from the PCA based on the single-subject dynamic change of gene expression within pathways gene set enabled a rough classification of the 2 clinical phenotypes (asymptomatic vs symptomatic; Figure 2 , right). Moreover, the bivariate relationship between the first and third components reveals a striking pattern of phenotypic separation (Figure 2, left). In light of these results, the single-subject scores of the 20 most relevant pathway dynamic expressions produced by the chi-squared test based on feature selection were used to train a classifier (Materials and Methods section). The random forest classifier achieved the best AUC of the ROC curve (86.7%) using leave-one-out cross-validation of the training set (see Materials and Methods section; classifier modeling; Supplementary Table S2).

Figure 2.

Metrics derived from responsive pathways discriminate asymptomatic from symptomatic subjects in the training set of in vivo HRV-stimulation data. In panel A, principal component analysis was conducted using responsive gene sets derived from each subject in the training set (Figure 1B and Materials and Methods section). The scatter plot on the left illustrates the bivariate relationship of the first and third principal components for each patient’s ternary-represented N-of-1-pathways scores derived from paired samples of PBMCs collected at baseline and after HRV exposure. Each point represents a subject in the training set through a linear combination of pathway gene set–level scores that explain the maximal variation in the data (see Materials and Methods section for details on the ternary representation and PCA construction). The first and third principal components show 2 clusters emerging that separate asymptomatic from symptomatic individuals. Thus, N-of-1-pathways scores are associated with the phenotype of interest. On the right, side-by-side box plots display the first principal component scores among asymptomatic and symptomatic subjects, and this component alone significantly dichotomizes the 2 phenotypes (Mann-Whitney U test, P = .0069). Panel B lists the 20 Gene Ontology biological processes used as features in our classifier, organized according to broad biological function categories.

Classifier derived from pathways responsive to gene-by-environment interactions, under HRV stimulation, enables asthmatic clustering and a clinical risk assessment representation of asthma exacerbation for each subject

We applied the random forest classifier to the dynamic expression change of pathways stimulated in vitro on PBMC samples of the validation set, which showed >70% accuracy in identifying children with and without recurrent asthma exacerbations (Figure 3A ). We also tried 4 other classification methods to assess the robustness of the signal regardless of the type of classifier (Table 2). An unsupervised principal component analysis of these pathway scores also showed a natural split between these 2 asthmatic phenotypes (Supplementary Figure S1). To further explore the dynamic change in expression within the validation set, we used a heatmap displaying an unsupervised clustering of individuals to visualize transformed N-of-1-pathways scores (Supplementary Figure S2). Further, to assess whether the classifier was dependent on HRV stimulation or due to a general transcriptomic viral response, we developed 2 additional classifiers using (1) influenza A and (2) respiratory syncytial virus infection data (GSE17156; Zaas et al.) as our training sets. We employed the same methodology as described to identify classifier features and attempted to use these to classify our validation set of HRV-stimulated PBMCs from asthmatic patients. The best classifier for either alternative virus was ineffective with every validation set patient-predicted as nonrecurrent, thus yielding no predictive stratification (52% accuracy and 61% AUC for influenza; 48% accuracy, 85% AUC for respiratory syncytial virus; see details in Supplementary Table S3). Moreover, our pathway-level metrics of HRV response outperformed prediction in the validation set compared to conventional mRNA expression–derived classifiers (Supplementary Table S4). Indeed, gene-level classifiers, built in a similar fashion from the HRV-stimulated training set, achieved a best accuracy of only 52% and AUC of 41% on the validation set. Again, this poor AUC is explained by a lack of prediction stratification, as all patients in the validation set were predicted to have REs.

Figure 3.

Fully specified classifier derived from responsive pathways in training set performs accurately in independent asthmatic validation/virogram set to predict exacerbation status. Panel A shows receiver operating curve (ROC) of the G × E classifier conducted in the validation set (overall accuracy, 74%; sensitivity, 75%; specificity, 73%). In panel B, the star plots illustrate the level of response to HRV stimulation for each pathway in that patient. The classifier is designed using 20 pathways, with each radial line representing the score of a pathway. The area above and below the gray zone represents upregulation and downregulation, respectively, of any given pathway (see Materials and Methods section for complete details). In panel C, each star plot represents a single subject, with label appearing above (eg, Subject 16 = SUB16). The star plot is located in the quadrant of the contingency table that represents the performance of the G × E classifier on predicting the clinical progression from a specific asthmatic patient’s data. This classifier prediction applied to the HRV-stimulation assay recapitulated the clinical progression in 17 out of the 23 asthmatic patients; 6 were misclassified: 3 false positives (SUBM1, SUB09, SUB23) and 3 false negatives (SUB2, SUB23, and SUBM2). One can straightforwardly identify that innate immunity (mauve) is upregulated in every asthmatic patient and does not contribute to the classification in the validation set. On the other hand, observed upregulation in acquired immunity (orange) pathways could be used to correctly classify 19 subjects.

Table 2.

Performance of 5 different pathway response classifiers built on N-of-1-pathway scores in the training set (in vivo) and validated without bias on N-of-1-pathway scores in the validation set (in vitro). The Random Forest (bold) classifier showed the highest scores and was chosen for all other analyses.

Classifier	#Features (GO-BP terms)	Classification performance in the validation set
		Accuracy	Sensitivity	Specificity	Precision	AUC	TP	FP	FN	TN
Random Forest	20	73.9	75	72.7	75	71.2	9	3	3	8
Naïve Bayes	20	73.9	75	72.7	75	67.4	9	3	3	8
Decision Tree	20	65.2	75	54.5	64.3	64.8	9	5	3	6
Support Vector Machine	20	69.6	75	63.6	69.2	69.3	9	4	3	7
Nearest Neighbor	20	69.6	66.7	72.7	72.7	69.7	8	3	4	8

The random forest and naïve Bayes classifiers showed the highest metrics on the validation set, and random forest was selected for further optimization.

TP: true positive; FP: false positive; FN: false negative; TN: true negative

We assessed the optimal number of pathway features to be included in our random forest classifier using the chi-squared method for feature selection. This showed robust performance by 20 to 30 features (Table 3). Guided by parsimony heuristics, we selected 20 features.

Table 3.

Our classifier performed optimally across all metrics using 20 gene-set (GO-BP term) features

Classifier	#Features	Classification performance in the validation set
		Accur	Sens.	Spec.	Prec.	AUC	TP	FP	FN	TN
Random Forest	10	47.8	0	100	NA	68.2	0	0	12	11
Random Forest	20	73.9	75.0	72.7	75.0	71.2	9	3	3	8
Random Forest	30	69.5	58.3	81.8	77.8	65.5	7	2	5	9

TP: true positive; FP: false positive; FN: false negative; TN: true negative

Focusing on each patient’s PBMC response to HRV stimulation, we developed a star plot representing personal profiles of these 20 pathways (Figure 3B) to report results as a visually interpretable risk assessment of asthma exacerbation. Comparing clinical records of individuals who experienced at least 2 exacerbations in a year to their retrospective molecular assessment of asthma exacerbation, 6 of the 23 individuals were misclassified (SUBM1, SUB09, SUB13, and SUB2, SUB23, and SUBM2; Figure 3C). Acquired immunity (orange), morphogenesis (blue), response to stimulus (green), and chromatin organization (yellow) responses were highly discriminatory between true positives and true negatives, while all of the asthmatic patients in the validation set showed substantial activation of innate immunity (mauve). Thus, while innate immunity–associated pathways are likely involved in common defense mechanisms to HRV and symptomatic vs asymptomatic infection of healthy subjects, it may be that asthmatic disease state alone predisposed all patients in our validation set to high innate activation. This is perhaps a limitation of our use of healthy patients in the training set, though it does highlight that our classifier was fully specified on the training set alone before being applied to the asthmatic patient data.

Nonetheless, despite limited additional discriminatory power from innate immunity–associated pathways, chromatin organization and acquired immunity clearly show a high propensity to be significantly upregulated in children who experienced RE relative to those who did not. As such, the star plot provides a visual tool for physicians to observe responsive mechanisms common and distinct between phenotypes, as well as individualized profiles.

DISCUSSION

In this study, we demonstrated the effectiveness of classifiers built using biological process pathways that are differentially responsive to clinically relevant stimulation (G × E classifier). Specifically, our case study uses HRV-stimulation response data from healthy patients to classify and predict the risk of exacerbation in asthmatic patients whose PBMCs were subjected to in vitro exposure of the same virus. Although the classifier was trained on healthy volunteers, we were able to accurately classify 70% of asthmatic patients as low risk (no exacerbation) or high risk (recurrent exacerbation) by using the same pathway gene sets, despite small numbers of patients in both the training and validation cohorts. A classifier trained on patient-specific dynamic responses at the pathway level (N-of-1-pathways: 74% accuracy, 71% AUC, 2-sided P = 0.039; Table 2 and Figure 3A) significantly outperformed a cohort-derived classifier trained on transcript-level variation between the sets of asymptomatic vs symptomatic patients at 1 static point (52% accuracy, 41% AUC; Results section). This increase in predictive power was achieved through 4 major mechanisms: (1) control of confounded environmental variables within each subject that impact global transcriptional activity (eg, preexisting habits, medications, and additional chronic diseases); (2) control of preexisting baseline differences in gene or pathway expression across patients that may introduce bias or lead to spurious conclusions; (3) reduced statistical dimension through aggregation of responsive transcripts into enriched gene sets with common biological properties prior to cross-group comparisons; and (4) a focus on analysis of clinically relevant elicited responses within each patient rather than total variability. Within a given clinical class, between-patient heterogeneity in specific transcript expression resulting from unique genetic architecture can also be subsumed into pathway scores if the summative effect on a functional pathway remains in common.

Unsurprisingly, these results suggest that dynamic changes in the personal transcriptome are more predictive than static measurements. This agrees with what is observed in clinical practice, where provocation tests are routinely shown as more correlated to the patient’s response to therapy than unprovoked phenotypes (eg, bronchial histamine challenge test with spirometry is superior to peak flow measurements). The stimulation assay we outline here corresponds to such a provocation assay, revealing a dynamic phenotype at the transcriptome level.

Along with the improved classifier performance, the N-of-1-pathways strategy provides an additional benefit over conventional classifiers, as pathways in each patient are individually scored. This permits researchers or clinicians to meaningfully examine each patient’s individual responsivity independently (Figure 3C). As our knowledge of the mechanistic involvement of each pathway grows, these star plots could be used in aggregate not only to stratify patients according to their classified risk, but also to suggest a more effective precision medicine approach to disease management for any given individual.

It is of great interest that several of the GO biological processes identified were highly relevant for the viral challenge under study. Acquired and innate immunity, for example, are critical elements in acute responses to viruses. Moreover, type I interferons (ie, α and β interferons) play a major role in activating antiviral responses. Activation of innate immunity in particular is also a hallmark of hypersensitivity in asthmatic patients,, so it is both expected and reassuring that our classifier included multiple GO terms related to this process. However, that same hallmark of hypersensitivity in asthma reduced the discriminatory power of innate immunity processes in classifying already asthmatic patients as exacerbating vs non-exacerbating (Figure 3). Not many classifiers of asthmatic exacerbation have been published, but when we compare the biological processes revealed by our classifier to one identified by Bjornsdottir et al. (2011), we see both similarities and differences. The classifier in their study used a larger dataset of 118 adult asthmatic patients, with RNA samples obtained during asthma-quiescent vs mid-exacerbation in the same cohort of patients, using the natural course of each patient’s disease progression rather than a provoked stimulus response assay. Their asthmatic exacerbation classifier, which was run through Ingenuity Pathway Analysis, revealed entirely immune-mediated responses, including Toll-like receptor signaling, T cell receptor signaling, B cell receptor signaling, interferon signaling, interferon responsive factors, and IL15 signaling. While we could not do a direct comparison of GO pathways or underlying transcript participation in classification due to the proprietary nature of Ingenuity Pathway Analysis, these appear to show substantial convergence with the GO pathways of acquired and innate immunity that we observed (Supplementary Table S5). Notably, however, while innate immunity–associated genes were significantly represented in their study, our study found that activation of the innate immune response was a characteristic of all asthmatic patients exposed to HRV, with little predictive discriminatory power to classify exacerbation from non-exacerbation. In addition, the classifier built by Bjornsdottir did not substantially prioritize transcripts or pathways associated with morphogenesis, chromatin organization, or generalized response to stimulus, which were highly discriminatory in our analysis. We found that these may provide insights into additional exacerbation pathology–associated processes that occur in the subset of predisposed patients who are exposed to triggering stimuli.

Another notable result of this study was that it underscores the specific clinical relevance of HRV infection to asthmatic exacerbation, as had been repeatedly described in clinical and epidemiological studies. While observational associations could be biased by the ubiquity of HRV (the common cold) relative to other respiratory infections, we found that the transcriptional response of PBMCs to HRV exposure was predictive of the outcome in our asthma cohort, while a classifier based on influenza A–provoked responses was not. This also highlights the importance of choosing a highly clinically relevant stimulus during design of the bioassay component in order to optimize the chance of success.

Though the training set data were the same, the classifier developed here using single-subject metrics of responsive pathways (N-of-1-pathways scores) diverges substantially from the classifier produced by Zaas et al. Their classifier was developed by first selecting a common gene signature indicative of acute respiratory infection across 3 distinct viral challenges (HRV, respiratory syncytial virus, and influenza A). This gene signature was then used to differentiate influenza infection from bacterial infection or asymptomatic status in an independent clinical dataset. On the other hand, our classifier was designed to predict a different phenotype – exacerbated asthma in patients using responsiveness of their PBMCs to HRV in an in vitro assay (virogram). It remains to be seen if the common gene signature of viral infection derived by Zaas et al. is predictive of asthmatic exacerbation, and this could be determined in future studies.

Although previous studies attempting to create asthmatic exacerbation risk classifiers used differential expression of pathway gene sets to classify clinical phenotypes,, these did not use a matched sample design (eg, stimulated vs unstimulated) and thus were not able to assess individual dynamic changes as features and/or true pathway-level signals with underlying transcript heterogeneity. We showed that developing an unbiased G × E classifier based on the genomic profiles of healthy individuals experiencing rhinovirus stimulation enabled risk assessment of exacerbation in children with asthma. Classification of asthma exacerbation was successfully demonstrated using the dynamic change of expression in pathways of single subjects, whereas conventional methods of gene-level profiling failed.

LIMITATIONS AND FUTURE STUDIES

One limitation of the current study is the small cohort size in both our training (healthy) and validation (asthmatic) sets. Designing more accurate classifiers on a larger and/or asthmatic cohort may improve patient stratification. Accurate classifiers using N-of-1-pathways analytics are powered to identify dynamic changes in gene expression patterns arising from 2 samples of a single subject and offer additional insights on individual responsive pathways.,, These improved genomic signals at the subject level can be compiled for a cross-patient statistical analysis to uncover common pathways and genes as well as allow patient risk stratification. In contrast, many conventional transcriptome analyses operating on single and static measurement profiles utilize statistics that generally require larger cohorts to combat cross-patient variability and random noise in the raw expression data, unrelated to the measure of interest.

As this study was conducted as a proof-of-concept, the number of subjects was limited and did not allow for designing a classifier from groups of children with asthma. The current accuracy of the classifier is sufficient for a proof-of-concept, but was generated from symptomatic vs asymptomatic non-asthmatic subjects and then applied a classifier to assess propensity for exacerbations in children with asthma. Further, the PBMC samples were obtained after the asthma exacerbations occurred, so it is not possible in this study to determine whether these responsive pathways lead to or result from recurrent asthma exacerbations. In future prospective studies, we plan to generate both training and validation cohorts from groups of children with asthma, which should increase the accuracy of the classifier and its clinical utility. We are also improving the N-of-1-pathways framework to identify better responsive transcripts in each patient, which may provide greater insight into the molecular mechanism(s) underpinning each pathway of the classifier. While expression-level classifiers implicitly harbor co-transcribed and co-regulated transcripts, the GO terms of our G × E classifier comprise more obvious overlap of mRNAs (Table 4). We will investigate methods to identify more independent classifier features in future studies. Per standard operating procedures, PBMC RNA is extracted from PBMCs as 1 homogenized sample without counts of cell-type subpopulations, thus the RNA expression changes between 2 samples confound changes in expressed transcripts with changes to cellular composition of the sample. In future studies, deconvoluting the signal according to cell type could be addressed for cell type–specific transcription, using expression change and proportion of cell type change between paired samples.

CONCLUSION

Given the importance of monitoring the dynamic gene expression driven by genome-by-environment interplay that influences the direction of biological processes and the course of a disease, improvements are needed in the design of advanced experimental and analytics assays intended for individual patients’ point of care and disease assessment in order to tailor treatment. Our proposed framework addresses this challenge by differing from traditional classifiers that employ mRNA expression as the features. In particular, we introduce the use of the N-of-1-pathways framework to first identify the responsive pathways induced gene-by-environment interaction via a paired stimulated vs unstimulated assay, where the stimulation is relevant to the clinical phenotype of interest. These responsive pathways are then utilized as features in classification algorithms to predict the phenotype of interest. In the context of asthma, at least 48% of pediatric admissions for exacerbated asthma are due to HRV, and, as such, our in vitro HRV-stimulated assay was designed ab initio for a single-subject analytic framework and is seen to identify children prone to asthma exacerbations. The type of stimulated PBMC assay proposed could scale to predict responses to therapy of a large number of clinical conditions that are immune-mediated, or for which PBMC is a justified target tissue. For example, chemical or biological ligands rather than intact viruses could be used to stimulate specific pathways to study drug-targeting receptors. N-of-1-pathways analytics provides an informative, quantitative, and visual characterization of subject-specific response to stimulation, influenced by the interaction of genomic variability and stimulus exposure, and is, in principle, applicable not only to the transcriptome but also to other omics measures (eg, proteome). When applied to easily accessible tissue in a minimally invasive way, these types of in vitro stimulation assays and associated analytics could potentially supplant more invasive in vivo prognosis procedures that also assay responses to the environment in the emerging field of precision medicine.

Funding

FDM and YAL are supported in part by the Arizona Health Sciences Center and the BIO5 Institute. This study was funded in part by the following grants: K22 LM008308‐04, 5U10HL064307, the University of Arizona Cancer Center (P30CA023074), and the Center for Biomedical Informatics and Biostatistics of the University of Arizona Health Sciences.

Competing interests

None.

Contributors

MJH, DJJ, and FDM conceived the bioassay experiments; VG, YAL, IA, AB, JL, HL, and DS conceived and conducted the computational biology and high-throughput experiments; VG, AB, IA, YAL, GS, JB, FDM, DJJ, MJH, JL, HL, DS, and CK contributed to Materials and Methods; VG, IA, YAL, AB, AB, JB, FDM, CK, QL, and YS contributed to figures, tables, and additional files; VG, IA, YAL, JL, JB, DS, MJH, FDM, QL, and YS analyzed the results; GS, JB, YAL, IA, CK, VG, JL, DS, DJJ, and FDM wrote and revised the manuscript.

Data-sharing

All software/programs are provided at http://Lussierlab.org/publications/N-of-1-pathways.

SUPPLEMENTARY MATERIAL

Supplementary material is available at Journal of the American Medical Informatics Association online.

Click here for additional data file.