We collected the expression profile datasets GSE111176 (22, 23), GSE77103 (24) and GSE38322 (25) of ASD from the GEO database by R package “GEOquery” (26). The dataset GSE111176 was contributed by Vahid H Gazestani et al., which contains information on the expression profiles of leukocytes derived from 119 ASD patients and 126 controls. One validation set GSE77103, which includes the expression patterns of peripheral blood mononuclear cells obtained from 4 healthy controls and 4 ASD patients, was submitted by Inoue R et al. Another validation set GSE38322, contains the expression patterns of brain tissues sourced from 18 controls and 18 ASD patients, was contributed by the Matthew Ginsberg et al. All samples were included in this analysis. The GPL10558 (Illumina HumanHT-12 V4.0 expression bead chip) was used for the GSE11176 dataset and the GSE38322 and the GPL17077 (Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381 (Probe Name version)) was used for the GSE77103 dataset. During data processing, the “removeBatchEffect” function of the “limma” (27) package was used to remove the batch effect of the sub-datasets of the GSE111176. And the expression profiles of three datasets were normalized by the “normalizeBetweenArrays” function of the “limma” package. Finally, Gene probes were annotated with official gene symbol, and mean values were taken if multiple gene probes matched to the same gene.

From the Molecular Signatures Database V7.0 (MSigDB) (28), we retrieved the endoplasmic reticulum stress-related gene sets “GOBP NEGATIVE REGULATION OF RESPONSE TO ENDOPLASMIC RETICULUM STRESS” and “GOBP POSITIVE REGULATION OF RESPONSE TO ENDOPLASMIC RETICULUM STRESS,” which includes a total of 256 genes.

The “limma” package of the R software was applied for differential gene analysis. The screening criteria for the DEGs of normal and ASD patients were |logFC| > 0 and an adjusted value of p<0.05. As for the ASD patterns, the threshold was set as |logFC| > 0.2 and pvalue<0.05. The R packages “ComplexHeatmap” (29) and “ggplot2” (30) were used to plot the findings of the differential analysis. The R package “RCircos” (31) is used to display information on the chromosomal location of differential genes.

Two machine learning (random forest (RF) and artificial neural network (ANN)) approaches were applied to construct the ASD classifier. Random forest regression is a machine-learning algorithm that takes an ensemble learning approach for prediction. It is made up of various decision trees, each trained on random features. ANN consists of the input layer, the hidden layer and the output layer. The number of neurons in the input layer indicates the number of features being evaluated; while the neurons in the output layer are the dependent variables; Each neuron in the hidden and output layers is linked to all neurons in the previous layer by the corresponding numerical weights. To develop the predictive signature, the differentially expressed endoplasmic reticulum stress genes were subjected to random forest analysis using the R package “randomForest” (32), and genes with gene importance scores greater than 4 were selected for the further neural network created using the “NeuralNetTools” (33) and “neuralnet” packages (34). Using the R package pROC (35), the ROC curve was drawn and the area under the curve was determined to assess the distinguishing performance of the signature. Additionally, the external validation set GSE77103 confirms the classifier’s diagnostic effectiveness.

Weighted Gene Correlation Network Analysis (WGCNA) (36) is mainly used to identify co-expressed gene modules and examine the connection between gene networks and phenotypes. For computational efficiency, the top 25% of variant genes based on an analysis of variance were selected for WGCNA analysis using the “WGCNA” package (36). The soft threshold was calculated by the “picksoftthreshold” function and scale-free networks were built. Subsequently, the adjacency matrix was transformed into a topological overlap matrix (TOM). The hierarchical clustering and dynamic tree-cut method were used to identify gene modules. Co-expression modules were defined using a minimum module size of 50 genes and by merging modules with a module eigengene dissimilarity below 0.2. Finally, the correlations between modules and clinical features were determined by spearman correlation analysis.

The varied ASD patterns were identified using the ConsensusClusterPlus (37) R package, and the clustering was carried out using the pam method with a sampling proportion of 0.8, and it was employed for 100 iterations to assure clustering stability. The clustering score for the cumulative distribution function (CDF) curve was calculated to estimate the ideal number of clustering. The consensus clustering’s reliability was confirmed by Principal Component Analysis (PCA).

Based on the differentially expressed genes, Gene Ontology (GO) (38) terms, which provide context for cellular component (CC), molecular function (MF), and biological process (BP), as well as the Kyoto Encylopaedia of Genes and Genomes (KEGG) analyses (39), were employed to assess the biological changes of ASD patterns. We used the R package “clusterprofiler” (40) for GO and KEGG analysis and the R package “pathview” (41) to visualize the significantly distinct KEGG pathways between ASD patterns. Terms with a false discovery rate (FDR) < 0.05 were considered statistically significant. Gene Set Enrichment Analysis (GSEA) (42) is a method that determines whether a priori-defined set of genes show statistically significant expression variations between two biological states. From the MSigDB database, we have gathered the gene set “c2 cp all v7.0 symbols.gmt.” The GSEA analysis was performed using the R package “clusterprofiler” and it was determined that the |normalized enrichment score| (|NES|) > 1 and FDR < 0.25 were significant results.

Master regulator analysis (MRA) (43) is an algorithm used to infer transcription factors (TFs) controlling the transition between the two phenotypes and the maintenance of the latter phenotype. In this study, we applied a recently developed coexpression-based gene network inference and interrogation tool, corto (44), to determine the master regulator genes of the transition from normal control to ASD subtypes. The collection of human transcription factors was used as the potential centroids for the model input, along with the gene expression matrix of normal controls and two ASD subclusters. And then, a list of master regulators sorted by normalized enrichment scores (NES) was obtained through the “mra” function of “corto” package. The NES is positive if the centroid network is higher in the sample vs. the mean of the dataset, negative if lower.

The analytical tool CIBERSORTx (45), developed by Newman et al., uses gene expression data to perform cell-type deconvolution and offers an estimate of the abundances of member cell types in a mixed cell population. The online tool CIBERSORTX (https://cibersortx.stanford.edu/) was applied to estimate immune cell infiltration in ASD patterns. Single sample gene set enrichment analysis (ssGSEA) was performed according to previous studies, which was used to determine enrichment scores (ES) for each coupling of a sample and immune reaction gene sets in the ImmPort database.¹ The R package “GSVA” (46)was used for the ssGSEA analysis and “ggplot2” for visualization. The Wilcox test was used to assess the enrichment scores representing immunocyte abundance and immune response activity between ASD patterns.

NetworkAnalyst² (47) is an online visualization and analysis platform which can perform differential gene expression analysis, protein–protein interaction analysis and integrated analysis of multiple datasets. With the use of this platform, we were able to obtain chromatin immunoprecipitation data from the ENCODE database (48) for the construction of TF-gene network and experimentally verified data for the construction of miRNA-gene interaction network from the miRTarBase database (49). Afterward, the Cytoscape software (50) was utilized for visualization.

The Connectivity Map database (Cmap) (51) focuses on the interaction between compounds, genes and disease states. It contains induced gene expression profiles of more than 25,200 perturbagens of which ~19,800 are small molecules, 314 biologics and ~ 5,075 genes with altered function by shRNA, cDNA, and/or CRISPR. These perturbagens were assayed in different cell lines to produce around 1.3 million individual gene expression profiles corresponding to ~473,000 gene expression signatures (52). Each ASD subcluster’s differentially expressed genes were submitted to the Camp database, and the connectivity score—which measures how closely compounds-induced transcript modifications resemble user-input differential gene alterations—was obtained. Additionally, the mechanism of action (MoA) of the compounds was analyzed.

R programming (version 4.0.2, available at https://www.rproject.org) was used for all data calculations and statistical analysis. The Spearman test was utilized for correlation analysis. p < 0.05 was considered statistically significant.

The ASD expression profiles were downloaded from the GEO database, and the mean values of gene expression for each sample were essentially the same after normalization (Supplementary Figure S1). We have estimated the ER stress score for each sample in the dataset using the ssGSEA method. The Wilcox test results highlighted the significant upregulation in the ER stress score between the ASD and normal samples (Figure 1A), which indicates that ER stress has an important role in the development and progression of ASD. Further, the expression of the ER stress genes was evaluated in the two groups, and 37 genes with differential expression were found, including 30 up-regulated genes and 7 down-regulated genes (Figure 1B). Figure 1C displays the chromosomal locations of the ER stress genes that are differentially expressed.

To investigate differentially expressed ER stress regulators’ contribution to the diagnosis of ASD, a series of bioinformatic algorithms were employed. We used a random forest to filter the variables, which showed the lowest average error rate at 99 trees (Figure 2A), and then we analyzed the importance of the variables by calculating the Gini index of a mean decrease (Figure 2B). Seven variables with the mean decrease Gini index greater than 4 were selected for subsequent analysis, and they were UBE4B, PIK3R1, MAP3K5, DNAJB14, USP13, P4HB, TMEM117, respectively. The above-mentioned seven genes were used as input variables into the artificial neural network, a five-neuron layer that extracts and fuses features from the input layer was used to determine the final binary classification result using well-trained weight coefficients and bias (Figure 2C). ROC analysis was implemented to determine the value of the classifier in diagnosing ASD, the results revealed that this classifier was effective at differentiating between healthy people and ASD patients (Figure 2D). Additionally, we verified its effectiveness and dependability in two independent validation set GSE77103 and GSE38322, both of which had AUC value of 0.75 (Figure 2E) and 0.759 (Figure 2F), respectively.

Given the previous analysis suggested that ER stress may be involved in the development of ASD. We further analyzed the co-expression genes associated with ERS scores in patients with ASD. A dendrogram of ASD patients in GSE11176 with ER stress score was clustered using the average linkage method and Pearson’s correlation method (Figure 3A). To obtain a signed network fulfilling the scale-free topology, the soft-thresholding power parameter was set to 7 (at this point scale-free R2 > 0.85) (Figure 3B). A total of 12 modules that have sizes between 65 and 774 genes were identified in the average hierarchical clustering and dynamic tree clipping (Figure 3C). To find the most correlated modules with ERS scores in ASD, we performed spearman correlation analysis (Figure 3D), and the results showed that the turquoise module appeared to have the highest association with ER stress scores (cor = −0.42, p = 2e-06), while magenta (cor = 0.26, p = 0.005), pink (cor = 0.3, p = 9e-04), brown (cor = 0.29, p = 0.002) modules showed a significant positive correlation with ERS scores (Figure 3E). We therefore selected the turquoise module for further analysis, which included 774 genes. A significant correlation existed between the module membership (MM) and gene significance (GS) of the turquoise modules (Figure 3F).

The overlapped results of the turquoise module and differentially expressed ER stress gene was considered to be the hub ER stress genes of ASD. As a result, five genes were obtained, namely SEC16A, EXTL3, DNAJB14, SGTA and UFM1 (Figure 4A). Further, they were subjected to the NetworkAnalyst web to collect the miRNA-gene and TF-gene interactions. The miRNA-gene network consists of 12 nodes and 14 interaction pairs, and DNAJB14 was regulated by 4 miRNAs in this network (Figure 4B). In addition, the TF-gene network includes 11 nodes and 17 interaction pairs (Figure 4C). The transcription factor SP2 may be closely associated with ASD because it regulated four hub ER stress regulators in this network. These 7 miRNAs and 6 TF regulate more than one hub ER stress gene of the network, which indicates high interaction of the miRNA/TF with hub ER stress genes.

To study the sub-clusters in ASD, the unsupervised consensus cluster analysis was carried out on ASD patients based on the expression of the hub ER stress regulators. There were found to be two distinct ASD modification subclusters (Figures 5A-C), with 45 samples in cluster 2 and 74 samples in cluster 1. Tsne analysis also verified that ASD patients can be categorized into two patterns (Figure 5D). We also looked at the hub ER stress gene expression levels across the different subtypes, and the results showed that UFM1 and DNAJB14 were up-regulated in subcluster 1, while SEC16A and EXTL3 were up-regulated in subcluster 2 (Figures 5E,F), indicating that the various ASD subclusters may have various transcriptome or other characteristics.

We performed differential analysis between ASD subclusters and obtained a total of 167 DEGs, of which 127 genes in ASD cluster 2 and 40 genes in ASD cluster 1 were expressed at lower levels (Figures 6A,B). Further enrichment analysis was carried out using these DEGs. GO enrichment analysis showed that ATP synthesis coupled electron transport, protein targeting to ER and phosphatidylinositol metabolic process were involved in the biological process subsection; Cellular component exhibits significant involvement of mitochondrial related complex in DEGs; Molecular function subsection data indicate the structural constituent of ribosome terms were most enriched (Figure 6C). Moreover, pathways of neurodegeneration-multiple diseases, oxidative phosphorylation, primary immunodeficiency interacted with the most number of genes according to the KEGG pathway database (Figure 6C). We also conducted the GSEA analysis, the results showed that REACTOME TRANSLATION (Figure 7A), REACTOME REGULATION OF EXPRESSION OF SLITS AND ROBOS (Figure 7B), REACTOME SELENOAMINO ACID METABOLISM (Figure 7C) and REACTOME RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY (Figure 7D) pathway were significantly enriched in ASD cluster2, while KEGG LEISHMANIA INFECTION (Figure 7E), (Figure 7F), WP EBOLA VIRUS PATHWAY ON HOST (Figure 7G), KEGG SYSTEMIC LUPUS ERYTHEMATOSUS (Figure 7H) and PID FAS PATHWAY (Figure 7I) were significantly enriched in ASD cluster1. Next, we focused on the master regulator genes significant enriched in ASD clusters. The findings demonstrated that the transition from normal control to ASD cluster 1 was supported by the repression of ATF3, UHRF1, BATF2, and the activation of the BCL6 (Figure 8A). While inhibition of AHR, BZLF1, and CREB1 assisted the transition from normal control to ASD cluster 2 (Figure 8B). Overall, based on DEGs, functional enrichment and the master regulatory analysis, we highlighted two distinct biological phenotypes of ASD, which validated the existence of heterogeneity of ERS processes in ASD.

To examine the variations in the immune microenvironment, we compared the immune infiltrating cells, immune reaction gene sets and HLA gene expression between distinct subtypes. Analysis of the abundance of 22 immune microenvironment infiltrating cells revealed that ASD subtype 2 had a higher amount of plasma cell infiltration while the other 21 immune cells had a similar abundance (Figure 9A). Regarding the immune reaction process, the activity of the synonyms B cell receptor (BCR) signaling pathway and interleukin receptor activity were both increased in ASD subcluster2, but interferon receptor activity was decreased in subcluster2 (Figure 9B). Therefore, we speculated that a more active immune system in ASD subcluster2. In addition, the expression of different HLA was also different between ER stress patterns, HLA-DRA and HLA-E were downregulated in ASD subtype2 (Figure 9C). These results confirmed again that ER stress modification plays a crucial regulatory role in forming different immune microenvironments in ASD patients.

The CMap database was employed to find prospective compounds that can effectively reverse the differential gene expression of different ASD subclusters. Differences between Subcluster 1 and Subcluster 2 were also apparent in compound prediction (Figure 10A). We found a total of 118 compounds that were significantly enriched in the ASD subcluster1, and 18 compounds that were enriched in the ASD cluster2. Among them, the PKC inhibitor BRD-K09991945 that targets Glycogen synthase kinase 3β (GSK3B) might have a therapeutic effect on both clusters, whereas the PI3K inhibitor BGT-226 and Radical formation stimulant temoporfin might worsen the situation of different ASD subclusters, respectively (Supplementary Table S1). Furthermore, we conducted the MoA analysis of the potential compounds. The result of four compounds (fexofenadine, VUF-5681, JNJ-10191584 and amodiaquine) shared the MoA of histamine receptor antagonist, another three compounds (imperatorin, CGP-60474 and BRD-K19136521) shared the MoA of the CDK inhibitor in ASD subcluster1 (Figure 10B). Compounds that target ASD subcluster2 enriched the MoA of Adrenergic receptor antagonists, Cytochrome P450 inhibitor and other substances (Figure 10C).