We used a dataset consisting of 732 samples including 366 risk genes of ASD and 366 disease genes (non-ASD) expressed across different developmental stages and time points of 26 brain structures. ASD-Risk was incorporated with a feature selection algorithm IBCGA to select a small set of temporospatial features associated with the risk genes of ASD. ASD-Risk obtained a 10-CV accuracy, sensitivity, specificity, Mathew correlation coefficient (MCC), AUC, and test accuracy of 81.83%, 0.84, 0.79, 0.63, 0.84, and 72.27%, respectively.

The major objective of ASD-Risk is to identify the temporospatial regions that are associated with the risk genes of ASD and accurately predict the risk genes of ASD and non-ASD. ASD-Risk identified 19 temporospatial regions and time points of brain structures and effectively distinguished the risk genes of ASD and non-ASD. Furthermore, prediction performance of ASD-Risk was compared with some standard machine learning methods of Waikato Environment for Knowledge Analysis (Weka) datamining software (Frank et al., 2004). Machine learning methods, including random forest (RF), logistic model tree (LMT), sequential minimal optimization (SMO), and simple logistic method, were used for performance comparison. The performance of ASD-Risk achieved a 10-CV accuracy, sensitivity, specificity, and MCC of 81.83%, 0.84, 0.79, and 0.63, respectively, whereas RF achieved 72%, 0.69, 0.66, and 0.44, respectively, LMT achieved 73%, 0.81, 0.71, and 0.45, respectively; SMO achieved 74%, 0.77, 0.76, and 0.48, respectively, and simple logistic method achieved 73%, 0.81, 0.71, and 0.46, respectively. The prediction performance of ASD-Risk is better than that of some machine learning methods. ASD-Risk on a full training dataset obtained a 10-CV, sensitivity, specificity, AUC, and MCC of 80.05%, 0.81, 0.78, 0.81, and 0.60, respectively. The performance comparison results are shown in Table 1. ASD-Risk identified 19 temporospatial regions of the developmental time points, which can efficiently distinguish risk genes of ASD from non-ASD genes and are listed in Table 2.

Additionally, the prediction performance of ASD-Risk was compared with the prediction model proposed by Cogill and Wang (2016). ASD-Risk achieved a mean training accuracy, mean specificity, and mean sensitivity of 86.03 ± 0.86, 0.92 ± 0.01, and 0.52 ± 0.01, respectively, while the SVM model proposed by Cogill and Wang achieved a mean training accuracy, mean specificity, and mean sensitivity of 76.7, 77.2, and 74.4%, respectively. ASD-Risk performed better than the previous method in terms of mean training accuracy and specificity.

We ranked the identified temporospatial features according to their contribution to the ASD risk gene prediction using main effect difference (MED) analysis (Tung and Ho, 2008). According to the MED analysis, the feature with the highest rank contributed more toward the risk gene prediction. The top 10 ranked temporospatial regions of the brain structures include the posteroventral parietal cortex at 13 post-conception weeks (pcw), the primary visual cortex at 8 years, the posterior superior temporal cortex at 16 pcw, the striatum at 13 pcw, the orbital frontal cortex at 40 years, the anterior cingulate (medial prefrontal) cortex at 8 pcw, the dorsal thalamus at 12 pcw, the amygdaloid complex at 8 years, the primary auditory cortex at 8 years, and the hippocampus at 11 years. The identified temporospatial features, corresponding ranks, and MED scores are given in Table 2. The selected features were described by the time points and the brain structure where the samples were collected.

Furthermore, we analyzed the importance of the top 10 ranked temporospatial features in ASD as follows.

(1) Posteroventral (inferior) parietal cortex—13 pcw: The age of developing embryo or fetus is often described in terms of pcw. A study on autistic children and their families reported that regional perfusion variations are observed in the parietal cortex of autistic children and their first-degree family members (Degirmenci et al., 2008). Raphael et al. reported that mutations in chromodomain helicase DNA-binding protein 8 (CHD8) were associated with ASD diagnosis; CHD8 is expressed in the fetal brain region of the posteroventral (inferior) parietal cortex (Bernier et al., 2014). A genome-wide association study observed a significant enrichment in the expressions of individual genes PANX1, SLC25A12, and PANX2 at quantitative trait loci in the parietal cortex that is implicated in autism disorder (Davis et al., 2012). Furthermore, we analyzed the genetic changes in brain regions and developmental stages that are linked to ASD using genome-wide prediction of the autism associated gene database (Krishnan et al., 2016). Gene expression signatures specific to the posteroventral (inferior) parietal cortex at early mid-fetal were Pumilio RNA-binding family member 2, guanine nucleotide binding protein (G protein), alpha activating activity polypeptide O, Nipped-B homolog (Drosophila), bromodomain PHD finger transcription factor, and K (lysine) acetyltransferase 6A.

(2) Primary visual cortex (striate cortex area V1/17)—8 years: Casanova et al. (2006) investigated a set of postmortem brains of individuals with ASD and observed a greater number of cells per minicolumn in the primary sensory and visual cortexes when compared to the brain samples of controls. A genomic analysis study on the coexpression gene network of ASD observed that ASD genes regulate different biological functions during human cortical development (Parikshak et al., 2013). Protein-altering rare de novo variation-affecting genes with high connectivity in the primary visual cortex were JMJD1C, RBM27, PPM1D, CNOT6, and MLL3 (Parikshak et al., 2013). Genome-wide prediction analysis results revealed that five ASD genes specific to the primary visual cortex at middle late childhood were ATPase, Ca2+ transporting, plasma membrane 2, sema domain, transmembrane domain (TM), and cytoplasmic domain, (semaphorin) 6D, Synaptosomal-associated protein, 25 kDa, zinc finger protein 148, and proteasome (prosome, macropain) subunit, alpha type 1.

(3) Posterior superior temporal cortex (area S1, area 3, 1, 2)—16 pcw: The posterior superior temporal cortex is part of the temporal lobe including Broca’s motor speech area and Wernicke’s area. The important role of the posterior superior temporal cortex is to integrate lexical-semantic and syntactic information during sentence comprehension (Friederici et al., 2009). A systematic analysis of the expression of molecular markers in postmortem brain samples from children with autism observed that focal disruptions of acritical laminar architecture were identified in the posterior superior temporal cortex of young children with autism (Stoner et al., 2014). The predicted ASD genes that are enriched in the superior temporal cortex at early fetal include pumilio RNA-binding family member 2, Nipped-B homolog (Drosophila), bromodomain PHD finger transcription factor, K(lysine) acetyltransferase 6A, and B-cell CLL/lymphoma 11A (zinc finger protein).

(4) Striatum-13pcw: Fuccillo and co-authors investigated multiple genetic mouse models of ASD to identify the abnormalities in striatal circuits that constitute a common pathophysiological mechanism in the development of autism-related behaviors, and reported that striatal dysfunction is intimately associated with the etiology and pathophysiology of ASD (Fuccillo, 2016). A gene expression analysis study observed a highly restricted pattern of expression in contactin associated protein-like 2, a member of the neurexin family in the striatum and frontal and anterior temporal lobes (Alarcón et al., 2008). Araujo et al. (2015) reported that fork head box p1 protein levels substantially reduced by approximately 50% in the stratum relevant to ASD. The preferential gene mutations in ASD were often reported in the striatum region (Shohat et al., 2017). The ASD gene signature of stratum regions consist of the sema domain, transmembrane domain (TM), cytoplasmic domain, (semaphorin) 6D, pumilio RNA-binding family member 1, guanine nucleotide binding protein (G protein), beta polypeptide 1, RING1 and YY1 binding protein, K (lysine) acetyltransferase 6A, and runt-related transcription factor 1, translocated to 1 (cyclin D-related).

(5) Orbital frontal cortex—40 years: Hu et al. (2015) reported that the expression of retinoic acid-related orphan receptor alpha was highly correlated with the NLGN1 gene in the orbital frontal cortex of female adults with ASD. Genome-wide prediction analysis results reported that the genes specific to the orbital frontal cortex-early mid-fetal 2 were spen homolog, transcriptional regulator (Drosophila), runt-related transcription factor 1, translocated to 1 (cyclin D-related), pumilio RNA-binding family member 1, K(lysine) acetyltransferase 6A, and AF4/FMR2 family, member 2.

(6) Anterior (rostral) cingulate (medial prefrontal) cortex—8 pcw: Chandley and co-authors reported that the expression alterations of genes SLC1A1, GRIN1, GRIP1, and GRM8 were observed in the neocortical pyramidal neurons of the anterior cingulate cortex. The reduced NTRK2 expression was observed in the anterior cingulate cortex of individuals with ASD (Chandley et al., 2015). The predicted genes specific to the medial prefrontal cortex-neonatal-early infancy were arginine-glutamic acid dipeptide (RE) repeats, latrophilin 1, catenin (cadherin-associated protein), delta 2, chondroitin sulfate proteoglycan 5 (neuroglycan C), and Rho GDP dissociation inhibitor (GDI) alpha.

(7) Dorsal thalamus—12 pcw: A gene expression analysis study on ASD reported that higher expression levels of CNTNAP2 were observed in the dorsal thalamus of a 19-week fetal brain (Alarcón et al., 2008). A high level expression pattern of autism susceptibility candidate 2 was detected in dorsal thalamus regions in patients with ASD (Bedogni et al., 2010). The predicted genes specific to the mediodorsal nucleus of the thalamus-early mid-fetal were synapsin II, RING1 and YY1 binding protein, neurexin 1, ankyrin 2, and neuronal and AF4/FMR2 family member 2.

(8) Amygdaloid complex—8 years: A microarray study reported that Homer1a is significantly upregulated in the amygdala and altered the function of ASD-related proteins such as metabotropic glutamate receptors and Shank3 (Banerjee et al., 2016). A magnetic resonance spectroscopic study demonstrated the higher concentrations of glutamate/glutamine observed in the amygdala-hippocampal region of individuals with ASD (Page et al., 2006). The predicted ASD genes specific to the amygdala complex-middle late childhood were synapsin II, slit homolog 1 (Drosophila), pleckstrin, and Sec7 domain containing 3,5-hydroxytryptamine (serotonin) receptor 2C, G protein-coupled and glutamate receptor, ionotropic, and AMPA 2.

(9) Primary auditory cortex (core)—8 years: Xiong et al. (2012) reported that autism candidate genes, phosphatase and tensin homolog, were found to be deleted on chromosome 10, and these genes were implicated in the primary auditory cortex of mouse models. The predicted ASD genes specific to the primary auditory cortex-young adulthood were trinucleotide repeat containing 6B, neurofascin, fragile X mental retardation 1, the connector enhancer of the kinase suppressor of Ras 2, and calmodulin binding transcription activator 1.

(10) Hippocampus (hippocampal formation)—11 years: A comparative gene expression analysis study on mouse models reported that two genes, BTBR and En2, were differentially expressed in the hippocampal region of ASD (Provenzano et al., 2016). A strong association between hippocampus and ASD in mouse models has previously been reported (Nadler et al., 2006). ASD genes predicted to be in the hippocampus region-middle late childhood were seizure-related 6 homolog (mouse)-like, neurobeachin, glutamate receptor, ionotropic, AMPA 2, dipeptidyl-peptidase 6, and doublecortin-like kinase 1. The genome-wide predictions of ASD-associated genes specific to different brain regions and time points are shown in Figure 2. The risk gene expression levels at different brain regions that are linked to ASD can be accessed from the genome-wide prediction of the autism-associated gene database (Krishnan et al., 2016).

Additionally, we employed the feature knockout analysis to investigate the individual feature contribution to the prediction accuracy. The feature knockout analysis revealed that removal of each feature reduced the prediction performance on an average of 6.64 ± 0.44. The performance difference showed that the identified features were potential candidates to distinguish the risk genes of ASD and non-ASD genes. The feature knockout analysis results are shown in Table 2.

Furthermore, we compared the relative gene expression levels between the risk genes of ASD and non-ASD genes. The significant difference was observed between the expression levels of the risk genes of ASD and non-ASD genes. Differences in the expression level of genes at temporospatial regions between the ASD and non-ASD are shown in Figure 3 using box plot analysis.

We downloaded the dataset consisting of 336 positive ASD risk gene and 1762 non-ASD disease genes from Cogill and Wang (2016). These data were generated across 13 developmental stages from 8 weeks post-conception to 40 years of age from 26 brain structures. The expressions of RNA-sequencing were read in the units of reads per kilobase of transcript per million mapped reads (RPKM), and aligned using the GENCODE consortium’s annotation release v10 (Harrow et al., 2012).

Cogill and Wang studied the brain gene expression data with 336 positive ASD risk gene and 1,762 non-ASD disease genes for prioritizing the autism risk gene candidates. The positive instances were derived by picking up the top 85% of the genes based upon the expression variance within the BrainSpan dataset as ASD risk genes compiled from the Simons Foundation Autism Research Initiative Gene database (Xu et al., 2012; De Rubeis et al., 2014). The genes associated with diseases but unrelated to the intellectual disability (ID) have been used as negative controls. We followed the same procedure to divide the positive and negative instances. The list of gene IDs and their expression values for the temporospatial timepoints and the risk gene information can be accessed from Cogill and Wang’s (2016) study. In the dataset, genes were instances and temporospatial regions, and time points were features for the training and test sets. In this work, we keep 336 positive ASD risk genes and randomly selected 336 non-ASD genes as a new balance dataset.

Support vector machines (SVMs) are statistical learning algorithms that are explicitly used in solving many biological problems (Vapnik, 1999; Srinivasulu et al., 2015). In this study, we incorporated an optimal feature selection algorithm IBCGA and SVM to build the ASD-Risk. An SVM works implicitly in the feature space by computing only the corresponding kernel K(x_i, x_j) between any two objects x_i and x_j:

where Φ(x) is a mapping function.

IBCGA was used to solve bi-objective combinatorial problems. The IBCGA uses an intelligent evolutionary algorithm (Ho et al., 2004a), which is good at deriving an optimized SVM with feature selection. The IBCGA has been successfully applied in solving several biological problems (Yerukala Sathipati et al., 2016, 2019; Yerukala Sathipati and Ho, 2017, 2018, 2020, 2021).

We used common genetic algorithm (GA) terms “GA-gene” and “GA-chromosome.” In this problem, “GA-chromosome” contains 524 binary genes. Two 4-bit “GA-genes” were utilized for tuning the C and γ of the SVM. So, this method encodes the parameter C in the 2^–7 to 2^–8 interval and 16 values of γ. Normalized and digitalized gene expressions were used as the input in the SVM classifier. Gene expressions corresponding to temporospatial time points were considered as features. Parameter tuning of ASD-Risk is as follows; the candidate feature range selected by the IBCGA is r_begin = 30 and r_end = 10. The steps involved in the IBCGA are as follows.

This work used the following equations to measure the performance evaluation.

where TP is true positive; TN is true negative; FP is false positive; FN is false negative; and MCC is the Matthews correlation coefficient.

To identify the lncRNA expression within genomic regions, Weka data mining software was used. Weka can implement all major learning techniques for classification and regression methods (Frank et al., 2004). We used the SVM, sequential minimal optimization (SMO), random forest, logistic model tree (LMT), and simple logistic method to distinguish the risk genes of ASD and non-ASD genes.