|Home | About | Journals | Submit | Contact Us | Français|
Gliomas are the most common type of primary brain tumors in adults and a significant cause of cancer-related mortality. Defining glioma subtypes based on objective genetic and molecular signatures may allow for a more rational, patient-specific approach to therapy in the future. Classifications based on gene expression data have been attempted in the past with varying success and with only some concordance between studies, possibly due to inherent bias that can be introduced through the use of analytic methodologies that make a priori selection of genes before classification. To overcome this potential source of bias, we have applied two unsupervised machine learning methods to genome-wide gene expression profiles of 159 gliomas, thereby establishing a robust glioma classification model relying only on the molecular data. The model predicts for two major groups of gliomas (oligodendroglioma-rich and glioblastoma-rich groups) separable into six hierarchically nested subtypes. We then identified six sets of classifiers that can be used to assign any given glioma to the corresponding subtype and validated these classifiers using both internal (189 additional independent samples) and two external data sets (341 patients). Application of the classification system to the external glioma data sets allowed us to identify previously unrecognized prognostic groups within previously published data and within The Cancer Genome Atlas glioblastoma samples and the different biological pathways associated with the different glioma subtypes offering a potential clue to the pathogenesis and possibly therapeutic targets for tumors within each subtype.
Primary brain tumors are an important cause of cancer mortality in adults and children in the United States (1). The molecular and genetic heterogeneity of gliomas undoubtedly contributes to the varied and often suboptimal response to treatment that is usually based on standard pathologic diagnoses (2, 3). Glioma diagnosis has been historically based on examining the cellular morphology of the tumor to assess its presumed cell of origin (astrocytic, oligodendroglial, ependymal) and surrogate markers of tumor aggressiveness (necrosis, nuclear pleomorphism, mitoses) to determine the tumor grade (4). Glioblastomas are the most common and aggressive gliomas and are thought to arise de novo (primary glioblastoma) or through the malignant transformation of lower-grade astrocytic and oligodendroglial tumors (secondary glioblastoma; ref. 5). Although some genetic aberrations and clinical characteristics (i.e., age) have been associated with each type of glioblastoma, currently, there are few distinguishable differences in the histopathology or prognosis associated with primary and secondary glioblastomas (6-9).
Histopathologic diagnoses are by nature subjective (10), and the diagnosis of glioma subtypes has historically been associated with significant intraobserver variability. Even when the diagnosis of a distinct subtype of glioma can be agreed upon (i.e., gliosarcoma, small cell glioblastoma), the molecular, genetic, or clinical relevance of such designations remains obscure. Although standard tumor morphologic observations and low throughput genetic studies have revealed some molecular characteristics relevant to particular histologic subtypes, the study of the molecular features of gliomas has only recently come to the forefront with the advent of high-throughput microarray technology (6, 8).
Gene expression profiles provides a transcriptomic snapshot of a biological phenotype and offers the opportunity for quantitative, reproducible evaluation of individual tumor biology (8, 11). Consequently, data derived from genomic-scale gene expression profiling facilitate the characterization of intertumoral variations and similarities (12, 13). Using this approach, several groups have recently attempted to identify glioma subtypes associated with particular molecular features. Although important steps forward, the findings of these studies have been limited by several methodologic constraints, including (a) incomplete coverage of whole-genome expression due to the usage of limited or outdated legacy microarray platforms (2, 3, 14, 15); (b) limited number of samples studied and/or incomplete coverage of the spectrum of glioma subtypes and grades (11, 13, 16, 17); and (c) group stratification using a priori knowledge (traditional histopathologic classification) or use of subjective thresholds when using objective clinical features as the class defining factor (i.e., survival; refs. 2, 3, 11, 13-17).
We believe that further refinement of a rational biology-based classification scheme would optimally be constructed by using the molecular data, obtained from a large number of samples with different histopathologic subtypes to define the classes according to underlying cellular processes without any preconceived investigator biases. To this end, we have now identified a novel group of glioma subtypes based on the most current, full-coverage gene expression profiles available without any a priori class prediction or bias. We have additionally determined classifier sets from the consensus subtypes and extensively validated them in three large independent test data sets (13) to assess the potential for clinical application of these classifiers in a biologically meaningful glioma classification system.
All tumor specimens used in this study were obtained from patients undergoing surgical treatment in several different institutions and hospitals following written consent in accordance with the appropriate clinical protocols (NABTC 01-07).
In total, 471 specimens were profiled via the HG-U133 Plus 2.0 array. The samples were provided as snap frozen sections of areas immediately adjacent to the region used for the histopathologic diagnosis. Each tumor was assigned a WHO glioma histopathologic subtype by the neuropathologist at the tissue contributing institution, and then the pathology slides were rereviewed by two NIH neuropathologists independently (C.O., M.Q.) who were blinded to the original designated diagnosis. The train set used for class discovery and classifier definition consisted of 159 samples. The results obtained were validated in a nonoverlapping, independent internal test set containing 189 samples, including all WHO histologic types and grades, and further validated in two independent external data sets containing 76 published high-grade gliomas and 265 glioblastoma samples taken from The Cancer Genome Atlas (TCGA) data base, respectively (see details in Supplementary Materials and Methods).
Approximately 50 to 80 mg of tissue from each tumor was used to extract total RNA using the Trizol reagent (Invitrogen) following the manufacturer’s instructions. The quality of RNA obtained was verified with the Bioanalyzer System (ref. 18; Agilent Technologies) using the RNA Pico Chips. Six micrograms of RNA were processed for hybridization on the Genechip Human Genome U133 Plus 2.0 Expression arrays (ref. 19; Affymetrix, Inc.). After hybridization, the chips were processed using Fluidics Station 450, High-Resolution Microarray Scanner 3000, and GCOS Workstation Version 1.3. Detailed procedures on the array hybridization and file processing can be found in the supplementary data.
Using specifically designed filters, we generated six independent probeset subsets for glioma stratification containing different numbers of gene probesets according to the variables described in the supplementary data. Two unsupervised machine learning algorithms, k-mean clustering and nonnegative matrix factorization (NMF), were then applied to each of the probeset subsets separately to identify the underlying classes following the workflow illustrated in Supplementary Fig. S1. The k-mean partition was done using Partek version 6.3 and NMF, using Matlab functions implemented by Brunet and colleagues (20). For the purpose of classifier identification and class prediction, six prediction data sets with 15,553 unique genes were initially created by retaining the probeset with the maximum signal intensity. Prediction analysis of microarray (PAM), a supervised machine learning method (21), was then applied to each prediction set to identify classifiers for the glioma subtypes (Supplementary Fig. S1). A detailed explanation of the methodology can be found in Supplementary Materials and Methods.
Kaplan-Meier survival analysis was used to estimate the survival distributions, and the log-rank test was used to assess the statistical significance between stratified survival groups using Prism 4.0.
Gene set enrichment analysis (GSEA) was used to identify up-regulated expression signatures associated with main types and subtypes (22). All 15 pairwise comparisons of the six subtypes were performed in our GSEA analysis by mapping all 1,687 c2 curated gene sets in MsigDB v 2.0 to the ranked gene expression profiles. The enrichment scores were calculated by walking down the ordered list, and the statistical significance of nominal P values of the enrichment scores was estimated using Kolmogorov-Smirnov statistics by constructing a cumulative null distribution with 1,000 permutations (22).
Two unsupervised machine learning algorithms, k-mean clustering and NMF, were applied to six independent probeset subsets of the 159 glioma specimens profiled via HG-U133 Plus 2.0 array (Supplementary Fig. S1; Supplementary Tables S1 and S2). After the model selection, samples were separated into two main consensus classes (Fig. 1A): one with 69 samples was designated as the O main type (containing the majority of oligodendroglial tumors) and the other, containing 67 samples, as the G main type (enriched with grade IV glioblastoma tumors) in line with their dominant sample composition (Table 1). These two main types were further stratified by applying the algorithms to each class separately. Hence, the O main type separated into two subtypes, designated as OA and OB (Fig. 1B), whereas the G main type separated into two, three, and four subtypes, respectively (Supplementary Fig. S2A). We found that the tumor samples in the four subtypes of the G main type were hierarchically nested under the two upper level subtypes. Thus, the upper level subtypes were designated as GA and GB and the lower level subtypes nested within them were designated as GA1/GA2 and GB1/GB2 (Supplementary Fig. S2B), respectively. Samples with low percentage agreement scores (<25%) were removed for clarity (Supplementary Fig. S3); however, when assessed using the resulting classifiers, the removed samples were assigned to the types and the subtypes that the class discovery algorithms derived (data not shown).
Survival analysis of the two main types shows that the patients in the O main type (median survival, 1,257 d) survive significantly longer than those in the G main type (median survival, 348 d, P < 0.0001; Fig. 1C). Similarly, patients in the OA subtype have a significantly longer survival than those in the OB subtype, suggesting that our unsupervised classification scheme indeed stratify tumors with similar clinical properties (Fig. 1D). Furthermore, there is a strong intragroup association between the histopathology and grade of the tumors within each of the major types and subtypes. The O main type is significantly enriched for a mixture of WHO grade II or grade III tumors (59 of 67), whereas the G main type is predominantly composed of grade IV (GBM) tumors (55 of 69; Table 1 and Supplementary Table S3). In agreement with the known association between patient age and histologic grade, patients with tumors that fall in the O main type are significantly younger than those in the G main type (median age, 43 versus 57 years, respectively; P < 0.0001; Table 1).
To apply this objective classification scheme to clinical applications, we proceeded to identify classifiers that can reliably assign an unknown tumor to the defined subtypes. A supervised machine learning method, PAM, was applied to all the prediction data sets resulting in six sets of independent classifiers consisting of 33 to 352 unique genes. These classifiers produced >92% prediction accuracy when assigning samples to the matching subtypes in a 10-fold cross-validation (Supplementary Table S4). All sets of the classifiers were differentially expressed between the pairwise subtypes across the train set and the validation data sets. One set of these classifier sets for predicting the O and the G main types is shown as an example in Table 2; complete classifier sets can be found in Supplementary Table S5.
To ensure that the classification scheme derived from our unsupervised machine learning system are truly representative of the expression signatures in the different tumor classes, we validated our scheme using an independent, nonoverlapping data set. In this analysis, a test set containing 189 specimens (Supplementary Tables S1 and S2) was stratified into two main types using the same methodology as in the train set (Supplementary Fig. S1), leading to one main type containing 78 samples and the other 87 samples (Fig. 2A, middle). We were able to further separate the smaller main type into two and three consensus subtypes (Fig. 2B, middle) and the larger one with 87 samples into two, three, and four consensus subtypes (Fig. 2C, middle). In agreement with the outcome in the train set, a hierarchically nested relationship was detected in the test set (Supplementary Fig. S3E). The reproducibility of the hierarchically nested subtypes in our test set attests to the reliability of our classification scheme derived from the train set.
To validate the classifier sets identified from the train set (Fig. 2A–C, left), we applied them to the recently stratified classes in test set using PCA analysis and these classifiers from the train set robustly project the intrinsic data variations of the matching subtypes in the test set (Fig. 2A–C, right). The centroids (fold changes of a given set of classifiers between pairwise subtypes) can be viewed as a prototypical expression pattern of the matching subtypes; the tumor subtypes were then identified according to the resemblance of the classifier expression to the centroids (13). Moreover, the classifier validation using a hierarchical clustering analysis on the stratified classes from the test set indicates clear separations along the subtypes (data not shown).
To ensure that these results are a representation of underlying biological processes and not just a statistical artifact of large data sets, we performed a mathematical validation by classifying six random data sets that were carefully designed with the same dimensions and equivalent scale of variations as our classification probeset subsets using NMF. Our results indicate that none of the six random data sets converged to consensus matrices as did the train set or the test set in our classification (data not shown) and their cophenetic correlation coefficients were much lower (0.45–0.73) than our real probeset subsets (0.90–1).
To determine if the classification scheme is independent of the data set used and to determine if the scheme correlates with one of the most comprehensive malignant glioma classification schemes published to date, we used our classifiers to stratify all 76 glioma tumors (GSE4271 data set) from the recent article by Phillips and colleagues (ref. 13; Supplementary Fig. S4). There were several significant differences between our data and the Phillips data. First, the Phillips data was obtained using the HG-U133AB platform, which contains less probesets and genes than the HG-U133 Plus 2.0 array we used. Second, the Phillips data set was comprised only of astrocytic tumors with a high-grade histology. Surprisingly, despite these differences, our classifiers separated the tumors in GSE4271 data set into the six subclasses. Fifty-three of our 54 classifiers for the O-G main types are available on the HG-U133AB array. Hierarchical clustering of all specimens in the GSE4271 data set using these 53 classifiers revealed two dominant tumor clusters, GSE4271-O and GSE4271-G (Fig. 3A). GSE4271-O and GSE4271-G were further assigned according with the centroids of the classifiers (Supplementary Fig. S5). Inspection of the samples in the two main types indicates that the GSE4271-O type is solely composed of tumors that fall into a class that Phillips and colleagues have called “proneural” secondary to the fact that it is characterized by a signature that is enriched for genes involved in neurogenesis (13). By contrast, the GSE4271-G type consists of samples from Phillips’ clusters designated as either “proliferative” or “mesenchymal,” as well as a small number of proneural tumors (Fig. 3A). Further stratification of the GSE4271-O type resolved two distinct subtypes and further stratification of the GSE4271-G type resolved four subtypes using HC (Supplementary Fig. S4A). After the identities of all the subtypes were assigned in terms of the centroids of the classifiers (Supplementary Fig. S5), we found that medium survival of GSE4271-O type (244 weeks) was significantly longer than GSE4271-G type (70 weeks, P < 0.0008; Fig. 3C, left). Possibly more impressively, this unsupervised classification system allowed us to resolve Phillip’s best single prognosis group (proneural) into two distinct subtypes, with the median survival of GSE4271-OA cases (445 weeks) being significantly longer than that in GSE4271-OB (203 weeks, P < 0.016; Fig. 3C, center). No separation in survival was detected among the four subtypes of GSE4271-G (data not shown), consistent with what we found in our own data sets.
As a final test of our classification scheme, we classified the expression data of the GBM collection generated by TCGA project (23). Two major differences exist between our train/test data sets and TCGA data set. First, TCGA data set is profiled on HT_HG-U133A (22,268 probesets) rather than on the HG-U133 Plus 2.0 array (54,000 probesets) we used resulting in a smaller gene representation in the TCGA data sets compared with ours (18,400 versus 47,000 transcripts). Second, all of the tumors profiled in TCGA data set are glioblastomas, whereas our original train set and test set represent a wide spectrum of glioma types and grades. Even so, we were able to classify the 265 GBM samples from TCGA data set using our classifiers, stratifying all TCGA GBM samples into two main classes: TCGA-O and TCGA-G. As we found with our other data sets, the TCGA-O and TCGA-G groups could be further separated into two and four subclasses, respectively (Fig. 3B and Supplementary Fig. S4B). The prototypical expression patterns of all the corresponding classes following their hierarchical nested relationship are identical to those subtypes from all the other data sets we used (Supplementary Fig. S5). Significantly, we identified a small group of GBM samples (37 patients) that grouped into the oligo-enriched group. Importantly, Kaplan-Meier analysis of these samples indicated that these patients had significantly prolonged survival compared with the other GBM patients in the TCGA data set (median survival, 561 versus 327 days, P < 0.0003; Fig. 3C, right).
GSEA detects cellular pathways and/or mechanism that are differentially expressed between two phenotypes in functional units of genes (22). To begin to elucidate some of the basic underlying biological differences between the subtypes of gliomas identified through our analyses, we performed pairwise comparisons between types and subtypes that were located at the same level in the nested hierarchical structure (Supplementary Table S6). A summary of the results obtained is presented in Fig. 4. When comparing the two main types, we found that the O main type exhibited enhanced activities of exogenous hormone stimulated growth and PAR1 signaling activity whereas up-regulated profiles associated with the G main type included a large number of cell cycle/mitotic pathways, hypoxia, tumor necrosis, and nuclear factor-κB (NF-κB) pathway signaling (Fig. 4; Supplementary Table S7), consistent with a more malignant and aggressive tumor type; a supposition supported by the poorer survival of the G main type.
In the case of the OA/OB comparison, we found that three gene sets (ARF pathway, TEL pathway, and programmed cell death signaling) were up-regulated signatures in the OA subtype (Fig. 4). These pathways are suggestive of tumors with relatively intact cell cycle regulation, maintenance of chromosome ends, and apoptosis, consistent with a less malignant phenotype, as supported by the longer survival of patients in this subtype (Fig. 4). By contrast, the OB subtype is notable for genes involved in plateletderived growth factor (PDGF) and epidermal growth factor (EGF) signaling, a common feature of more aggressive tumors (Supplementary Table S8).
The expression profiles in GA subtype are dominated by pathways involved in mitosis as indicated by 13 up-regulated gene sets in cell cycle regulation, DNA replication, and proliferation (Supplementary Table S9), making the GA subtype one of high proliferative potential. Interestingly, the proliferative subtype of the GSE4271 data set maps overwhelmingly to the GA subtype, indicating concordance between the two different schemes. Meanwhile, the GB subtype is notable for overexpressed signatures of the NF-κB pathway, tumor necrosis factor pathway, transforming growth factor-β signaling pathway, interleukin (IL-1, IL-10, IL-22, IL-6) signaling pathways, and up-regulated vascular endothelial growth factor signatures (Supplementary Table S9), all consistent with a mesenchymal phenotype, in agreement with the observation that the majority of the tumors designated as mesenchymal in the GSE4271 data set map to the GB subtype. Further comparisons of the functional annotations between the GA1/GA2 and GB1/GB2 subtypes are detailed in Supplementary Tables S10 and S11.
The need for an objective, biological-based classification of gliomas is exemplified by the high rate of divergent diagnoses, the inexact prognostic capabilities, and poor therapeutic predictive properties of the current histopathologic classification schemes (10). Here, we report the development of an unbiased, gene expression–based, histology-independent glioma classification system by applying unsupervised machine learning algorithms to the gene expression profiles using the largest collection of gliomas published to-date, covering the main WHO histologic types and grades. Our analyses show that gliomas can be separated into two main types, O and G. The O main type is further divided into OA and OB subtypes, and the G can be divided into four subtypes, showing hierarchically nested relationships designated as GA1, GA2, GB1, and GB2 (Supplementary Fig. S2B).
Our classification system differs from previously reported ones in several important aspects. First, the classes were discovered based on the expression profiles of all glioma histologic types and grades without any a priori exclusion of the samples based on some clinical variable. Second, we did not select the genes to be analyzed from a manually curated list of genes thought to be relevant but rather to include all informative genes to represent the entire transcriptomic profiles of individual glioma patients. Third, we selected the subtype-defining classifiers based solely on their statistical ability to separate the classes rather than our preconceived notion of their functionality. That is, we elected to allow the computational analyses determine the appropriate classification scheme based totally on the biology of a wide spectrum of gliomas, as manifested by their genomic-scale transcriptomic profiles, in the hope that the biology would translate to a clinically relevant classification. By respecting the concordant molecular signature behavior of the tumors in the analyses, we have achieved a truly unsupervised solution. The difference in how this analysis was performed compared with prior studies is not trivial because the classification scheme we have constructed is based purely on observed gene expression profiles, without predefined genetic, pathologic, or clinical assumptions, and thereby represents a glioma classification scheme based solely on unbiased biological data. The power and reliability of the classification system described here rest not only on the purely statistical and objective way in which the analyses were performed but also on the fact that all six glioma subtypes and their corresponding classifiers were derived from a train set, confirmed and validated using an in-house and totally independently generated test set larger than the train set, and further validated using two independent, externally generated data sets.
Despite the use for the first time of a purely unsurpervised methodology, our results are generally consistent with and extend the findings from previously published smaller glioma classification studies (11, 13, 24). For example, our O main type is very similar to GSE4271 proneural class and our G main type extensively overlaps with the GSE4271 mesenchymal and GSE4271 proliferative classes. More specifically, the GSE4271 proliferative class largely overlaps with our GA subtype, both being represented by up-regulated cell proliferation–associated genes, whereas the mesenchymal class largely overlaps with our GB subtype, both characterized by invasive and mesenchymal tissue–associated genes. Our subtypes, however, further refine this classification system as exemplified by the fact that the OA and OB classification separates the GSE4271 proneural tumors into two different subtypes with significantly different overall survivals. Furthermore, our analyses show that the GSE4271 proliferative and GSE4271 mesenchymal classes can be broken into four different subtypes based on the differences in their classifier signatures.
It should be noted that the lack of survival difference between the most aggressive tumors that fell into the GBM-rich subtypes should not be seen as a weakness in or failure of our classification system. In contrast to diseases, such as breast cancer and lymphoma, that tend to have much greater variability in the natural history of the disease (i.e., survival), the vast majority of glioblastoma patients have a rather uniform and poor survival. Our classification system is consistent with those by Phillips and others and with data from large clinical trials conducted over the last two decades that consistently show a small subtype of GBM patients with better-than-expected survival (as represented by the GBMs in our O main type) but with a majority of GBM patients having uniformly short survival. Thus, although survival can be a correlation of biology, biology is not necessarily a correlation of survival and is the reason that we intentionally chose to perform an unsupervised analysis that was not based on survival, as has been done so often in the past. Certainly, patients with aggressive tumors with vastly different biologies can have roughly equivalent survivals (i.e., GBM, pancreatic cancer, metastatic lung cancer). In the case of tumors within the brain, survival may, in particular, be more a surrogate of tumor growth rate and lack of therapeutic responsiveness than biology given the physiologic constraints of a mass growing within the closed compartment of the cranium and the vital and sensitive nature of the underlying brain tissue. This is likely why nearly all malignant tumors within the brain, whether they be primary tumors (gliomas, medulloblastomas) or secondary tumors (metastatic lung cancer, melanoma), tend to have a uniformly poor and homogenous survival (6–12 months) if the specific tumor is not inherently sensitive to therapy. Nevertheless, the accumulating molecular data clearly shows that, despite survival homogeneity, GBMs are heterogeneous at the molecular and genetic level (11, 16, 25). Our classification scheme defines and categorizes this heterogeneous biology that will be crucial for designing clinical trials of molecularly target agents that are enriched for patients most likely to respond to therapy and ultimately for the practice of patient-specific or “personalized medicine.” Thus, the clinical utility of these GBM-rich subtypes will only be realized with the acquisition and the analyses of much more corollary clinical data.
We entered into this project with the bias that the standard WHO classification system would be very poorly representative of the underlying tumor biology. For the most part, we found this to be true, for although the WHO classification of GBM versus non-GBM pathology was largely upheld in our two major expression groups (96% of all GBMs in GBM-rich group versus 75% of non-GBMs in the oligo-rich group), tumors designated as WHO grade 2 or grade 3 astrocytomas, oligodendrogliomas, or mixed gliomas randomly distributed between the expression subgroups whether that designation was based on the original pathologic diagnosis (home institution) or by our central pathology review (data not shown). The fact that there was nearly a 30% to 40% discrepancy in the designation of grade and glioma subtype between the original and central pathology review of non-GBM tumors (data not shown) reinforces the subjective nature of the currently used classification system and testifies to the potential problem of using the WHO system to group specific patients into particular biological strata. By contrast, our classification system, derived and based purely on computational analyses of gene expression, consistently groups tumors into one of six groups across four very divergent data sets constituting over 700 gliomas.
The classifier gene sets described here showed >92% prediction accuracy in a 10-fold cross-validation (Supplementary Table S4). Additionally, the classifiers were successfully used to predict and to assign the derived subtypes in a large independent test set (Fig. 2) and further to stratify two external data sets into six subtypes (Fig. 3 and Supplementary Fig. S4). The robustness of our gene classifiers suggested a potential for a useful clinical tool for glioma diagnosis once more extensive validation is undertaken. Indeed, our classification system allowed us to identify a subgroup of GBMs in the TCGA data base that fell into the oligo-rich group and had markedly better survival than the remainder of the group, demonstrating the potential power of this biologically based classification system.
Finally, our GSEA analysis points to potential functional properties of the different subtypes identified in this analysis. For example, the survival advantages of OA subtype might in part be explained by an intact p53 regulatory pathway as represented by the activity of the ARF pathway and by a tendency toward genomic stability through the maintenance of chromosome ends, as suggested by activity of the Tel pathway. In contrast, the protooncogene signaling of the PDGF pathway (known to be aberrantly regulated in a significant percentage of astrocytomas) may confer higher proliferative properties to tumors in the OB subtype. Needless to say, the true functional significance of pathway activation within tumors in each of these subtypes remains to be elucidated through biological studies. Nevertheless, this analysis, as well as others like it, begins to build a new framework by which basic and clinical scientists can investigate the biological, functional, and clinical significance of these novel molecular classes with the hope of ultimately deriving a tumor classification system that will have both biological and therapeutic predictive value.
Grant support: Intramural Research Program of NIH National Cancer Institute Center for Cancer Research.
J.C. Zenklusen and H.A. Fine contributed equally.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.