Construction of Autophagy-Related Gene Classifier for Early Diagnosis, Prognosis and Predicting Immune Microenvironment Features in Sepsis by Machine Learning Algorithms

Figures & data

Table 1 Dataset Included in the Study

Accession	Cohort Description	Timing of Gene Expression Profiling	Country	Normal/Control Sample	Mortality/Sepsis Sample
GSE54514	Sepsis	Day1/2/3/4/5 of ICU admission	Australia	36	31/127
GSE57065	Septic shock	0h/24h/48h of ICU admission	France	25	–/82
GSE63042	Sepsis/SIRS	Day of enrollment upon presentation to the ED	United States	23	28/106
GSE65682	Sepsis due to CAP and HAP +AS	On ICU admission	Netherlands and UK	42	48/231
GSE69063	Sepsis	0h/1h/3h of ICU admission	Australia	33	–/57
GSE69528	Sepsis due to CAP	On ICU admission	USA	55	-/83
GSE95233	Septic shock	Day1/3 of ICU admission	France	22	34/102
GSE106878	Septic shock	0h/24h of ICU admission	Germany	–	26/94
GSE131761	Septic shock	On ICU admission	Spain	15	–/114
GSE154918	Septic shock	Day 1 of ICU admission	Germany	40	–/53
E-MTAB-4421	Septic shock	On ICU admission	UK	–	56/265
E-MTAB-4451	Sepsis due to CAP	On ICU admission	UK	–	57/114
E-MTAB-7581	Septic shock	At enrollment	UK	–	48/176

Abbreviations: CAP, community-acquired pneumonia; HAP, hospital-acquired pneumonia; AS, abdominal sepsis; SIRS, systemic inflammatory response syndrome; ICU, intensive care medicine; ED, emergency room.

Figure 1 The illustrations for this study. The overall protocol utilized in the current study to comprehensively characterize of costimulatory molecule gene for diagnosis, prognosis and recognition of immune microenvironment features in sepsis.

Abbreviations: ARGs, autophagy-related genes; DEARGs, differentially expressed autophagy-related genes.

Figure 2 Identification of differentially expressed ARGs in sepsis and Functional enrichment analysis. (A) An UpSet diagram exhibited the interactions result of GEO datasets difference analysis and ARG cohorts. (B–I) Principal component analysis for the expression profiles of DEARGs to distinguish sepsis patients from healthy control/uninfected control/control patients in multi-transcriptome cohorts. (B) GSE54514 datasets; (C) GSE57065 datasets; (D) GSE65682 datasets; (E) GSE69063 datasets; (F) GSE69528 datasets; (G) GSE95233 datasets; (H) GSE131761 datasets; (I) GSE154918 datasets. (J and K) Bubble plots illustrating functional enrichment analysis using DEARGs. (J) The top 15 significant terms of Gene Ontology (GO), including the molecular function (MF), cellular component (CC), and biological process (BP). (K) The top 15 significant pathways in Kyoto Encyclopedia of Genes and Genomes (KEGG) enriched.

Figure 3 Identification of prognosis-related ARGs (sepsis-specific ARGs) and construction of a ARG classifier based on machine learning algorithms (including Lasso, SVM, ANN). (A–D) Modified Lasso was used to identify candidate ARGs with 10-fold cross-validation in multi-transcriptome datasets. The Y-axis shows mean-square error and the X-axis is Log (λ). Dotted vertical lines represent minimum and 1 standard error values of λ. The genes selected at minimum standard error values of λ were finally used for further analysis. (A) GSE54514 datasets (N=20); (B) GSE63042 datasets (N=18); (C) GSE65682 datasets (N=12); (D) GSE95233 datasets (N=24). (E–H) Modified Lasso was used to identify candidate ARGs with 10-fold cross-validation in multi-transcriptome datasets. SVM algorithm was applied to screen candidate ARGs. The red dots indicated the lowest error rate and the highest precision when genes are this number. (E) GSE54514 datasets (N=25); (F) GSE63042 datasets (N=22); (G) GSE65682 datasets (N=40); (H) GSE95233 datasets (N=16). (I) An UpSet diagram exhibited the interactions result of Lasso and SVM in multi-transcriptome datasets. (J) The visualization of artificial neural network (ANN). The neural network contains 8 input layers, 4 hidden layers, and 2 output layers.

Abbreviatons: I, input layers; H, hidden layers; O, output layers.

Figure 4 The diagnostic efficacy of ARG classifier in multi-transcriptome datasets. (A–H) ROC curves analyzed the diagnostic accuracy of ARG classifier in multiple datasets and time points. (A) GSE54514 datasets; (B) GSE57065 datasets; (C) GSE65682 datasets; (D) GSE69063 datasets; (E) GSE69528 datasets; (F) GSE95233 datasets; (G) GSE131761 datasets; (H) GSE154918 datasets. (I and J) Comparison of the ARG scores between sepsis and SIRS in GSE63042 datasets (I), uncomplicated Infection in GSE154918 datasets (J). p value was calculated using the Wilcoxon test. (K) Comparison of the ARG scores between sepsis, healthy control and anaphylaxis in GSE69063 datasets. p value was calculated using the Kruskal–Wallis test. (L) Comparison of the ARG scores between sepsis, healthy control and head trauma in GSE69063 datasets. p value was calculated using the Kruskal–Wallis test. (M–P) Principal component analysis based on hub ARGs to distinguish sepsis from SIRS (M), uninfected control (N), anaphylaxis (O), and head trauma (P).

Abbreviation: SIRS, systemic inflammatory response syndrome.

Table 2 Univariable and Multivariable Logistic Regression Analysis for Prediction of Survival in GEO and ArrayExpress Databases

Dataset	Factors	Subgroup	Univariable Analysis		Multivariable Analysis
GSE54514			OR (95% CI)	P	OR (95% CI)	P
	Age		1.02 (1.02–1.08)	0.004*	1.04 (0.99–1.09)	0.057
	APACHE II		1.17 (1.08–1.26)	< 0.001*	1.20 (1.07–1.34)	0.001*
	ARG score		54.71 (7.15–418.96)	< 0.001*	79.67 (9.16–693.03)	< 0.001*
GSE63042
	ARG score		5.56 (2.02–15.27)	0.001*	5.56 (2.02–15.27)	0.001*
GSE65682
	Age		1.02 (0.99–1.04)	0.114	NA	NA
	Sex	Female	1
		Male	1.42 (0.75–2.67)	0.284	NA	NA
	MARS	1–2	1
		3–4	0.93 (0.67–1.31)	0.688	NA	NA
	DM	No	1
		Yes	1.63 (0.71–3.74)	0.248	NA	NA
	ARG score		4.53 (2.17–9.43)	< 0.001*	4.40 (2.10–9.20)	< 0.001*
GSE95233
	Age		0.99 (0.95–1.03)	0.657	NA	NA
	Sex	Female	1
		Male	0.65 (0.35–1.18)	0.156	NA	NA
	ARG score		12.13 (4.13–35.59)	< 0.001*	12.14 (4.13–35.69)	< 0.001*
GSE106878
	Age		1.05 (1.01–1.09)	0.019*	1.04 (0.99–1.08)	0.107
	ARG score		5.06 (1.80–14.23)	0.002*	6.54 (3.06–13.97)	0.009*
E-MTAB-4421
	Age		1.05 (1.02–1.07)	< 0.001*	1.05 (1.02–1.08)	< 0.001*
	Sex	Female	1
		Male	1.16 (0.64–2.10)	0.620	NA	NA
	SRS	1
		2	1.77 (0.98–3.20)	0.060	NA	NA
	ARG score		6.45 (3.09–13.48)	< 0.001*	6.54 (3.06–13.97)	< 0.001*
E-MTAB-4451
	Age		1.00 (0.98–1.03)	0.934	NA	NA
	Sex	Female	1
		Male	0.95 (0.40–2.28)	0.913	NA	NA
	SRS	1
		2	2.70 (1.18–6.19)	0.019	1.35 (0.50–3.70)	0.555
	ARG score		4.90 (2.15–11.15)	< 0.001*	4.35 (1.87–10.15)	0.001*
E-MTAB-7581
	Age		1.02 (0.99–1.04)	0.091	NA	NA
	Sex	Female	1
		Male	1.39 (0.71–2.72)	0.343	NA	NA
	SRS	1
		2	1.31 (0.68–2.55)	0.423	NA	NA
	APACHE II		1.08 (1.04–1.13)	0.001*	1.09 (1.04–1.15)	0.001*
	ARG score		5.08 (2.37–10.86)	< 0.001*	5.01 (2.28–11.14)	< 0.001*

Notes: These variables were eliminated in the multivariate logistic regression model, so the HR and P values were not available. *P < 0.05.

Abbreviations: OR, odds ratio; CI, confidence intervals; MARS, the Molecular Diagnosis and Risk Stratification of Sepsis; DM, diabetes mellitus; NA, not available.

Figure 5 The prognostic efficacy of ARG classifier in multi-transcriptome datasets. (A) An Lollipop chart showed the AUC values of prognostic prediction in multiple datasets. (B and K) Correlation between ARG score and APACHE-II in GSE54514 datasets (B), and E-MTAB-7581 datasets (K). Correlation coefficient and p value were calculated by Spearman correlation analysis. (C, D and F) The distribution of mortality rate between low ARG and high ARG subgroups in GSE54514 datasets (C), GSE63042 datasets (D), and GSE95233 datasets (F). p value was calculated using the Chi-square test. (E) Comparison of the ARG scores between uncomplicated sepsis, severe sepsis and septic shock in GSE63042 datasets. p value was calculated using the Kruskal–Wallis test. (G, I and J) The Sankey diagram depicts the flow from the two subgroups (low ARG and high ARG) to different phenotypes and survival outcomes, in which the width of the flow rate is proportional to the patient number. (G) E-MTAB-4421 datasets; (I) E-MTAB-4451 datasets; (J) GSE65682 datasets. p value was calculated using the Chi-square test. (H and L) The heatmap depicted the distribution of clinical characteristics arranged by the increasing ARG score in GSE106878 datasets (H), and E-MTAB-7581 datasets (L). Comparison of the different clinical characteristics between low ARG and high ARG subgroups. The asterisks indicate a statistically significant p-value calculated using the Chi-square test (*p < 0.05; **p < 0.01; ***p < 0.001).

Abbreviation: APACHE-II, Acute Physiology and Chronic Health Evaluation II.

Figure 6 Clinical usefulness of ARG classifier. (A–H) Multiple ROC analysis was performed to compare the prognostic performance of the ARG classifier against clinical features in multiple cohorts. (A) GSE54514 datasets; (B) GSE63042 datasets; (C) GSE65682 datasets; (D) GSE95233 datasets; (E) GSEGSE106878 datasets; (F) E-MTAB-4421 datasets; (G) E-MTAB-4451 datasets; (H) E-MTAB-7581 datasets. (I–P) Decision curve analysis was applied to evaluate the clinical usefulness of ARG classifier against clinical features in multiple cohorts. The Y-axis represents the net benefit. The black line represents the hypothesis that no patients die. The X-axis represents the threshold probability. The threshold probability is where the expected benefit of treatment is equal to the expected benefit of avoiding treatment. (I) GSE54514 datasets; (J) GSE63042 datasets; (K) GSE65682 datasets; (L) GSE95233 datasets; (M) GSEGSE106878 datasets; (N) E-MTAB-4421 datasets; (O) E-MTAB-4451 datasets; (P) E-MTAB-7581 datasets.

Figure 7 Analyzing the correlation between ARG classifier/hub ARGs and infiltrating immune cells in GEO datasets. (A) Comparison of infiltrating immune cells between ARG different subgroups based on CIBERSORTx tool. p value was calculated using the Wilcoxon test. (B) Comparison of infiltrating immune cells between ARG different subgroups based on ssGSEA algorithms. p value was calculated using the Wilcoxon test. (C) Correlation between ARG classifier/hub ARGs and immune cells. Correlation coefficient and p value were calculated by Spearman correlation analysis.

Figure 8 Analyzing the correlation between ARG classifier/hub ARGs and infiltrating immune cells in ArrayExpress datasets. (A) Comparison of infiltrating immune cells between ARG different subgroups based on CIBERSORTx tool. p value was calculated using the Wilcoxon test. (B) Comparison of infiltrating immune cells between ARG different subgroups based on ssGSEA algorithms. p value was calculated using the Wilcoxon test. (C) Correlation between ARG classifier/hub ARGs and immune cells. Correlation coefficient and p value were calculated by Spearman correlation analysis.

Figure 9 Analyzing the correlation between ARG classifier/hub ARGs and biological pathways in GEO and ArrayExpress datasets based on Gene set variation analysis (GSVA). (A and C) Comparison of biological pathways between ARG different subgroups in GEO datasets (A), and ArrayExpress datasets (C). (B and D) Correlation between ARG classifier/hub ARGs and biological pathways in GEO datasets (B), and ArrayExpress datasets (D). Correlation coefficient and p value were calculated by Spearman correlation analysis.

Figure 10 Analyzing the correlation between ARG classifier/hub ARGs and cytokines expression levels in GEO and ArrayExpress datasets. (A and C) Differential expression of cytokines between ARG different subgroups in GEO datasets (A), and ArrayExpress datasets (C). The asterisks indicate a significant statistical p value calculated using the Wilcoxon test (*p < 0.05; **p < 0.01; ***p < 0.001). (B and D) Correlation between ARG classifier/hub ARGs and cytokines in GEO datasets (B), and ArrayExpress datasets (D). Correlation coefficient and p value were calculated by Spearman correlation analysis.

Figure 11 Correlation between Hub ARGs, immune-related genes and HLA-related genes and Preliminary experimental validation. (A and B) Correlation between Hub ARGs, immune-related genes and HLA-related genes in GEO datasets (A), and ArrayExpress datasets (B). Correlation coefficient and p value were calculated by Spearman correlation analysis (*p < 0.05; **p < 0.01; ***p < 0.001). (C) Cytoscape software depicted visually the gene-gene network regulation relation. Pink lines represent a significant positive correlation. Green lines represent a significant negative correlation. (D) Differential expression of hub ARGs between sepsis patients and control patients by qRT-PCR. The asterisks indicate a significant statistical p value calculated using the Wilcoxon test (*p < 0.05; **p < 0.01; ***p < 0.001). (E) ROC curves compared the diagnostic efficacy of ARG classifier, PCT and CRP based on clinical specimens.