Co-expression analysis identifies long noncoding RNA SNHG1 as a novel predictor for event-free survival in neuroblastoma

Abstract

Despite of the discovery of protein therapeutic targets and advancement in multimodal therapy, the survival chance of high-risk neuroblastoma (NB) patients is still less than 50%. MYCN amplification is a potent driver of NB, which exerts its oncogenic activity through either activating or inhibiting the transcription of target genes. Recently, long noncoding RNAs (lncRNAs) are reported to be altered in cancers including NB. However, lncRNAs that are altered by MYCN amplification and associated with outcome in high-risk NB patients are limitedly discovered. Herein, we examined the expression profiles of lncRNAs and protein-coding genes between MYCN amplified and MYCN non-amplified NB from microarray (n = 47) and RNA-seq datasets (n = 493). We identified 6 lncRNAs in common that were differentially expressed (adjusted P ≤ 0.05 and fold change ≥ 2) and subsequently validated by RT-qPCR. The co-expression analysis reveals lncRNA, SNHG1 and coding gene, TAF1D highly co-expressed in NB. Kaplan-Meier analysis shows that higher expression of SNHG1 is significantly associated with poor patient survival. Importantly, multivariate analysis confirms high expression of SNHG1 as an independent prognostic marker for event-free survival (EFS) (HR = 1.58, P = 2.36E-02). In conclusion, our study unveils that SNHG1 is up-regulated by MYCN amplification and could be a potential prognostic biomarker for high-risk NB intervention.

Keywords: SNHG1; co-expression study; event-free survival; long noncoding RNAs; neuroblastoma.

Figure 1. Heatmap of differentially expressed lncRNAs and biological functions of de-regulated coding genes in NB

(A) Heatmap showing the expression values of lncRNAs found differentially expressed in MYCN amplified compared with MYCN non-amplified patient samples. Each column indicates a MYCN amplified (grouped in red bar) or non-amplified (green) patients samples. Each row represents the lncRNA ordered by hierarchical clustering analysis. Expression value of each lncRNA was scaled across samples and represented in a blue-red color scale. (B) Functional enrichment map of differentially expressed coding genes. Nodes represent GO terms and edges represent genes shared between GO terms (kappa score threshold = 0.4). Edge color gradient represents kappa score. Overview GO terms are shown by the highest statistically enriched term in each group formed by fusion of the GO terms sharing similar genes. Clusters with orange color are enriched functions for up-regulated genes while blue clusters are for down-regulated genes. Node color gradient represents the proportion of up (orange) or down-regulated (blue) genes associated with the term. A node (GO term) with equal proportion of up and down-regulated genes is represented in grey color.

Figure 2. Expression profile of potential lncRNAs

(A) Boxplot showing six lncRNAs expression profiles in MYCN amplified (n = 92) and MYCN non-amplified (n = 401). (B) Bar charts of expression levels measured by RT-qPCR in NB cell lines. The expression levels were normalized to endogenous GAPDH and relative to SK-N-F1 and displayed in log2 scale. MYCN and TAF1D are also shown for comparison.

Figure 3. Correlation of lncRNA and coding gene expression profiles

Fisher's Z-transformed score of Spearman's correlation coeffiecient (SCC) between lncRNAs (n = 13) and coding genes (n = 591) in (A) MYCN amplified and (B) MYCN non-amplified neuroblastoma tumor samples. Dendrograms of coding genes (rows) are displayed.

Figure 4. Co-expression network of lncRNAs and protein-coding genes in MYCN amplified neuroblastoma

Nodes represent the lncRNAs and coding genes whereas edges represent the z-scores of expression correlation between lncRNAs and coding genes. Red nodes represent up-regulated lncRNAs and blue ones represent down-regulated lncRNAs. Yellow nodes represent protein-coding genes. Node size and edge width are proportional to the degree of a node and z-score, respectively. Dashed edges indicate that the connected lncRNA and coding gene pairs also appeared in the MYCN non-amplified network.

Figure 5. Positive correlation between SNHG1 and TAF1D expression

(A) Scatter plot of SNHG1 and TAF1D expression levels in neuroblastoma patients measured by RNA-seq (n = 493) shows positive correlation. (B) The positive correlation was confirmed in NB cell lines via RT-qPCR (r = 0.75).

Figure 6. SNHG1 is a prognostic marker for neuroblastoma

(A, B) Boxplots showing the normalized log2RPM expression values of SNHG1 in different risk groups and stages in a NB cohort (n = 493). The P-values presented were determined by Mann-Whitney-Wilcoxon test (A) and Dunn's multiple comparison test (B). (C, D) Scatter plots showing the correlation between MYCN and SNHG1 in MYCN non-amplified (n = 401) and MYCN amplified patients (n = 92). SCC and the corresponding P-values are displayed. (E) Bar chart showing the ordered expression levels across 16 normal human tissues, based on the RNA-seq data from the Illumina Body Map project. (F) Bar chart showing the ordered expression levels of SNHG1 per survival status of the patient. Here, blue and red bars represent patients did not die and died of disease, respectively.

Figure 7. Kaplan-Meier survival analysis for NB patients

The Kaplan-Meier plots for (A) event-free survival (EFS) and (B) overall survival (OS) of low-expression versus high-expression groups based on the median SNHG1 expression level of GSE62564 patients (n = 493). (C) Kaplan-Meier curves of low-expression (n = 401) versus high-expression group (n = 92) based on the ordered expression of SNHG1. (D) Kaplan-Meier plot curve for EFS of GSE16476 patients (n = 88) (E) Kaplan-Meier plot for EFS of low-expression versus high-expression group based on the median expression of TAF1D (n = 493). (F) Kaplan-Meier plot for EFS of combinatorial low and high expression of both SNHG1 and TAF1D (n = 493). The P-values were obtained by log-rank (Mantel-Haenszel) test.