Study Design

Participants are drawn from the Longitudinal Youth-at-Risk Study (LYRIKS), a prospective observational study on youths susceptible to psychosis (Lee et al., 2013). LYRIKS participants are drawn from mental health care and community-based services and from various educational institutions in Singapore. Eligibility requirements mandate a narrow age range (14–29 years), no existing antipsychotic treatment, no other psychotic disorder or neurological disease, and no illicit substance use history. UHR status is ascertained using CAARMS, and positive identification is made if subjects exhibit attenuated psychotic symptoms (APS), brief limited intermittent psychotic symptoms (BLIPS), or general vulnerability (Yung et al., 2002). Healthy controls in this study are participants who did not fulfill UHR criteria and had no psychiatric disorders when evaluated on the Structured Clinical Interview for DSM–IV Axis I Disorders (First, Spitzer, Gibbon, & Williams, 1996). Details on LYRIKS’s study methodology are obtainable from prior publications (Lee et al., 2013; Mitter, Nah, Bong, Lee, & Chong, 2014). The dataset comprises 55 UHR subjects and 28 healthy controls (Table 1). Gender and ethnic ratios are balanced between subjects (UHRs) and healthy controls (non-UHRs) to minimize factor proportion-imbalance disparities (Patil, Bachant-Winner, Haibe-Kains, & Leek, 2015).

Gene Expression Measurement

Peripheral blood is drawn immediately following assessment into a Tempus Blood RNA tube (Applied Biosystems, Foster City, CA) and is stored at − 80°C until RNA extraction. Total blood RNA is extracted using a Tempus Spin RNA isolation kit (Applied Biosystems) and amplified using an Illumina TotalPrep RNA amplification kit (Ambion, Austin, TX). mRNA expression profiles are assessed on Illumina HumanHT-12 v4 Expression BeadChip arrays. Experimental quality controls are performed prior to RNA amplification and before RNA hybridization to ensure that RNA concentrations are ≥100 ng/ml, A₂₆₀/A₂₈₀≥2.0, and clearly defined ribosomal peaks were seen on agarose gel. All procedures adhered to design protocols as per manufacturer recommendations. The readings of the beads in the array are analyzed by Illumina iScan following hybridization.

Microarray Quality Control

Probe intensities on the Illumina Microarray BeadChips are summarized using GenomeStudio (Illumina, San Diego, CA). To ensure read quality, sample-based and probe-based quality control (QC) is performed. In sample QC, samples with signal-to-noise ratio below 10 or samples with intensity of negative control probes higher than control probes are considered unreliable and discarded. In probe QC, probes with detection p values below 0.05 or below background are considered undetectable and discarded. A total of 84 subjects (56 UHRs + 28 non-UHRs) across 18,029 gene probes meets QC standards.

Normalization Methods

Following background subtraction, probe-to-gene name mapping, and following median centering of expression, four normalization approaches are applied on the quantified data: Log-conversion (None), quantile normalization (Quantile), gene fuzzy scoring (GFS; Belorkar & Wong, 2016), and surrogate variable analysis (SVA; Leek & Storey, 2007).

Strictly speaking, None is not true normalization; we use this primarily as a negative control and to facilitate parametric statistics. We use the natural log, such that for each observed measurement x, the transformed value is log_ex, where e = 2.718. None reduces the effect of high orders of magnitude (because e ∼ 3, the data range is reduced by three orders) and improves the distribution symmetry. However, None is usually insufficient to render samples cross comparable and is usually accompanied by a second data transformation, typically z scaling, linear interpolation, or quantile normalization.

Quantile is a gold-standard data transformation procedure and reportedly performs well in stabilizing variance in Illumina BeadChip HT-12 (Schmid et al., 2010). Although we could include linear interpolation and z scaling, these are already known to be inferior to Quantile anyway and detract from more interesting new methods that explicitly deal with heterogeneity. GFS is an unsupervised signal-boosting transformation (Belorkar & Wong, 2016; Wang et al., 2016). In GFS, the log-transformed expression matrix is transformed by weighting individual genes per sample based on expression ranks. GFS uses two thresholds, θ₁ and θ₂. Features with ranks above θ₁ are assigned a weight of 1, features with ranks between θ₁ and θ₂ are assigned an interpolated weight between 1 and 0, and features with ranks below θ₂ are weighted as 0. Let r(g_i,p_j) be the rank of a biological feature g_i in patient p_j and q(p_j,θ) be the rank corresponding to the upper θth level of feature ranks in p_j. The GFS score s(g_i,p_j) assigned to feature g_i for patient p_j is determined by the function

As arbitrarily defined thresholds, θ₁ and θ₂ can take on different values, for example, the default settings in GFS set θ₁ to 5% and θ₂ to 15%. However, in their evaluations, Belorkar and Wong (2016) stated that varying θ₁ between 5% and 10% and θ₂ between 15% and 20% produces similar results. Comparing GFS against standard normalization techniques—for example, mean-scaling, z, and quantile normalization—GFS consistently gives better class discrimination, is robust against batch effects, provides improved power even when sampling at small sample sizes, and facilitates reproducible selection of biologically relevant features.

SVA is applied on the log-transformed expression matrix. In contrast to GFS, SVA is a supervised method—that is, it requires specification of class labels, UHR subject and UHR control, a priori, meaning that the algorithm is aware of class information—and is therefore autobiased toward class effects (Leek & Storey, 2007). It assumes that consistent sources of variation nonassociated with the class factor are likely associated with some unknown heterogeneity factor. Having identified the part of the data not associated with class variation, it then estimates heterogeneity factors via singular value decomposition. These irrelevant variations (expressed as surrogate variables) are then removed via regression. The beauty of SVA is that we may isolate UHR-specific signals against a backdrop of various confounding factors that we cannot possibly control for, especially genetic factors. Although it is extremely powerful, SVA can be complex to use (Jaffe et al., 2015).

Principal Component Analysis (PCA)

PCA is a linear summarization technique that collapses high-dimensional data (e.g., thousands of genes) into a lower number of orthogonal dimensions, or principal components (PCs; Giuliani, 2017). We use PCA to observe whether individual samples group by class consistently, without feature selection and regardless of normalization method.

We also use the PCA-generated loading matrix to determine the association of individual PCs against each factor (class, gender, and ethnicity). Assuming that a factor may have more than two levels, association is evaluated using the nonparametric Kruskal–Wallis (KW) test (p ≤ 0.05).

Feature Selection and Multiple-Test Correction

The F test is used to evaluate a differential expression for each gene, following Benjamini–Hochberg (BH) multiple-test correction.

The F statistic is expressed as

Classifier Training and Cross-Validation

Functional Analysis Based on Networks

GeneMANIA (http://www.genemania.org/) is used to evaluate functional relationships among differential genes (Montojo, Zuberi, Rodriguez, Bader, & Morris, 2014). It features an adaptive weighting method to combine various network data, such that after the networks are combined, the differential genes interact maximally with each other, while interacting as little as possible with genes not in the list. This allows GeneMANIA to learn which networks mediate the underlying functional relationship among the differential genes but also permits the inclusion of additional genes strongly associated with the differential set.

Enriched gene ontology annotations associated with the finalized network are also reported, and we used these for checking whether the associated network, induced by the set of differential genes, has strong neurological associations.

Ultra-high Risk (UHR) Blood-Derived Gene Expressions Are Strongly Distinct Between Data Normalization Techniques

Previously, it was unknown whether UHR subjects possessed distinct but conserved gene expression profiles. Plotting individual samples using the first 3 PCs without prior feature selection suggests that they do (Figure 1A).

The existence of strong class-distinguishing signals without a priori feature selection is good evidence that UHRs form a distinct class. However, this may be a consequence of data-processing bias. Therefore we rechecked using various data-processing and normalization techniques: None, Quantile, GFS, and SVA. The first two techniques are common, GFS is a signal-boosting transformation, and SVA deals with unwanted heterogeneity.

While the scatterplots present a consistent view that UHR subjects are distinct (Figure 1A), inter- and intracluster variances differ. In None, interclass variability is lower relative to intraclass variability. Quantile, GFS, and SVA all produce significant improvement, demonstrating the importance of normalization. Note that good class separation is expected in SVA anyway.

PCs are high-dimensional projections based on the expression patterns across thousands of genes such that each holds information regarding a conserved pattern of expressional change correlated with various factors (e.g., class, age, gender, ethnicity, and even biological pathways; Giuliani, 2017; Goh & Wong, 2016c). Each gene-based variable loads differently onto each PC, and its loading can be determined from the PCA loading matrix.

We may also investigate PCs for evidence of bias (Figure 1B): A disproportionate proportion of variance in PC1 suggests technical bias, for example, batch effect (Giuliani, 2017; Goh & Wong, 2016c) or spurious correlations (Giuliani, 2017; Goh & Wong, 2016c). Expectedly, None holds a high proportion of total variance in PC1, while total variance is distributed more evenly in the other normalization methods.

The intersample distances change as more PCs are considered, and there is no reason why subsequent PCs (after PC3) should be ignored. Each PC uniquely encapsulates a distinct portion of total variance and is mutually orthogonal. We may deploy the PCs more effectively by quantifying the association of each PC against known factors (class, gender, ethnicity). If sample class is associated significantly among the top PCs, then this supports the scatterplot analysis. Conversely, we may also reverse engineer interesting PCs and identify which genes load strongly onto it (Goh & Wong, 2016c).

To test association, the KW test is used, accounting for scenarios where a factor has more than two levels. If a PC exhibits strong differential behavior between factor levels (p ≤ 0.05), then we assume that this factor can explain this PC. As a rule, for good normalization methods, the top PCs should be significantly associated with class, followed by other important factors.

Table 2 is consistent with the scatterplots. To improve visual impact, the KW p values are summarized to two decimal places, and boldface indicates a significance level below 0.05. It should be noted that in hypothesis-based statistical testing, the magnitude of the p value is actually inconsequential; of importance is only whether it falls below the threshold (Goodman, 1992). Although all PCs can be tested, the top 10 are adequate to make our point regarding factor rankings.

Class effect ranks highly. Interestingly, in Quantile, the top two PCs correlate strongly with gender and ethnicity, while class is relegated to lower PCs. This does not mean that a class effect is absent in the first two PCs; rather, these PCs exhibit a stronger association with other factors. GFS brings class, gender, and ethnicity to the fore (PCs 1–3). Note that GFS is unsupervised and robust against heterogeneity and technical bias (Belorkar & Wong, 2016; Wang et al., 2016), so this is strong evidence that UHRs are genotypically distinct.

In SVA, various factors besides class can be considered. Because we only supplied class, ethnicity and gender effects are suppressed. It is known that psychosis is subclassifiable via these factors (Bresnahan et al., 2007; Canuso & Pandina, 2007). Our test results agree with expectations: Class association ranks highly in SVA (PC2), while gender and ethnicity are relegated to PCs 8 and 9, respectively.

In None, Quantile, and GFS, UHR subjects and controls are consistently separated (Figure 1), although intra- and intercluster distances vary. PC-based factor association suggests that class ranks highly. The unsupervised approach of GFS greatly improves signal-to-noise ratio, and relevant factors rank highly among the top PCs following GFS transformation.

Statistical Feature Selection Depends Strongly on the Data Normalization Technique

We expect statistical feature selection to provide disparate information regarding differential genes (Figure 1A). Biological relevance is inferred from statistical testing (differential expression), but there are many caveats. While different statistical feature-selection methods can yield different results (Christin et al., 2013; Langley & Mayr, 2015), normalization is also important—and less well explored.

If normalization has limited impact on feature selection, then given a single feature-selection method, selected features from differently normalized datasets should be similar. Genes are selected using the F test/BH correction (p ≤ 0.01). The F test is used to ensure cross-comparability of the other methods with SVA as it reduces overall variability, while the degrees of freedom, given any statistical test, remain the same. While this amplifies effect size (such that if a gene is truly differential, it will be more easily detected), it also increases the false-positive rate (a nondifferential gene is also likely to be reported as significant). Thus SVA requires an appropriate correction (provided within the SVA package on the F test).

Based on p value distributions (extreme right skew), it is clear that None is problematic (Figure 2A). A standard cutoff at 0.05 introduces many statistically significant features. After raising the cutoff to 0.01, ∼6,000 genes are still declared differential (likely most are false positives). Note that the issue is unresolvable via p value ranking (and selecting the top n genes), because rank instability increases as p values get smaller (Wang et al., 2016).

The other methods’ p value distributions are more reasonable, suggesting that only a few genes are differential. Quantile completely reverses the p value distribution (Figure 2A). GFS is extremely stringent, while SVA predicts slightly more features than Quantile.

None, Quantile, GFS, and SVA select a total of 5,877, 256, 5, and 556 significant genes, respectively (p ≤ 0.01; Figure 2B). Among these, only one gene (MAGEB16) is shared. Quantile, GFS, and SVA exhibit deeper overlaps with each other. In particular, GFS shares 80% (four out of five) of its selected features with Quantile and SVA, but not with None. And these overlapping genes are important: Their associated p values (based on SVA) are far lower (and therefore more significant) than nonoverlapping genes (Figure 2C). Thus we conclude that normalization has strong downstream implications for feature selection.

We are concerned with which approach produces cleaner data. If noise were eradicated, then even without feature selection, the best approach would produce highly accurate models anyway. To do this, we turn to cross-validation without feature selection. Both GFS and SVA reduce heterogeneity, particularly for GFS where the median cross-validation accuracy is highest at 88% (Figure 2D). This confirms that GFS and SVA remove unwanted confounding variation, unlike None and Quantile. Therefore we conclude that GFS is the best normalization technique and is more likely to provide a useful signature.

Gene Fuzzy Scoring Gene Signature Is Functionally Relevant and Beats Randomly Generated Signatures

GFS selected five genes (IGSF1, LOC653712, LRRTM2, MAGEB16, and PCSK1; p ≤ 0.01). The stringent cutoff is used for limiting the high number of selected features due to log-conversion. For phenotypic correlation, we may bolster sensitivity by relaxing the p value cutoff to 0.05 (additional genes are CDH11, CXORF55, CYR61, NOVA2, NTRK2, PARVA, and RBMY1J).

We use the term signature liberally: We do not mean it in the sense that it has valid diagnostic capabilities in the greater UHR population but rather use it to indicate a set of genes being evaluated for class relevance and predictive power, within the confines of our study design. Although small, the GFS signature can distinguish UHR subjects and controls (hierarchical clustering; Euclidean distance; Ward’s linkage; Figure 3A); that is, we may replicate the strong PC-based interclass segregation (Figure 1A) with only 12 genes.

To determine functional interrelationships, we supplied the GFS signature to GeneMANIA (Montojo et al., 2014; Vlasblom et al., 2015) and searched for functional links based on coexpression, co-localization, and shared functional domains (Figure 3B). We also recovered additional implicated genes (in gray) strongly associated with the signature. The induced gene network (including both signature and implicated genes) is strongly associated with neurological functions, including genes associated with behavior, signal release, and neuron projection.

GFS provides the highest cross-validation accuracies (Figure 2D), suggesting that it removed a high amount of noise. Not all genes are useful though; we may simplify model building and explanation if we isolate only differential genes (via statistical feature selection) to form signatures. Unfortunately, this is not straightforward. Venet et al. (2011) demonstrated that inferred signatures across breast cancer studies seldom generalize and do not outperform random signatures. In fact, the larger the gene signature is, the more likely it is that randomized gene sets will outperform it. Although the breast cancer signatures in themselves are predictive, they offer no more information than any random assortment of genes, meaning that they cannot be used to reveal mechanism or biological insight. This occurs because a large fraction of genes are autocorrelated with cancer but play no role in disease progression. It is unclear if this is also true for UHR.

We repeated cross-validation with feature selection using the F test/BH correction (p ≤ 0.01). It is interesting that feature selection with GFS reduces classifier accuracy (Figure 3C). This means that although the top five GFS genes have some predictive power, UHRs are quite heterogeneous. Moreover, despite the drop in classifier accuracy, GFS signatures do much better than random signatures (Figure 3C).

Feature selection is important. Certainly the F test failed to preserve the original high accuracy given all features. However, we also removed much noise, such that randomly picked gene signatures do no better. In comparison, feature selection with Quantile, SVA, and None appears to preserve, or improve, classifier accuracy, but the random signatures accompanying these also do well (Goh et al., 2017, Supplementary Figure 1). Inferred signatures derived from these normalization methods are likely less meaningful.

It is worthwhile to evaluate feature-selection stability during cross-validation. Goh et al. (2017, Supplementary Figure 2A; compare Figure 2D) suggested that None/SVA/Quantile-normalized data still contain some noise and thus that feature selection helps a little, whereas the GFS-normalized data have less noise, so feature selection does not help as much. Moreover, the feature-selection stability across all normalization is poor, although GFS is still superior in this respect (Goh et al., 2017, Supplementary Figure 2B). Moreover, the gene reproducibility index (see Goh et al., 2017, Supplementary Methods) for GFS is the highest and approximately five times more reproducible than SVA (Goh et al., 2017, Supplementary Figure 2B).

Limitations of the Current Study and Follow-Up Investigations

The present study demonstrates that UHR subjects possesses distinct blood-based signatures, but we do not claim generalizability. Ostensibly, the sample size is limited. However, given careful subject selection to facilitate cross-comparability and minimize background (Qin et al., 2014), UHR subjects can be consistently distinguished from UHR controls, but normalization can present disparate views on the underlying biology.

Good normalization methods should minimize or eliminate technical bias and increase signal-to-noise ratio (i.e., amplify signal, reduce noise). We should not rely on feature-selection methods to work properly without proper data cleaning. As it turns out, most feature-selection methods are not resistant to technical bias or noise (Goh & Wong, 2016c).

Commonly used normalization methods fall short, while newer approaches, designed for dealing with heterogeneity, are superior. The unsupervised approach of GFS is promising. In follow-up studies, we aim to determine if GFS signatures are generalizable in larger cohorts and across a wide variety of independent studies on the same phenotype.

This is the first study to investigate a genetic basis underlying UHRs. We are positive that specific blood-based gene expression signatures from UHR subjects can be elucidated, which can in turn facilitate accurate early detection and therapeutic intervention.

We are also aware that issues pertaining to blood–brain gene expression correlations are not addressed (Cai et al., 2010; Jasinska et al., 2009). We avoid overinterpretation of mechanism but merely point out that the GFS blood signature is phenotypically coherent. In subsequent studies, we aim to relate blood- and brain-derived UHR gene signatures.

Although we focus mostly on the downstream effects of normalization in functional analysis, we also introduce some useful analytical techniques. PCA and normalization are discussed in subsequent sections. Here we comment on our cross-validation setup where we randomly sampled from UHR subject and control, performed a standard feature selection, trained a classifier, and obtained a cross-validation accuracy. This is standard procedure, but it does not tell us whether this accuracy is meaningful in that the selected features are phenotypically relevant in an exclusive sense (Goh & Wong, 2016b). This can be counterchecked by randomly sampling an equal number of features and checking if the derived cross-validation accuracies do better. An empirical p value for the cross-validation accuracy is the number of times a random signature does better than the observed. This is a powerful approach, because it deals specifically with the issue of random signature superiority (Venet et al., 2011): If the signature is relevant, it must be exclusive, and it must perform better than any random assortment of genes. Intuitively, a signature passing this check is more likely to be generalizable.

Finally, when study designs are imbalanced, it is important never to derive statistics for each gene across samples, as this would lead to test-set bias during classifier building (Patil et al., 2015). Test-set bias occurs when the classification of any subject or tissue in a dataset is dependent on the composition of the other subjects or tissues as the expression of each protein or gene is normalized with respect to the entire dataset (Patil et al., 2015). Any gene signatures derived based on any normalization strategy that has dependency across subjects or tissues of the test set will be irreproducible when the population changes composition or size. Thus readers should be careful when encountering any normalization method that adjusts across multiple subjects or tissues rather than information strictly within a subject or tissue itself (Goh & Wong, 2016a). The problem is easily avoided simply by not normalizing across samples but rather within samples, as we have done here. However, the unsolved problem is that we have less information regarding variance in UHR controls, and this can indeed undermine our analysis methods. But as mentioned earlier, we do not claim generalizability but work within the confines of this study.

Normalization Is Vital and Contributes Toward Reproducible Signatures

It is clear that log-conversion (None) alone is insufficient and may produce false effects. The oft-used quantile normalization (Quantile) provides considerable improvement but may not be good enough for practical applications. Finally, the highly sophisticated technical variability removal method, SVA, is powerful, but as a supervised method, it is susceptible to the veracity of the class labels: If wrong, for example, owing to misdiagnosis, then SVA can make mistakes. Finally, we introduce the unsupervised signal-boosting transformation GFS, which provides the highest data-cleaning utility (Figure 2D).

Although GFS’s performance suffers when combined with a feature-selection approach, it is the only method that can strongly outperform random signatures (Figure 3C; Goh et al., 2017, Supplementary Figure 1). This is vital, because when applied on random subsets of data, classifiers built with GFS are “meaningful.” If classifiers built with random genes are accurate as well, this means that many genes in the dataset are noncausally correlated with the phenotype, such that any inferred signature has no real meaning anyway. Although the other approaches generally produce higher classifier accuracies following feature selection, only the GFS-derived classifiers outperform the randomly built classifiers. This provides greater confidence that the five-gene GFS signature inferred from the full dataset is not associated with the UHR phenotype simply by chance.

PCA Can Be Used More Effectively

PCA scatterplots are commonly used for determining the presence of batch effects and checking sample separation before and after feature selection. Typically, analysis is limited to PCs 1–3, with the assumption that most variance is contained therein (the remaining variance being uninteresting). This is not always the case where, following certain normalization approaches, the distribution of variance among PCs becomes more balanced (Figure 1B). Although each subsequent PC contains a lower proportion of total variance, it is not negligible, and there is no valid reason to discard any of them. Conversely, a high proportion of variance in PC1 suggests high technical bias (Goh & Wong, 2016c). It is important to deal with bias, instead of blindly accepting the top three PCs.

We may use PCs for conducting correlation checks with known factors: If technical batches are known, we may determine which PC is associated with batch. If correlated with PC1, batch effect dominates, and skews feature selection. Similarly, examining PC correlations with other factors determines if they contribute strongly and meaningfully to the model. In particular, we expect gender and racial effects to be present. Indeed, these are correlated with top PCs, except in SVA. (This is expected, as we did not supply this to SVA, and gender/racial effects were thus suppressed.)

Detailed PC examination reveals the influence of various factors, but we need to figure out which genes are responsible. This is possible (but we did not deploy this technique explicitly here): Suppose PC1 is strongly associated with batch; we can go back to the PCA loading matrix and isolate genes that load strongly onto PC1, thus isolating the batch effect–susceptible genes (Goh & Wong, 2016c). If PC1 is strongly and uniquely associated with class, then we can extract the “class differential” genes from there. This provides useful information, in addition to statistical feature selection.

Certainly PCA is a very versatile technique, and it can be used more creatively, to powerful effect.

UHR Phenotype Class Is Inferable via Blood-Based Signatures

Although our sample size is moderate, we are able to consistently observe strong differences between UHRs and non-UHRs, even without feature selection. Examination of the gene signature derived from GFS suggests association with known psychosis genes. Many of these genes are implicated in psychosis and schizophrenia: PCSK1 is reported to be dysregulated (Hokama et al., 2014). LRRTM2 is a maternally suppressed gene that is associated paternally with handedness and schizophrenia (Francks et al., 2007). NTRK2 together with BDNF is associated with paranoid schizophrenia (Lin et al., 2013). PARVA is associated with cognitive control, the loss of which is a core trait in schizophrenia (Lewis, Curley, Glausier, & Volk, 2012). Although some genes (e.g., MAGEB16) are reported consistently, they are immune associated and do not seem consistent with the UHR phenotype. However, the immune system and schizophrenia have a long history of intricate association (Horvath & Mirnics, 2014; Jenkins, 2013; Muller & Schwarz, 2010; Vetlugina, Lobacheva, Semke, Nikitina, & Bokhan, 2016). Although it is useful to be able to phenotypically relate gene signatures, doing so does not offer nonquantifiable evidence. Therefore additional means to demonstrate signature specificity and power are necessary.

Taken together, the results suggest that blood-based diagnosis is useful for characterizing UHRs.