DNA Microarray Analysis in Screening Features of Genes Involved in Spinal Cord Injury.

BACKGROUND Spinal cord injury (SCI) is the most critical complication of spinal injury. We aimed to identify differentially expressed genes (DEGs) and to find associated pathways that may function as targets for SCI prognosis and therapy. MATERIAL AND METHODS Seven gene microarray expression profiles, downloaded from the GEO database (ID: GSE33886), were used to screen the DEGs of leg tissue and to compare these between SCI patients and corresponding normal specimens. Then, GO enrichment analysis was performed on these selected DEGs. Afterwards, interactions among these DEGs were analyzed by String database and then a PPI network was constructed to obtain topology character and modules in the PPI network. Finally, roles of the critical proteins in the pathway were explained by comparing the enrichment results of the genes in sub-modules and all the DEGs. RESULTS A total of 113 DEGs were determined. We found that 21 up-regulated genes were enriched in 7 biological processes, while 9 down-regulated genes were significantly enriched in 4 KEGG pathways. The PPI network was constructed, including 40 interacting genes and 73 interactions. Three obvious function modules were identified by exploring the PPI network, and ACTC1 was identified as the critical protein in the 3 enriched signal pathways. However, no obvious difference was found in the signal pathway in which both the 11 genes in module 1 and all 113 DEGs participated. CONCLUSIONS Core proteins in the signal pathway associated with spinal cord injury may serve as potential prognostic and predictive markers for the diagnosis and treatment of spinal cord injury in clinical applications.

Background

Spinal cord injury (SCI) is the most critical complications of spinal injury [1]. It usually leads to serious dysfunctions of limbs below the injured segments [2]. Spinal cord injury harms the physical and mental condition of patients and causes substantial economic burden to society. Therefore, the prevention and treatment of spinal cord injury has become an important topic in medicine. Scientists are investigating various avenues for treatment of spinal cord injury; however, aside from methylprednisolone, none of these developments have reached even limited use in the clinical care of human spinal cord injury patients in the USA [3,4].

Microarray data analysis features high throughput and high sensitivity, which has made it possible to test the expression changes of the whole genome [5]. There have been many reports about spinal cord injury based on gene expression profiling [3,4,6]. Therefore, the development of microarray analysis has provided new insights into diagnosis and treatment of SCI [7-9]. Toshiya, using complementary DNA microarray analysis, identified 3 up-regulated and 7 down-regulated genes that may play a role in response to tissue damage or repair following SCI [10]. Jason identified targets at the mRNA level for SCI research and possible therapeutic intervention [8]. However, most of these reports only concentrated on several individual genes, lacking overall data. In addition, the DEGs of interest were rarely the same in different reports, and results of analysis of chips varied due to many factors, such as samples, platforms, and analysis methods. Thus, many experts have made great efforts to obtain reliable results from chips by use of meta-analysis; however, many essential issues, such as annotation and comparison of cross-platforms, could not be solved appropriately.

In the present study we identified 113 differentially expressed miRNAs between normal and spinal cord injury patients. Then we constructed a PPI network by use of the String database to explore 3 function modules to determine the biological processes and signal pathways they were involved in. We found no significant differences in signal pathway between genes participating in function module 1 and all 113 DEGs. GO functional and KEGG pathway analyses were conducted on up-/down-regulated expressed DEGs. Eventually, we found that a critical protein, ACTC1, could function as a prognostic and predictive marker in clinical treatment of SCI.

Material and Methods

Data source

The microarray expression profiles of spinal cord injury were extracted from the GEO (Gene Expression Omnibus, ) database under the accession number of GSE33886 [11,12]. A total of 7 specimens divided into 2 groups were available for the analysis, including 4 samples from normal vastus lateralis leg tissues and 3 specimens from leg tissues of patients with spinal cord injury. The data platform used was Affymetrix Human Gene 1.0 ST Array.

Data preprocessing

The original data are preprocessed by RMA function in the Affymetrix package [13] of R language. The original CEL files were converted into probe expression measures, and the probe-level data were converted into gene names by an annotation package based on the Bioconductor platform. Probes matching more than 1 gene were eliminated and the average value was used for probes matching the same single gene.

Screening of differentially expressed genes

DEGs that were significantly differentially expressed among the mesenchymal stem cell samples from SCI patients and the matched normal tissues were screened out by use of the LIMMA package [14] in R language. P=0.05 was used as the threshold of screening differentially expressed genes. Fold change >2 was used as the threshold to determine the significance of gene expression difference. We used cluster analysis and the corresponding cluster figure to ensure that the screened genes perfectly expressed the differences between samples of SCI patients and the matched normal tissues. The DAVID (Database for Annotation Visualization and Integrated Discovery) database [15] was used to perform the GO functional and KEGG pathway analysis. Significantly enriched GO terms and KEGG pathways with FDR (Fold Discovery Rate) <0.05 were screened out.

Protein-protein interaction (PPI) network construction

Protein-protein interactions regulating differentially expressed genes were obtained by use of String () [16], a database of known and predicted protein-protein interactions. Interactions with scores higher than 0.4 were screened out to construct the PPI network.

Network analysis

Topology characters of the PPI network, such as network nodes contribution, the shortest path, the average assemble coefficient, and centrality, were analyzed by use of Network Analyzer [17] in Cytoscape software. Then, the regulation relationship of each gene was obtained and KEGG pathway enrichment analysis was performed on the nodes in the network.

Network module analysis

Genes of each module were analyzed by ClusterONE [18] in Cytoscape software. P-value less than 0.05 and nodes larger than 5 were used as the cut-off criteria. Comparison analysis of functional enrichment between selected modules and all the DEGs was performed.

Results

Microarray data of 19 433 genes from 7 samples were obtained after data preprocessing. Data before and after normalization are listed in Figure 1.

Figure 1

Data before and after normalization. The horizontal axis represents the name of samples, while the vertical axis represents the expression value. The black lines stand for median, which can be used to identify the degree of standardization. In this figure, it can be seen that the black lines were almost on the same line, indicating an excellent degree of standardization.

DEGs screening

A total of 113 DEGs were identified between 4 normal and 3 spinal cord injury specimens, including 54 up-regulated and 59 down-regulated genes, accounting for 47.78% and 52.22%, respectively, of all the DEGs (Figure 2, Table 1).

Figure 2

Hierarchical clustering dendrogram of gene expression. The horizontal axis at the bottom represents the name of samples and the vertical axis on the left side represents the degree of gene clustering. The vertical axis on the right side represents the name of genes and the horizontal axis at the top represents the degree of clustering of samples. The red color stands for down-regulated while the green color stands for up-regulated. Generally, the samples can be divided into 2 clusters: the cluster of spinal cord injury samples and the cluster of normal tissue samples. Similarly, the genes can be divided into 2 clusters: the cluster of down-regulated genes in spinal cord injury tissue cells and the cluster of up-regulated genes in spinal cord injury tissue cells.

Table 1

Top10 significantly enriched up- and down-regulated DEGs.

Gene ID	Log2 (fold change)	q-value	Expression change
ACTC1	3.319406483	0.033053615	UP
ARHGAP36	2.464939379	0.03687478	UP
UNC13C	2.339913333	0.025644411	UP
MIR1-1	2.276586053	0.011153346	UP
HCN1	2.261911782	0.010131561	UP
MED12L	2.108777782	0.013091823	UP
PLCE1	1.924863506	0.010131561	UP
DNAH11	1.874092112	0.005566248	UP
MYH8	1.858627221	0.013827859	UP
TEP1	1.849621413	0.012077379	UP
TECRL	−5.03171	0.005566	Down
MYL3	−4.7519	0.010132	Down
MYH7	−4.65842	0.011153	Down
TPM3	−4.22697	0.010132	Down
TNNT1	−4.05069	0.020847	Down
LGR5	−3.54974	0.00687	Down
ATP2A2	−3.30509	0.022695	Down
MYOZ2	−3.17909	0.032815	Down
IDI2	−2.82965	0.006899	Down
PDLIM1	−2.76184	0.000293	Down

GO enrichment and KEGG pathway analysis

There were no significantly enriched GO terms or pathways among the up-regulated genes (Table 2). Down-regulated genes were enriched in more than 1 GO term and KEGG pathways. Among the up-regulated genes, 21 genes were enriched in 7 biological processes, including fatty acid metabolism, system regulation, and ventricle morphology. We found 16 down-regulated genes enriched in pathways composed by cells, including sarcomere, actin cytoskeleton, and muscle fibril. Nine down-regulated genes were significantly enriched in 4 KEGG pathways, including myocardial contraction, hypertrophic cardiomyopathy, dilated cardiomyopathy, and PPAR signal pathway.

Table 2

Top 10 GO terms and KEGG pathways enrichment results of up-regulated genes.

Category	Term	Count	FDR
BP	GO: 0006631~fatty acid metabolic process	7	0.012487908
BP	GO: 0044057~regulation of system process	8	0.007831171
BP	GO: 0055010~ventricular cardiac muscle morphogenesis	4	0.006895875
BP	GO: 0055008~cardiac muscle tissue morphogenesis	4	0.009712116
BP	GO: 0060415~muscle tissue morphogenesis	4	0.009712116
BP	GO: 0006941~striated muscle contraction	4	0.038267294
BP	GO: 0055114~oxidation reduction	9	0.040984753
CC	GO: 0030017~sarcomere	6	0.00253617
CC	GO: 0015629~actin cytoskeleton	8	0.001978778
CC	GO: 0030016~myofibril	6	0.001545624
CC	GO: 0044449~contractile fiber part	6	0.001263462
CC	GO: 0043292~contractile fiber	6	0.001403307
CC	GO: 0016459~myosin complex	4	0.030578677
CC	GO: 0032982~myosin filament	3	0.035582676
CC	GO: 0031966~mitochondrial membrane	7	0.032787333
CC	GO: 0005859~muscle myosin complex	3	0.030764112
CC	GO: 0016460~myosin II complex	3	0.033479545
KEGG	hsa04260: Cardiac muscle contraction	5	0.022085116
KEGG	hsa05410: Hypertrophic cardiomyopathy (HCM)	5	0.015368817
KEGG	hsa05414: Dilated cardiomyopathy	5	0.013827022
KEGG	hsa03320: PPAR signaling pathway	4	0.047785825

PPI network analysis

A total of 73 protein-protein interaction pairs with reliability scores higher than 0.4 were screened out, each with 40 codes, accounting for 35.39% of all DEGs. The PPI network of DEGs appeared to have high cluster properties (Figure 3). Most genes with strong interactions were significantly down-regulated (Figure 4).

Figure 3

PCA of the 20 most obvious DEGs. The horizontal axis represents the scores of first principal components and the vertical axis represents the scores of second principal components of samples. In the first principal components, 77.6% of variances of 20 samples was explained, while in the second principal component, 21.058% of variances of 20 samples was explained. In total, the resolution degree of variances was 98.658%. The start letter C_ represent normal tissues (total of 4), while others are spinal cord injury tissue.

Figure 4

Network nodes distribution. (A) represents degree distribution, (B) represents shortest pathways distribution, (C) represents average clustering coefficient, and (D) represent closeness centrality.

GO Enrichment and KEGG pathway analysis in sub-modules

Three modules were selected out with the threshold of node >5 and p-value <0.05 by ClusterONE plug in Cytoscape (Table 3, Figure 5). Enrichment analysis was performed on the 3 modules. Eight genes were enriched in 29 biological pathways in the first module, while 10 genes were enriched in 18 pathways related with molecule functions. Five genes were enriched in 5 KEGG pathways (Table 4). Several pathways associated with biological process were enriched, but no obvious signaling pathway or molecular function was enriched in the second or the third module (Table 5). Twelve biological processes, 9 cell constitution, 2 molecular function, and 3 KEGG pathways were enriched in 3 function modules and all DEGs (Figure 6), except for 5 biological processes that were enriched in the DEGs but not in the function module. This result revealed that there was no obvious difference in signaling pathways between the 11 genes that participated in module 1 and the 113 differentially expressed genes.

Table 3

The PPI network modules of DEGs.

Ranking	Numbers of nodes	Numbers of sides	Density	p-value of significance
1	11	32	0.582	3.908e+5
2	7	13	0.619	0.016
3	6	9	0.6	0.04

Figure 5

The PPI networks of DEGs. Red color represents up-regulated genes, green color represents down-regulated genes, and connecting lines represent the relationships.

Table 4

GO and KEGG enrichment analysis of the first module.

Category	Term	Count	FDR
BP	GO: 0006936~muscle contraction	7	5.94E-09
BP	GO: 0003012~muscle system process	7	5.25E-09
BP	GO: 0030048~actin filament-based movement	4	9.09E-06
BP	GO: 0055008~cardiac muscle tissue morphogenesis	4	1.12E-05
BP	GO: 0060415~muscle tissue morphogenesis	4	1.12E-05
BP	GO: 0006941~striated muscle contraction	4	4.64E-05
BP	GO: 0030705~cytoskeleton-dependent intracellular transport	4	5.62E-05
BP	GO: 0048738~cardiac muscle tissue development	4	6.37E-05
BP	GO: 0003007~heart morphogenesis	4	1.18E-04
BP	GO: 0033275~actin-myosin filament sliding	3	1.39E-04
CC	GO: 0015629~actin cytoskeleton	9	1.47E-11
CC	GO: 0032982~myosin filament	6	1.23E-11
CC	GO: 0016459~myosin complex	7	1.81E-11
CC	GO: 0030016~myofibril	7	3.69E-10
CC	GO: 0044449~contractile fiber part	7	3.30E-10
CC	GO: 0043292~contractile fiber	7	4.17E-10
CC	GO: 0005859~muscle myosin complex	5	3.75E-09
CC	GO: 0016460~myosin II complex	5	4.95E-09
CC	GO: 0030017~sarcomere	6	1.57E-08
CC	GO: 0005856~cytoskeleton	9	7.12E-07
MF	GO: 0008307~structural constituent of muscle	6	2.97E-09
MF	GO: 0003774~motor activity	6	7.56E-07
MF	GO: 0008092~cytoskeletal protein binding	7	8.76E-06
MF	GO: 0003779~actin binding	6	2.36E-05
MF	GO: 0005524~ATP binding	8	1.84E-04
MF	GO: 0032559~adenyl ribonucleotide binding	8	1.68E-04
MF	GO: 0030554~adenyl nucleotide binding	8	2.03E-04
MF	GO: 0001883~purine nucleoside binding	8	1.97E-04
MF	GO: 0001882~nucleoside binding	8	1.83E-04
MF	GO: 0005198~structural molecule activity	6	2.41E-04
KEGG	hsa05416: Viral myocarditis	5	2.34E-05
KEGG	hsa04530: Tight junction	5	1.49E-04
	hsa04260: Cardiac muscle contraction	4	6.13E-04
KEGG	hsa05410: Hypertrophic cardiomyopathy (HCM)	4	5.94E-04
KEGG	hsa05414: Dilated cardiomyopathy	4	6.01E-04

The TOP 10 result of every item, according to the p-value of significance, were selected.

Table 5

The enrichment results of biological pathways of modules 2 and 3.

Category	Term	Count	FDR
Cluster2	GO: 0006732~coenzyme metabolic process	3	0.036877789
Cluster2	GO: 0051186~cofactor metabolic process	3	0.029982387
Cluster3	GO: 0045859~regulation of protein kinase activity	4	0.016195912
Cluster3	GO: 0043549~regulation of kinase activity	4	0.009001935
Cluster3	GO: 0051338~regulation of transferase activity	4	0.006794111
Cluster3	GO: 0042325~regulation of phosphorylation	4	0.009964111
Cluster3	GO: 0051174~regulation of phosphorus metabolic process	4	0.008983695
Cluster3	GO: 0019220~regulation of phosphate metabolic process	4	0.008983695
Cluster3	GO: 0042127~regulation of cell proliferation	4	0.031166907

Figure 6

Three function modules and the GO enrichment of all DEGs. The yellow color stands for the transcription factors, red color represents up-regulated genes, and green color represents down-regulated genes.

Discussion

In this study we identified 113 differentially expressed genes from the most significantly altered genes using data from 2 datasets on the same platform of the GEO dataset, which can effectively avoid differences resulting from many factors, such as different samples, platforms, and analysis methods. Much of the previous research concentrated on the several genes they selected instead of the overall data. Thus, there are few common genes in the reported data. Although much effort was made, such as using meta-analysis to obtain more reliable results, critical problems could not be solved, such as annotation and comparison of trans-platforms. Our study effectually overcomes the above defects by using 2 datasets on the same platform, thus enhancing the reliability of our results.

As 2-gene interaction is more sensitive and accurate in clinical prognosis and diagnosis of SCI compared to use of a single gene, we selected 40 interacting genes by String method to construct a PPI network for better analysis of data, which turned out to be a network with biological functionality. The PPI network of DEGs in our study appeared to be highly clustered, which is the essential character of a biological network, reflecting a high degree of modularization. We divided the network into 3 separate modules and performed functional enrichment analysis on these 3 modules. Module 1 was proven to be the module with specific functions. The PPI network we constructed in the research was larger and more reliable than the previous ones.

Three function modules – cardiac muscle contraction (hsa05410), hypertrophic cardiomyopathy (HCM) (hsa05414), and dilated cardiomyopathy (hsa04260) – were explored using ClusterONE. More importantly, the critical protein in the 3 pathways – ACTC1 – was found to be up-regulated, while myosin heavy-chain 6 (MYH6), myosin heavy-chain 7 (MYH7), and myosin light-chain 3 (MYL3) were down-regulated. In addition, the availability and rationality of these functions were supported by the fact that there was no obvious difference between the signal pathway participation of genes in module 1 and those of all the DEGs, which make our results more reliable. Our results suggest a potential application of target genes and the proteins they express as prognostic and predictive markers in development and improvement of SCI therapies.

Previous research groups have identified several pathways associated with spinal cord injury. STAT3 (signal transducer and activator of transcription 3) was reported to be involved in scar formation after spinal cord injury by regulating astrocyte activity [19], which can in turn restrict the spread of inflammatory cells after SCI. In another study, the Rho signaling pathway was reported to be a potential target for therapeutic interventions after spinal cord injury [20]. However, these studies concentrated primarily on pathways after spinal cord injury, and there has been little research into the mechanisms that regulate the progression and development of SCI.

In this study, through comparing functions between module 1 and all the DEGs, we obtained 3 pathways – the myocardium contraction pathway, hypertrophic cardiomyopathy pathway, and the dilated cardiomyopathy pathway – that may function as predictors of clinical outcome in SCI patients. Evidence from previous research showed that SCI has many common complications, and cardiac dysfunctions were among the most important complications following SCI. Cardiovascular disturbances were proven to be the leading causes of morbidity and mortality in both acute and chronic stages of SCI [21,22]. Therefore, we suggest that the 3 significantly enriched pathways found this study may be associated with the occurrence of spinal cord injury because cardiovascular deconditioning may lead to loss of skeletal muscle pumping activity [23].

ACTC1 protein, the critical protein in these 3 pathways, was significantly up-regulated, while myosin heavy-chain 6 (MYH6 [20]), myosin heavy-chain 7 (MYH7 [21]), and myosin light-chain 3 (MYL3 [22]) were down-regulated.

These 3 signal pathways have a homo-regulator, ACTC1. Therefore, it is important to study the regulation mechanism of ACTC1 to help understand its role in predicting targets in the injury process in SCI patients [24]. The up-regulation of ACTC1 can affect the myofilament contraction in the myocardium contraction pathway, and the state of striated muscle in hypertrophic cardiomyopathy pathway and the dilated cardiomyopathy pathway. In fact, in the last decade, ACTC1 [25] has increasingly been recognized as an important property of muscle state [26]. It mediates a variety of essential biological functions in all eukaryotic cells, providing a structural framework around which cells shape themselves [27]. Our results in the present study suggest its momentous role in the cytoskeleton and a potential role as a prognostic factor and predictor in clinical management of SCI patients. In addition, we discovered that DEGs were significantly enriched in cytoskeleton-dependant intracellular transport (GO: 0030705) and actin filament-based movement (GO: 0030048), confirming the role of actin in SCI development and progression. We also that found myosin heavy-chain 6 (MYH6 [28]), myosin heavy-chain 7 (MYH7 [29]), and myosin light-chain 3 (MYL3 [30]) were down-regulated. Recent studies have discovered that the regulation of MHC isoform expression involves a complex interaction of multiple control mechanisms, including myogenin and calcineurin, and indicated that other intracellular signaling pathways were also likely to contribute. The enriched pathway in the present study – muscle tissue morphogenesis (GO: 0060415) – is likely involved in the progression of SCI by participating in this process. Additionally, it has been reported that mutations in MYH6 can cause a spectrum of phenotypes, ranging from HCM (hypertrophic cardiomyopathy to DCM (dilated cardiomyopathy) [31], which is correlated with the pathways that were found to be enriched in our study – the hypertrophic cardiomyopathy pathway and the dilated cardiomyopathy pathway. Thus, we hypothesize that down-regulation of MYHC is related with SCI and has a potential to function as a new marker in clinical diagnosis and development of therapies for SCI.

Conclusions

Three function modules were explored using ClusterONE and 3 enriched signaling pathways were found: the myocardium contraction pathway, the hypertrophic cardiomyopathy pathway, and the dilated cardiomyopathy pathway. We found no differences among the signal pathways between genes participating in module 1 and all the DEGs, confirming its potential role as a prognostic and predictive indicator for use in SCI patients. Moreover, critical proteins, such as ACTC, MYH6, MYH7, and MYL3, in the enriched pathways may also provide useful information for better understanding of genetic changes in SCI. Some of these crucial genes might be utilized in diagnosis and prognosis of SCI.