We used an autism data set from the SSC [30]. The SSC data consist of 2600 simplex families, each of which has 1 child affected with ASD (a proband), unaffected parents, and at least one unaffected sibling. The data consist of 3931 individuals whose exome sequences are available (Figure 1A), and 2249 samples of these individuals are labeled as diseased group (ASD). From the SSC data set, we can query the specific variables for exome variants (Figure 1A), and the variants are in the variant call format (VCF). There are more than 1.5 million variants in the data set, which has the genotype information along with read depth, allele depth, and genotype quality. In the VCF data, VCF_GT represents the genotype quality, encoded as allele values separated by “/,” such as “0/1” and “2/3”, where 0 represents the reference allele, 1 for the first allele listed in the alternate allele, 2 for the second allele listed in the alternate allele, and so on. Thus, VCF_GT can be “0/0”, “0/1”, “2/0”, “1/2”, and so on. The read depth is denoted as VCF_DP, and the conditional genotype quality is denoted as VCF_CQ. We mainly used the information of VCF_GT, VCF_DP, and VCF_CQ in the SSC data for this study. The Harvard Medical School Research Ethics Committee approved this study.

For all the variants, we have their unphased genotype information using the format of VCF_GT. To make the data processable for deep learning models, we encoded the VCF_GT data by creating categorical values to represent different types of genotype [27]. Specifically, 0 denotes that both allele values are reference alleles, 1 represents one allele value is a reference allele and the other one is the alternate allele, and 2 represents both are alternate alleles. For example, 0/0 is made as 0, 1/0 is made as 1, 1/1 is made as 2, and so on. Therefore, the variants consist of 3 categories: 0, 1, and 2. We used VCF_DP (read depth at a position for a sample) and VCF_CQ (conditional genotype quality) as a filter to control the genotype information quality (Figure 1A). We extracted the genotype information for each variant that has a read depth no less than 10 and genotype quality no less than 90 [31]. Therefore, the genotypes of read depth less than 10 (VCF_DP < 10) or genotype quality less than 90 (VCF_CQ < 90) were excluded. Since we only explored common variants in our study, we removed all the variants with occurrence frequency less than 1% over the whole data set, resulting in 153,347 variants selected as common variants after the genotype quality filters (Figure 1A). We used these common variants for our study. After selecting the common variants, the SSC samples were partitioned into 2 sets based on random sampling of individuals into a training set (80%) and a hold-out test set (20%). There was no overlap of individuals across the 2 partitions. The test set was only used after model fitting to assess performance.

As the number of variants was too large to apply deep learning models directly, to construct the features for the deep learning models, we used feature selection to reduce variant dimension (Figure 1B). Feature selection is one of the core concepts in machine learning that hugely impacts the performance of a model [32-35]. The data features that are used to train machine learning models have a huge influence on the performance that we can achieve. Therefore, our hypothesis is that not all variables contribute to the predictive performance of the models we built. Variant selection is the process wherein we automatically select those features that contribute most to our prediction accuracy and are considered as contributory variants to ASD diagnosis. Therefore, significant common variation selection was applied because variant selection is the process of removing redundant or irrelevant features from the original data set to reduce overfitting. To this end, we analyzed the importance scores of the common variants that are related with ASD development mechanisms in the training data set. For each individual, a 153,347-dimensional vector was constructed, corresponding to 153,347 common variants identified from the data preprocessing. Chi-square test was applied to evaluate the importance of each variant to distinguish the class in order to select the most significant common variants. Given the training data D, we estimated the following quantity for each variant and ranked them by their scores:

where N is the observed frequency in D and E is the expected frequency, e_t takes the value 1 if the training data contains term t and 0; otherwise, e_c takes the value 1 if the training data is in class c and 0 otherwise. For each variant, a corresponding high score indicates that the null hypothesis H₀ of independence (meaning the individual’s category has no influence over the term’s frequency) should be rejected and the occurrence of the variant and class are dependent. In this case, we select the variant for the ASD diagnosis prediction. We used the implementation from scikit-learn [36] for “Chi-Square Feature Selection” with default settings.

By calculating the chi-square scores for all the variants, we can rank the variants by the chi-square scores and then choose the top ranked variants as significant variants for model training. Figure 2A lists the variant importance of the high scoring 20 features (variants) that are selected via the chi-square test. In Figure 2A, while the Y-axis corresponds to variant IDs of the variants, the X-axis corresponds to relative importance, which is calculated using the chi-square score. We selected the top 100 most significant variants as inputs to train a deep learning classifier. Therefore, the number of input contributory variants to our classifier after selection was 100 for ASD prediction. In order to analyze whether the variation data can be divided into 2 clusters representing control and ASD cases, the first 2 groups of data were obtained using t-distributed stochastic neighbor embedding (t-SNE) as an unsupervised learning approach. The visualization of clusters for the top 100 variants using t-SNE in both case group and control group is shown in Figure 2B. From the visualization, accurate genetic classification of control group versus ASD is possible using 100 common variants determined to be highly significant. Therefore, for each individual in the training set, a 100-dimensional input vector was constructed corresponding to 100 selected significant common variants for training a deep learning model.

The overall framework of the proposed DeepAutism (Figure 1) consists of 3 components, namely, data preprocessing, variant selection, and neural network classifier. Figure 1C illustrates the convolutional neural network (CNN) architecture. We used Keras and TensorFlow version 2.0 for constructing and training the CNN model. We used a block of two 1D convolutional layers, followed by a max-pooling layer to generate feature maps that contain only the most important features. The max-pooling layer is followed by a dropout layer to avoid overfitting the data. Then, the learned feature maps are combined using a fully connected layer. The final layer contains a sigmoid function to produce probabilities of output from 0 to 1, with the diseased group belonging to class 1 and the control set belonging to class 0. All the parameters, including the weights and biases of hidden layers, are learned through backpropagation [37]. The detailed network topology used in our CNN architecture is shown in Multimedia Appendix 1.

For training, DeepAutism uses a set of selected common variants (top 100 significant common variants) to estimate the probabilities of an individual belonging to control case or autism. For a set of variants v from a testing individual, DeepAutism computes a probability p(v) using 4 states:

The sigmoid function is used for computing probabilities of a set of variants v belonging to either control group or autism group, and the produced probabilities are from 0 to 1, with the control set belonging to class 0 and the ASD group belonging to class 1. The convolution stage (convf) scans a set of filters as feature maps across the variants. Each neuron consists of a rectified linear unit (ReLU) activation function to introduce nonlinearity between 2 neural networks. The pooling layer only picks the maximum values from the convolved feature maps. Since the variant data are categorical values, one-hot encoding is employed to ensure that the DeepAutism model is unbiased and does not favor one genotype over the others. The DeepAutism is then trained using mini-batch gradient descent by backpropagation algorithm [37]. The performance was evaluated using the area under the receiver operating characteristic curve (AUC). We also used the most common procedure for evaluating classifiers for ASD prediction, including accuracy, sensitivity (recall), specificity (precision), F1-score, and false discovery rate, in which the lower value indicates better performance to evaluate the classifiers for ASD diagnosis.

Apart from CNNs, we also employed conventional machine learning techniques to evaluate the effectiveness of DeepAutism for classifying autism diagnosis. The conventional machine learning models that we compared were random forests, logistic regression, and Naive Bayes. We used the same training and test data (with the selected 100 common variants) for the conventional machine learning models as used for the DeepAutism model, aiming to evaluate whether the CNN model outperforms other machine learning classifiers. To evaluate whether the selected top 100 common variants are significant for ASD diagnosis, we also compared the chi-square–based variant selection method with random variant selection by using the same training and test data sets. We randomly selected 100 common variants as inputs that were fed into both DeepAutism and conventional machine learning models to compare the changes in their performance.

Statistical analyses focused on the selected top 100 common variants, which most significantly contributed to the classifiers of ASD. Of the 100 common variants within our classifier, 66% are exonic mutations and 23% are intronic mutations, while small proportions are splicing mutations or from an untranslated region. Within the 66% exonic mutations, about half are synonymous single nucleotide variants and about half are nonsynonymous single nucleotide variants. It is important to point out that the selected contributory common variants were significantly enriched on chromosome X while chromosome Y is also discriminatory in identifying individuals with ASD from individuals without ASD.

A number of variants were populated by the same genes. Related to the contributory common variants, the statistically significant genes were ARSD, MAGEB16, and MXRA5. There are 18 common variants in the ARSD gene. ARSD is a protein-coding gene and is located within a cluster of similar arylsulfatase genes on chromosome X, while a related pseudogene has been identified in the pseudoautosomal region of chromosome Y. Variants rs209372, rs2109135, and rs1047248 in 3 genes, namely, NRK, TLR8, and MAGEA4, respectively, have the highest scores in determining an individual’s classification as with ASD or with no ASD.

After the training phase was over, we picked the same common variations from the test data for each individual. We used the rest of the 787 samples for testing. Based on the trained DeepAutism model, each test individual was predicted the probabilities of belonging to the control group or the diseased group. The deep learning model was extremely accurate in classification of the holdout test set with an AUC of 0.955 (Figure 3A). Figure 3B describes the performance of the DeepAutism classifier on the test data. DeepAutism predicted ASD in 423 samples out of 456 samples with ASD.

Apart from deep learning, we also employed Naive Bayes, logistic regression, support vector machine, random forest, and deep neural network classifiers to compare the prediction of ASD diagnosis. We applied five-fold cross-validation to evaluate the selected significant common variants. Our classifier performed better than the conventional machine learning techniques in terms of AUC, accuracy, specificity, sensitivity, and F1-score. As shown in Table 1, accuracy was 0.886 in the case of DeepAutism, followed by 0.808 for random forest in the same test data set for ASD diagnosis prediction. DeepAutism also yielded the best sensitivity of 0.881 for prediction of ASD and best specificity of 0.893 for non-ASD prediction. The false positive (discriminatory) rate is minimum for DeepAutism with 7% compared with other machine learning techniques. These results are shown in Table 1.

We assessed the classification performance by using randomly picked 100 common variants as inputs to train classifiers. We used the same training and test data as in the above experiment. As shown in Table 2, when the classifiers classify ASD using randomly selected common variants, all the classifiers achieved reduced performance compared to using selected significant common variants. For instance, the AUC and accuracy of DeepAutism dramatically dropped from 0.955 to 0.670 and from 0.885 to 0.689, respectively. The random 100 common variants yielded accuracy of 0.454 and 0.583 using Naive Bayes and logistic regression classifiers, respectively, which is like random guessing. This revealed that the random 100 common variants are not discriminative in distinguishing ASD diagnosis. These results suggest that variant selection is important for identifying significant common variants that are more correlated and significant in improving the classification accuracy.