Integrated Analysis of lncRNA-miRNA-mRNA ceRNA Network and the Potential Prognosis Indicators in Sarcomas

BMC Medical Genomics 2021 AI 8 Explanations View Original
Original Paper (PDF)

Unable to display PDF. Download it here or view on PMC.

Plain-English Explanations
Pages 1-2
Why Noncoding RNA Networks Matter for Sarcoma Prognosis

Sarcomas are rare malignant tumors of mesenchymal origin, meaning they arise from connective tissues including bone, cartilage, fat, muscle, and fibrous tissue. The WHO classification recognizes more than 100 distinct sarcoma subtypes, broadly divided into bone sarcomas (such as osteosarcoma, Ewing sarcoma, and chondrosarcoma) and soft tissue sarcomas (such as synovial sarcoma, liposarcoma, and rhabdomyosarcoma). Despite aggressive multimodal treatment including surgical resection, chemotherapy, and radiotherapy, five-year overall survival rates remain poor for many subtypes, particularly in metastatic and recurrent disease. Novel molecular biomarkers are urgently needed to stratify patients by prognosis and identify new therapeutic targets.

The ceRNA hypothesis: Competitive endogenous RNAs (ceRNAs) represent a newly appreciated regulatory layer in cancer biology. The ceRNA framework proposes that messenger RNAs (mRNAs), long noncoding RNAs (lncRNAs), and other RNA species compete for a shared pool of microRNAs (miRNAs) by harboring complementary miRNA-binding sites called microRNA response elements (MREs). When lncRNAs are upregulated in tumor cells, they can act as molecular "sponges," sequestering specific miRNAs and preventing those miRNAs from binding to and repressing their target mRNAs. The net effect is increased expression of oncogenic mRNAs that would otherwise be silenced. Conversely, downregulated lncRNAs free up miRNAs to suppress their target mRNAs more aggressively, potentially reducing expression of tumor-suppressive genes.

Prior ceRNA work in other cancers: The ceRNA paradigm has been applied productively in lung adenocarcinoma, breast cancer, gastric cancer, esophageal cancer, cholangiocarcinoma, ovarian cancer, and endometrial cancer, yielding networks that identify prognostically relevant RNA species and potential drug targets. In sarcoma specifically, individual lncRNA studies have described roles for NEAT1, XIST, SOX2OT, and H19 in Ewing sarcoma and osteosarcoma, but no comprehensive integrative ceRNA network analysis spanning multiple GEO datasets had been reported at the time of this study.

This 2021 paper published in BMC Medical Genomics addresses that gap by constructing a genome-wide lncRNA-miRNA-mRNA ceRNA network for sarcomas using publicly available GEO microarray datasets, then validating the prognostic relevance of network components against TCGA survival data. The goal is to generate a curated list of candidate prognostic biomarkers and potential therapeutic targets grounded in the sarcoma-specific regulatory landscape.

TL;DR: Sarcomas encompass 100+ WHO subtypes with poor prognosis in metastatic disease. This study applies the ceRNA framework, where lncRNAs sponge miRNAs to regulate mRNA expression, to build a comprehensive lncRNA-miRNA-mRNA interaction network in sarcoma using GEO microarray data, then identifies prognostically relevant RNAs against TCGA survival outcomes.
Pages 2-3
Dataset Selection, Differential Expression Analysis, and the Bioinformatics Pipeline

The study began with a systematic search of the Gene Expression Omnibus (GEO) database for high-throughput gene or miRNA expression profiling datasets from sarcoma patients and normal tissue controls. Four candidate datasets were initially identified: GSE55625, GSE31045, GSE17674, and GSE18546. After quality filtering, two datasets were excluded because they contained multiple zero and negative expression values that would compromise statistical analysis. The two retained datasets were GSE17674, based on the Affymetrix Human Genome U133 Plus 2.0 Array platform, and GSE18546.

Dataset composition: GSE17674 contained expression profiles from 32 Ewing sarcoma tumor samples and 18 normal muscle tissue samples. This dataset served as the primary source for identifying differentially expressed mRNAs (DEGs) and lncRNAs (DELs). Because the Affymetrix U133 Plus 2.0 Array was originally designed to detect protein-coding transcripts, an additional lncRNA annotation step was required. Transcription sequences beginning with NM and XM prefixes in the RefSeq database were designated as the mRNA reference database, while noncoding RNAs with accession numbers beginning with NR were defined as lncRNAs. Any probes that could not be mapped to Ensembl gene IDs were discarded, yielding a cleaned probe-to-lncRNA mapping. GSE18546 included 10 synovial sarcoma samples, 5 Ewing sarcoma samples, and 5 normal muscle samples, and served as the source for identifying differentially expressed miRNAs (DEMs).

Differential expression analysis: DEGs, DELs, and DEMs between sarcoma and normal muscle samples were computed using the limma package in R version 3.3.2, a linear models framework designed for microarray and RNA-seq data. The significance thresholds applied were P less than 0.05, false discovery rate (FDR) less than 0.05, and fold change greater than 3 for DEGs and DELs. For DEMs, the same P value and FDR thresholds were used but with a fold change cutoff of greater than 2. Hierarchical clustering was applied using EPCLUST to confirm sample groupings and visualize expression patterns across all detected differentially expressed RNA species.

Functional enrichment: Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed on all 3,415 DEGs to identify overrepresented biological processes, molecular functions, and metabolic pathways. Fisher's exact test was used to assess statistical enrichment. The results of GO and KEGG analyses were then used to filter the DEG list, retaining only those genes assigned to significantly enriched biological categories, which yielded 1,296 intersecting DEGs for downstream ceRNA network construction.

TL;DR: Two GEO datasets were analyzed: GSE17674 (32 Ewing sarcoma vs. 18 normal muscle) for DEGs and DELs, and GSE18546 (10 synovial sarcoma + 5 Ewing sarcoma vs. 5 normal) for DEMs. Limma in R identified differentially expressed RNAs at P < 0.05, FDR < 0.05, and fold change > 3 (mRNA/lncRNA) or > 2 (miRNA). GO and KEGG enrichment filtered 3,415 DEGs down to 1,296 for network construction.
Pages 3-4
Scale and Composition of Differentially Expressed RNAs in Sarcoma

Applying the defined statistical thresholds to the GSE17674 dataset produced a large set of differentially expressed coding transcripts. A total of 3,415 mRNAs were significantly dysregulated in sarcoma relative to normal muscle, of which 2,554 were upregulated and 861 were downregulated. The predominance of upregulated transcripts reflects the widespread activation of proliferative, cell cycle, and oncogenic gene programs in these mesenchymal tumors. Alongside these protein-coding changes, 338 lncRNAs were differentially expressed: 234 upregulated and 104 downregulated.

miRNA differential expression: The miRNA analysis of GSE18546 was conducted separately for two sarcoma subtypes because the dataset contained samples from both Ewing sarcoma and synovial sarcoma. In Ewing sarcoma versus normal muscle, 52 miRNAs were differentially expressed: 39 upregulated and 13 downregulated. In synovial sarcoma versus normal muscle, 145 miRNAs were differentially expressed: 109 upregulated and 36 downregulated. The study then applied a Venn diagram intersection across both DEM sets to identify shared miRNAs, yielding 26 upregulated DEMs and 10 downregulated DEMs common to both subtypes. These 36 shared miRNAs were used for all downstream ceRNA network analyses, on the reasoning that they represent a conserved miRNA dysregulation signature across at least two sarcoma subtypes.

Functional enrichment of DEGs: GO and KEGG analyses of the 3,415 DEGs revealed that upregulated genes were most significantly enriched in transcription-related processes. The top three enriched upregulated GO terms were transcription, DNA-templated (GO:0006351; P = 2.64 x 10-104), cell division (GO:0051301; P = 2.57 x 10-82), and regulation of transcription, DNA-templated (GO:0006355; P = 7.99 x 10-75). Downregulated DEGs were significantly enriched in muscle contraction and striated muscle tissue development processes, consistent with the normal muscle origin of the control tissue and the loss of muscle-specific differentiation programs in sarcoma transformation. The final filtered set of 1,296 DEGs retained after GO/KEGG intersection formed the pool of candidate mRNAs for ceRNA network integration.

Heatmap visualization confirmed clear separation between sarcoma tumor and normal muscle samples across mRNA, lncRNA, and miRNA expression profiles, validating the robustness of the differential expression analysis and the quality of both input datasets.

TL;DR: GSE17674 yielded 3,415 DEGs (2,554 up, 861 down) and 338 DELs (234 up, 104 down). Across Ewing and synovial sarcoma subtypes in GSE18546, 36 shared miRNAs were identified (26 up, 10 down). Top upregulated GO terms: transcription (P = 2.64 x 10-104) and cell division (P = 2.57 x 10-82). GO/KEGG filtering produced 1,296 final DEGs for network input.
Pages 4-6
Building the lncRNA-miRNA-mRNA Interaction Network

Constructing the ceRNA network required two parallel prediction steps: identifying which of the 36 DEMs target which of the 1,296 DEGs, and identifying which lncRNAs from the 338 DELs are targeted by those same miRNAs. For miRNA-mRNA target prediction, three complementary computational algorithms were applied in combination: miRanda, TargetScan, and miRWalk. Each algorithm uses different sequence-based features and scoring systems to predict miRNA binding to mRNA 3' untranslated regions (3' UTRs). Only target gene predictions that overlapped with genes already identified as significant in the GO/KEGG enrichment analyses were retained, adding a biological plausibility filter to reduce false positive predictions. This step yielded 448 miRNA-mRNA interactor pairs, encompassing 34 miRNAs and 269 mRNAs.

lncRNA target prediction: For lncRNA-miRNA interactions, two algorithms were employed: miRanda and PITA (Probability of Interaction by Target Accessibility). PITA specifically accounts for the accessibility of target sites within the RNA secondary structure, which is a determinant of miRNA binding efficiency. Predicted lncRNA-miRNA pairs were filtered to include only those where the lncRNA was already identified as a DEL in the sarcoma versus normal comparison. This step produced 454 lncRNA-miRNA interaction pairs, covering 36 miRNAs and 117 lncRNAs.

ceRNA selection logic: The ceRNA theory requires that lncRNAs and mRNAs should be regulated in an inverse direction relative to their shared miRNA sponge. Specifically, when a miRNA is upregulated, its competing lncRNA and mRNA targets should be downregulated (because the miRNA is actively repressing them). Conversely, when a miRNA is downregulated, its competing lncRNA and mRNA targets should be upregulated (because the miRNA is no longer effectively repressing them). Applying this inverse expression constraint to the predicted interaction pairs filtered the network down to biologically coherent ceRNA triplets. The final ceRNA network contained 1,440 total interactions, incorporating 29 miRNAs, 69 lncRNAs, and 113 mRNAs. The network was visualized using Cytoscape version 2.8.2, with rectangles representing miRNAs, circles representing mRNAs, and triangles representing lncRNAs, and red and green coloring indicating upregulated and downregulated nodes, respectively.

Network scale in context: The 1,440-interaction network represents a substantial regulatory architecture, with each miRNA node connecting to multiple lncRNA and mRNA partners. This topology is consistent with ceRNA networks described in other cancer types, where a relatively small number of miRNA "hubs" coordinate the regulation of dozens of coding and noncoding RNA species simultaneously.

TL;DR: miRNA-mRNA pairs were predicted using miRanda, TargetScan, and miRWalk, yielding 448 pairs (34 miRNAs, 269 mRNAs). lncRNA-miRNA pairs from miRanda and PITA produced 454 pairs (36 miRNAs, 117 lncRNAs). After applying the inverse-expression ceRNA constraint, the final network contained 1,440 interactions across 29 miRNAs, 69 lncRNAs, and 113 mRNAs, visualized in Cytoscape.
Pages 6-7
PPI Network Analysis Identifies Six Key Hub Proteins

To further prioritize the 113 mRNAs included in the ceRNA network, the authors performed protein-protein interaction (PPI) network analysis using the STRING database version 11.0. STRING integrates experimental evidence, co-expression data, genomic context, and curated pathway information to assign confidence scores to predicted protein interactions. Only protein pairs with a combined STRING interaction score above 0.4 (the medium-confidence threshold) were included in the PPI network. Cytoscape was again used for visualization, with nodes colored red for upregulated proteins and green for downregulated proteins.

Hub protein identification: The PPI network analysis identified six hub proteins representing the most highly connected nodes in the interaction graph. Among upregulated hub proteins, the three most significant were IGF1 (insulin-like growth factor 1), PRKCB (protein kinase C beta), and GNAI3 (G protein subunit alpha i3). Among downregulated hub proteins, the three most significant were AR (androgen receptor), CYCS (cytochrome c, somatic), and PPP1CB (protein phosphatase 1 catalytic subunit beta).

Biological roles of hub proteins: IGF1 is a well-characterized growth factor that activates the PI3K-AKT and MAPK signaling cascades, both of which are frequently dysregulated in sarcomas. IGF1 receptor (IGF1R) has been an active therapeutic target in Ewing sarcoma, with multiple clinical trials testing IGF1R inhibitors, though with limited success in unselected patient populations. PRKCB is a serine-threonine kinase involved in cell proliferation, survival, and angiogenesis signaling. GNAI3 encodes a G-protein alpha subunit that couples to G protein-coupled receptors and modulates downstream signaling cascades. The downregulation of CYCS is notable because cytochrome c is a central mediator of intrinsic apoptosis, and its reduced expression may contribute to the apoptotic resistance characteristic of aggressive sarcomas.

PPP1CB encodes a catalytic subunit of protein phosphatase 1, a serine-threonine phosphatase with broad regulatory functions in cell cycle control, DNA damage response, and RNA splicing. AR, though primarily studied in prostate cancer, has been implicated in sarcoma biology, particularly in leiomyosarcoma. The identification of these six hub proteins as network nodes provides prioritized targets for functional validation studies and potential therapeutic intervention in sarcoma.

TL;DR: PPI analysis via STRING (score > 0.4) identified six hub proteins. Upregulated hubs: IGF1 (a PI3K-AKT/MAPK activator and validated Ewing sarcoma target), PRKCB (proliferation and angiogenesis kinase), and GNAI3 (G-protein signaling). Downregulated hubs: AR (androgen receptor), CYCS (apoptosis mediator), and PPP1CB (phosphatase 1 catalytic subunit).
Pages 7-9
Twelve ceRNA Network Components Significantly Predict Sarcoma Overall Survival

Survival analysis was performed by cross-referencing the ceRNA network components against high-throughput expression data with clinical outcome profiles from The Cancer Genome Atlas (TCGA) sarcoma cohort. Log-rank tests were used to compare overall survival between patients with high versus low expression of each RNA species in the network, with P less than 0.05 considered statistically significant. Of all mRNAs, miRNAs, and lncRNAs included in the ceRNA network, 12 RNA species were significantly associated with overall survival of sarcoma patients: seven mRNAs, four miRNAs, and one lncRNA.

Survival-associated mRNAs: High expression of SMARCC1 (SWI/SNF related matrix associated actin dependent regulator of chromatin subfamily C member 1), SRSF10 (serine and arginine rich splicing factor 10), PRPF38A (pre-mRNA processing factor 38A), JARID2 (Jumonji and AT-rich interaction domain containing 2), and GNAI3 were each significantly associated with shorter overall survival in sarcoma patients (all P values below 0.05). Conversely, high expression of ARF3 (ADP-ribosylation factor 3) and PRKCB were associated with shorter overall survival through lower expression patterns relative to better-prognosis patients.

Survival-associated miRNAs: Four miRNAs in the ceRNA network were significantly associated with survival outcomes. High expression of miR-301a-3p (P less than 0.0001), miR-106b-5p (P = 0.0046), miR-130b-3p (P = 0.0128), and miR-423-3p (P = 0.0147) were each associated with shorter overall survival. The miR-301a-3p association was the most statistically significant among all network components examined. Prior literature supports a functional role for miR-301a in Ewing sarcoma: overexpression of miR-301a has been demonstrated in Ewing sarcoma cell lines, and transfection with anti-miR-301a inhibited proliferation and cell cycle progression in those cells.

Survival-associated lncRNA: LINC01296 was the only lncRNA in the ceRNA network that reached statistical significance for overall survival. High expression of LINC01296 was associated with shorter overall survival in sarcoma patients. LINC01296 has been reported in other cancer types as a promoter of tumor progression, but this study provides one of the first demonstrations of its prognostic relevance specifically in sarcoma.

TL;DR: Log-rank tests against TCGA sarcoma data identified 12 prognostic RNAs: 7 mRNAs (including SMARCC1, SRSF10, JARID2, GNAI3), 4 miRNAs (miR-301a-3p at P < 0.0001, miR-106b-5p at P = 0.0046, miR-130b-3p at P = 0.0128, miR-423-3p at P = 0.0147), and 1 lncRNA (LINC01296). High expression of all 12 was linked to shorter overall survival.
Pages 9-10
Molecular Roles of Key ceRNA Nodes in Sarcoma Biology

The study discusses the biological plausibility of the identified prognostic biomarkers within known sarcoma molecular pathways. SMARCC1, a core subunit of the SWI/SNF chromatin remodeling complex, is involved in epigenetic regulation of gene expression through ATP-dependent repositioning of nucleosomes. SWI/SNF subunits are among the most frequently mutated genes in human cancer, and dysregulation of chromatin remodeling has been specifically implicated in synovial sarcoma, where the SS18-SSX fusion oncogene directly disrupts SWI/SNF complex function. The association of high SMARCC1 expression with poor survival in the TCGA sarcoma cohort adds a quantitative prognostic dimension to the known qualitative role of SWI/SNF in sarcoma pathogenesis.

JARID2 and epigenetic regulation: JARID2 is a Polycomb group protein and a regulatory subunit of the Polycomb repressive complex 2 (PRC2), which catalyzes trimethylation of histone H3 at lysine 27 (H3K27me3), a repressive chromatin mark. JARID2 modulates PRC2 activity and target gene selection. PRC2 has been implicated in sarcoma aggressiveness, and its catalytic subunit EZH2 has been actively investigated as a therapeutic target in sarcoma subtypes including SMARCB1-deficient malignant rhabdoid tumors and epithelioid sarcomas. The identification of JARID2 as a survival-associated node in the ceRNA network suggests that PRC2-dependent epigenetic silencing may be a broader driver of poor prognosis across sarcoma types.

ceRNA interactions involving miR-301a-3p: Within the ceRNA network, miR-301a-3p participates in three documented lncRNA sponge interactions: NEAT1/miR-301a-3p, XIST/miR-301a-3p, and another lncRNA-miRNA pair. Both NEAT1 and XIST are large, well-characterized nuclear lncRNAs with documented oncogenic and tumor-suppressive roles depending on context. The authors propose that in sarcoma, these lncRNAs sequester miR-301a-3p, preventing it from suppressing its mRNA targets and thereby indirectly promoting tumor growth. The combination of statistical significance at P less than 0.0001 for overall survival and the functional literature on miR-301a in Ewing sarcoma makes this node a particularly credible candidate for future experimental validation.

miR-423-3p and miR-130b-3p: These two miRNAs have documented roles in other cancer types as regulators of apoptosis and proliferation signaling pathways, including through the PI3K-AKT and TGF-beta pathways. Their presence in the sarcoma ceRNA network at statistically significant survival thresholds extends their potential relevance into the sarcoma context and motivates in vitro validation studies using sarcoma cell lines to confirm their functional contributions to tumor biology.

TL;DR: Key prognostic nodes connect to established sarcoma biology: SMARCC1 is a SWI/SNF subunit disrupted by the SS18-SSX fusion oncogene in synovial sarcoma; JARID2 modulates PRC2 epigenetic silencing (related to EZH2 targeting in epithelioid sarcoma); miR-301a-3p is sponged by NEAT1 and XIST lncRNAs and experimentally validated in Ewing sarcoma cell lines.
Pages 10-11
Study Constraints and the Path to Clinical Translation

Subtype coverage: The authors explicitly acknowledge that sarcomas are an extraordinarily heterogeneous disease group spanning more than 100 distinct subtypes, yet the ceRNA network in this study was constructed using data from only two subtypes: Ewing sarcoma and synovial sarcoma. The GEO datasets with acceptable data quality (GSE17674 and GSE18546) happened to contain only these two histological types, and the differential expression profiles are therefore most directly applicable to Ewing sarcoma and synovial sarcoma rather than sarcomas broadly. Extending this analysis to additional subtypes such as osteosarcoma, rhabdomyosarcoma, liposarcoma, leiomyosarcoma, and undifferentiated pleomorphic sarcoma would require larger and more diverse expression datasets. Many sarcoma subtypes are rare enough that suitable public GEO datasets simply do not exist yet at sufficient sample sizes.

Complexity of ceRNA regulation: The ceRNA framework assumes that changes in lncRNA abundance meaningfully alter miRNA availability for target mRNAs. However, ceRNA crosstalk is quantitatively sensitive to the relative abundance of all three RNA species involved. If a miRNA is expressed at very high levels, sponging by a single lncRNA may have negligible functional impact. Additionally, the subcellular compartmentalization of ceRNA components matters: a nuclear lncRNA cannot effectively sponge a cytoplasmic miRNA. These context-dependent constraints mean that not all predicted ceRNA interactions will function equivalently in living tumor cells, and computational predictions require experimental validation to confirm functional relevance.

Purely bioinformatic study: All findings in this paper are derived from reanalysis of publicly available microarray data and cross-referencing with TCGA survival information. No in vitro experiments (cell line knockdown, overexpression, or reporter assays) or in vivo experiments (xenograft models) were performed to validate the predicted interactions or confirm the causal relationships between ceRNA network components and sarcoma behavior. The survival associations identified via log-rank tests are observational correlations that are consistent with the ceRNA hypothesis but do not establish causality.

Future directions: The authors identify several important next steps. Validation of the 12 survival-associated RNA species as functional biomarkers using independent sarcoma patient cohorts with complete clinical annotation would strengthen confidence in their prognostic utility. Functional experiments using CRISPR-based knockout or RNA interference targeting key miRNAs (especially miR-301a-3p), lncRNAs (especially LINC01296, NEAT1, and XIST), and mRNA hub nodes in sarcoma cell lines would confirm their contribution to proliferation, invasion, and drug resistance phenotypes. Integration of additional omics layers including copy number alterations, DNA methylation, and mutational profiles could refine the ceRNA network and reveal synthetic lethal dependencies exploitable for therapy.

TL;DR: Key limitations include subtype restriction to only Ewing and synovial sarcoma, the theoretical complexity of ceRNA abundance thresholds and compartmentalization, and the purely computational design with no experimental validation. Future work requires independent cohort validation, functional cell-line experiments for miR-301a-3p and LINC01296, and multi-omic integration across a broader range of sarcoma histologies.