We devised a comprehensive study flowchart ( Figure 1 ). The gene expression matrix of ASD patients involved in this study comes from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database created by the National Biotechnology Information Center. After screening, three data sets were included, namely GSE6575, GSE26415 and GSE42133. Subsequently, we first filter the genes recorded in the three data sets and correct the background to ensure that the data will not be disturbed. Subsequently, we logarithmically transformed the data and further normalized it. Finally, the three data sets are properly merged and batch corrected by using the Combat method in the “sva” package. GWAS data acquisition for ASD from the Psychiatric Genomic Consortium (https://www.med.unc.edu/pgc/results-and-downloads).

In this study, many algorithms were used to calculate the immune infiltration. Firstly, XCELL software package was used to calculate the gene expression profile of ASD and quantify their internal immunity and relative abundance of stromal cells. Subsequent algorithms for evaluating the scores and relative abundance of immune cells come from EPIC (11), ssGSEA (12), quanTIseq (13), TIMER (14), CIBERSORT (15), MCPCounter, XCELL and ESTIMATE.

We utilized the WGCNA software package to construct the WGCNA network with the aim of identifying gene modules associated with two subclasses of ASD and the clinical features of ASD patients (16). To ascertain the optimal threshold, additional scale-free topological criteria are applied. The cluster tree analysis encompasses over 50 modules, distinguished based on diverse colors. The R software package “Cluster Profiler” is utilized to conduct gene ontology (GO) and Kyoto Gene and Genome Encyclopedia (KEGG) analyses, elucidating the positions and pathways of pivotal genes within the cyan module.

In this study, we used the following R packages (openxlsx, seqinr, plyr, randomForestSRC, glmnet, plsRglm, gbm, caret, mboost, e1071, BART, MASS, snowfall, xgboost, ComplexHeatmap, RColorBrewer, pROC) and harnessed the capabilities of 12 machine learning algorithms. These machines learning included Lasso, Ridge, Enet, Stepglm, SVM, glmBoost, LDA, plsRglm, RandomForest, GBM, XGBoost, NaiveBayes. Importantly, several algorithms, such as RF and Lasso, inherently offer feature selection capabilities. These feature selection capabilities are instrumental in refining our model, allowing it to focus on the most predictive attributes. Therefore, we mainly choose feature variables based on RF or LASSO (17).

In ML framework, the RF model was implemented via the randomForestSRC package. RF had two parameters ntree and rf_nodesize, where ntree represented the number of trees in the forest and rf_nodesize was the number of randomly selected variables for splitting at each node. In this study, we set ntree to 1000 and rf_nodesize to 5 for RF model construction. The Enet, Lasso, and Ridge were implemented via the glmnet package. We utilized the built-in cv.glmnet function from the glmnet package to perform 10-fold cross-validation. This approach was employed to accurately determine the optimal value of the lambda parameter, which is critical for regularization in the Elastic Net model. Additionally, the alpha parameter was dynamically set within the framework to distinguish between the Lasso (alpha = 1) and Ridge (alpha = 0) penalties. The Stepglm was implemented via stats package. A stepwise algorithm using the AIC (Akaike information criterion) was applied, and the direction mode of stepwise search was set to “both”, “backward”, and “forward”, respectively. The SVM was implemented via e1071 package. This approach allows for the construction of a hyperplane or set of hyperplanes in a high-dimensional space, which can be used for classification, regression, or other tasks. Particularly, the SVM implementation in e1071 is versatile and supports various SVM kernels. The GBM was implemented using the gbm package. The gbm function facilitates the construction of the GBM model, which is particularly adept at handling complex nonlinear relationships within high-dimensional datasets. In this analysis, this function was employed to identify the optimal number of trees by minimizing the cross-validation error through an iterative process. The glmboost model was implemented using the mboost package. It employs cross-validation to determine the optimal number of boosting iterations, enhancing model performance while preventing overfitting. The function supports two modes: returning the final model for predictions (“Model” mode) or extracting significant predictors (“Variable” mode), thus facilitating both predictive accuracy and interpretability of influential variables. LDA is implemented in your function using the caret package, emphasizing the use of training data, labels, and cross-validation to improve model performance. The plsRglm was implemented using the plsRglm package. The plsRglm mainly conducts cross-validation with ten folds to optimize model parameters and validate the model’s predictive performance. This approach aims to enhance model accuracy while addressing high-dimensionality and multicollinearity in the dataset. The XGBoost was implemented using the xgboost package. This function performs 5-fold cross-validation to determine the optimal number of boosting rounds by minimizing the test log loss, ensuring a balanced approach between prediction accuracy and model simplicity.The NaiveBayes was implemented using the e1071 package. It constructs the model by combining feature data with the target variable, ensuring the target is treated as a categorical outcome, thus facilitating both predictive modeling and feature analysis. Finally, in ML framework, we employ a novel approach by combining two distinct machine learning algorithms: one for variable selection and the other for model construction. This methodology leverages the strengths of different algorithms to enhance the predictive performance and interpretability of the resulting models.

We utilized GNOVA to evaluate the impact of single nucleotide polymorphisms (SNPs) on the genetic inheritance of ASD and neutrophils. This involved conducting regression analysis on the z-statistics produced from two separate studies on these traits. The studies utilized LD scores precomputed using 1000 Genomes European data (18, 19). Through this methodology, we estimated both the heritability (h2) of ASD attributed to SNPs and the overall genetic correlation (rg) between ASD and neutrophils.

First, we defined the set of SNPs as those located within a particular gene using the VEGAS annotation file. Then, the P values of the SNPs within each gene are weighted averaged to obtain a gene-level P value, and the P values are simultaneously transformed into Z-statistics. Finally, the newly determined Z-statistic was subjected to a pleiotropy test using the PLACO method. PLACO is an innovative method for detecting pleiotropy at the level of SNPs using the concept of composite null hypothesis from high-dimensional mediation analysis (20). Previous simulations and variance-component-based mediation analyses under the composite null hypothesis have suggested the potential use of this method to assess validity at the gene level (21, 22). Consequently, we used it to identify polymorphic associations at the gene level. To mitigate the impact of excessively large effects, SNPs with extreme Z2 (>80) values were excluded. PLACO assumes three sub-null scenarios for each SNP studied using the composite null hypothesis of pleiotropy: (i) H00: The SNP is not associated with either disease. (ii) H01: The SNP has an effect only on the first disease. (iii) H02: The gene has an effect only on the second disease. (is H1: the SNP effect on both diseases, which represents a pleiotropic relationship.

A pairwise cross-trait meta-analysis was performed utilizing Cross Phenotype Association (CPASSOC). CPASSOC integrates effect estimates and standard errors obtained from GWAS summary statistics. This approach aims to examine the hypothesis of association between a single nucleotide polymorphism (SNP) and two traits (23). For this analysis, we utilized the heterogeneous version of cross-phenotype association (SHet). SHet employs a fixed-effect model weighted by sample size, rendering it more potent for detecting heterogeneous effects when they exist among studies (24).

We performed differential expression analysis and gene set enrichment analysis with FUMA for the pleiotropic genes analyzed by PLACO and CPASSOC employing FUMA as the tool (25). expression trait loci (EQTL) analysis will from GTEx and other database were collected for 53 different tissues, with a final consideration of 22,146 genes. The resulting gene sets with adjusted P ≤ 0.05 were reported as significant findings.

In order to further study the expression level of the selected target genes in ASD and normal people, we collected new blood samples (10ml, EDTA anticoagulation) from ASD children and normal children, with 5 people in each group. Subsequently, the collected blood was stored at 4 degrees, and serum was obtained by low-speed centrifugation. The people included in this study are volunteers, who voluntarily provide blood samples after being fully explained, and pass the examination and approval of the ethics Committee of the hospital.

In this research, we included 5 ASD and 5 healthy individuals from Xiangyang No.1 People’s Hospital for qPCR analysis. The gathered blood was subjected to centrifugation at a speed of 3000 revolutions per minute for a duration of 10 minutes. Subsequently, the upper portion of the resulting serum was preserved for mRNA extraction. The Trizol reagent (Life Technologies, USA) was employed to extract the total RNA from the serum. Afterwards, cDNA synthesis was conducted using the RevertAid First Strand cDNA synthesis kit (Fermentas, Canada) with the isolated total RNA. PCR amplification was executed using the QuantiTect SYBR Green PCR kit (Qiagen, Inc) utilizing the ABI Prism 7000 sequence detection system (Applied Biosystem, CA, USA). The cycling conditions were as follows: an initial denaturation step at 95°C for 10 minutes, followed by 40 cycles of denaturation at 95°C for 15 seconds and annealing/extension at 60°C for 1 minute. The mRNA expression level was determined employing the 2^-△△Ct method, and subsequently normalized to the expression level of β-actin. All experiments were repeated five times, and the primer sequences are listed in Table 1 .

R language (version 4.2.0) is used for statistical analysis. Wilcoxon test was used for comparison between groups. P-value < 0.05 were considered as potential associations.

The flowchart systematically described our study ( Figure 1 ). The essence of using the ESTIMATE algorithm to determine the subgroup characteristics of ASD is to achieve the goal by scoring immunity and cell matrix. Interestingly, there was a significant difference in the immune score between the two groups, while the ASD group showed a higher matrix score ( Figure 2A ). In order to further clarify immune infiltration and understand the immune characterization, we quantified the abundance of 22 kinds of immune cells in the microenvironment ( Figure 2B ). Several immune checkpoint genes found in the ASD group showed higher levels, suggesting their potential as immunotherapy targets in the future, including TLR9, TNFRSF4, CTLA4, TNFRSF9, PDCD1 and PDCD1LG2 ( Figure 2C ). In addition, ASD showed different expression levels in six immune cell subgroups, including Plasma cells, CD8 cells, activated CD4 cells, resting NK cells, Macrophages M0 and Macrophages M1. It is worth noting that, unlike the significantly increased gene expression, macrophages of immune cell subgroup M1 showed significantly high infiltration, but M2 showed significantly low infiltration risk. Notably, unlike the significantly increased gene expression, macrophages of immune cell subgroup M1 showed significantly high infiltration, but M2 showed significantly low infiltration risk.

We further applied WGCNA depth analysis to the merged data set to extract the clinical phenotype related to ASD and immune cells. When the threshold is 4, the scale-free network and connectivity show the greatest compatibility ( Figure 3A ). Further clustering tree algorithm divides ASD-related genes into six gene modules, and each module is represented by a special color ( Figure 3B ). Among them, cyan modules contain the most genes, and have the most significant correlation with the poor prognosis of ASD (R=0.49). Specifically, CD8 cells (R=0.97), endocrine cells (R = 0.93) and fibers (r = 0.89) showed a strong positive correlation ( Figure 3C ). Therefore, the cyan module was selected as the hub module, and the hub gene was extracted from it by using the selection criteria cor.MM>0.7 and cor.GS>0.4.

In order to further clarify the diagnostic and therapeutic effects of feature genes on ASD, we combined common machine learning algorithms in pairs and evaluated the diagnostic ability of each feature gene in predicting ASD progression in the internal dataset through calibration curve and subject operating characteristic (ROC) curve analysis. The AUC values of pairwise combination in the training dataset ranged from 0.649 to 0.822, and we chose the most appropriate RF+Ridge (AUC=0.811) as the final model ( Figure 4 ; Supplementary Table 1 ). Subsequently, we compared the expression of key genes in wild-type and mutant ASD, and tracked the changes in key genes according to the different ages of the patients. The single-cell sequencing results indicate that the transcriptome fragments of the 47 core genes exhibit a consistent trend ( Figure 5 ).

Genome-wide single nucleotide polymorphism (SNP)-based heritability (eri was estimated to for ASD and Neutrophils. Then, we found a significant positive genetic correlation between them ( r ^ g = 0.123, P = 0.0003) by cross-trait genetic correlation analysis. This link implies a possible common genetic cause between them. Hence, it is necessary to further investigate the genetic mechanism.

Subsequently, pleiotropy analysis we sued by PLACO and CPASSOC. Therefore, PLACO analysis identified 3187 statistically significant regions (P < 5E-08) ( Supplementary Table S2 ), with ASD and Neutrophils. Similarly, we found that 13501 genes were statistically significant (P < 5E-08) ( Supplementary Table S3 ), with ASD and Neutrophils. In the end, we found that PLACO and CPASSOC identified a total of 16 potential pleiotropic regions ( Supplementary Table 4 ).

First, we performed gene enrichment analysis using FUMA for the 16 pleiotropic regions identified. And the EQTL results indicated that the differentially expressed genes were predominantly enriched in tissues such as pancreas, liver, heart, blood, brain, and muscle based on the analysis of expression levels ( Supplementary Table 5 ).

The comparison of serum mRNA expression levels between patients with ASD and the normal population shows that, β-actin is basically similar, while the mRNA levels of IGF2R, UBN1 and CFLAR in plasma of children with ASD showed a significant downward trend, and the difference was significant (P<0.05), suggesting that the surge of transcriptome had an impact on the development of nervous system ( Figure 6 ).