Diagnosis of Kawasaki Disease Using a Minimal Whole-Blood Gene Expression Signature

Key Points Question Can Kawasaki disease be accurately diagnosed on the basis of the pattern of host gene expression in whole blood? Findings In this case-control study of 606 children (404 in the discovery cohort; 202 in the validation cohort), a 13-transcript signature was identified that accurately discriminated Kawasaki disease from comparator febrile diseases in discovery and validation cohorts. Meaning A diagnostic blood test based on measurement of small numbers of host gene transcripts might enable early discrimination of Kawasaki disease from other infectious and inflammatory conditions.


Pathogen diagnosis
Viral diagnostics were undertaken on nasopharyngeal aspirates using immunofluorescence (RSV, adenovirus, parainfluenza virus, influenza A+B) and nested PCR (RSV, adenovirus, parainfluenza 1-4, influenza A+B, bocavirus, metapneumovirus, rhinovirus/enterovirus). Bacterial cultures included blood, CSF, urine and tissue sites. Pneumococcal antigen was measured in blood and urine, and bacterial DNA was detected by meningococcal and pneumococcal PCR.

Diagnostic process in febrile controls
Patients had a diagnostic work-up as directed by the clinical team, including blood count, blood chemistry, Creactive protein (CRP), blood urine and throat swab cultures; cerebrospinal fluid analysis and chest radiographs were performed where appropriate. Multiplex PCR was used to detect common respiratory viruses in nasopharyngeal aspirates or throat swabs, and common viruses in blood. Once the results of all investigations were available, patients were assigned to diagnostic groups using predefined criteria (Figure 1), as follows: Bacterial infection: Patients assigned to the bacterial pathogen group had a bacterial pathogen (gram-positive coccus or gram-negative bacillus) identified by culture or by molecular techniques in a sample from a sterile site (blood, CSF, pleural space, joint, urine), and a clinical syndrome in keeping with the identified bacterial species. This group included patients with and without viral co-infection. Children diagnosed with other bacterial infections (for instance mycoplasma, pertussis, mycobacteria) were not included in this group. No threshold for inflammatory markers was set for this group, as identification of bacteria in a sterile-site sample was taken as conclusive evidence for a confirmed bacterial infection. Uncertain bacterial or viral infection: When children with an acute febrile illness and features of infection could not be assigned confidently to one of the above groups, they were labelled as 'Uncertain Bacterial or Viral'. Children in this group had inconclusive features of bacterial or viral infection, negative microbiological findings or absent virological investigations, a syndrome inconsistent with their microbiological findings, inflammatory markers inconsistent with other clinical features of their illness, or insufficient clinical data for confident coding in another group. Patients in this group did not have bacterial infection detected at a sterile site, and some patients did have detectable virus.
Other inflammatory syndromes: a) Henoch-Schönlein purpura (HSP) was diagnosed in children presenting with palpable purpura, typically over the buttocks and extensor surfaces in association with abdominal pain, arthralgia or renal abnormalities (hematuria and proteinuria); b) Juvenile idiopathic arthritis (JIA) was defined according to International League of Associations for Rheumatology 1 . Patients with JIA included i) treatmentnaïve and ii) active-exacerbation/smouldering. © Wright VJ. JAMA Pediatrics.

Microarray pre-processing -The Discovery Dataset
Background subtraction and robust spline normalisation (RSN) were applied to the raw expression data using the R package lumi 2 . Sample outliers were assessed by Principal Component Analysis (PCA). One sample from a Kawasaki patient, was a clear outlier on PC1 and was removed from the analysis (eFigure 1).
The samples in the discovery dataset were randomly assigned to ten different folds conditional on equal numbers of each comparator group (KD Kawasaki Disease, DB Definite Bacterial, DV Definite Viral, U infections of uncertain bacterial or viral aetiology, JIA juvenile idiopathic arthritis, HSP Henoch-Schönlein purpura, HC healthy controls). Two folds (20%) were reserved as the test set and the remaining eight folds made up the training set. As a diagnostic test for KD would be of most value early in the course of the illness, we developed our signature using only samples from patients at 7 or fewer days of fever in the discovery cohort.

Microarray pre-processing -The Validation Dataset
The validation dataset was constructed by merging two gene-expression datasets: one with acute and convalescent Kawasaki samples 3 and one with bacterial and viral infections 4 . All convalescent samples had ESR (erythrocyte sedimentation rate) levels less than 40mm/hr and all acute samples were taken within ten days of onset of illness. Background subtraction and RSN normalisation were applied to the two datasets separately in the R package lumi 2 . At this stage, there were differences between the cohorts. This is evident from a PCA plot which shows that PC1 clearly distinguishes samples by batch (eFig. 2a). We therefore employed the ComBat 5 method to remove batch effects. Two binary covariates were passed to ComBat which assigned samples to three groups -healthy, KD and other diseases. The Kawasaki convalescent samples were assigned as healthy. A PCA after ComBat shows samples from both batches overlap on a plot of PC1 against PC2 with no significant batch effects (eFig. 2b). T-tests for difference in Principal Component values between healthy controls and convalescent KD patients were non-significant (PC1 p=0.32; PC2 p=0.98).

Model estimation
Before model estimation probes were pre-filtered to identify robustly expressed transcripts with log2 fold change ≥1 between the relevant disease groups. This was implemented by selecting probes meeting all of the following criteria in the training data: 1. Probes measured on both V3 and V4 Illumina Beadchips 2. Robustly expressed transcripts: for each probe, we calculated the proportion of samples in each comparator group for which the detection threshold p-value<0.01, and selected those probes for which this proportion was > 80% in at least one disease group 3. The majority of Kawasaki patients were recruited in UCSD. To ensure probe selection was not biased by batch effects emanating from UCSD, we excluded probes which showed association with recruitment at UCSD at p<0.05 in a linear model conditional on age in months and all disease groups which also included non-KD patients recruited from UCSD (DV, U, KD and HSP) 4. log2 fold change (conditional on age) was calculated between Kawasaki and each other comparator group; we took forward those probes with |log2 fold change|>1 for at least one of these comparisons The functions lmFit and eBayes in the R package limma 6 were used to calculate probe association statistic used in steps (3) and (4) above.

Discovery using Parallel Regularised Regression Model Search (PReMS)
We used PReMS, to derive a parsimonious gene-expression signature, which balances small transcript number with accurate discrimination 7 . The PReMS method used the pre-filtered transcripts that were robustly expressed with log2 fold change ≥1 between groups. PReMS uses cross-validation and selects the optimal model size as the one which minimises the out-of-sample log-likelihood. In the analysis presented in this paper, 20 crossvalidation folds were used.

Calculating model accuracy
The area under ROC curves, and corresponding confidence intervals of the models' application to the test and validation datasets were calculated using the R package pROC 8 .
Results for each patient were summarised as a Disease Risk Score (DRS) to determine the accuracy of classification by the 13-transcript signature, and the optimal threshold-cut-off for classification as KD or not KD, based on training set data, was determined according to Youden's J statistic by the point in the ROC curve © Wright VJ. JAMA Pediatrics.
that maximizes the distance to the identity line (maximum of (sensitivities + specificities)) 9 . The same threshold was used in accuracy calculations for the validation data.
Confidence intervals (CI) for sensitivity and specificity were calculated using Jeffrey's method. Jeffrey's method is derived from a Bayesian perspective in which the underlying proportion of interest is assigned the non-informative Jeffrey's reference prior: Beta( ½, ½ ) 10