Research Article - Biomedical Research (2018) Volume 29, Issue 4
Accepted date: November 26, 2017
Objective: The aim was to research Spinal Cord Injury (SCI) related genes and functions for postoperative recovery and complications and prevention.
Methods: E-GEOD-18179 was used to screen Differently Expressed Genes (DEGs) between SCI samples and controls. Adjacency matrix was constructed as gene co-expression network. Sub-network was identified based on spearman correlation coefficient. Assortativity coefficient was calculated to represent the relationship among nodes. Gene Ontology (GO) functional enrichment was processed and guilt-byassociation was undertaken.
Results: Total 129 DEGs were obtained and Adjacency matrix with 129 × 129 genes was constructed. There into, TYROBP and MIA were with weight of 0.9695. Total 129 genes were enriched in 127 GO terms, such as molecular function and catalytic activity. K terms and related AUC values were obtained. The guilt-by-association node degrees AUCs were generated. In addition, the optimal list genes were obtained, such as JUN and PTK2B. Finally, optimal GO terms were obtained, such as regulation of signal transduction, positive regulation of macromolecule metabolic process and cellular component organization.
Conclusion: Screened genes including JUN, PTK2B, CCND1 and ADAM8 might be potential genes for SCI postoperative recovery and complications and prevention.
Spinal cord injury, Guilt-by-association, Differently expressed genes, Functions, Adjacency matrix
Spine is a group of nerves which is protected by the individual vertebrae of the spine . The main function of spin cord is to send signals from the brain to other regions . Spinal Cord Injury (SCI) can be classified into complete injury and incomplete injury, which was with different symptoms, such as pain, numbness, paralysis to incontinence and even quadriplegia . Moreover, various complications may occur, including muscle atrophy, pressure sores, infections and respiratory problems.
In 2006, Schmitt et al.  found that different degree of recovery was showed in different strains of rats, and difference in functional recovery was closely related with gene expression after SCI. In Wistar rats with high-level SCI, α1-adrenergic receptors gene were a contributing factor for hypertension . Moreover, EPH Receptor A4 (Epha4) was confirmed to be related with deprenyl, which was found to protect nervous system by regulating the expression of p75NTR and TrkB after SCI .
Though various genes were researched, recovery and complication were not well controlled. The changes of genes after spinal cord injury were complexly and closely interacted. Gene Ontology was an effective method to describe the function of genes. In addition, Gene Ontology (GO) function was introduced to gene expression network, which could functionally classify the network. In order to research SCI related genes and functions with integrality; guilt-by-association was used in this study to increase the statistical reliability of GO terms. It might lay a foundation for SCI treatment and complication prevention.
Microarray data acquiring and preprocessing
Raw microarray profile of E-GEOD-18179 including 32 spinal cord injury samples and 30 controls were obtained from Array Express database with the platform of A-AFFY-6-Affymetrix GeneChip Murine Genome U74Av2 (MG_U74Av2). Random Multiple Access (RMA) algorithm was processed for raw data preprocessing with 3 steps: background correction, replicated probes removing and standardization . After mapping probes to genes, gene expression profile was obtained.
Differently expressed genes screening
The expression levels of genes in gene expression profile were calculated by limma package, and examined by T-test and Fisher exact test. LmFit function was used for linear fitting, expected Bayesian estimation (eBayes) was processed for statistics and FDR was corrected p value. The threshold of DEGs was |logFC|>0 and p<0.05.
Co-expression network of DEGs
Adjacency matrix was constructed as gene co-expression network. Spearman correlation coefficient between two genes was calculated. Based on the strength of relationships, adjacency matrix was constructed. Nodes of sub-network were identified with the threshold of spearman correlation coefficient>0.8 and degree>1.
Density of node degree
Based on the occurrence of degree, R package was used to show the density of node degree. Assortativity coefficient was a kind of Pearson correlation coefficient based on degree. The results of assortativity coefficient represented the relationship among these nodes. Positive value represented synergetic relationship of nodes with same degrees, while negative value represented a link between nodes with different degrees. Normally, the value of assortativity coefficient ranged from -1 to 1.1 represented assortativity pattern, 0 represented non-assortativity pattern and -1 represented negative relationship.
GO functional enrichment analysis
A total of 18421 GO terms with 18402 genes were downloaded from GO consortium. Screened DEGs were enriched in GO terms, and enriched GO terms with more than 20 genes were screened.
Inference analysis of network correlation
Gene i in gene co-expression network was selected and all genes except i were enriched in K terms. MF values of each gene i were calculated by following formula:
Numink represented numbers of neighbor genes in K terms and Numoutk represented numbers of neighbor genes which did not exist in K terms.
Three times folding method was used to calculate AUC value of K terms, then average AUC was calculated and evaluated for each GO function.
Based on emergent frequency of genes, cool read was processed to evaluate these predicted functions. If a fine AUC was obtained, the GO term was a good candidate with differently expressed genes. Finally, the most excellent gene list was generated. Then, the most excellent gene list was used to predict GO terms and AUCs were obtained. GO terms with AUC more than 0.6 were screened.
Acquisition of gene expression profile and DEGs screening
After mapping probes to genes, gene expression profile with 8713 genes was obtained. Based on the threshold of DEGs, a total of 129 DEGs were obtained, including Small Proline Rich Protein 1A (SPRR1A), Small Proline Rich Protein 1B (SPRR1B), Activating Transcription Factor 3 (ATF3), Desmoplakin (DSP) and Adenylate Cyclase Activating Polypeptide 1 (ADCYAP1).
Co-expression network of DEGs
Adjacency matrix with 129 × 129 genes was obtained. Based on the threshold of spearman correlation coefficient>0.8 and degree>1, a total of 84 correlations were obtained. There into, TYRO Protein Tyrosine Kinase Binding Protein (TYROBP) and Melanoma Inhibitory Activity (MIA) were with weight of 0.9695. Small Proline Rich Protein 1B (SPRR1B) and PDZ and LIM Domain 1 (PDLIM1) were with weight of 0.9760.
Density of node degree
Based on the occurrence of degree, density of node degree was shown in Figure 1. The degree of node ranged from 54 to 72. The median of these degrees was 64.7674. In addition, the assortativity coefficient of this density of node degree was 0.3422, which was confirmed that the subnetwork was with assortativity pattern.
GO functional enrichment analysis
A total of 129 genes were enriched in 127 GO terms (Table 1). Top 5 GO terms were molecular function, catalytic activity, binding, protein binding and cellular component. Various genes were enriched in molecular function, such as Activating Transcription Factor 3 (ATF3), Desmoplakin (DSP), Adenylate Cyclase Activating Polypeptide 1 (ADCYAP1) and Jun Proto-Oncogene, AP-1 Transcription Factor Subunit (JUN). Similarly, different genes in catalytic activity included Tubulin Beta 6 Class V (TUBB6), ATP Binding Cassette Subfamily A Member 1 (ABCA1) and CCR4-NOT Transcription Complex Subunit 4 (CNOT4).
Table 1: Top 10 GO terms enriched by 129 genes.
Relevant inference and cool reading inference
Based on relevant inference analysis, K terms and related AUC values were obtained. As shown in Table 2, top 3 K terms were regulation of cell differentiation (GO: 0045595, degree=66.697, AUC=0.5767), positive regulation of macromolecule metabolic process (GO: 0010604, degree=66.6754, AUC=0.6539) and regulation of macromolecule metabolic process (GO: 0060255, degree=66.5429, AUC=0.6694). The guilt-by-association node degree AUCs were generated (Figure 2). The average of AUC was 0.52198. After cool read of K functions, cool read node degree AUCs were obtained (Figure 3). The average of cool read node degree AUCs was 0.5488. In addition, the optimal list genes were obtained, such as JUN (MF score=0.0338), PTK2B (Protein Tyrosine Kinase 2 Beta, MF. score=0.0322), CCND1 (Cyclin D1, MF. score=0.0323) and ADAM8 (ADAM Metallopeptidase Domain 8, MF. score=0.0319). Finally, optimal GO terms were obtained. Therein, top 3 optimal GO terms with higher AUC were regulation of signal transduction, positive regulation of macromolecule metabolic process and cellular component organization (Table 3).
|GO. term||Node. degree||AUC||Node. degree. AUC|
Table 2: Top 10 K terms by relevant inference analysis.
|GO: 0009966||Regulation of signal transduction||Biological_process||0.653968||0.680476||Ptk2b, Edn1 Adm, Adcyap1, Sfrp1, Mllt3, AI464131, Sox11, Tnfrsf19, Gadd45a , Map4k2, Adam8, Atp6ap1, Crlf1 Fas, Ccnd1, Lgals3, Scg2, Rdh11, Igfbp2, Psap, Jun, Abca1|
|GO: 0010604||Positive regulation of macromolecule metabolic process||Biological_process||0.691506||0.663036||Edn1, Ccl2, Ptk2b, Cst3, Jun, Sox11, Hoxa7, Atf6, Adcyap1, Atf3 , Ccl3 , Crem, Foxo4, Sfrp1, Ank3, Star , Cdkn1a, Gadd45a, Map4k2, Adam8, Ccnd1, Crlf1, Fas, E330034, G19Rik|
|GO: 0016043||Cellular component organization||Biological_process||0.603535||0.662069||Asf1a, Pom121, Jun, Mrps25, Gadd45a, Ttc17, Arpc5, Ptk2b, Rhou, Sfrp1, Dsp, Ccl2, Maea, Abca1, Stxbp1, Topbp1, Ank3, Ctss, Tubb6, Atp6, ap1, Fas, Mmp9, Gmnn, Adm, Plxna3, Adcyap1, Edn1, Adam8, Acan, Cst3, Lgals3, Mia|
|GO: 0050789||Regulation of biological process||Biological_process||0.600944||0.610714||Hoxa7, Cdkn1a, Jun, Edn1, Ccl2 , Scg2, Ccnd1, Mmp12 , Sox11, Sfrp1, Igfbp2, Adcyap1, Adm, Atf3, Crlf1, Cst3, Ptk2b, Mmp9, Foxo4, Lgals3, Neu1|
Table 3: Top 5 optimal GO terms.
In most previous studies, microarray studies failed to research sufficiently the biological variability. Thereby, the validation rate was used within a system in this study to identify accurate genes and their interactional functions. In 2014, Ego network was firstly reported , which could screen potential biomarkers and lay the deeper foundation for breast cancer. Importantly, this method increased the prediction accuracy of gene functions with higher statistical confidence [9,10]. In addition, relevant inference and cool reading inference using make the functional enrichment prediction was with statistical confidence. Genes in same GO function were with preferential connection [10,11]. In this study, several genes including JUN, PTK2B, CCND1 and ADAM8 were screened and their related functions were identified.
JUN encoded a protein which was highly similar to the viral protein, and was confirmed to interact with specific target DNA sequences to regulate gene expression . After cervical spinal cord injury, JUN expressed highest at 7 days, especially in axotomized medullary neurons . Moreover, JUN could be re-induced by axonal regeneration after chronic spinal cord injury . Vinit et al.  also confirmed that after SCI, there were difference between molecular cell body response of spared neurons and axotomized neurons because of the regulation of JUN. In other diseases, the expression level of JUN and its’ related pathways were confirmed to be closely associated with neuronal apoptosis and virus infection . In addition, JUN N-terminal kinase-dependent was found to play important role in respiratory . Besides, PTK2B was screened as a potential key gene of SCI. As we all known, PTK2B encoded a cytoplasmic protein tyrosine kinase, which was involved in calcium-induced regulation of ion channels and activation of the map kinase signaling pathway. Previous study showed that the common variant in PTK2B was related with late-onset Alzheimer’s disease . Liu et al.  researched ion channel blockers and found that ion channel blockers were effective on SCI treatment but with complex effects. Interestingly, the inhibitor of acid-sensing ion channels could significantly decrease neurological deficits following SCI . In larval lamprey, SCI could induce changes of ion channels of reticulospinal neurons . Thereby, JUN and PTK2B might be potential targets for SCI postoperative recovery and complications and prevention.
Two other genes, CCND1 and ADAM8 were screened as optimal list genes in this study. CCND1 encoded a protein which belongs to the highly conserved cyclin family, whose members were characterized by a dramatic periodicity in protein abundance throughout the cell cycle . In previous study, CCND1 was found to be targeted by mir-195, which suppressed metastasis of osteosarcoma cells. In addition, et al. found that mir-19a was related with various genes of chondrogenic differentiation, such as CCND1 . In addition, CCND1 genomic amplification was occurred during malignant transformation of a cell tumor of bone . Similarly, ADAM8 encoded a protein which might be involved in cell adhesion during neurodegeneration . In adult SCI mice models, ADAM8 was upregulated in endothelial cells, and also related with angiogenesis . Thereby, both CCND1 and ADAM8 might be important genes for SCI postoperative recovery and complications and prevention.
In this study, these screened genes enriched in optimal GO terms, including regulation of signal transduction, positive regulation of macromolecule metabolic process and cellular component organization. Quan et al.  found that acupuncture could activate signal transduction and further promote the repair of SCI. After SCI, signal transduction of erythropoietin receptor was also confirmed to inhibit apoptosis of neuron cells . Besides, Zhao et al.  also found that multitude gene signal transduction was contribution to the development of neural stem cells after SCI. Though various macromolecules can’t be used with the limitation of across blood spinal cord barrier, the functions of these macromolecules for SCI recovery couldn’t be ignored . In rat models, convective delivery of macromolecules was used as an effective treatment method for injury or spinal cord disease . It was worth noting that Fehlings et al.  have obtained remarkable advances in recruiting cell therapy for SCI. Besides, cellular transplantation could replace the lost elements, and further regulate cell functions after SCI .
In this study, relevant inference and cool reading inference using make the functional enrichment prediction was with statistical confidence. Screened genes including JUN, PTK2B, CCND1 and ADAM8 might be potential targets for SCI postoperative recovery and complications and prevention. Various GO terms, such as regulation of signal transduction, positive regulation of macromolecule metabolic process and cellular component organization, were important in these processes.