1.
Introduction
Spinal cord injury (SCI) is serious damage to the central nervous system, which once meant a lifetime of wheelchair and medication [1]. SCI not only brings physical and psychological pain to patients but also has become a major burden on society [2,3]. With advances in neuroscience, advanced interventions offer the possibility of neural regeneration and functional recovery. However, traditional treatment methods cannot completely repair SCI.
The development of precision medicine offers the opportunity to optimize personalized treatment plans based on individual genetic characteristics and pathogenesis, and blood and tissue-based bioassays are routinely used in treatment planning [4]. It is now known that in the process of gene expression, mRNA transcripts have undergone multiple interactions with various factors, and have undergone extensive reshaping [5]. For example, the combination of cis-acting factors and trans-acting elements before transcription [6], and ceRNA competition regulation mechanism [7], and mRNA N6 methylation modification [8] in the post-transcriptional level.
The genome-wide distribution of N6-methyladenosine (m6A) was unclear until 2012, but m6A is the most common RNA modification of eukaryotic mRNA [9]. M6A modification requires the participation of regulatory factors. These regulatory factors are usually defined as three categories including the m6A methyltransferases called "writers", the demethylases called "erasers", and m6A-binding proteins called "readers" [10]. Moreover, the perturbation of m6A mediated by these regulators has been shown to dysregulate the physiological mechanism of cells and lead to the occurrence of diseases [11].
In this study, we analyzed the gene expression profile of peripheral blood leukocytes from patients with acute SCI. Through network integration analysis, we discovered the pathogenesis and potential therapeutic targets mediated by m6A regulators.
2.
Method
2.1. Data collection
The RNA-seq profiles in peripheral white blood cells of spinal cord injury (SCI) samples (n = 38) and healthy individuals (HC, n = 10) without a history of central nervous system (CNS) pathology were collected from the gene expression omnibus [12] (GEO, https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE151371. We have collected 20 m6A regulators including 11 readers, 7 writers, and 2 erasers from the previous study [13]. The transcription factor of human was collected from AnimalTFDB 3.0 database [14] (http://bioinfo.life.hust.edu.cn/AnimalTFDB/). The data of the regulatory relationship between transcription factors and target genes were collected from TRRUST [15] (https://www.grnpedia.org/trrust/) and ORTI [16] (http://orti.sydney.edu.au/about.html). The relationship pairs between lncRNA and mRNA in the ceRNA competition mechanism were obtained from the starBase database [17] (http://starbase.sysu.edu.cn/). Further, downloaded the GRCh38 v29 version of the human genome annotation data (including the information of long non-coding RNA (lncRNA)) from GENCODE [18] (https://www.gencodegenes.org/) for identifying lncRNA. The gene sets of KEGG signal pathways were collected from the Molecular Signatures Database [19] (MSigDB, http://software.broadinstitute.org/gsea/msigdb). Data for drug-target gene relationships were also collected from the DrugBank database [20] (www.drugbank.ca).
2.2. Statistical analysis of RNA-seq profiles
First, the R package DEseq2 (v1.28.1) is used to standardize the count data and statistically test the differentially expressed genes between the HC and SCI two sets of samples. DESeq2 [21], which is a method for differential analysis of count data, uses the shrinkage estimation for dispersion and fold change to improve the stability and interpretability of the estimation. Differentially expressed genes (DEGs) were identified through the cutoff of the adjusted p-value as 0.05, and the differentially expressed m6A regulators were defined as m6A-DEGs.
2.3. Construction of the lncRNA-mRNA regulation network
Pearson correlation analysis was used in the standardized gene expression matrix to identify other genes that are significantly related to m6A-DEGs. We defined that the gene pairs with p-value < 0.01 and correlation coefficient |R| > 0.7 are significantly correlated. There were multiple lncRNAs in genes co-expressed with differentially expressed m6a regulators (m6A-lncRNAs), which play an important role in post-transcriptional regulation. Using the ceRNA collected from the starBase database, we identified m6A-lncRNAs involved in the ceRNA regulatory mechanism. Since microRNA (miRNA) expression information was not available in the RNA-seq profiles collected in this study, a direct measure of the regulatory relationship between lncRNA and mRNA. We have identified that each mRNA is significantly positively correlated with m6A-lncRNAs in expression (p-value < 0.01 and |R| > 0.7) and has been experimentally confirmed to compete with the m6A-lncRNA for the same miRNA in previous studies. Further, Cytoscape (v3.7.0) was also used to visualize the relationship of m6A-lncRNA-mRNA and the NetwokerAnalyzer tool was used to calculate the topological properties of the network.
2.4. Feature screening to identify spinal cord injury (SCI) patients
Since lncRNA plays an important role in the pathogenesis of SCI, we have used the R packages of boruta (v7.0.0) algorithm [22] to screen features from m6A-lncRNA involved in the ceRNA regulatory mechanism, where the p-value was set to less than 0.01 and the default was to run up to 100 times. The results of the Boruta algorithm divide the features into three categories: Tentative temporary/pending features (not enough to accept or reject), Confirmed features, and Rejected features. The importance of tentative features is so close to the best genetic attributes that boruta cannot make a strong confidence decision on the default number of random forest runs. For this case, we used the TentativeRoughFix function under the Boruta package to judge the importance of Tentative features. Further, we randomly divide the sample into two groups according to the ratio of 1:1, namely the train set and the test set. Two algorithms, support vector machines (SVM) and decision trees, were used to build classifiers to identify SCI samples. The plotting of the receiver operating characteristic curve (ROC) and the calculation of area under curve (AUC) was used to evaluate the effectiveness of the classifier.
2.5. Construction of the transcriptional regulatory network
We define m6A regulators whose gene expression fold change is greater than 1.5 between HC and SCI samples as true-m6A-DEGs. Using the human transcription factor (TF) data obtained from the AnimalTFDB database, we separately extracted TFs (including positive and negative TFs) from the gene sets related to each true-m6A-DEGs. We required that the target genes of the TF are expressed in the same direction as this TF during the change in expression with m6A-DEGs and has been demonstrated in previous studies. Further, Cytoscape [23] was used to visualize the relationship of TF-target genes and the NetwokerAnalyzer tool was used to calculate the topological properties of the network.
2.6. Spinal nerve repair-related signaling pathway research and potential drug prediction
First, we extracted 25 signal transduction-related pathways from the KEGG pathways downloaded from the MisigDB database. We calculated the enrichment scores of these signaling pathways in the SCI samples using the R package of the gene set variation analysis [24] (GSVA, v1.36.3) algorithm. Based on the enrichment score matrix, we have used the R package of the ConsensusClusterPlus (v1.52.0) [25] algorithm to cluster the SCI samples. In previous studies, it was found that MAPK, NOTCH, MTOR, and WNT signaling pathways are essential in nerve repair [26,27,28,29]. Next, we screened out genes related to m6A regulators (p-value < 0.01 and |R| > 0.7) from the genes involved in these four signal pathways. For these genes, we retained key genes that were significantly differentially expressed between HC and SCI samples and between clusters. Then, analysis of variance (ANOVA) was used and the cutoff of p-value was 0.05. Finally, we identified drugs that target these key genes and used cytoscape to visualize the key gene-drug relationship.
2.7. Statistical analysis
All statistical analyses and graph generation were performed in R (version 4.0.2, https://www.r-project.org/).
3.
Results
3.1. Identified six m6A regulators specifically expressed in spinal cord injury
The m6A regulators can regulate gene expression at the post-transcriptional level through epigenetic modification [13,30]. To explore the role of m6A regulators in SCI, we first identified genes that were significantly differentially expressed between SCI samples and HC samples. A total of 1710 up-regulated and 1884 down-regulated genes were recognized (Figure 1A). In these results, six m6A-DEGs including FTO, METTL14, RBMX, YTHDF2, YTHDC2, and HNRNPA2B1 were included in the significantly down-regulated genes (Figure 1B). In the previous study, FTO is considered an eraser, and METTL14 is considered a writer, RBMX, YTHDF2, YTHDC2, and HNRNPA2B1 are considered readers (Figure 1C). Through hierarchical clustering of all samples with the six m6A-DEGs as the characteristics, we found that there are obvious differences between SCI samples and HC samples (Figure 1C). Moreover, we found that SCI samples have different types, which may be determined by the origin of SCI, revealing individual differences in SCI. To explore the relationship between m6A-DEGs and other functional genes, Pearson's algorithm was used to calculate the correlation between the six m6A-DEGs and all other genes. Then, 6081 genes were identified to be significantly related to m6A-DEGs. Further, The R package clusterprofiler [31] was used to perform GO function enrichment and KEGG pathway analysis for the 6081 genes. We found that genes co-expressed with m6A-DEGs are mainly enriched in functions related to RNA transport and modification and protein processing (Figure 1D, E), indicating that the imbalance of the six m6A-DEGs in SCI plays an important role in the imbalance of physiological mechanisms.
3.2. The three key lncRNAs regulating the progress of SCI
LncRNA is an emerging biomarker and has emerged as significant player in almost every level of gene function and regulation [32], which played an important role in the development of disease [33]. To explore the regulatory mechanism of lncRNA related to m6A regulators in SCI, we extracted 184 m6A-lncRNAs from 6081 genes related to m6A regulators based on the lncRNA annotation information collected from GENCODE. Further, we found that only 92 m6A-lncRNAs are involved in the ceRNA regulation mechanism based on previous research results. Combined with the lncRNA-RNA interaction pairs collected from starbse [17], finally, 808 putative m6A-lncRNA-mRNA regulatory relationship pairs are identified (Figure 2A), which contains 53 lncRNA and 635 mRNA (Figure 2B). We found that the top three m6A-lncRNAs in degree including MIRLET7BHG, LINC00963, and PSMA3-AS1, which have ceRNA regulatory relationships with 77, 71, and 70 mRNAs (Figure 2C-E), respectively. Previous studies have shown that IMIRLET7BHG can regulate the activity of the Hedgehog signaling pathway through the ceRNA mechanism [34] and this signal plays an important role in the regulation of nerve cell fate, axon growth, and guidance [35]. LINC00963 and PSMA3-AS1 were also involved in the occurrence and development of a variety of complex diseases through the ceRNA regulatory mechanism [36,37]. Moreover, 635 mRNAs in the m6A-lncRNA-mRNA network were used for KEGG pathway enrichment analysis, and 21 significantly enriched signal pathways were identified (Figure 2F). We found that these protein-coding genes regulated by m6A-lncRNA are significantly enriched in the Neurotrophin signaling pathway and MAPK signaling pathway, indicating that the ceRNA regulation mechanism plays an important role in the development of nerve injury. Besides, we found that AKT2/3 and PIK3R1 are involved in the signal transduction of multiple signal pathways (Figure 2G), which may be potential therapeutic targets for SCI.
3.3. Machine learning-assisted identification of SCI patients
Since the diagnosis and treatment of SCI are still suboptimal, we constructed classifiers to identify patients with SCI based on m6A-lncRNAs. To obtain accurate signatures to construct the classifier, we first performed the boruta algorithm to screen features from 184 m6A-lncRNAs. After a rough judgment on the importance of tentative features, we finally rejected 153 m6A-lncRNAs and accepted the remaining 31 m6A-lncRNAs for classifier construction (Figure 3A). We extracted the expression profiles of the 31 features for clustering and found that there are obvious differences between HC and SCI samples (Figure 3B), indicating that the features are appropriate. Further, the samples were randomly divided into two groups (the training set and the test set) in the ratio of 1:1. Two machine learning algorithms, support vector machine (SVM) and decision tree, were used to construct the classifier. The decision tree algorithm emphasized that m6A-lncRNA MKNK1-AS1 is an important marker in the peripheral blood of SCI (Figure 3C), which will play an important role in the diagnosis of the severity for SCI. Moreover, the classifiers were used to validate the test set and the ROC curve was used to measure the reliability of the prediction results and the stability of the classifiers. We found that the AUC value of the ROC curve of the prediction result for the SVM classifier was 0.972 (Figure 3D), and the AUC value of the ROC curve of the prediction result for the decision tree classifier was 0.944 (Figure 3E). In comparison, the predictive performance of the classifier constructed using the SVM algorithm is better, which will serve the clinical diagnosis of SCI.
3.4. m6A regulators-TF-target genes-biological pathways cascade responses
TFs represent the convergence point of multiple signaling pathways in eukaryotic cells and can be used as potential therapeutic targets [38,39]. To explore the physiological mechanism of TF regulated by m6A regulators in SCI, we first identified TFs that were co-expressed with m6A regulators FTO and RBMX that were significantly differentially expressed between HC and SCI samples (p-value < 0.05 and the cutoff of foldchange is 1.5). Then, 48 TFs (32 positive and 16 negative) were identified to be significantly associated with FTO, and 302 TFs (219 positive and 83 negative) were identified to be significantly associated with RBMX. Further, we required that target genes regulated by TFs are also significantly related to the expression of FTO and RBMX. We identified eight target genes in the genes related to FTO that are regulated by four TFs (Figure 4A) and identified 1067 target genes in the genes related to RBMX that are regulated by 30 TFs. Four TFs related to FTO (ATC2, NFE2, STAT3, and STAT5B) were also specifically expressed in SCI (Figure 4B), and their regulated target genes next were used for GO function enrichment and KEGG pathway analysis. We found that the genes are significantly enriched in the production of skeletal muscle cells and the recruitment and differentiation of immune cells (Figure 4C and D). Additionally, we have identified the top four TFs of degree (ETS1, E2F4, CEBPA, and SPI1) in the transcriptional regulatory network constructed with genes related to RBMX were specifically expressed in SCI (Figure 4E). We found that the target genes in the network are significantly enriched in the functions of RNA transport and processing and protein synthesis (Figure 4F, G). Taken together, TFs-target genes related to FTO and RBMX were significantly enriched in the NF-kappa B signaling pathway related to spinal cord nerve injury [40].
3.5. Research on signal pathways related to spinal nerve repair reveals potential therapeutic drugs
The treatment of patients with SCI has not improved significantly in recent years, so it is important to explore the classification and drug repositioning of SCI. We first used 25 signal transduction pathways collected from the MisigDB database as functional features, and the GSVA algorithm was used to calculate the enrichment scores of SCI samples in the signal pathways. Based on the enrichment score matrix, the consistency clustering algorithm was used to cluster the SCI samples and the samples were clustered into four clusters (Figure 5A). We found that cluster 3 has an overall low signaling pathway activity (Figure 5B) and that such patients may have more severe neurological damage. In previous studies, the MAPK, NOTCH, MTOR, and WNT signaling pathways were found to be essential in nerve repair. The enrichment scores of the four signaling pathways in each cluster were extracted and compared. We found significant inter-cluster differences in the activity of the four signaling pathways, with cluster 1 possessing the highest pathway activity, followed by clusters 4, 2, 3 in that order (Figure 5C). To explore potential drugs for the treatment of SCI, we identified genes related to m6A-DEGs from the four pathways and require them to be significantly differentially expressed between clusters. Finally, 15 key genes related to five m6a-DEGs were identified. Combined with the drug-target genes data collected from Drugbank, we have identified 8 key genes that can be used as target genes for 21 drugs. We found that Fostamatinib is the drug with the highest degree (Figure 5D), indicating that Fostamatinib may play an important role in regulating the activity of spinal cord nerve repair-related signaling pathways and can be used as a potential drug for the treatment of SCI.
4.
Discussion
In this study, we emphasized the key factors that play an important role in SCI through the network integration analysis of N6 methylation-related genes. Six m6A regulators (FTO, METTL14, RBMX, YTHDF2, YTHDC2 and HNRNPA2B1) were identified and expressed specifically in SCI, and FTO and RBMX were significantly down-regulated. Through m6A-lncRNA-mRNA network analysis, three lncRNAs (MIRLET7BHG, LINC00963, and PSMA3-AS1) were emphasized to participate in the N6 methylation mechanism and regulate multiple signal pathways related to nerve injury. Moreover, we also used these m6A-lncRNAs to construct an accurate classifier for the identification of SCI patients to assist in clinical diagnosis. The transcriptional regulatory network was constructed to reveal the key regulators involved in m6A regulation. Based on the identified genes involved in the repair of nerve damage, we conducted potential drug discovery to serve the treatment of SCI.
Similar to DNA methylation modification, m6A RNA modification can be added by writing enzyme and removed by erasing enzyme, as well as the combination of reader enzyme to improve the translation efficiency of target RNA [41,42]. Therefore, the dysregulation of the expression of m6A regulators will affect the coding of genes and thus affect the physiological functions of cells. Chen's group has reported that FTO as N6-methyladenosine RNA demethylase exerts carcinogenic effects in acute myeloid leukemia[43]. Similarly, we found that FTO and RBMX play an important role in SCI.
With the rapid development of high-throughput technology, multi-omics big data analysis has become popular.
AKT2/3 and PIK3R1, co-expressed with m6A regulators, were found to be involved in multiple signaling pathways related to nerve injury, which can be used as potential therapeutic targets for the treatment of SCI. Target genes co-expressed with m6A regulators FTO and RBMX were found to be significantly enriched in the NF-kappa B signaling pathway. Studies have shown that NF-kappa B is a key regulator of the inflammatory process of reactive glial cells and its transgenic inhibition can reduce pain and inflammation after peripheral nerve injury [44].
With the rapid development of high-throughput technology, multi-omics big data analysis has become popular. It is expensive and time-consuming to identify new disease-biomarker associations based on traditional one-by-one experimental research. Therefore, the use of bioinformatics strategies in research can help reveal valuable results. We believe that our analysis will contribute to the development of clinical diagnosis and treatment strategies for SCI. Additionally, compared to identifying the molecular characteristics of the lesion tissue, exploring disease-related biomarkers in peripheral blood has unique advantages for clinical diagnosis and treatment.
In summary, we revealed the role of m6A regulators in SCI by network integration analysis and identify potential therapeutic targets and drugs. FTO and RBMX, which are significantly down-regulated in SCI, are emphasized to play an important role in the repair of nerve injury. Additionally, we have constructed a classifier to identify and assist the clinical diagnosis of SCI patients.
Acknowledgments
The authors gratefully thank the GEO for providing data for this work. S. W. and X. X. performed analysis and wrote the manuscript. C. L. and J. J collected the data. C. C. designed the study. All authors read and approved the final manuscript.
Funding
The work was supported by the National Natural Science Foundation of China (81572109), Wuxi Health Committee Research Grants for Top Talent Support Program (2020) and Nanjing University of Chinese Medicine Research Grant (XZR2020075).
DownLoad:
