The gene expression profiles of 63 septic shock blood samples and 337 normal whole blood samples were downloaded from the Gene Expression Omnibus (GEO¹) (Braga et al., 2019) and the Genotype-Tissue Expression Portal (GTEx, 2015²), respectively. The gene expression profiles were retrieved in formats of HTseq-counts and were normalized to Transcripts Per Million (TPM) based on gene length. Clinical data including demographics (age and gender) and blood collection time points (T1: within 16 h of ICU admission; T2: 48 h after study enrollment; T3: on day 7 from ICU admission or before discharge from the ICU) were also extracted from the database.

Samples and patients with incomplete clinical information were excluded, and conformers who met the inclusion and exclusion criteria were extracted for further analysis. In order to reduce the error caused by the different length and sequencing depth of each single gene in the sequencing process, firstly, raw counts (RCs) of RNA-seq data obtained from GTEx and GSE131411 were quantified, respectively, which were then standardized using voom function in limma (Linear Models for Microarray Analysis) package. To eliminate the batch effect, these two batches of RNA-seq data were corrected using the normalizeBetweenArrays function in the limma package and then merged for differential expression analysis (Smyth, 2004).

Differential expression genes (DEGs) between 63 septic shock blood samples and 337 normal whole blood samples were determined using the Edge R package (Smyth, 2004). The standard to define the DEGs was | log₂ Fold Change (FC)| > 1.0 along with False Discovery Rate (FDR) < 0.05. Besides, the Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) items of enrichment analysis were performed to explore the biological processes (BPs), cellular components (CCs), molecular functions (MFs), and KEGG pathways that septic shock-specific DEGs enriched.

A total of 2,498 immune-related genes (IRGs) were retrieved from the ImmPort database³ and Molecular Signatures Database (MSigDB) v7.1,⁴ respectively (Liberzon et al., 2015). The intersection of the 2,498 IRGs and DEGs was filtered by further non-parametric test (Kruskal–Wallis H test) to identify SSSIRGs.

The “CIBERSORT, version x” algorithm was applied to estimate the fraction of 22 types of immune cell in each whole blood sample (Newman et al., 2019). The non-parametric test (Mann–Whitney U test) was conducted to illuminate the different immune cell components between septic shock blood samples and normal whole blood samples. In addition, the fraction of 29 immunomodulating gene sets was quantified using single-sample Gene Set Enrichment Analysis (ssGSEA) (Barbie et al., 2009).

Correlation analysis was applied to SSSIRGs, 22 immune cells, and 29 immune gene sets. The pairwise interactions between SSSIRGs and immune cells (| correlation coefficient| > 0.500, p-value < 0.001) along with SSSIRGs and immune gene sets (| correlation coefficient| > 0.600, p < 0.001) were identified as immune components potentially regulated by SSSIRGs.

Fifty hallmark gene sets (signaling pathways) and the 318 TFs were retrieved from the MSigDB v7.1 (see text footnote 4) and Cistrome database,⁵ respectively (Liberzon et al., 2015; Zheng et al., 2019). Differential expression and correlation analyses were conducted to determine upstream TFs, while hallmark pathways of SSSIRGs were absolutely quantified by Gene Set Variation Analysis (GSVA) (Hanzelmann et al., 2013). Next, the correlation analysis was applied to key SSSIRGs, TFs, and aforementioned signaling pathways. Pairwise interactions between TFs and SSSIRGs (| correlation coefficient| > 0.900, p < 0.001), SSSIRGs, and the 50 hallmark pathways (| correlation coefficient| > 0.600, p < 0.001) were extracted for further analysis. All above significant interaction pairs were integrated to construct the regulation network among SSSIRGs and their upstream TFs and downstream pathways/immune cells/immune gene sets.

In order to better understand and validate the relationship between SSSIRGs and prognosis of patients with septic shock, patients were categorized into three groups (T1–T3) based on blood sample collection time points (T1: within 16 h of ICU admission; T2: 48 h after study enrollment; T3: on day 7 from ICU admission or before discharge from the ICU). Further differential expression analysis based on different time points was carried out, respectively, using the edge R package.

In addition, gene sets over-representation analysis (GSORA) was performed to evaluate the fraction of genes in tested gene clusters (Pomyen et al., 2015). In this study, 54 septic shock-related gene sets downloaded from MSigDB were divided in nine subclusters on the basis of function features, which included C1–C8 and the hallmark pathway cluster (H) (Liberzon et al., 2015). In addition, ORA was carried out to determine the enrichment level of these gene sets.⁶

Based on the key SSSIRGs obtained previously and different lengths of hospital stay (LOS) of septic shock patients, non-parametric test and clinical correlation analysis were performed, and the intersection of the results from these two analyses which were statistically significant was extracted for further analysis. Specifically, LOS served as an important sign in clinical judgment of the severity of septic shock.

Here, we applied the Connectivity Map (CMap, build 02) to find potential inhibitors that may target key SSSIRGs. In total, 6,100 gene expression cases covering 1,309 drugs were obtained from the CMap database⁷ (Lamb et al., 2006), that is, a potential drug might correspond to various gene expression cases. Genes in each case were ranked via taking differential expression values between drug-untreated and drug-treated cell lines, and 6,100 gene lists associated with drugs were then generated. On the basis of identified key SSSIRGs involved in septic shock and 6,100 drug-related cases, we conducted a non-parametric test to explore the relationship between drugs and septic shock.

Information of targeting compounds is available in the mechanism of actions (MoAs)⁸ that includes transcriptional responses of various human cell lines to perturbagens, structural formulas, and various targets. On the basis of MoA, compounds which may target to septic shock-related SSSIRGs in this study were extracted.

Chromatin Immunoprecipitation sequencing (ChIP-seq) and Assay for Targeting Accessible-Chromatin with high-throughput sequencing (ATAC-seq) data were used to validated the regulation mechanism of the network. ChIP-seq and ATAC-seq data were obtained from Cistrome database (see text footnote 5) and GSE139099,⁹ respectively (Zheng et al., 2019; Grubert et al., 2020). The WashU Epigenome Browser and Gviz package were used to visualize the binding peaks (Hahne and Ivanek, 2016; Li et al., 2019).

Many studies supported that COVID-19 severe forms share clinical and laboratory aspects with various pathologies such as hemophagocytic lymphohistiocytosis, sepsis, or cytokine release syndrome. Fatal COVID-19 patients, as would be expected with septic shock patients, produce inflammatory reaction and cytokine storm (Bost et al., 2021). Therefore, single-cell RNA sequencing (scRNA-seq) data of peripheral blood mononuclear cells (PBMCs) from fatal COVID-19 patients were obtained from GSE150728 to determine the specific cellular location of key SSSIRGs and cell cycle regulation in this study, which can further explain the molecular pathogenesis of septic shock spatiotemporally.

The preprocessing of scRNA-seq data of seven COVID patients and six healthy samples was conducted using 10x Genomics Chromium¹⁰ (Chen et al., 2020). After demultiplexing, the sequencing results were divided into two paired-end reads fastq files that were then trimmed to eliminate the template switch oligo (TSO) sequence and polyA tail sequence. In addition, clean reads were aligned with the GRCh38 (Version: 100) genome assembly, which were quantified using the Cell Ranger Software (Version 1.0.0).¹¹

The quantitative gene expression matrices (the row names of matrices were genes and col names were barcodes) were analyzed with Seurat pipeline (Version: 3.2.2) for further analysis (Butler et al., 2018). Only cells with less than 10% mitochondrial gene mapped and expressing more than 100,000 transcripts were extracted for further analysis. In addition, genes expressed in over three single cells were included in follow-up analyses. After the completion of quality control (QC), all samples were integrated into one Seurat object with the function of “IntegrateData,” which were then scaled and normalized with the function of “ScaleData.” The “vst” method was utilized to determine the top 1,500 variable genes. To reduce model dimensionality, principal component analysis (PCA) was initially conducted, and the top 20 PCs were incorporated as input file for Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) analysis. UMAP plots illustrating cell subclusters were constructed using the “DimPlot” and “RunUMAP” function.

Genes with significantly differential expression patterns from the top 1,500 variable genes were defined as DEGs with the “wilcox” method using the “FindAllMarkers” function.

To identify the cell type of each unsupervised cluster, DEGs in all subclusters were applied as potential references, which were combined with known specific cell surface biomarkers obtained from CellMarker¹² (Zhang et al., 2019) for a comprehensive annotation of cell type. Considering the variable gene expression patterns, a specific cellular annotation method was utilized in this study. In addition, the cellular feature plots, dot plots, violin plots, and heat maps were constructed to show the marker genes of each cell type using SCANPY (Version: 1.7.1) and the Seurat R package (Version: 3.2.2) under the environment of Python 3.6 (Butler et al., 2018; Wolf et al., 2018).

To elucidate the key cellular communication patterns and ligand-receptor pairs among various different cell types in peripheral blood, cellular communication analysis was conducted using the iTALK R package (Version: 0.1.0)¹³ (Wang et al., 2019b). Firstly, because of the scarcity of musculus resources in present mainstream cellular communication algorithms, top 1,500 variable genes were transformed into human genes using the biomaRt package (Version: 2.46.0) in a certain homologous degree with the “getLDS” function (Durinck et al., 2009). Secondly, the normalized expression matrix of these genes was incorporated into the iTALK object with the “rawParse” function. Finally, top 200 ligand-receptor pairs were shown by ligand-receptor plots and iTALK networks.

In this study, all statistics analysis processes were performed by the R software (version 3.6.3¹⁴). In descriptive statistics, mean ± SD was utilized for continuous variables in normal distribution. In addition, when encountering continuous variables in abnormal distribution, the median was utilized. A two-sided p < 0.05 was suggested to have statistical significance.

The flowchart of the present study is illustrated in Figure 1. A total of 14,843 genes (9,402 down-regulated genes and 5,441 up-regulated genes) were identified as DEGs between 63 septic shock blood samples and 337 normal whole blood samples (Figures 2A,B). The most important GO items of BPs, CCs, and MFs for DEGs were regulation of GTPase activity, actin cytoskeleton, and nucleoside-triphosphatase regulator activity (Figure 2C). The most significant KEGG pathway where DEGs mainly enriched was MAPK signaling pathway (Figure 2D).

A total of 839 SSSIRGs (666 down-regulated genes and 173 up-regulated genes) were extracted from the intersection of 14,843 DEGs and 2,498 IRGs in septic shock blood, and the expression level of these SSSIRGs were shown in the heat map (Figure 3A) and volcano plot (Figure 3B).

The CIBERSORT algorithm was used to estimate the fraction of 22 types of immune cell in each blood sample (Figure 3C). The results of the Mann–Whitney U test suggested that B cells naïve (p < 0.001), B cells memory (p < 0.050), Plasma cells (p < 0.010), T cells CD8 (p < 0.001), T cells CD4 naïve (p < 0.001), T cells CD4 memory activated (p < 0.001), T cells regulatory (Tregs) (p < 0.001), T cells gamma delta (p < 0.001), NK cells resting (p < 0.001), NK cells activated (p < 0.010), Monocytes (p < 0.001), Mast cells activated (p < 0.001), Dendritic cells activated (p < 0.001), Mast cells resting (p < 0.001), Macrophages M2 (p < 0.001), and Neutrophils (p < 0.010) were differentially enriched between the septic shock and normal whole blood samples (Figure 3D).

Furthermore, heat map and volcano plot were used to illustrate the DEG analysis results of the TFs (Figures 4A,B). Absolute quantification of hallmark pathways was illustrated by the heat map as well as volcano plot (Figures 4C,D). In addition, the t-score of each signaling pathway was shown in the bar plot (Figure 4E). Further, expression levels of 29 immune gene sets in normal whole blood and septic shock blood acquired by ssGSEA were displayed in the heat map (Figure 4F). Only 93 of 839 IRGs were not only differentially expressed between septic shock and normal whole blood samples but also significantly associated with blood collection time points (T1–T3). Thus, they were defined as the SSSIRGs (Figure 4G).

Finally, three SSSIRGs Toll-like receptor 8 [TLR8, protein phosphatase 3, catalytic subunit, alpha isozyme (PPP3CA), and Kirsten rat sarcoma viral oncogene homolog (KRAS)] were screened by the co-expression analysis with 28 co-expression interaction pairs including four TF-SSSIRG interaction pairs, nine SSSIRG-immune-cell interaction pairs, four SSSIRG-immune-gene-set interaction pairs, and 11 SSSIRG-pathway interaction pairs (Figure 5A), based on which a regulation network was constructed. Besides, all correlation patterns of SSSIRGs with both the upstream TFs and the downstream signaling pathways/immune cells/immune gene sets were illustrated by a heat map (Figure 5B).

Importantly, vascular endothelial zinc finger 1 (VEZF1) may upregulate the expression level of PPP3CA (R = 0.530, p < 0.001), which may further increase the activities of T cells gamma delta (R = 0.390, p < 0.001) and HLA (R = 0.620, p < 0.001) while suppressing the activity of hallmark TNFα signaling via NFκB (R = −0.580, p < 0.001); VEZF1 may enhance the expression of KRAS (R = 0.480, p < 0.001), which may further activate the activities of hallmark late estrogen response (R = 0.590, p < 0.001), T cells gamma delta (R = 0.370, p < 0.001), and HLA (R = 0.580, p < 0.001); VEZF1 may stimulate the expression of TLR8 (R = 0.480, p < 0.001), which may inhibit the activities of hallmark TNFα signaling via NFκB (R = −0.510, p < 0.001) and macrophages (R = −0.490, p < 0.001) while increasing the activity of T cells gamma delta (R = 0.210, p < 0.001).

Significantly, gene differential expression patterns between normal samples and septic shock patients in different time point groups (T1–T3) are illustrated by the heat maps and volcano plots, respectively (Figures 6A–C). Importantly, the key SSSIRGs identified previously (TLR8, KRAS, and PPP3CA) were all significantly differentially expressed in the analysis of three time point groups, which was consistent with the previous results. Additionally, the ORA suggested that the DEGs were belonging to several gene sets of immune response and inflammation, including immune response mediated by circulating immunoglobulin, humoral immune response gene set, and hallmark coagulation (Figures 6A–C).

The Venn plot showed there were 12,302 DEGs that were significantly differential expressed between normal samples and all time point groups. The box plots of non-parametric test (Kruskal–Wallis H test) showed that the three identified key SSSIRGs (TLR8, KRAS, and PPP3CA) were significantly related to the length of hospital stay (LOS) (p < 0.050) (Figure 6D). Besides, the CMap analysis suggested that 15 bioactive small molecules including alpha-estradiol, camptothecin, daunorubicin, ebselen, hexestrol, irinotecan, lysergol, menadione, podophyllotoxin, sanguinarine, streptozocin, thioguanosine, trichostatin A, viomycin, and vorinostat were identified as significant inhibitors (p < 0.001). Importantly, podophyllotoxin with the highest specificity was considered as the best compound targeting to key SSSIRGs (Figure 6E).

The transcriptional regulation mechanisms between four key TFs (CTNNB1, XRN2, VEZF1, and NR3C1) and three key IRGs (TLR8, PPP3CA, and KRAS) were verified using ChIP-seq and ATAC-seq data. Binding peaks of three key IRGs and four key TFs in ATAC-seq data of 24 different cell lines from a variety of tissues are shown in Figures 7, 8, respectively. Additionally, transcriptional regulation mechanisms of four TF-IRG pairs were validated by ChIP-seq data available from the Cistrome database (Figure 9).

In total, 10 main cell types were identified using the UMAP analysis of unsupervised clustering (Figure 10A). DEGs from the top 1,500 variable genes were identified in all single cells. Cleveland’s dot plot was used to show the expression levels of key SSSIRGs, which were TLR8, PPP3CA, KRAS, and CD8A (Figure 10B). In addition, the proportions of 10 cell types identified previously in each sample were illustrated by the bar plot (Figure 10B). The heat map illustrated expression levels of the top 5 markers in 10 aforementioned subclusters (Figure 10C). Further, the bar plot was constructed to illustrate the fractions of these cell types (Figure 10D). Moreover, analysis of cell cycle demonstrated that the major cell types were in the G2M and S phases (Figure 10E).

The cellular feature plots of key SSIRGs (CD44, MKI67, TLR8, PPP3CA, and KRAS) and TFs (EGR1, CTNNB1, XRN2, VEZF1, and NR3C1) were used to show the expression levels and specific cellular locations of these key biomarkers in different immune cells (Figure 11A). In addition, Figure 11B illustrates the expression levels of three key SSSIRGs (TLR8, KRAS, and PPP3CA) in 10 immune cell clusters. Importantly, TLR8 was significantly highly expressed in Macrophage–Monocyte, Myelocyte, and Neutrophil; KRAS was significantly highly expressed in Myelocyte, CD8+ T cell, and Macrophage–Monocyte; PPP3CA was significantly highly expressed in CD8+ T cell, Myelocyte, and Macrophage–Monocyte (Figure 11B). Intersected cellular communication ligand-receptor plots showed the specific mechanism of intercellular signal transduction. Importantly, there was significantly strong intercellular communication among these cell clusters (Figure 11C).