To collect all existing datasets that used children and adults with ASD’s peripheral blood for transcriptome studies, we searched the datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the condition terms “Autism spectrum disorder,” “Autism,” “Autistic disorder,” “ASD,” the tissue term “blood,” and the organism term “Homo Sapiens.” The inclusion criteria were as follows: the datasets contained gene expression profiling by microarray or high-throughput RNA sequencing, the datasets included ASD and typically developing (TD) samples, and 20 or more samples were in the dataset. Thus, six eligible GEO Series (GSE) datasets were adopted, including GSE6575 (25), GSE18123 (26), GSE42133 (27), GSE111175 (28), GSE26415 (29), and GSE89594 (30). All transcriptomic datasets are listed in Supplementary Table 1 . The GSE18123 dataset consisted of two RNA sequencing profiles from different GEO Platforms (GPL), so there were five datasets from children and two from adults. Among the five children’s datasets, four were used as “training” datasets, and one was used as the validation dataset. Figure 1 illustrates the workflow of this study. The R package “limma” (v3.50.3) was used to analyze each gene expression matrix with the threshold of p < 0.05 (31).

The biological functions and corresponding pathways of the identified DEGs were determined through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the R package “clusterProfiler” (32) (v4.2.2) with a significance level of p_adjusted < 0.05. The online databases GENEMANIA (https://genemania.org/) (33) and Metascape (http://metascape.org/gp/index.html#/main/step1) (34) were used to predict gene function. Gene Set Enrichment Analysis (GSEA) was conducted in R software using the hallmarks gene set “h.all.v7.5.1.symbols.gm” (35). The results were visualized with the R packages “enrichplot” (v1.14.2) and “ggplot2” (v3.3.6).

The CIBERSORT analytical tool (https://cibersortx.stanford.edu/) was used to evaluate the abundance and differences of 22 types of immune cells between the ASD and TD group (36), including different kinds of B cells, dendritic cells, macrophages, mast cells, NK cells, T cells, eosinophils, monocytes, neutrophils, and plasma cells. The results were analyzed using the R packages “ggplot2” and “Kruskal.test”. For the training dataset, LASSO logistic regression was carried out to further select significantly different immune cells in combination with Kruskal.test. The Wilcoxon test and receiver operating characteristic (ROC) curves were calculated using R-based tools to assess the classifier performance for the LASSO regression model (v4.1.0).

A PPI network was constructed among the consistently changing DEGs using the STRING database (https://string-db.org) (37). The CytoHubba module in Cytoscape (v 3.7.1) was used to score the top 10 node genes using 12 different algorithms, namely MCC (Maximum Clique Centrality), DMNC (Density of Maximum Neighborhood Component), MNC (Maximum Neighborhood Component), Degree, EPC (Edge Percolated Component), BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, and Clustering Coefficient. The hub genes generated from each algorithm were shown by the R package “UpSet” and validated using GSE111175.

To validate the most significant immune cells in children, the whole blood samples from 30 children with ASD and 30 age and sex-matched TD children were collected at the 3rd Affiliated Hospital of Zhengzhou University. All participants provided informed consent to participate in the study, which was approved by the Medical Ethics Committee of the 3rd Affiliated Hospital of Zhengzhou University (Ethical number 2020–56). All blood samples were processed within 8 h for flow cytometry with a mix of antibodies according to the manufacturer’s instructions. The antibodies purchased from BD Biosciences included CD45 (HI30, Cat# 564105) CD14 (M5E2, Cat# 561712), CD16 (3G8, Cat# 563692), and HLA-DR (G46-6, Cat# 560896). The data were analyzed using FlowJo software (TreeStar) and presented using the UMAP method.

The association of the hub genes with the infiltration levels of immune cells was explored and analyzed in R, and the correlation results were visualized using the “ggplot2” and “ggstatsplot” packages.

The GSE111175 dataset was used to assess the predictive effectiveness according to the validation group between the children with ASD and the children in the TD group using ROC curve analysis. The results were visualized with the “pROC”, “ggplot2”, and “ComplexHeatmap” packages.

The starBase (predicted program ≥ 2) (38) and miRDB (39) databases were used to predict the upstream miRNAs of hub mRNAs by taking intersections. Cytoscape (v 3.7.1) was used to visualize the mRNA-miRNA network.

The DSigDB database (http://tanlab.ucdenver.edu/DSigDB) was used to predict potential agents for hub genes related to ASD, with agents–hub genes ≥3 as the condition for the agents screen. The results were displayed using a Sankey diagram based on the “ggalluvial” (v0.12.3) package.

All data from the bioinformatics analyses were analyzed in R, and the data from routine blood tests were analyzed using SPSS software version 23 (IBM Corp). All data are presented as the mean ± standard deviation.

We conducted the analysis using the four children’s datasets and two adults’ datasets. For the children’s datasets, the “RRA” and “Batch correction” methods were used to get accurate DEGs from multiple children’s datasets. For the “Batch” dataset using the same platform, including GSE6575-GPL570 and GSE18123-GPL570, batch effects were removed using the ComBat function of the R package “SVA” (40). The 3D principal component analysis plots indicated that the processed data were more reliable after batch effect removal ( Figures 2A, B ). For “RRA” dataset using different platforms, including GSE18123-GPL6244 and GSE42133-GPL10558, the batch effect was corrected using the R package “RobustRankAggreg” (RRA) (41) ( Figure 2C ). As analyzed by “limma”, the final list of DEGs used for further analysis was generated as shown by the Venn diagram in Figure 2D . There were a total of 113 DEGs, of which 95 were consistently changing DEGs between these two datasets in the childhood group, including 57 co-upregulated and 38 co-downregulated genes ( Supplementary Table 2 and Figure 2D ). For the two adult datasets, DEGs were extracted as the intersection between GSE26415 and GSE89594 as shown in Supplementary Figure 1 . A total of 25 consistently changing DEGs were screened between the GSE26415 and GSE89594 adult datasets, including 17 co-upregulated and 8 co-downregulated genes ( Supplementary Table 3 and Supplementary Figure 1A ).

To reveal the biological processes and signaling pathways behind the DEGs in children and adults, we performed a systematic analysis, including GO, KEGG, GSEA, and DisGeNET database prediction. For the childhood group, the 95 consistently changing DEGs were analyzed with Metascape and GSEA. In total, biological functions related to neuronal development (axon guidance, nervous system development, and Roundabout (ROBO) receptors), immunity (adaptive immune system, IFN-alpha/gamma response, complement, and IL2/Stat5 signaling), and cell development and metabolic processes (mTORC1 signaling, mitotic spindle formation, KRAS signaling, heme metabolism, and protein secretion) ( Table 1 ; Figure 3A ) were identified.

Interestingly, the DisGeNET prediction suggested that the DEGs were mainly associated with erythrocyte-related parameters and disorders such as mean corpuscular hemoglobin (MCH), extramedullary hematopoiesis function, red cell distribution width (RDW), mean corpuscular hemoglobin concentration (MCHC), and anemia ( Figure 3B ; Supplementary Table 4 ). Indeed, these results were confirmed and validated when using peripheral blood samples collected from children with ASD. We observed significant differences between the ASD and TD groups for red blood cell count (RBC), hematocrit (HCT), MCH, coefficient of variation of RDW (RDW-CV), and MCHC ( Table 2 ). Moreover, GO enrichment analysis of the DEGs suggested biological functions associated with erythrocyte differentiation ( Supplementary Figure 2 ).

It is noteworthy that when analyzing the functional enrichment of DEGs using the “RRA” and “Batch” method separately, we found that the B cell receptor (BCR) signaling pathway was significantly enriched in both ( Supplementary Figures 3A, B ). Similarly, when analyzing all 113 DEGs instead of the consistently changing 95 DEGs between the two methods using the Metascape database, we also found the BCR signaling pathway to be a significant pathway ( Supplementary Figure 3C ). In addition, several DEGs, including PIR-B (LILRB1), IgA, CD81, and VAV, were found to be involved in the BCR signaling pathway (data not shown).

For the adult group, Metascape analysis of the 25 consistently changing DEGs suggested that deubiquitination, regulation of hormone levels, cell morphogenesis involved in neuron differentiation, and positive regulation of catabolic process were the most enriched terms ( Supplementary Figure 1B ). GENEMANIA analysis showed that the TGF-β/SMAD signaling pathway was a significantly altered signaling pathway ( Supplementary Figure 1C ), and some metabolic disorders and immune disease-related processes were significantly enriched according to the DisGeNET analysis, such as hyperinsulinism, scleroderma, juvenile arthritis, inflammation, and celiac disease ( Figure 3C ). In addition, GSEA analysis showed that oxidative phosphorylation, MYC targets v1, E2F targets, and protein secretion were suppressed, while heme metabolism and KRAS signaling DN were activated, all of which were the opposite of what was seen in the childhood group. It is also noteworthy that estrogen response was among the biological processes enriched in the adult ASD datasets compared to children ( Figure 3D ).

Overall, the biological process analysis using DEGs showed that immune function, neuronal development, and metabolic disorders were common in both children and adults with ASD. However, there were obvious differences between the two groups. For example, erythrocyte differentiation was involved in children with ASD, while estrogen response affected the course of adults with ASD.

The analysis above clearly showed that immune function is among the most significantly affected biological functions in both children and adults with ASD. In addition, immune responses and inflammation have been suggested to play a crucial role in childhood ASD progression, so we conducted a comprehensive analysis of the immune cell profiles in the selected datasets. The “Batch” dataset for children (101 ASD children and 45 TD children) and the GSE26415 dataset for adults (21 ASD adults and 21 TD adults) were used because of the high quality of their sequencing data and the availability of the CIBERSORT algorithm.

Using the CIBERSORT algorithm, 22 immune cell types were identified in each sample from each group ( Supplementary Figure 4 ). Among these, naive B cells, monocytes, neutrophils, resting NK cells, naive CD4⁺ T cells, and CD8⁺ T cells were the main immune cells in children, while monocytes, neutrophils, memory-activated CD4⁺ T cells, and CD8⁺ T cells were similarly enriched in adults. In addition, there were significant differences between the ASD and TD groups regarding resting/activated dendritic cells, M0/M2 macrophages, and monocytes in children and M0 macrophages, resting mast cells, and resting/activated NK cells in adults ( Figure 4A ).

LASSO regression analysis can more accurately identify significant differences in immune cells. As shown in Supplementary Table 5 and Figure 4B , the LASSO regression model was found to have good classifier performance according to Wilcoxon’s test and area under the ROC curve (AUC) analysis. When further analyzing the data in children with ASD, 12 types of immune cells were extracted ( Figure 4C ). Among these, four cell types (monocytes, M2 macrophages, and resting/activated dendritic cells) from two algorithms (Kruskal.test and LASSO logistic regression) were identified.

It has been suggested that monocytes play a major role in ASD development in children (42). To validate the results obtained from the integration analysis above, monocyte phenotypes were further identified and confirmed using peripheral blood samples. We collected a total of 60 peripheral whole blood samples, including 30 ASD and 30 age and sex-matched TD children, and performed flow cytometry analysis ( Figure 5A ). In the UMAP analysis, each of the monocyte subtype populations, including intermediate monocytes, classical monocytes, and nonclassical monocytes (ncMos), were shown as distinct cell clusters ( Figure 5B ). The flow cytometry analysis showed that among all analyzed monocyte subtypes, there were significant differences in total monocytes (p < 0.05) and ncMos (p < 0.001) ( Figure 5C ).

Hub genes usually refer to genes with a connectivity degree greater than 10 in the genetic interaction network, and these play important roles in biological systems. Here, we used the “Batch” and “RRA” datasets as the training dataset and GSE111175 as the validation dataset for children with ASD. Notably, due to the integration and consistency of the data in the “Batch” and “RRA” datasets, we used the expression data from the “Batch” dataset to represent the expression levels of the hub genes in the training dataset.

The PPI network used to identify the hub genes was obtained when 95 consistently changing DEGs in children with ASD were analyzed by STRING ( Supplementary Figure 5 ). Using the CytoHubba function, there were a total of 15 hub genes ( Figure 6A ), including AHSP, ALAS2, SELENBP1, AP2A1, BCL2L1, CD3G, MAP1LC3A, CD69, DCAF12, EPB42, GMPR, IGF2R, LILRB1, PSMC4, and SLC4A1, and the expression levels of the hub genes were determined in the training ( Figure 6B ) and validation ( Figure 6C ) datasets in children with ASD and TD. Clearly, the validation dataset confirmed only ALAS2, SELENBP1, CD3G, MAP1LC3A, CD69, DCAF12, IGF2R, PSMC4, and SLC4A1 as significantly differently expressed hub genes between children with ASD and TD and not the genes AHSP, AP2A1, BCL2L1, EPB42, GMPR, and LILRB.

Furthermore, as shown in Figure 7A , ALAS2, CD3G, DCAF12, IGF2R, PSMC4, SELENBP1, and SLC4A1 were significantly and consistently altered hub genes in the training and validation datasets, which suggested that these genes were critical to ASD and might be used as predictors for ASD in children. An AUC more than 0.60 was considered to be relatively good predictive accuracy. The AUC values of PSMC4 and ALAS2 were 0.633 (95% CI 0.534–0.732) and 0.651 (95% CI 0.559–0.744) in the training group and 0.677 (95% CI 0.574–0.779) and 0.625 (95% CI 0.517–0.733) in the validation group, respectively. In the logistic regression model the combined AUC values of PSMC4 and ALAS2 reached 0.668 (95% CI 0.575–0.762) and 0.729 (95% CI 0.630–0.829), respectively, in both datasets ( Table 3 ; Figure 7B ). In addition, we also undertook diagnostic performance analyses of the remaining five possible marker genes, including single and multiple gene combination diagnostic analyses. However, we did not find any additional predictors in the two datasets when considering higher diagnostic efficiency and lower cost. Therefore, PSMC4 and ALAS2 might act as predictors with moderate strength (43).

To further estimate the sensitivity of using these hub genes as indicators for ASD severity, we extracted the Autism Diagnostic Observation Schedule Social Affect (ADOS-SA) deficit scores that were positively correlated with ASD severity from the GSE111175 dataset. The children with ASD were divided into the following three groups according to severity: mild (5 to 11), medium (12 to 15), and high (16 to 21). The heatmap showed the expression of 15 hub genes in the three groups ( Figure 7C ), and we used scatter plots to show gene expression patterns among the three groups to make it more intuitively understandable. In particular, LIlRB1 could distinguish between high and mild severity ASD (p < 0.05), and CD69 had relatively lower expression in the medium severity group compared to the mild and high severity groups (p < 0.05) ( Figure 7D ). These results indicate that LIlRB1 and CD69 can be used as indicators for ASD severity in children.

To investigate the correlations between hub genes and immune cell profiles in children with ASD, we performed Spearman correlation analysis using the “Batch” dataset. First, we explored the correlation between immune cells in children with ASD and TD, and we found that monocytes were positively associated with M0/M1 macrophages but negatively related to M2 macrophages in children with ASD, whereas M0 macrophages were negatively related to M2 macrophages in TD children ( Figure 8A ). Second, the correlation between 15 hub genes and 4 significantly differentially expressed immune cells ( Figure 8B ) showed that activated dendritic cells were positively correlated with AHSP, LILRB1, and CD69 (r > 0.20, p < 0.05) but negatively correlated with PSMC4 (r < –0.20, p < 0.05). In addition, M2 macrophages showed a positive correlation with AP2A1, MAP1LC3A, EPB42, BCL2L1, GMPR, SLC4A1, and AHSP (r > 0.20, p < 0.05) and a negative correlation with CD3G (r < –0.20, p < 0.05) in children with ASD ( Figure 8C ). Lastly, the graphs of the linear correlation indicated that monocytes were significantly and positively related to LIlRB1 (r > 0.40, p < 0.001) ( Figure 8D ).

To gain further insight into the regulatory mechanisms and to identify potential agents affecting ASD in children, we next tried to identify the miRNAs and agents that influence identified 15 hub genes.

A total of 118 candidate miRNAs were predicted to target the 7 identified hub genes. There were 10 miRNAs targeting AP2A1, 18 miRNAs targeting BCL2L1, 33 miRNAs targeting CD69, 37 miRNAs targeting DCAF2, 5 miRNAs targeting GMPR, 33 miRNAs targeting IGF2R, and 2 miRNAs targeting LIlRB1 ( Figure 9A ). We then summarized those miRNAs that regulated more than one hub gene in Table 4 . Of these, hsa-miR-342-3p and hsa-miR-1321 interacted most intensively with robust hub genes (target hub genes = 3).

In addition, using the DSigDB database, the 23 agents affecting ASD were shown in a Sankey diagram with agents–hub genes ≥ 3, and the agents-disease associations were also displayed, including positive, negative, and unknown correlations with ASD status ( Figure 9B ).