Sub-pathway analysis for severe burns injury patients: Identification of potential key lncRNAs by analyzing lncRNA-mRNA profile.

The aim of the study was to identify key long non-coding RNAs (lncRNA) and related subpathways following severe burn injuries and research their functions. The miRNA-mRNA and lncRNA-miRNA interactions were downloaded from starBase v2.0 database. In addition, mRNA-miRNA interactions were obtained from TarBase, mirTarBase, mir2Disease, miRecords (V4.0) databases. The relationships of lncRNA-miRNA-mRNA were constructed. Genes of expression profiling were intersected with mRNA and lncRNA in lncRNA-mRNA interaction. Screened mRNAs were enriched into various pathways and screened lncRNAs were embedded into candidate pathways. Wallenius approximation methods were used to calculate the false discovery rate value of each sub-pathway. Based on the results of significant sub-pathways, the related lncRNA-mRNA network was constructed. A total of 18,081 genes were obtained. The lncRNA-mRNA intersections including 835 lncRNAs, 1,749 mRNAs and 7,693 interacting pairs were constructed. The enriched mRNAs were further enriched into various candidate pathways such as ribosome biogenesis in eukaryotes. Several sub-pathways were screened, including ribosome biogenesis in eukaryotes and MAPK signaling pathway. The network of pathway-lncRNA-mRNA was constructed. Hub-genes were identified, including C14orf169 and YLPM1. Several hub-lncRNAs were obtained, including PRKAG2 antisense RNA 1 and LEF1 antisense RNA 1. Several hub-lncRNAs including C14orf169, YLPM1, TTTY15, and PCBP1-AS1 were screened. The sub-pathways regulated by these lncRNAs were identified, and functions were predicted.

Introduction

Burn injury is a common traumatic injury with considerable mortality and morbidity, requiring a long period of hospitalization and rehabilitation (1). Severe injury may induce a series of physiological changes, such as the change of metabolic rate, core temperature and various substrate cycling (2). In addition, both complications and burn wound of burn injury influence the physical and psychological health (3).

In previous studies, the expression of various genes was found to be altered after severe burns. Noronha et al (4) found that several genes, including C-X-C motif chemokine ligand 8, interleukin 6, tumor necrosis factor and major histocompatibility complex, Class I, E played a contributory role in wound infection induced by burns. In addition, Padfield et al identified various genes that encode chemokines, complement, oxidative-stress and immune functions following severe burn trauma (5). Long non-coding RNAs (lncRNA) were defined as non-protein coding transcripts longer than 200 nucleotides (6). Initially they were regarded as ‘noise’ with no biological function in the transcription of the genome. However, lncRNAs were confirmed to regulate gene expression in epigenetics, transcription and at the post-transcriptional level (7). In addition, lncRNAs were closely associated with the occurrence, development and prevention of disease by participating in chromatin modification, transcriptional activation, transcriptional interference and intranuclear transport (8). Haijun et al found a differential expression of lncRNAs compared to controls in skeletal muscles of rats with burn injuries (9). However, lncRNA and related functions and pathways in human have not been studied.

Identification of sub-pathways regulated by lncRNA are essential for monitoring of the occurrence and development of disease, and also help to research the functions of lncRNAs. In this study, we suggested a new method for identification of function regulated by lncRNA. The expression profiles of lncRNA and mRNA were integrated. Moreover, the regulation role of miRNA-target mRNA and topological analysis of pathways were processed to identify sub-pathways regulated by lncRNA and predict the function of lncRNAs in severe burn injury patients.

Materials and methods

Samples

The expression profiling of E-GEOD-37069 with 553 injury samples and 37 controls, were downloaded from ArrayExpress Archive (http://www.ebi.ac.uk/arrayexpress/), which was based on the platform of [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA, USA). These 553 injury samples were gathered from 244 severe burns patients over time. The patients had burns covering >20% of the total surface and were admitted at 96 h. The raw profiles were downloaded for further analysis.

Data preprocessing and genes screening

Based on the database of hgu133plus2, the raw data were preprocessed. Genes were obtained by combing various corresponding probes, and the expression value of this gene was assigned by the mean value of these probes.

Construction of candidate lncRNA-mRNA interaction

miRNA-mRNA interactions and lncRNA-miRNA interactions were downloaded from starBase v2.0 database. Moreover, some mRNA-miRNA interactions were obtained from various databases including TarBase, mirTarBase, mir2Disease, and miRecords (V4.0). According to shared miRNA in the aforementioned interactions, the relationships of lncRNA-miRNA-mRNA were constructed, and candidate lncRNA-mRNA regulatory relationships were obtained. The lncRNA-mRNA regulatory relationships were screened using the criteria: i) Hypergeometric functions were used to assess the significance of miRNA that have a regulatory relationship with both lncRNA and mRNA. The Benjamini and Hochberg (BH) method was used to adjust the P-value. The lncRNA-miRNA intersections with false discovery rate (FDR) <0.05 were screened. ii) Top 20% lncRNA-mRNA intersections with higher Jaccard coefficient were retained. The formula of hypergeometric functions was shown as:

Construction of co-expression lncGenePairs

The genes of expression profiling were intersected with mRNA and lncRNA in a lncRNA-mRNA interaction, respectively. Based on these matched expression profiles of lncRNA and mRNA, the Pearson's correlation coefficient was used to assess candidate lncRNA-mRNA intersections (10), and Fisher Z transformation was processed to calculate the P-value of lncRNA-mRNA-intersections (11). Furthermore, the P-value was adjusted as per the FDR using the BH method, and lncRNA-mRNA-intersections with FDR <0.05 were screened with significance (12).

Sub-pathway analysis

The lncRNA-regulated sub-pathway comprised two steps: Mapping interested mRNAs and lncRNAs into linked pathways, and identifying lncRNA-regulated sub-pathways (Fig. 1).

Figure 1.

Analysis process of identification of lncRNA-regulated sub-pathways. LncRNA, long non-coding RNAs.

Pathway enrichment analysis of screened mRNA

Based on the Kyoto Encyclopedia of Genes and Genomes database, the screened mRNAs were enriched into various pathways. Fisher's exact test was used to identify significant pathways. Then, the P-value was adjusted into FDR by the BH method. The pathways with FDR <0.01 were regarded as candidate different pathways.

LncRNAs embed into enriched pathways

The screened lncRNAs were embedded into candidate pathways according to the results of lncGenePairs. Then, these pathways enriched lncRNA and related genes were obtained.

Identification of sub-pathways regulated by lncRNA

The screened lncRNA and related mRNA were regarded as flag nodes. Based on the situation of lncRNA embedding, the shortest route between two flag nodes was calculated. Once the molecular number between two flag nodes ≥1, the two flag nodes were combined to one node. In addition, the nodes in pathways were also calculated, and pathways with ≥ 8 nodes were regarded as candidate sub-pathways.

Assessment of sub-pathways and identification of hub-lncRNA

The Wallenius approximation methods were used to calculate the FDR value of each sub-pathway. In this process, the pathway weight was calculated using the formula:

While the significance was evaluated using the formula of P-value = F (x, m1, m2, n, weighti), where x, represents the number of lncGenePairs that were enriched into this hub-pathway; n, was the total number of mRNA; m1, represents the number of mRNA enriched into this hub-pathway; m2, the number of mRNA that could be annotated to lncGenePairs; w, was the weight of this hub-subway; P and G were mRNA numbers in this sub-pathway; G and L were mRNA numbers regulated by lncRNA; β=1. In addition, the sub-pathway with FDR <0.01 was screened. Based on the results of significant sub-pathways, the related lncRNA-mRNA network was constructed. The lncRNAs with a degree more than the average degree were regarded as hub-lncRNAs.

Results

Genes screening and obtaining of lncRNA-mRNA interaction

After data preprocessing and combination, 18,081 genes were obtained. In addition, lncRNA-mRNA intersections including 835 lncRNAs, 1,749 mRNAs and 7,693 interacting pairs were constructed by several databases and sharing miRNAs.

Taking intersection between screened genes and lncRNA-mRNA, 1,641 mRNA expression profiles and 125 lncRNA expression profiles were obtained. After Fisher Z transformation, 656 lncGenePairs including 101 lncRNAs, and 503 mRNAs were constructed with P<0.05. The top 5 lncGenePairs were EMX2 opposite strand/antisense RNA-ACTN3 (P=4.80×1027), TBX5-AS1-ADAMTS16 (P=3.79×1018), IQCH-AS1-ANKRD37 (P=0.02), MFI2-AS1-APLP1 (P=1.16×1026) and HOXA-AS2-ARL16 (P=0.04) (Table I).

Table I.

The top 10 hub-lncRNA genes.

Lnc	Gene	corValue	P-value
EMX2OS	ACTN3	4.23×10−1	4.80×1027
TBX5-AS1	ADAMTS16	3.47×10−1	3.79×1018
IQCH-AS1	ANKRD37	9.48×10−2	2.12×10−2
MFI2-AS1	APLP1	4.20×10−1	1.16×1026
HOXA-AS2	ARL16	8.5×10−2	3.8×10−2
HOXA-AS2	ASB9	9.97×10−2	1.53×10−2
FOXN3-AS1	ASPSCR1	3.31×10−1	1.22×1016
LBX2-AS1	ASTN1	1.77×10−1	1.43×1005
LEF1-AS1	BBIP1	2.97×10−2	1.71×1013
TDRG1	BNC1	2.55×10−2	6×1010

The enriched mRNA were enriched into various candidate pathways including ribosome biogenesis in eukaryotes (FDR =8.19×105), neurotrophin signaling pathway (FDR =8.07×105) and thyroid cancer (FDR =7.99×106) (Table II).

Table II.

The top 5 pathways enriched by candidate mRNAs.

Pathway_id	Pathway_P-value	P.adjust (FDR)
hsa03008: Ribosome biogenesis in eukaryotes	8.19×105	8.19×105
hsa04722: Neurotrophin signaling pathway	8.07×105	8.07×105
hsa05216: Thyroid cancer	7.99×106	7.99×106
hsa05222: Small cell lung cancer	7.59×106	7.59×106
hsa05145: Toxoplasmosis	7.15×105	7.15×105

FDR, false discovery rate.

After identification and assessment, several sub-pathways were screened, including ribosome biogenesis in eukaryotes, and the MAPK and ErbB signaling pathways (Table III). In the top 3 pathways, various lncRNA and related genes were involved. For example, YLP motif containing 1 (YLPM1) in ribosome biogenesis pathway was closely associated with heat repeat containing 1 (HEATR1), WD repeat domain 3 and spermatogenesis associated 5 (Fig. 2). Besides, HOXA transcript antisense RNA, myeloid-specific 1 in MAPK signaling pathway was regulated by AKT serine/threonine kinase 3 (AKT3), ribosomal protein S6 kinase A5 and mitogen-activated protein kinase kinase 6 (Fig. 3). Furthermore, YLPMI in the ErbB signaling pathway was associated with MAP2K1 and AKT3 (Fig. 4).

Table III.

The top 5 hub-pathways.

Pathway_id	Pathway_name	Molecule ratio (m2/x)	BgRatio (m1/n)	Weight	P-value	FDR
03008_1	Ribosome biogenesis in eukaryotes	15/489	19/26232	1.341037	<0.001	0
04010_1	MAPK signaling pathway	38/489	51/26232	1.387023	<0.001	0
04012_1	ErbB signaling pathway	18/489	25/26232	1.556393	<0.001	0
04066_2	HIF-1 signaling pathway	20/489	25/26232	1.395929	<0.001	0
04110_1	Cell cycle	26/489	35/26232	1.428843	<0.001	0

FDR, false discovery rate.

Figure 2.

Pathway of ribosome biogenesis in eukaryotes related to lncRNAs and genes. The blue and yellow nodes represent lncRNAs and genes, respectively. The dark and light nodes represent nodes in and out of pathways, respectively. Edges represent regulated relationship between lncRNAs and genes. LncRNA, long non-coding RNAs.

Figure 3.

MAPK signaling pathway-related lncRNAs and genes. The blue and yellow nodes represent lncRNAs and genes, respectively. The dark and light nodes represent nodes in and out of pathways, respectively. Edges represent regulated relationship between lncRNAs and genes. LncRNA, long non-coding RNAs.

Figure 4.

ErbB signaling pathway-related lncRNAs and genes. The blue and yellow nodes represent lncRNAs and genes, respectively. The dark and light nodes represent nodes in and out of pathways, respectively. Edges represent regulated relationship between lncRNAs and genes. LncRNA, long non-coding RNAs.

Identification of hub-lncRNA

Based on the above results, the network of pathway-lncRNA-mRNA was constructed (Fig. 5). Through topological analyses of this network, the hub-genes were identified, including chromosome 14 open reading frame 169 (C14orf169), (YLPM1), testis-specific transcript, Y-linked 15 (TTTY15), PCBP1 antisense RNA 1 (PCBP1-AS1) and deleted in lymphocytic leukemia 2. More importantly, several hub-lncRNAs were obtained, including PRKAG2 antisense RNA 1, LEF1 antisense RNA 1, EMX2 opposite strand/antisense RNA (EMX2OS) and BOLA3 antisense RNA 1.

Figure 5.

The network of pathway-lncRNA-mRNA. The red, yellow, light blue, dark blue and green nodes represent hublnc, lnc, mRNA, mRNA and path, respectively. The size of node represent degree of node. Edges represented regulated relationship among hublnc, lnc, mRNA, mRNA& and path.

Discussion

Increasing evidence has shown that lncRNAs were important factors for regulating gene expression (13). However, the molecular mechanism in burn injuries including lncRNA functions and regulatory genes remain unknown. In the present study, an effective method for identification of lncRNA-regulated functions was presented. LncRNA-regulated subpathways were identified, including ribosome biogenesis in eukaryotes, as well as MAPK and ErbB signaling pathways. More importantly, several hub-lncRNAs were identified, including C14orf169, YLPM1, TTTY15, and PCBP1-AS1. In addition, the sub-pathways regulated by the lncRNAs were identified, and functions of these lncRNAs were predicted.

C14orf169 was found to regulate various genes in this study, including eukaryotic translation initiation factor 1 (EIF1) and PPP2R5C. Thus EIF1 was mainly present in nucleus and cytosol. In a previous study, EIF1 was confirmed to be related to the generation of stable 40S preinitiation complex (13). In addition, the eukaryotic translation initiation factor contributed to muscle regeneration by participating in the skeletal muscle stem cell differentiation (14). Besides, PPP2R5C was found to be overexpressed in myeloid leukemia cells, and closely related to leukemic cell proliferation and differentiation (15). In the present study, PPP2R5C was enriched into the PI3K-Akt and Wnt signaling pathways. Consistent with the study of Shi et al (16), anti-fibrotic effects by activation of the PIK3-Akt signaling pathway was beneficial for the reduction and prevention of scar formation. Furthermore, alteration of the Wnt signaling pathway could induce muscle fibrosis, and even lead to failure in muscle generation and tissue disorganization (17). Thus, C14orf169 served as a potential key lncRNA for the recovery of severe burn injury and post-operative complications by regulating EIF1 and PPP2R5C.

Moreover, YLPM1 was screened as a hub-lncRNA in the present study and resulted in the regulation of various genes, such as RPL10L, PPT1, HEATR1 and RPL9. Both RPL10L and RPL9 encoded ribosomal proteins, and were mainly present in cytosol and nucleus. Kim et al (18) found that RPL9 interacted with Ftn-2, and played an important role in the maintenance of iron homeostasis. In addition, ribosomal protein phosphorylation was confirmed to interact with myostatin, and further affect the mature muscle (19). In mouce models, knockout PPT1 induced neuronal ceroid lipofuscinosis (20). In several diseases, HEATR1 was confirmed to regulate functional cytotoxic T lymphocytes and enhance cell immunity (21,22). In the present study, HEATR1 was also found to enrich into the hub-pathway of ribosome biogenesis in eukaryotes. Based on this information, we inferred that YLPM1 was an important lncRNA after severe burn injuries by regulating RPL10L, PPT1, HEATR1 and RPL9.

TTTY15 has been widely researched and found to have functions in the pathogenesis of prostate cancer, but it is less known regarding other diseases. In the present study, the gene was found to regulate the expression of different genes including KCTD13, UFL1 and ZNF322. In zebrafish models, KCTD13 was found to be closely associated with neurodevelopmental phenotype (23). In addition, KCTD13 interacted with Bacurd1, and further influenced the dendritic maturation and long-term positioning of cerebral neurons (24). It is well known that UFL1 is mainly involved in ufmylation of ubiquitin fold modifier 1 targets, which was necessary for erythroid differentiation (25,26). Besides, ZNF322 is a crucial gene that regulates transcriptional activation in MAPK signaling pathways (27). Consequently, TTTY15 was important for recovery of severe burn injury by regulating KCTD13, UFL1 and ZNF322.

PCBP1-AS1, which was closer to the gene of PCBP1, resulted in the regulation of the expression of various genes including ZAP70, SOCS5 and SIRT1. Moreover, these regulated genes were involved in the T-cell receptor and Jak-STAT signaling pathways. In 2005, Wu et al (21) found that ZAP70 activated the generation of T-cell receptor microclusters and sustained T-cell activation, and participated in the immune response. Furthermore, SOCS5 was confirmed to be a specific inhibitor of IL-4 signaling, which was important for eosinophil infiltration (28). Of note, Schug et al (29) found that SIRT1 had a beneficial role in the treatment of inflammation and other associated diseases. Thus, PCBP1-AS1 may be involved in the mechanism following severe burn injury by regulating, the immune and inflammation process.

In conclusion, by integrating the expression profiles of lncRNA and mRNA, several hub-lncRNAs including C14orf169, YLPM1, TTTY15, and PCBP1-AS1 were screened. In addition, the sub-pathways regulated by these lncRNAs were identified, and functions of these lncRNAs were predicted.

Competing interests

The authors declare that they have no competing interests.