Normalization and Principal Component Analysis of Three-Day Different Sequencing Data for Spinal Cord Injury

To identify transcriptomic changes that occur following spinal cord injury, our goal to compare transcriptomic expression profiles between spinal cord injured and healthy mice, with the goal of discovering novel targets that may contribute to explaining the pathological progression of spinal cord injury and exploring potential therapeutic interventions. Therefore, we used the two datasets GSE132242 and GSE166967 in the GEO database to obtain the original microarray sequencing data of seven mice each in the Sham group and the SCI3d group. Regarding data processing, we first use the GEO query, limma, and affy packages in the R language to obtain the original data in the data set, and then perform single-gene sequencing target correspondence, group correspondence, dimensionality reduction, and normalization processing on the original data (Fig. 1A, B). Subsequently, to verify the independence of each of the 14 samples, we used PCA analysis for data verification and found that the independence of the 14 samples in the two databases was good, and the independence and correlation between group differences met our subsequent analysis requirements (Fig. 1C, D).

Fig. 1

Normalization and principal component analysis of three-day different sequencing data for spinal cord injury. (A) Normalization of 8 samples in the GSE132242 dataset. (B) Normalization of 6 samples in the GSE166967 dataset. (C) PCA distribution plot of 8 samples in the GSE132242 dataset. (D) PCA distribution plot of 6 samples in the GSE166967 dataset

Different Genes in Mice Exhibiting Spinal Cord Injury 3 Days Old by Screening

Use the target mapping and the expression table corresponding to the data set to search for different genes and determine the up-regulated and down-regulated genes after 3 days of spinal cord injury for subsequent in-depth verification. Therefore, we used the pheatmap package in the R language to organize the two data sets and visualize all the different genes (Fig. 2A, B). Due to the large number of different genes in the data set, we then sorted them according to their statistical differences, using pheatmap to display the specific entries of the top 50 genes in the difference (Fig. 2C, D), some of these genes are related to inflammation, tissue fibrosis, cytokine regulation.

Fig. 2

Differential genes in mice exhibiting spinal cord injury 3 days old by screening. (A) Gene expression heatmap of 8 samples in the GSE132242 dataset. (B) Gene expression heatmap of 6 samples in the GSE166967 dataset. (C) The heat map of the top 50 gene expression of 8 samples in the GSE132242 dataset. (D) The heat map of the top 50 gene expression of 6 samples in the GSE166967 dataset

The Cross-Validation of the Two DATA sets Obtained the Differential Genes with Strong Reliability

We have analyzed the gene expression heatmap obtained by sorting the raw data. We discovered that, even though there was overlap among the factors with the most significant differences in the two datasets, they were not the same. Consequently, we opted for a strategy of mutual verification between the two separate datasets to pinpoint more dependable genes that are differentially expressed. We use the R program package ggplot2 to analyze the data set this time and consider that the condition for setting differential genes is that P. Value less than 0.05 is a differential gene, logFC greater than 1.5 is an ascending gene, and logFC is less than -1.5 is a descending gene. Generate volcano maps (Fig. 3A, B). It can be seen from the figure that in the GSE132242 dataset, 858 genes were upregulated after SCI (red), and 32 genes were downregulated after SCI (blue); In the GSE166967 dataset, 1041 genes were upregulated after SCI (blue), and SCI The post-down-regulated genes are 150 (red). Different genes obtained but it is also according to our conditional criteria. 3 days after spinal cord injury, most of the up-regulated genes, less down-regulated genes. Next, we used the Venn diagram to cross-validate the different genes in the two databases and found that there were overlapping genes in the differential genes in the two databases, which indicated that in the two different batches of sequencing, these genes were all present at 3 days after spinal cord injury Changes occur. In subsequent experiments, we will focus on the crossover genes in the two datasets (691 up-regulated genes and 13 down-regulated genes after SCI) (Fig. 3C, D).

Fig. 3

The cross-validation of the two data sets obtained the differential genes with strong reliability. (A) Volcano map of all differential gene expression in the GSE132242 dataset. (B) Volcano map of all differential gene expression in the GSE166967 dataset. (C) Venn diagrams of upregulated genes in GSE132242 and GSE166967 datasets. (D) Venn diagrams of downregulated genes in GSE132242 and GSE166967 datasets

GO Enrichment Analysis of SCI3d Differential Genes Found that the Spinal Cord Microenvironment was Still in a State of Immune Disorder

Subsequently, GO analysis was performed on the obtained data set containing 691 up-regulated genes and 13 down-regulated genes. Results revealed two highly enriched pathways: myeloid leukocyte activation (GO: 0002274) and positive regulation of cytokine production (GO:0001819) (Fig. 4A). Next, the GO enrichment results were analyzed from the aspects of cellular location (Cellular Component, CC), molecular function (Molecular Function, MF), and biological process (Biological Process, BP) (Fig. 4B, C, D). The results showed that the top three pathways with the highest enrichment of cellular component (CC) factors were all located in the extracellular region (GO: 0005576, GO:0044421, GO:0005615). The most enriched molecular function (MF) factors were interaction (GO: 0005488) and protein interaction (GO: 0005515). The most abundant biological process (BP) factors are positive regulation of biological processes (GO: 0048518), positive regulation of cellular processes (GO: 0048522), and stress response (GO: 0006950). In the secondary classification of GO, we further explored the number of differential genes in each category (Fig. 4E). Among biological processes, the highest number of enriched factors was cellular processes (GO:0009987). Among the cellular components, it is the cell (GO:0005623). Among the molecular functions, the most common is the binding activity (GO:0005488). In summary, the spinal cord injury microenvironment is in a state of severe immune disorder and the interaction and regulation of multiple factors.

Fig. 4

GO enrichment analysis of SCI3d differential genes found that the spinal cord microenvironment was still in a state of immune disorder. (A) GO enrichment analysis circle diagram of differential genes in two datasets. (B) GO enrichment analysis of differential genes in the two data sets Cellular ComponentTop25. (C) GO enrichment analysis of differential genes in the two data sets Molecular Function Top25. (D) GO enrichment analysis of differential genes in the two datasets Biological Process Top25. (E) Two-level classification diagram of GO enrichment analysis of differential genes in the two datasets

KEGG Enrichment Analysis Found that Various Chemokines and Metabolic Processes Play an Important Role in SCI3d of Spinal Cord Injury

To comprehend the common physiological consequences exhibited by differential genes, we sought to determine whether they are concentrated and enriched in certain pathways. Therefore, we conducted the KEGG enrichment on the obtained 691 up-regulated genes and 13 down-regulated genes. KEGG enrichment results indicated that association with organic systems, cellular processes, and environmental information processing (Fig. 5A). Pathways of the TOP 5 enriched factors are cytokine-cytokine receptor interactions (42), phagosomes (28), tuberculosis (27), osteoclast differentiation (25), proteoglycans in cancer (25) (Fig. 5B). Additionally, pathways enriched in differential genes also include the classic TNF signaling pathway, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, IL-17 signaling pathway. The results of enrichment were then annotated by KEGG pathways. Three days after SCI, differential genes were enriched in six major pathways: metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases (Fig. 5C). According to the enrichment results the KEGG, in SCI3d, inflammatory reactions and some necrotic reactions are still occurring in large quantities, abundant simultaneously, the system mobilizing numerous cell chemokines to address adverse reactions that occur and heal the damage through metabolic processes. Cells provide energy. It is noteworthy that certain factors rich in proteoglycans, lectins, phagosomes, and other pathways may serve as potential therapeutic targets.

Fig. 5

KEGG enrichment analysis found that various chemokines and metabolic processes play an important role in SCI3d of spinal cord injury. (A) Circle diagram of KEGG enrichment analysis of differential genes in two datasets. (B) KEGG enrichment analysis of top25 pathways of differential genes in the two datasets. (C) Pathway annotation map of KEGG enrichment analysis of differential genes in the two datasets

Reactome Enrichment Analysis Revealed the Important Role of Neutrophils in SCI3d of Spinal Cord Injury

To further elucidate the changes in SCI3d, we used Reactome, a clinically friendly database that has emerged in recent years, to conduct an enrichment analysis. From the Reactome enrichment analysis, we can see that the top five enriched pathways are Immune System, Innate Immune System, Neutrophil degranulation, Cell Cycle, Cell Cycle, Mitotic. Among these pathways and the enriched Top 25 pathways, the disorder of the 3D immune microenvironment of spinal cord injury and the significant of chromatin regulation at the molecular level can be confirmed again, in which many neutrophil-related factors are recruited to damaged areas to potentially against the adverse state of inflammation and gliosis scarring (Fig. 6A, B).

Fig. 6

Reactome enrichment analysis revealed the important role of neutrophils in SCI3d of spinal cord injury. (A) Circle diagram of Reactome enrichment analysis of differential genes in two datasets. (B) Reactome enrichment analysis of top25 pathways of differential genes in the two datasets

The Hub Gene UHRF1 was Calculated by Differential Protein Interaction Score 3 days After Spinal Cord Injury

Faced with such many differential proteins, we expected to identify a subset of the most important roles among them. Therefore, we predicted the interactions of all differential factors using the online website https://cn.string-db.org/ and exported and imported the results into Cytoscape 3.9.1 for central gene screening (Fig. 7A). We first screened the clusters with the highest scores among all the differential factor networks using the MCODE plugin (yellow group in the figure), and then displayed the clusters individually (Fig. 7B), followed by scoring cluster 1 with multiple algorithms using the cytoHubba plugin, showing the TOP10 pivotal factors selected by the MCC, DMNC, MNC, Degree algorithms (Fig. 7C–F). Next, we took the intersection of multiple algorithms (Fig. 7G). Seven differential factors that appeared in all the multiple intersections were selected for subsequent study, which were Tpx2, Kif11, Ncapg, Bub1b, Ttk, Uhrf1, Cdk1, and then an online tool http://genemania.org/ was used to generate a central gene map of the TOP 7 and their interactions with the predicted proteins (Fig. 7H). In this prediction result, we found that UHRF1 seems to be different from all the other factors, he does not appear to have the same biological function as the other 6 factors, but can interact with them, so that we became interested. In summary, the seven factors with the highest calculated scores in the protein interaction network 3 days after spinal cord injury were Tpx2, Kif11, Ncapg, Bub1b, Ttk, Uhrf1, Cdk1, and that Uhrf1 may independently play other biological roles.

Fig. 7

The hub gene UHRF1 was calculated by differential protein interaction score 3 days after spinal cord injury. (A) All differential protein interaction plots, with the group of proteins with the highest cluster score shown in yellow. (B) Enlarged view of the interaction of the group of proteins with the highest cluster score. (C–F) Four different algorithms MCC, DMNC, MNC, Degre calculated the TOP 10 factors for cluster 1. (G) Intersection of TOP 10 factors of four different algorithms. (H) Protein interactions predicted for the top 7 factors in the intersection set, with colors representing different biological functions

Inhibition of UHRF1 Improves Motor Function in Spinal Cord Injured mice

Next, we evaluated the Hub gene expression counts derived from the algorithm in the dataset (Fig. 8A–G). Further confirmation indicates that in the spinal cord injury model, the predicted mRNA expression and protein interaction scores for these seven factors are significantly different. Subsequently, we will validate the mRNA obtained from the spinal cord tissue sequencing results of C57 mice. We have chosen the specific inhibitor NSC232003 of UHRF1 and our team has been conducting intraperitoneal injection of ZnG for many years. The experimental evaluation of mouse motor function using footprint testing (Fig. 8H) and BMS score (Fig. 8I) confirmed that intervention with UHRF1 after injury improved mouse motor function. Meanwhile, we validated the expression and role of UHRF1 in neuronal cells using an in vitro neuronal cell line PC12 (Fig. 8J). We also validated the expression of UHRF1 after spinal cord injury and the role of inhibitors in mouse spinal cord tissue slices (Fig. 8K). The above experimental results indicate that three days after spinal cord injury, UHRF1 co localizes in neuronal cells. Inhibiting UHRF1 can improve spinal cord injury to a certain extent and ultimately promote the recovery of motor function in spinal cord injury mice.

Fig. 8

Inhibition of UHRF1 improves motor function in spinal cord injured mice. (A–G) 7 factors of Tpx2, Kif11, Ncapg, Bub1b, Ttk, Uhrf1, Cdk1 are expressed in transcriptome sequencing (n = 7, *P = 0.0104, **** < 0.0001). (H) Footprints of C57 mice in each group. (I) BMS scores of C57 mice in each group (n = 12, *P = 0.232, **P = 0.0091, ****P < 0.0001). (J) Immunofluorescence staining of UHRF1 in PC12 cell line. (K) Immunofluorescence staining of UHRF1 in C57 mouse tissues. (L) Immunofluorescence statistical analysis of PC12 cells in each group (*P = 0.0152, **P = 0.0027, ***P = 0.0002). (M) Immunofluorescence statistical analysis of each group of spinal cord tissue sections (*P = 0.0428, **P = 0.0040, ****P < 0.0001)