Development of a blood-based gene expression algorithm for assessment of obstructive coronary artery disease in non-diabetic patients

Patient selection and clinical methods

All patients were clinically referred for angiography and angiograms were performed based on local, institutional protocols. The first microarray cohort of 198 subjects (88 cases and 110 controls) was derived from the Duke University CATHGEN registry, a retrospective blood repository, enrolled between August 2004 and November, 2005 [16]. For CATHGEN patients, clinical angiographic interpretation defined cases as ≥75% maximum stenosis in one major vessel or ≥50% in two vessels and controls as <25% stenosis in all major vessels. Clinical inclusion and exclusion criteria were described previously and included both diabetic and non-diabetic patients [13]. All CATHGEN patients gave written informed consent and the study protocol was approved by the Duke University IRB.

The second microarray cohort of 210 subjects (105 case: control pairs, matched for age and sex) and the RT-PCR algorithm development cohort (210 cases and 430 controls) were part of PREDICT, a multi-center US study of patients referred for coronary angiography (http://www.clinicaltrials.gov, NCT00500617). For PREDICT patients, core laboratory QCA reads (Cardiovascular Research Foundation New York) were used for case: control classification. Cases had ≥50% stenosis in at least one major coronary vessel and controls <50% stenosis in all major vessels.

Subjects from PREDICT were eligible if they had a history of chest pain, suspected anginal-equivalent symptoms, or a high risk of CAD with no known prior MI, revascularization, or CAD. Detailed inclusion/exclusion criteria have been described [6]. Diabetic status was defined by clinical identification, blood glucose (non-fasting ≥200 or fasting ≥126), rorhemoglobin A1c, (≥6.5), or diabetic medication prescription. Complete blood counts with differentials were obtained for all patients. PREDICT patients gave written informed consent, and the study protocol was approved by the Institutional Review Boards.

Blood collection, RNA purification and RT-PCR

Whole blood samples were collected in PAXgene^® tubes prior to coronary angiography, according to the manufacturer's instructions, then frozen at -20°C. For the CATHGEN samples RNA was purified as described (PreAnalytix, Franklin Lakes, NJ), followed by quantitative analysis (Ribogreen, Molecular Probes, Eugene, OR). For the PREDICT samples an automated method using the Agencourt RNAdvance system was employed.

Correlation between gene expression and cell type distributions

Correlations with complete blood counts and database gene expression analysis (SymAtlas, http://biogps.gnf.org) were used to identify highly cell-type selective genes. In addition, whole blood cell fractionation by density centrifugation or through positive antibody selection followed by RT-PCR was performed on specific cell fractions (see Additional file 1).

Statistical methods

All statistical methods were performed using the R software package.

Microarray methods

Microarray samples were labeled and hybridized to 41K Human Whole Genome Arrays (Agilent, PN #G4112A) using the manufacturer's protocol. For PREDICT microarrays all matched pairs were labeled and hybridized together to minimize microarray batch effects. Microarray data sets have been deposited in GEO (GSE 20686).

Normalization

Agilent processed signal values for array normalization were scaled to a trimmed mean of 100 and then log2 transformed. Standard array QC metrics (percent present, pair-wise correlation, and signal intensity) were used for quality assessment, resulting in 3 of 198 CATHGEN and 12 of 210 PREDICT samples being excluded.

Array analysis

For the CATHGEN array, logistic regression (unadjusted and sex/age adjusted) was used to assess gene expression association with case: control status. For the PREDICT array, given the paired design, conditional logistic regression was used. False discovery rates were used to account for multiple comparisons. BINGO was used to assess enrichment of gene ontology terms in the set of 655 genes[17]. A hyper-geometric test was used to identify overrepresented terms; results were corrected for multiple testing using Benjamini & Hochberg False Discovery Rate (FDR) correction.

Gene selection

Genes for RT-PCR were selected based on statistical significance, gene ontology pathway analysis, and literature support. Hierarchical clustering based on gene: gene correlations ensured that RT-PCR genes represented multiple clusters. Normalization genes were selected based on low variance, moderate to high expression, and no significant association with case: control status, sex, age, or cell counts. Cell-type genes were selected based on known literature or correlation to known cell-type specific markers.

PCR methods

Amplicon design, cDNA synthesis, and RT-PCR were performed as previously described [6, 13]. All PCR reactions were run in triplicate and median values used for analysis. Clinical/demographic factors were assessed for CAD association using univariate and multivariate logistic regression. Gene expression association with CAD and other clinical/demographic factors was assessed by robust logistic regression (unadjusted and sex/age adjusted) [13].

Algorithm development and validation

Hierarchical clustering was used to group genes using a correlation cutoff. Clusters were reduced to meta-genes[18] and normalization genes based on correlation structure, known biology, and cell count correlation. In general, a meta-gene was a set of 1-4 genes from a specific cluster, chosen to best represent the cluster center using a parsimonious number of genes. Genes within meta-genes were equally weighted with one exception (Additional File 1). For meta-gene pairs with high correlation and opposite disease regulation, ratio terms (differences on the log scale) were defined. Meta-genes independently associated with outcome were selected by the LASSO method, with sex by meta-gene interactions allowed during variable selection [19].

The final algorithm was fit using Ridge regression [20], where the outcome variable was case: control status and the predictors the LASSO-selected meta-genes and sex-specific age terms. Sex was a binary predictor, and age a linear predictor with separate slopes for men, women >60, and women <60 (the slope for women age < 60 was estimated to be approximately 0 and thus was set to 0 in the final algorithm). The LASSO was fit using the glmnet package in R and ridge regression was fit using the Design package in R; in both cases the shrinkage parameter lambda was estimated using 10-fold cross validation. Model performance was estimated using leave-one-out cross-validation.

Gene selection

Analysis of algorithm development cohort: clinical and gene expression factors

The most significant clinical/demographic factors for CAD association were age, sex, chest pain type and neutrophil count. Age and sex were independent risk factors for CAD (Table 3) and showed significant gene expression correlation. Chest pain type was also a significant independent risk factor (p = 0.0004), especially in men, but was gene expression independent. Neutrophil count was significantly correlated (positively or negatively) to expression of 93 of 113 RT-PCR genes, and was significantly associated with CAD in men (p = 0.049) but not women (p = 0.77). Neutrophil-associated genes showed both up and down regulation with CAD status, whereas lymphocyte-associated genes were generally down-regulated. There was significant gender-specific regulation of neutrophil correlated genes (men 40/42 genes up-regulated, women, 41/42 down-regulated) whereas lymphocyte gene down-regulation was gender independent.

Hierarchical clustering of the 113 PCR genes resulted in 18 correlated clusters (Figure 5, Table 2), a significant fraction of which could be mapped to cell-type specific gene expression groups, with finer correlation substructure within the lymphocyte and neutrophil associated genes. There were 3 lymphocyte subgroups representing T-cells (clusters 1,2,3), B-cells (cluster 3) and NK cells (cluster 12). Three neutrophil subgroups were also identified: previously described neutrophil genes (IL8RB, S100A8, S100A12, TXN, BCL2A1; cluster 13, 16); newly identified up-regulated neutrophil genes (CLEC4E, CASP5, TNFAIP6; cluster 16) and down-regulated neutrophil genes (KCNE3, TLR4, TNFRSF10C; clusters 13, 14) [13]. Cluster 8 appears to be eosinophil specific. The 26 genes in clusters 4-7 and 9-11 did not have clear cell-type association.

Algorithm derivation and performance

Based on gene expression correlation clustering and cell-type analyses, 15 meta-genes and 3 normalization genes were defined as inputs for model variable selection (Table 2, Figure 6). Selection by the LASSO method and further weight penalization by Ridge regression resulted in the final, locked algorithm, comprising 20 CAD-associated genes and 3 normalization genes in 6 meta-gene terms (Figure 6), where each term represents a ratio of meta-genes or meta-gene to normalization genes. The algorithm score is calculated as described (Additional file 1) and was defined as the predicted regression model value.

The estimated cross-validated algorithm AUC in ROC analysis in the PREDICT development set was 0.77 (95% CI 0.73 to 0.81) (Figure 7); prospective validation in an independent PREDICT validation set of 526 patients (192 cases, 334 controls) yielded an AUC of 0.70 (95% CI = 0.65 to 0.75) [6].

Biological significance of algorithm terms

The use of correlated meta-genes as building blocks for the algorithm is significantly reflective of gene expression cell-type specificity. The algorithm genes are expressed selectively in multiple types of circulating cells including neutrophils, NK cells, B and T-lymphocytes, [27], supporting roles for both adaptive and innate immune responses in atherosclerosis [10].

A role for neutrophils in both the early and later stages of atherogenesis has recently been suggested, especially in connection with hyperlipidemia [28, 29]. Algorithm term 1 is a ratio of neutrophil expressed meta-genes that are up and down regulated with CAD (Figure 6). This term may particularly reflect neutrophil apoptosis, as Caspase-5 is increased with CAD, whereas TNFRSF10C, an anti-apoptotic decoy receptor of TRAIL, is decreased [30]. Term 2 genes up-regulated with CAD are also expressed largely by neutrophils and likely reflect both innate immune activation, (S100A8 and S100A12), [31] and a cellular necrosis response (CLEC4E) [32]. S100A8 and S100A12 are up-regulated in chronic inflammatory conditions, including asthma, rheumatoid, and inflammatory arthritis, perhaps reflecting a more general pathophysiological signal, consistent with increased CAD in disorders such as rheumatoid arthritis [33, 34].

Term 2 is highly correlated with the signature previously identified by us [13] and includes the most significant gene from that work, S100A12. This term is normalized in a sex-specific manner, perhaps reflecting sex-specific differences in the significance of neutrophil counts in CAD and MI [26]. In men normalization to RPL28 which is strongly expressed in lymphocytes, reflects the neutrophil to lymphocyte ratio, which is prognostic for death or MI in a CAD population [14]. In women normalization to AQP9 and NCF4, two CAD insensitive neutrophil genes, permits assessment of neutrophil up-regulation of the S100s and CLEC4E, independent of neutrophil count.

Term 3 consists of 2 NK cell receptors, SLAMF7 and KLRC4, normalized to T-cell specific genes (TMC8 and CD3D). SLAMF7 may specifically activate NK cell function, while inhibiting B and T cells [35]. KLRC4 is also likely involved in NK cell activation [36]. NK cells have been associated with atherosclerosis in both mouse models and humans, and reduced lymphocyte counts associated with cardiac events [14, 37].

Term 4 is a gene expression based measure of the B/T-cell ratio. The roles of both T and B cells in atherosclerosis development are complex; mouse models have shown B cells to be both athero-protective and more recently, atherogenic [38–40]. In this study apparent upregulation of B-cell specific genes is correlated with CAD, perhaps supporting the latter. The last two terms, based on AF289562 (AF2) and TSPAN16 are genes of unknown function that will require further work to clarify their role in atherosclerosis.

Previous work by Sinnaeve and coworkers also examined peripheral blood gene expression in a coronary disease population [12]. As noted by these authors, there is little overlap between their genes and the signatures identified in our previous study [13] or this one. A number of significant differences in the study populations (restricted age range, combining two sex specific cohorts) in their study may have contributed to this. In addition, there are differences in both RNA isolation methodology and microarray platforms. Further work is needed to resolve these issues.

Algorithm development

For algorithm development, as described above, we used an approach that minimized the effect of any single gene by using meta-genes as building blocks [18, 41] Penalized stepwise logistic regression (LASSO) selected significant meta-genes from a 640 patient data set which greatly exceeded the number of candidate variables (15 meta-genes), reducing the likelihood of over-fitting. Further, in order to minimize over-weighting of individual terms, meta-gene coefficients were penalized using Ridge regression. An alternative approach would have been to use a combined two-step shrinkage method such as the elastic net [42]. Although correlations with specific cell types was a key observation, recent methodologies described for deconvoluting gene expression data sets from complex cell mixtures might have led to improved results [43].

The cross-validated model AUC was 0.77 (95% CI 0.73 to 0.81), suggesting that the algorithm score was a significant CAD predictor. A decrease to an AUC of 0.70, with overlapping confidence intervals (95% CI = 0.65 to 0.75), was observed in the independent validation set [6]. This decrease may reflect an over-optimistic cross-validation estimate, as we did not re-select terms during each iteration. Ultimately, the validation results provide the most informative measure of a model's prediction accuracy.

Limitations

Although this is one of the largest studies examining gene expression in peripheral blood in CAD patients and has yielded a specific algorithm for the assessment of CAD status, it has several limitations.

From a clinical perspective, diabetics and patients with known chronic inflammatory disorders were excluded. The differences observed between diabetics and non-diabetics with CAD could be due to differences in the molecular pathophysiology of the disease, medications, or some combination of the two. In addition, although race was not an independent risk factor after adjustment for age and sex, the number of minority patients was low, so conclusions with respect to them are significantly underpowered. The use of a dichotomous angiographic endpoint does not account for variations in disease burden or external remodeling, and is not a measure of ischemia. Finally, the contribution of atherosclerosis in other vascular beds is outside the scope of this study, but may be important in asymptomatic high-risk individuals.

From a cellular and gene expression perspective, the relative ease of obtaining peripheral blood cell RNA is counter-balanced by not directly interrogating changes in the diseased vascular wall. Another complementary approach could be to examine secreted proteins in the blood that might reflect endothelial or vascular dysfunction. Finally, given the chronic nature of atherosclerotic disease, it is likely the gene expression signature observed reflects a response to disease rather than the underlying cause.