Blood RNA biomarkers of coronary artery disease

Clinical parameters: A total of 113 patients were entered into the study, with 1 patient excluded post-consent due to undetected exclusion criteria of heart failure. Thus, 112 patients were available for analysis, and RNAseq was successfully completed on 96 patients. The following example analyzed the RNAseq result of the first 48 patients, as well as the overall clinical parameters of the entire cohort, in order to determine whether the time and expense of RNAseq was justified on the remaining patients.

Clinical prediction model: The accumulated clinical risk factors yielded a range of 5 to 12 points per patient. A receiver-operator curve (ROC) was computed using JROCFIT between the clinical risk prediction compared to a true positive for CAD based on 70% stenosis on angiogram. The fitted ROC curve yielded a C=0.663 with an accuracy of 60.7%, sensitivity of 57.7%, and specificity of 61.6%, for all 112 patients, and similar values of X, Y, Z for the 96 RNAseq+pts. These values are slightly lower than the published values for this model (C=0.71−0.77), possibly due to differences in the threshold for stenosis (50% there vs 70% here), and comparable to a Diamond-Forrester prediction model (C=0.64).

Analytical parameters: The yield of RNA and the number of reads per patient, and per group did not vary systematically. With 16 additional patients failing at various technical quality control steps in RNA sequencing, RNAseq was completed on 96 patients. Most failures occurred at 3 key steps: low yield from RNA purification or ribosomal depletion, inefficient cDNA synthesis, or low yield of reads from RNAseq. Because of the high loss during ribo-depletion (>90%), there was not always adequate RNA to repeat the RNAseq process.

Sources of variation in RNA yield. Patient blood samples collected with either Paxgene or Tempus RNA preservation tubes show a surprisingly large variation in the RNA yield, with Tempus generally producing higher yield. In the present Tempus-preserved samples, total nucleic acid yield (prior to DNAse treatment) ranged from (0.6-35 μg/3 ml whole blood, mean=10.6 μg). The correlations between RNA yield and any single blood cell count were quite weak, with the strongest correlation to absolute lymphocyte count (r=0.55 with N=112 (R 2 =0.31) ( FIG. 1 ). While one outlier seems to drive this correlation, it actually has a small impact on the correlation (r=0.45 w/o). This demonstrates that the white cell counts account for only about 30% of the variability in RNA yield, and suggest there may be important differences in the recovery of RNA between patients.

Initial identification of TRACs at midpoint: As shown in FIG. 2 , there is excellent linear quantitative ability of RNA levels over 25 log 2 orders of magnitude in RPKM without any artificial normalization of the data. To identify TRACs, the patients were divided into LOW (<20% stenosis), MID (20-69%), and HIGH CAD (>70%) groups. To identify differentially expressed transcripts, we compared the LOW vs. the HIGH CAD group, and required a >1.4-fold change and t test p<0.05 uncorrected, resulting in 238 transcripts, highlighted in dark gray in FIG. 2 . This combined fold-change/t test strategy has been established in large, multicenter control studies using spiked samples as a reliable approach to identify true differences. The 238 transcript list is presented in Table 1, wherein the relevant transcripts are identified with several different descriptors that are common in the art, especially the UCSC ID (ie uc002rtm.2), maintained by the University of California at Santa Cruz, and the RefSeq ID (ie NM_001113494), maintained by the US National Center for Biotechnology Information (NCBI), as well as the consensus Gene Symbol, and Gene Description.

Classification algorithm: Using ANOVA, we re-analyzed all 3 groups and identified 84 differentially expressed transcripts (Table 2). Using Genespring GX13 software, a Support Vector Machine (SVM) algorithm was built which was 100% sensitive and specific. Leave One Out Cross Validation (LOOCV) applied to the same samples found the model was 87.5% accurate at diagnosing CAD (33.3% is random). For comparison purposes, using the CardioDx transcripts (Elashoff, et al., 2011) to build a comparable SVM yielded only 36.7% accuracy (CAD 37.5%, LOW 45.1%, MID 10%, w/33%=random). This SVM of the CardioDx transcripts, however, does not include the algorithm they employed. For reference, classifying the 48 patients by gender, identified 34 transcripts, mainly from X and Y chromosomes, and yielded an SVM prediction model that was 100% accurate.

Example 2. Description of RNA Biomarkers from Transcriptome Alignment of RNAseq Data from 96 Patients

The entire 112 patient cohort was subjected to RNA sequencing using the Helicos (now SeqLL) platform. After RNA purification, ribosomal RNA depletion, and RNA sequencing, it was determined that 96 patient samples were suitable for further analysis.

Discovery and validation sets. Those 96 samples were randomized into a ‘Discovery Set’ of 45 patients and a ‘Validation set’ of 51 patients.

Discovery Set. The 45 patients were subdivided into groups according to the degree of CAD as determined upon the same day as blood was drawn for RNAseq. The groups are described above and are composed of LOW (0-20% stenosis), MID (21-69%) and HIGH (≥70%). To identify a relatively small set of biomarkers, with relatively high confidence, the LOW group was compared to the MID+HIGH groups via a combined T test (p<0.001 uncorrected, and a minimum fold change of 1.5) resulting in 59 transcripts (see Table 3).

A smaller list was constructed of just the top 5 changed transcripts based on both fold-change and a scientific analysis of the transcripts. A PLSD model built on only the 5 transcripts in the DISC group was 84% accurate (91% for LOW, 77% for MID+) (see Table 4).

Because the LOW group includes true ‘normals’ with essentially 0% stenosis and some patients with small amounts of CAD, we conducted another analysis on the Discovery group that re-arranged the patients by ‘normal’ (0%) vs ‘abnormal’ (>0%) and conducted a similar T test yielding 500 transcripts (see Table 5).

Overall prediction models built with the 96 patients. To gain the most powerful model, we included all samples and divided them as in the Discovery/Validation design into LOW vs MID/HIGH (MID+). This yielded 48 patients in each group. Using this design, transcripts were first filtered to exclude transcripts with RPKM<0.01 in 70% of samples, and then TRACs were identified by one-way ANOV of the LOW vs MID+ patients at an uncorrected p value of 0.001, yielding 198 transcripts (see Table 6). These 198 transcripts contain some duplicates that cover essentially the same RNA transcript, yielding 169 non-redundant, unique RNA transcripts.

A PLSD model build on these 198 transcripts was very accurate at a discrimination showing an overall accuracy of 98.9% (100% for LOW, 97.9% for HIGH). This remained fairly robust even with N-fold internal validation yielding overall accuracy of 80% (77% for LOW, 83% for MID+).

Filtering for higher level expression. By filtering the 198 transcripts for those which had a >20 th percentile level in both groups, the list was narrowed to 96 transcripts with relatively higher absolute expression than the 198 transcripts. This however, did not improve the predictive ability of the PLSD model built on it, with overall accuracy of 93% (92% for LOW, 94% for MID+), but still produced a quite powerful test, with fewer transcripts (see Table 7).

Identifying a small set of predictive transcripts. To narrow the list even further, the 198 transcripts predictive in the overall cohort were compared for overlap with the transcripts predictive in the DISC set using the LvM+ criteria (Table 3, 59 transcripts). Common to both lists were 14 transcripts that were used to build a new model for the overall cohort. A PLSD model remained powerful with an overall accuracy of 86% (87% for LOW, 85% for MID+) which dropped to 80% overall accuracy in a LOOCV (see Table 8).

Identifying the most comprehensive list of predictive transcripts. For both diagnostic and potential therapeutic uses, there is value in knowing the most comprehensive list of transcripts that might be employed. Thus, by relaxing the ANOVA stringency to p<0.005 (vs. 0.001 in the prior examples), a larger set of 1132 transcripts is identified (see Table 9).

Interpreting the TRAC signature: Both the midpoint and full analysis showed that almost all of the differentially expressed transcripts were down-regulated in the CAD patients, a pattern that rarely occurs in RNA expression analysis, where there is typically a balance between increased and decreased transcripts. Furthermore, the changes are essentially all of the same magnitude (mean=˜1.7 fold), while we would normally expect a range of magnitudes. After extensive analysis of the TRACs, the most likely explanation is that TRACs are markers of a particular cell type, and that cell type is reduced in abundance in persons with developing or existing CAD. An extensive literature documents that two relatively selective cell type can be altered in abundance in patients with CAD.

TRACs as markers of a modified population of circulating cells. One way to obtain the uniform decrease in the transcripts is if these transcripts are preferentially associated with a subset of cells, and that subset is reduced in the CAD patients. Certain transcripts identified in the midpoint analysis, such as CCL28, CD2BP2, CD3E, CD5, CD81, KDM2B, LEF1, and MIB2 have been previously associated with various stem and progenitor cell populations. CD81 is interesting because of its well-characterized role in induced pluripotent stem cells (iPSCs), mesenchymal stem cells (MSCs), and hematopoetic stem cells (HSCs).

There is a substantial literature that consistently reports reductions in circulating progenitor cell (CPC) populations in patients with stable CAD or preclinical atherosclerosis. Hypercholesterolemia, smoking and obesity are associated with reduced levels of circulating progenitors. Conversely, circulating endothelial progenitor cells (EPC) are increased in acute MI cases. However, it is unlikely that a decrease in EPC numbers, which are rare (<1%), could cause the substantial shift in RNA levels detected in whole blood. Rather, both EPC decreases and the RNA profile might be two aspects of some larger change in the cellular composition of blood in atherosclerosis.

TRACs as putative markers of Treg: A second potential explanation for TRACs is that there are known reductions in the T regulatory (Treg) subset of lymphocytes in atherosclerosis. An extensive literature documents reduced Treg abundance, and a relative imbalance in Treg vs T effector (Teff) cells in patients with CAD. To test for the potential changes in Treg, we examined the regulation of known Treg markers in the 3 CAD groups. As shown in FIG. 3 , FoxP3, and other markers such as CD4, CD25, ETS1, and RUNX1 mRNA expression levels show a stepwise decrease in expression from LOW, MID, to CAD. By comparison, the expression of an irrelevant marker such as the prostaglandin E receptor 3 (PTGER3) does not show this CAD-related trend.

Treg imbalance is sufficient to create the observed RNA expression pattern. To determine whether a reduction in Treg cell counts in blood would be sufficient in magnitude to produce the observed changes in RNA levels, 8 publications that reported Treg % in normal and CAD patients, such as unstable angina or acute coronary syndrome (ACS), were reviewed and the change in Treg % computed. The average reduction in Tregs, typically defined as CD4+CD25+CD127low by flow cytometry, was 30.25%, which would translate to a 1.43-fold difference in Treg RNA levels, assuming that these markers are essentially unique to Tregs. Thus, the 1.47-fold reduction in mRNA for the consensus FoxP3 marker is quite consistent with the reported reduction in Treg cell numbers in CAD.

There are not pre-specified pathways or gene ontologies specific to the Treg maturation and differentiation process. Recent review articles by Li and Rudensky, and others, were adapted to create a functional pathway of Treg signaling to compare with the pattern of RNA changes observed. A striking overlap was observed in the DGEs and known Treg signaling.

Biomarkers of CAD can also be defined by their biological relationship to the Treg pathways and thus we identified a series of markers with known relationship to Tregs differentiation and which varied systemically with CAD status. These include FOXP3, CTLA4, ZAP70 and others (see Table 10).

Using a gene list of just 3 transcripts (FoxP3, Runx1, IL2RA/CD25) associated with relative specificity for Treg cells, but not covered in other TRAC lists, a prediction model was built on the 96 patient cohort whereby 48 subjects are classified as LOW and 48 as MID or HIGH (MID+). Using only these 3 transcripts, a PLSD model showed 60% overall accuracy (67% for LOW, 54% for MID+). However, upon N-fold validation this drops to 46% overall accuracy (47% for LOW, 44% for MID+). This suggests that while the TRACs may be related to Treg numbers or function, the TRACs perform better at predicting CAD than known markers of Treg cells.

Other published and patented works have identified transcripts with predictive value for CAD based on Affymetrix array technology and PaxGene blood RNA preservation tubes (U.S. Pat. Nos. 9,122,777 and 8,914,240, the contents of which are incorporated by reference). Prior studies using microarray profiling of PaxGene preserved blood have led to one blood RNA-based diagnostic (Corus CAD by CardioDx). While the CardioDx test (C=0.745) improves slightly on a purely clinical risk assessment (C=0.732), the present RNAseq test proved to be significantly more accurate (C=0.870). Further, the CorusCAD test by CardioDx uses the age and gender of the subject in calculating their risk score, while the present TRAC method is solely based on gene expression levels.

To compare the methods more directly, we identified 17 best matched transcripts from those patents and publications, without knowledge of the expression levels or predictive values in the current dataset. The expression levels of those 17 CardioDx markers did not differ between Low and MID+ patients in our RNAseq dataset (see Table 11). Nonetheless, we attempted to build a classification model using them, and these 17 transcripts in a PLSD model had 73% overall accuracy (75% for LOW, 71% for MID+), but the accuracy fell to 63.5% in N-fold validation, and 50% would be random chance. Thus, while it is difficult to make an accurate head to head comparison, because they are fundamentally different methods, the present TRAC biomarkers are better at predicting CAD than existing RNA biomarkers.

Based upon the medium-sized TRAC lists, specifically Table 6 (198 transcripts, 169 unique), a smaller panel of transcripts was selected based upon specific criteria for their suitability in human diagnostics. These criteria related to the putative association of the transcript with Treg cells, a high level of expression, and relative independence from close gene family members. Using these criteria, 7 transcripts were selected: DLG1, DGKA, ICOSLG, IKZF4/Eos, SMYD3, TCF3, and TRIM28. The levels of these transcripts were extracted from the RNAseq data, expressed as a percentage (ratio) of the mean expression of that transcript in the entire cohort, to minimize overweighting of highly expressed transcripts, and then the normalized levels of the 7 TRACs were added to make a cumulative score, as shown in FIG. 4 (right panel). The classification of the patients into CAD or low CAD was highly accurate (80.2%) based upon this TRAC score alone, with an ROC area (AUC, C-stat) of 0.873. By comparison, using clinical predictors alone, ( FIG. 4 left panels), only 54% accuracy (with 50% as random chance) and a C-stat of 0.636 (with 0.500 as random chance) could be achieved.

As discussed, other inventions have described gene panels for CAD that differ completely from the present panels. Those studies employed a preservative for blood RNA called PaxGene, which is based on a cationic detergent, while the present studies employed Tempus preservative, which is based upon guanididium-like denaturants. To determine whether the preservative had an effect on the detection of the TRACs, blood from the same donors were drawn into either PaxGene or Tempus tubes, RNA was prepared by the manufacturer's directions, and then equal amounts of RNA were quantitated by droplet digital PCR, after reverse transcription with iScript (BioRad). The results, shown in FIG. 5 , clearly indicate that all of TRACs measured (7) were detected better in the Tempus preservative than in the Paxgene preservative by a range of 1.5 fold (TRIM28) to >7-fold (ICOSLG). Thus, the method of RNA preservation is important to the detection of these TRACs, and the use of differing preservatives may explain some differences between the results.

The accurate diagnosis of CAD via a blood-based test offers several major advantages in clinical medicine. First, it avoids invasive testing, such as cardiac catheterization, as well as allowing more judicious use of imaging resources, such at CT and MR angiography. Secondly, it can improve diagnosis of CAD in rural areas worldwide, where advanced imaging methods are unavailable. Finally, the individual biomarkers within the biomarker panel can serve both as therapeutic targets, and markers to monitor the efficacy of new or existing therapies, such as statins. For instance, Treg numbers have been shown to be responsive to statin therapy, and so it might possible to use TRACs to monitor therapies.

The connection between the immune system and atherosclerosis is well known. Blood components, especially monocytes/macrophages, neutrophils, lymphocytes, and platelets mechanistically contribute to the development of CAD. The present results are consistent with the extensive evidence that CAD is associated with changes in the Treg/Teff ratio, which are thought to be mechanistically related to atherosclerotic progression.

Mechanistically, the present results suggest that changes in the immune system may correlate with the presence of CAD, which has been a substantial line of investigation for many years, progressively building connections between lipid imbalances, inflammation, microbiome changes, and autoimmunity in atherosclerosis. There is a large and fairly consistent literature demonstrating changes in the Treg/Teff ratio in patients with CAD, and those observed cellular changes would be consistent in both direction and magnitude with the current changes in mRNA expression. One interpretation of the beneficial effects of statins is that in addition to lowering LDL cholesterol, statins can induce FoxP3+ Treg cells, via modulation of TGF-ß signaling. Beyond the reproducible clinical correlation, experimental manipulation of Treg levels in mouse models of atherosclerosis lends credence to a potentially causal relationship. Furthermore, it has been suggested that an immune modulatory approach may offer therapeutic potential for atherosclerosis. It is plausible that a stimulatory effect of statins on Treg number could contribute to their anti-atherogenic effect, an effect that has not yet been reported for PCSK9 inhibitors.

An additional relationship can be seen between Treg dysfunction and atherosclerosis through the well-known incidence of atherosclerosis in various autoimmune diseases, but most notably in the case of systemic lupus erythematosus (SLE). While the relationship between Treg and SLE is complex, there is a general consensus that deficient Treg numbers and/or function is one element of SLE, and thus, might also be a component of SLE-associated atherosclerosis. Likewise, psoriasis and psoriatic arthritis, which are associated with Treg imbalances, have well-established associations with atherosclerosis. The immune-CAD connection is further strengthened by an apparently causal relationship in immune-mediated transplant arteriosclerosis.

Potentially the strongest evidence for the immune-CAD connection derives from the proven efficacy of rapamycin, and related compounds which are antibiotics/immunosuppressants, to block restenosis. Rapamycin is known to increase Treg numbers and function at clinically relevant levels.

Many of the TRACs have known relationships with Treg function. Foremost, the forkhead-box transcription factor FoxP3 is considered the definitive marker for the Treg subset, and is thought to transcriptionally regulate a set of transcripts involved in Treg function. FoxP3 is itself epigenetically controlled by promoter demethylation to allow stable expression in Treg, and in turn regulates Treg-specific transcription via known promoter binding sites. Other studies indicate that two isoforms of diacylglyceral kinase (DGKA, DGKZ) have been implicated in T cell anergy, and DGKZ has been implicated in the generation of natural Tregs via modulation of the NFkB signaling through c-rel. The FoxP3 transcription factor is known to interact with Ets transcription factors, and prior work on aortic aneurysm transcript profiling demonstrated that this atherosclerotic condition was associated with coordinated changes in Ets transcription factor-dependent signaling in the aortic tissue. Zap70, the zeta chain associated protein, is preferentially phosphorylated in Treg cells during TCR engagement.

While the underlying cause of the Treg imbalance is currently unknown, there is increasing evidence that changes in the microbiome could alter the type of short-chain fatty acids released into the circulation, thereby modulating both GI and circulating Treg differentiation.

In addition to composing a diagnostic test, it is likely that TRACs can help to elucidate some of the underlying components of CAD, and possibly new targets for therapeutics, in particular, by modulating the Treg/Teff ratio. Through high-throughput screening, dozens of FDA compounds that stimulate Treg generation can be rapidly assessed.

The transcript profile identified here can be adapted to creating a screening test for CAD in asymptomatic individuals, where early detection would allow risk modification and drug therapy. Using microarrays, we've used a similar blood RNA approach to generate new insights into the origin of other cardiovascular diseases, such as adriamycin cardiotoxicity and aspirin resistance (AR), and made the innovative observation that AR is potentially due to an autoimmune antiplatelet syndrome.

There are certain limitations to the present studies: 1) Confounding variables: the TRAC signature could be related to a risk factor or drug treatment that differs between groups. However, based on the clinical variables collected, we cannot identify a co-morbid condition or prescription drug use that would differ sufficiently to create this effect, but it is difficult to completely rule it out. 2) Specificity: It is possible that the TRACs would detect disease in any artery, but even if the test was detecting atherosclerosis in other arteries, it would still have tremendous diagnostic value. The combination of chest pain and a positive TRAC profile in blood would be strong evidence that catheterization, or other follow-up diagnostics, is warranted. Likewise, numbness in the legs, combined with a positive TRAC test, would be sufficient to justify peripheral vascular ultrasound. The Cardiogram group observed that their predictors for CAD also predicted aortic atherosclerosis, suggesting that blood markers are likely to be general predictors of atherosclerosis. Prior work has suggested that the location of atherosclerosis is effectively an endophenotype associated with specific risk factors in a given patient, which one can speculate may be related to the particular autoantigens that are prevalent in particular vascular beds.