Performance of case-control rare copy number variation annotation in classification of autism

ASD sample collection

Samples were collected as part of the Autism Genome Project (AGP), an international consortium with more than fifty contributing sites in North America and Europe; informed consent was obtained from all participants. ASD was diagnosed based on the ADOS (Autism Diagnostic Observation Schedule) and ADI-R (Autism Diagnostic Interview - Revised) [26, 27]. Patients with karyotypic abnormalities, Fragile X syndrome or other genetic syndromes causing congenital malformations were excluded. Only samples from subjects with European ancestry were used for this analysis, as determined by multidimensional scaling analysis (MDS) of known SNPs (single nucleotide polymorphisms).

Unrelated controls were assembled from three large studies in which subjects had no previous psychiatric history: SAGE (Study of Addiction Genetics and Environment), Ontario Colorectal Cancer study, HABC (Health Aging and Body Composition) [28–31]. Controls were genotyped on the same array platforms (Illumina 1M single or duo arrays) as ASD subjects and parents, applying the same quality control procedures and CNV calling algorithms.

Part of the samples (stage 1) was analyzed and published in 2010 [8], whereas the remainder (stage 2) is being published in a separate paper [32].

CNV calling and clinical classification

Only samples meeting quality thresholds were used for CNV analysis. CNVs (of size 30 kb or greater) were detected using an analytical pipeline optimized for Illumina 1M arrays [8, 33]. All de novo CNVs were experimentally validated. Samples with copy number variation greater than 7.5 MB were excluded.

The clinical classification of de novo and inherited CNVs as pathogenic, uncertain or benign was established according to the American College of Medical Genetics guidelines [18], on the basis of genetic loci known to be implicated in ASD (124 genes and 55 loci) [16]; large and very rare CNVs were also classified as pathogenic, even if they were not reported before. Experimental CNV validation details can be found in [33].

CNV gene annotations and gene-set construction

Rare CNVs were mapped to genes whenever at least one transcript overlapped the CNV; transcript coordinates were based on RefSeq hg18 (the same build as for the BeadChip arrays). Similar performance results were obtained using the more stringent exonic mapping (i.e. requiring at least one exon to be overlapped by a CNV).

Curated neurally-relevant gene-sets captured (i) predicted haploinsufficiency [33], (ii) brain expression levels, based on the spatiotemporal Brainspan RNA-seq data-set [35], as well as brain expression specificity compared to other tissues, based on the Novartis Tissue Atlas microarrays [36], (iii) experimentally-determined synaptic complex membership [37] and regulation by FMR1, a key modulator of mRNA translation required for synaptic plasticity, based on two different methods [38, 39], (iv) implication in neurological or neuropsychiatric disease in humans (according to HPO, the Human Phenotype Ontology) [40], (v) implication in abnormal nervous system or abnormal behavior in mice, according to Mammalian Phenotype Ontology (MPO) annotations provided by Mouse Genome Informatics (MGI) [41], and (vi) neuronal or nervous system function based on Gene Ontology (GO) [42] annotations or pathway membership. A more detailed description is provided in the supplementary methods and results [see Additional file 1].

Cross-validation strategy

Stratified three-fold cross-validation was used to avoid overfitting [43]. Subjects with at least one rare CNV were randomly divided into three equal subsets, each with the same proportion of cases and controls. The union of two subsets was used to train the model, while the remaining subset was used as the test set, to assess the classification performance. This process was repeated three times without re-dividing the data-set, so that each subset was used once as test set and twice as training set. Twenty independent cross-validation iterations were performed to estimate the mean and standard deviation of the area under the curve (AUC), to model the classification stochasticity. Absence of overfitting was further assessed by replacing real classification features with randomized features based on gene identity permutation.

Classifier implementation details

For RF, we used the implementation provided by the R/CRAN package ‘randomForest’ version 4.6-7 [44], which follows the original algorithm proposed by Breiman 2001 [20]. We used default settings unless otherwise specified. For CF, we used the implementation provided by the R/CRAN package ‘party’ version 1.0-9 [45]. We used default settings unless otherwise specified. R 2.15.2 was used for all RF and CF analyses.

For the linear SVM, the libSVM package was used [46] in MATLAB R2013a. The cost parameter was kept at default as 1 and class weights were kept even. Each feature was independently normalized and rescaled to a 0-1 interval prior to being input into the classifier.

The Neural Network was built with two middle layers of 100 and 50 nodes each, a learning rate of 0.005 (which affects the rate of connection weight adaptation) with a 0.9 momentum (which affects the acceleration of connection weight adaptation). The network was trained through back-propagation, and unlike the SVM, it did not require feature normalization or scaling. The NN was also implemented in MATLAB R2013a.

Feature relevance metrics and feature selection

Two feature relevance metrics were utilized for RF: Mean Decrease Accuracy (MDA) and Mean Decrease Gini, as implemented in the randomForest package. MDA is calculated by permuting the value of the feature in OOB (out of bag) samples and comparing the accuracy of prediction before and after permutation tree by tree. Mean Decrease Gini is a measure of how each feature contributes to the homogeneity of the nodes and leaves.

MDA was also used for CF, as implemented by the party package function ‘varimp’, using the standard version (unless stated otherwise); in fact, the conditional version, which adjusts for correlations between predictor variables, is extremely computationally demanding.

For MDA-based selection, we ranked features based on MDA, and selected a predefined number of features. The procedure used for step-wise decorrelation and the detailed implementation of Minimum Redundancy Maximum Relevance Feature Selection (MRMR) [46] are described in detail in the supplementary methods and results [see Additional file 1].

Percentage of correctly classified ASD subjects

The percentage of correctly classified ASD subjects was calculated as the number of ASD subjects correctly predicted in at least 15 out of 20 iterations divided by the study total (1,892); the prediction probability cutoff was chosen ensuring that the percentage of correctly predicted ASD subjects without pathogenic or de novo CNVs exceeds the percentage of incorrectly predicted controls by more than 1.5x. Carriers of pathogenic and de novo CNV are easier to classify, thus the ratio described above is a more conservative estimate of the true positive to false positive ratio. Results did not change substantially when requiring subjects to be classified in at least 10 out of 20 iterations.

Feature construction

We analyzed 1,892 ASD subjects (1623 males and 270 females) and 2,342 platform-matched controls (1093 males and 1250 females) with at least one rare CNV (frequency 1% or less); all subjects are of Caucasian ethnicity. Rare CNVs were mapped to gene transcripts; de novo and inherited rare CNVs were labeled as pathogenic, uncertain or benign following clinical annotation guidelines. Univariate burden and pathway analysis, as well as details on the clinical CNV classification, can be found in a separate publication [32].

Classification features were constructed for every subject as gene counts for each gene-set, i.e. counting how many genes participating in a gene-set harbor a rare genic gain or loss; in particular, separate features were constructed for gains and losses. For CNV mapping, we decided to focus on transcripts of known genes as there is a wealth of information that can be used to predict their implication in autism, and their boundaries are well characterized in the human genome; upstream transcription start site regulatory motifs, or other non-coding sequence could be the object of another paper.

Only subjects harboring at least one rare genic CNV were used for classification, as features would be constantly zero for the other subjects, but all subjects were considered when reporting percentage “explained” statistics. This resulted in a subset of 1,570 ASD subjects (80.8%) and 1,916 controls (81.8%), of which 958 ASD subjects and 1,113 controls had at least one genic loss, whereas 1,132 ASD subjects and 1,363 controls had at least one genic gain. 78 ASD subjects harbored at least one genic de novo CNV (30 gains, 51 losses); 57 ASD subjects subjects harbored at least one de novo or inherited CNV recognized as “pathogenic” according to clinical significance annotation [16, 18] (21 gains, 37 losses). All subsets presented a gender composition similar to the full data-sets.

For different classifiers and parameter settings, we classified all ASD subjects, or only the ASD subjects carrying a pathogenic or de novo rare CNV. This was helpful to evaluate the performance of our classification approach for more obviously implicated genes and loci. In addition, we classified either all ASD subjects with any rare genic CNV, or only the carriers of rare genic losses (using only loss-based features) or rare genic gains (using only gain-based features). This was helpful to evaluate the predictive power of the two CNV types.

Classification results: RF and CF

For all classifiers, we used a robust cross-validation approach and tested randomized features to ensure the absence of significant overfitting (Figure 1).

Random Forest (RF) was our original choice for its resilience to overfitting and its capability to handle a large feature space [20]. However RF has been criticized for its positive selection bias towards features with more categorical values [22]; for this reason, Conditional Inference Forest (CF) has been suggested as an alternative to RF [23].

The performance of RF and CF was compared for the 20 curated neurally-relevant features plus the total gene count. We found CF to have slightly higher AUC than RF when classifying all subjects (all CNV types, loss-only, gain-only), whereas it had slightly lower AUC than RF for de novo and pathogenic subjects (Table 2). We inspected the feature relevance metrics for both classifiers, using the Mean Decrease Accuracy and Mean Decrease Gini for RF, and Mean Decrease Accuracy, with or without correlation adjustment, for CF. Both RF metric showed greater relevance for gain-based features when classifying all subjects, in contradiction with the fact that classification based uniquely on losses displayed a better AUC than the one based on gains. On the other hand, CF relevance metrics showed the opposite pattern, with loss-based features more relevant than gain-based features, in accordance with expectations (Figure 2). For this reason, we used CF in place of RF for the rest of the analyses.

Subject	Classifier	All CNV	Gain CNV	Loss CNV
All subjects	RandomForest	0.531±0.005	0.509±0.004	0.544±0.006
All subjects	CForest	0.533±0.004	0.513±0.005	0.546±0.003
De novo	RandomForest	0.805±0.012	0.769±0.024	0.840±0.010
De novo	CForest	0.787±0.008	0.732±0.013	0.846±0.011
Pathogenic	RandomForest	0.913±0.014	0.913±0.012	0.935±0.016
Pathogenic	CForest	0.880±0.012	0.897±0.008	0.922±0.030

CF classification using the 20 curated neurally-relevant features (gains and losses), together with the total gene count, resulted in 7.9% correctly classified ASD subjects and less than 3% incorrectly classified control; this result is reasonably close to the expected contribution of rare CNVs to ASD risk in 10% of the ASD subjects. Losses alone correctly classified 5.9% of the ASD subjects and less than 2.3% incorrectly classified controls; this suggests that gains, besides the ones already recognized to be pathogenic, have limited predictive power, a result in line with univariate burden analysis for this data-set [32].

A detailed analysis of feature relevance for all subjects, or only de novo and pathogenic, is presented in the supplementary methods and results [see Additional file 1].

Classification results for GO and pathways using different feature selection strategies

The set of features based on GO annotations and pathways presented the additional challenges of having many features, with a high degree of mutual overlap and presence of many features with limited or no classification relevance.

To address those issues, we performed classification based on GO and pathway features in two ways: (a) using all features, (b) embedding a feature selection step in the cross-validation procedure; in particular, we used (i) Mean Decrease Accuracy (MDA) based selection, (ii) MDA based selection with stepwise decorrelation, (iii) MRMR (Minimum Redundancy Maximum Relevance Feature Selection). For each procedure, we selected the top 20, top 15% and top 40% ranking features excluding the total gene count, and then added the total gene count. It is important to point out that feature selection was based on the feature relevance metrics calculated on the data subset used for training, and performed independently for every training set, to avoid any overfitting issues.

We assessed the classification performance of different feature selection strategies in comparison to classification without any extra feature selection step or performing manual feature selection based on previous knowledge of ASD biology.

GO without feature selection produced a suboptimal performance (the AUC was very similar to using the gene count alone). When classifying all subjects, the best results for GO were achieved by MDA, either by taking the top 20 features using MDA with decorrelation or by taking the top 15% using MDA without decorrelation; after decorrelation, the top 15% features had a lower performance, suggesting that many relevant yet highly correlated features are removed by decorrelation (see Additional file 1). The best feature selection strategy achieved a slightly better performance than the manually selected GO sub-set (1 sd unit or more), but still inferior to the 20 curated neurally-relevant features (Table 3).

The performance for pathway-based features was markedly worse, with the AUC very close to the total gene count, even when restricting to the manually-selected pathway subset and with the best feature selection strategy (top 20 features using MDA and decorrelation) (Table 3).

These results suggest that feature selection is able to manage the large number and redundancy of GO and pathway-based, although their information content for ASD classification based on rare CNV genes appears to be more limited than the curated neurally-relevant features. Pathway-based features may have particularly disappointing performance results because of the small size of the corresponding gene-sets. Based on these results, we preferred using only the 20 curated neurally-relevant features for the other analyses. It is important to consider that the 20 curated neurally-relevant features include two features derived from the gene-set union of the GO and pathways manually-selected on the basis of known ASD pathobiology.

Classification results using other classifiers and modifying CF parameters

Classification of all subjects using linear SVM or NN achieved lower or comparable AUC compared to CF [see Additional file 1]. This suggests CF is an optimal classification method for this problem.

Additional analysis demonstrated that CF performance does not vary substantially modifying model construction settings (inferential tests and p-value threshold for feature selection operating for tree construction) [see Additional file 1].

Finally, we compared the performance of our classifier based on the 20 neurally-relevant feature to an existing classifier designed to distinguish benign from pathogenic CNV for mental retardation (GeCCO [19]). We found our classifier to have a significantly better performance [see Additional file 1].

Prioritization of subjects

We used the classification results based on the 20 curated neurally-relevant features and total gene count to prioritize subjects with potentially interesting inherited rare losses; we focused on inherited rare variation because it is more difficult to assess its significance outside known ASD loci, and we selected rare losses as these display better classification performance than rare gains.

To ensure the classification results are robust to the algorithm’s stochasticity, we required a subject to be classified in at least 15 cross-validation iterations out of 20. We scanned the prediction probability cutoffs to maximize the number of correctly predicted ASD subjects without pathogenic or de novo CNVs with respect to incorrectly classified controls. At the ASD prediction probability cutoff of 0.52, the correctly classified ASD subjects represented 3.76% of all ASD subjects with at least one genic rare loss, which dropped to 1.25% (12 subjects) when removing subjects with pathogenic or de novo rare CNVs (TP*); this was still larger (about 3.5x) than 0.36% of incorrectly classified controls (FP, false positives), interpretable as a false discovery rate < 30% [see Additional file 1].

We manually inspected these 12 subjects, and found that one of them had a pathogenic loss missed by the clinical annotation, whereas the others had at least one VUS (Variant of Unknown Clinical Significance) each; we used a set of top-scoring loss-based features to identify specific genes that are more likely to contribute to ASD risk within the VUS variant [see Additional file 5].

Finally, we also noticed that prioritized subjects were significantly enriched in females compared to all autism subjects with at least one rare genic loss (Fisher’s Exact Test p-value 0.0066 and OR 2.7). A similar enrichment was reported when considering pathogenic CNVs further classified as “highly penetrant” based on clinical genetic literature [32].