Characterization of long non-coding RNA and messenger RNA profiles in laryngeal cancer by weighted gene co-expression network analysis

Construction of co-expression modules of LC

The 5000 genes (including 2500 lncRNA and 2500 mRNA expression data) were normalized by Limma package with Voom function. The auxiliary data was then removed and the expression data was transposed for further WGCNA analysis. The expression values of top 2500 lncRNAs in 99 LC samples (Supplementary Table 1) and top 2500 mRNAs in 99 LC samples (Supplementary Table 2) were used to develop co-expression modules with WGCNA package. The clinical characteristics information of these LC patients is listed in Supplementary Table 3, including age at initial pathologic diagnosis, history of smoking, history of alcohol consumption, intermediate dimension of tumor, lymph node count, neck lymph node dissection, pathologic N stage, radiation therapy, targeted molecular therapy, tumor status, success of follow-up treatment, and overall survival. The red line (cut height = 50) was the filter which we used to remove outlier samples in sampleTree. The TCGA-KU-A66S sample was excluded after removing outliers in the sample based on gene expression data (Figure 1A). The sample dendrogram and trait heatmap grouped the selected samples into the different clusters, and provided the distribution map of clinical trait data (Figure 1B). The independence and the average connectivity degree of the co-expression modules were decided by power value (β) and scale R² value. First, a set of soft-thresholding powers were plotted (Supplementary Table 4). When the power value was equal to 5, the scale R² was up to 0.87 (Figure 2A). Therefore, we define the adjacency matrix using soft thresholding with beta=5 to construct and identify distinct co-expression gene modules in LC samples. A cluster dendrogram of all selected genes was constructed based on a TOM-based dissimilarity measure. These identified co-expression modules were distributed in different colors (Figure 2B). The interactions of these co-expression modules were analyzed with the Pearson correlation coefficient (Figure 2C). The darker background indicated higher module correlation. Hierarchical clustering of module eigengenes summarizing the modules were found in the clustering analysis. Branches of the dendrogram (the meta-modules) were grouped together based on the correlation of eigengenes. Each row and column in the heatmap plot of topological overlap corresponded to one module eigengene labeled by a different color. Each module contains different gene clusters. In the heatmap, blue represented negative correlation, while red represented positive correlation. Squares of red along the diagonal were the meta-modules (Figure 2D).

Gene co-expression modules correspond to clinical traits

Principal component analysis of each module elected the first principal component as its eigengenes. Eigengene external traits (including age at initial pathologic diagnosis, history of smoking, history of alcohol consumption, intermediate dimension, lymph node count, neck lymph node dissection, pathologic N stage, radiation therapy, targeted molecular therapy, tumor status, success of follow-up treatment) were correlated with different co-expression modules and the most significant associations were identified (Figure 3). A heatmap of the correlation between module eigengenes and clinical traits of LC showed correlation coefficient (R) and significant difference (p value). The Module-trait relationships plot demonstrated that the co-expression Red module contained 121 genes, co-expression Yellow module contained 336 genes, co-expression Blue module contained 1123 genes, co-expression Pink module contained 73 genes, co-expression Black module contained 88 genes, co-expression Green module contained 301 genes, co-expression Brown module contained 372 genes, co-expression Turquoise module contained 1186 genes, and co-expression Grey module contained 1390 genes (Table 1). The hub genes of each module might be potential biomarkers for the specific clinical characteristics. The correlation analysis of gene co-expression module and clinical traits demonstrated that the co-expression Red module was significantly associated with documented alcohol history (R =0.23, p=0.02); the co-expression Yellow module was significantly associated with age at initial pathologic diagnosis(R =0.23, p=0.03), and the number of packs of cigarettes per year the patient had smoked (R=0.23, p=0.02); the co-expression Blue module was significantly associated with intermediate dimension of tumor (0.2 to 1.4 cm) (R=0.25, p=0.01), and success of follow-up treatment (R=0.27, p=0.007); the co-expression Green module was significantly associated with success of follow-up treatment (R=0.27, p=0.008); the co-expression Brown module was significantly associated with tumor status (R=0.27, p=0.008), and success of follow-up treatment (R=0.36, p=3e-04); the co-expression Turquoise module was significantly associated with pathologic N stage (R=0.27, p=0.008), and targeted molecular therapy (R=0.36, p=3e-04); and the co-expression Pink module was significantly associated with cigarette packs per year (R=-0.26, p=0.009), lymph node examined count (R=-0.24, p=0.02), tumor status (R=0.21, p=0.03), and success of follow-up treatment (R=0.42, p=1e-05). A scatterplot of gene significance (y-axis) vs. module membership (x-axis) was shown in the most significant module. Most interestingly, the scatterplot of gene significance (GS) vs. module membership (MM) was plotted in the co-expression Pink module. In modules related to a trait of interest, genes with high module membership often had high gene significance, implying that hub genes of the co-expression Pink module tend to be highly correlated with selected clinical characteristics. The results consistently revealed that MM in the Pink module significantly correlated with history of smoking (R=-0.28, p=0.016), lymph node count (R=-0.219, p=0.028), tumor status (R=0.263, p=0.003), and success of follow-up treatment (R=0.43, p=0.00015) (Figure 4A–4D). The correlation results in the co-expression Pink module showed consistency between module-trait relationships plot and the scatterplot of gene significance (GS) vs. module membership (MM) plot. Considering the correlation coefficient, p value, and consistency between module-trait relationships plot and the scatterplot, we chose the co-expression Pink module for further analysis.

Functional enrichment gene analysis in the co-expression Pink module

Based on the heatmap, we extracted all 73 of the co-expression Pink module genes (Supplementary Table 5) and hub genes (including IFIT2, XAF1, UBE2L6, IFITM3, HLA-C, CTSL, ARHGDIB, LGALS3BP, IFITM1, MLKL, SERPING1, TRIM21) (Supplementary Table 6). A total of 10 statistically significant signaling pathways were obtained by pathway analysis based on the co-expression Pink module (Table 2). The data revealed 10 distinct clusters which might possibly indicate that LC is potentially related to multiple signaling pathways (Figure 5). Cluster 1 contains cytokine signaling in the immune system, interferon signaling, immune system, interferon alpha/beta signaling, interferon alpha-beta signaling, and type II interferon signaling. Cluster 2 contains interferon gamma signaling, the herpes simplex infection pathway, and the Epstein-Barr virus infection pathway. Cluster 3 contains an endosomal/vacuolar pathway, antigen processing and presentation, antigen processing cross presentation, phagosome, and Ebola Virus pathway. Cluster 4 contains natural killer cell which mediate cytotoxicity, viral myocarditis, human cytomegalovirus infection, and human immunodeficiency virus 1 infection. Cluster 5 contains classical complement pathway, classical antibody-mediated complement activation, pertussis, complement activation, oxidative damage, systemic lupus erythematosus, initial triggering of complement, and creation of C4 and C2 activators. Cluster 6 contains regulated necrosis, RIPK1-mediated regulated necrosis, DNA damage response, and programmed cell death. Cluster 7 contains the human immune response to tuberculosis, B cell receptor signaling pathway, choline metabolism in cancer. Cluster 8 contains antigen presentation, autoimmune thyroid disease, allograft rejection, Graft-versus-host disease, Type-I diabetes mellitus, proteasome degradation, cellular senescence, viral carcinogenesis, human T-cell leukemia virus 1 infection, immunoregulatory interactions between a lymphoid and a non-lymphoid cell, and allograft rejection. Cluster 9 contains JAK-STAT molecular variation 1, and cytokine-cytokine receptor interaction. Cluster 10 contains rheumatoid arthritis and TNF signaling pathway. Most enriched genes of those identified pathways were cytokines, chemokines, histocompatibility antigens, complement factors, autoantigens and/or innate immune response molecules, all of which were thought to play a role in tumorigenesis, immunity, and gene regulation. In recent years, new molecules have been widely reported which are associated with cancer progression and immune escape and there are also strong correlations between antitumoral and immunomodulation effect [19]. The mechanisms that support evolution of the inflammatory environment and its relationship with neoplasms are unclear [20]. Tumor-associated inflammation is predictive of poor prognosis and associated with a variety of tumorigenic phenotypes, including angiogenesis, apoptosis, tumor proliferation and survival, autophagy, invasiveness, and metastasis [21]. Here, we identified an inflammation activation related pathway in the co-expression Pink module. The present understanding of this area is not yet comprehensive, but certain inflammatory pathways have emerged as important mediators of the crosstalk between tumor biology behavior and inflammation in tumors [22].

The 73 genes in the co-expression Pink module were uploaded to STRING for protein–protein interaction (PPI) analysis (Figure 6A). The combined scores of nodes ranged from 0.400 to 0.999. The PPI results also determined co-expression coefficient between proteins (Supplementary Table 7). High combined scores and co-expression coefficients indicated that there were interaction effects on spatial position and expression events between proteins. Some of interaction results showed both high combined scores (value > 0.9) and high co-expression coefficients (value > 0.8), such as RSAD2 and OASL, RSAD2 and XAF1, DDX60 and RSAD2, RSAD2 and IFIT2, IFITM2 and IFITM1, OAS1 and IFI6, C1R and C1S, IFITM2 and IFITM3 (Supplementary Table 7). Some interaction results showed high combined scores (value > 0.9), but low co- expression coefficients (value = 0), including TRIM5 and HLA-C, HLA-G and IFI6, HLA-G and IFITM3, IFITM2 and HLA-G, HLA-G and IFITM1, HLA-G and XAF1, HLA-G and OAS2, HLA-G and IFIT2, HLA-G and TRIM5, HLA-G and TRIM21, HLA-G and TRIM22, HLA-C and XAF1, TRIM22 and HLA-C. This indicated that interaction effects between proteins were not only at the level of expression, but also signaling cascade.

GO enrichment analysis of genes also revealed biological processes in the co-expression Pink module (Supplementary Table 8). The genes in the co-expression Pink module were mainly distributed in three parts, including protein trimerization, regulation of complement activation, and response to interferon. Among the biological processes, the responses to interferon also included many detailed functions. This includes: negative regulation of multi-organism process, regulation of symbiosis, encompassing mutualism through parasitism, response to type-I interferon, response to interferon-gamma, negative regulation of viral process, regulation of viral process, defense response to virus, response to interferon-alpha, response to interferon-beta, viral genome replication, cellular response to interferon-gamma, cellular response to type I interferon, regulation of viral life cycle, negative regulation of viral life cycle, negative regulation of viral transcription, regulation of viral transcription, interferon-gamma-mediated signaling pathway, type I interferon signaling pathway, regulation of viral genome replication, negative regulation of viral genome replication, viral entry into host cell, regulation of viral entry into host cell, and negative regulation of viral entry into host cell (Figure 6B). We found consistency was between KEGG and GO analysis, which indicated the function of tumor immunity.

Network analysis and survival-associated lncRNAs

A ceRNA network analysis was used to determine whether lncRNAs regulate the identified mRNAs in the co-expression Pink module through miRNAs. The ceRNA network analysis found three lncRNAs (CYTOR, HCP5, and DANCR) in the co-expression Pink module: 33 miRNAs (hsa-miR-106a-5p, hsa-miR-106b-5p, hsa-miR-128-3p, hsa-miR-135a-5p, hsa-miR-135b-5p, hsa-miR-139-5p, hsa-miR-140-5p, hsa-miR-144-3p, hsa-miR-17-5p, hsa-miR-186-5p, hsa-miR-203a, hsa-miR-20a-5p, hsa-miR-20b-5p, hsa-miR-214-3p, hsa-miR-216a-5p, hsa-miR-216a-5p, hsa-miR-22-3p, hsa-miR-27a-3p, hsa-miR-27b-3p, hsa-miR-299-3p, hsa-miR-29a-3p, hsa-miR-29b-3p, hsa-miR-29c-3p, hsa-miR-328-3p, hsa-miR-33a-5p, hsa-miR-33b-5p, hsa-miR-3619-5p, hsa-miR-496, hsa-miR-519d-3p, hsa-miR-653-5p, hsa-miR-758-3p, hsa-miR-761, hsa-miR-93-5p). 23 identified mRNAs (FLJ36031, SLC22A3, TNFSF10, TRIM22, UBE2L6, OAS2, IL15, DDX60, MLKL, SP100, C1S, CFI, OAS1, HLA-G, CFH, PLSCR1, CTSL1, PARP12, PARP10, RAC2, LGALS3BP, TRIM5, SIX5) in the co-expression Pink module were involved in a ceRNA network (Figure 7A).

LncRNA-RNA binding protein-mRNA network analysis was used to determine whether lncRNAs regulate identified mRNAs through RNA binding proteins. Through this analysis, we found in the co-expression Pink module 13 lncRNAs (CYTOR, CTD-2341M24.1, DANCR, RP11-661A12.4, HCP5, RP11-1398P2.1, RP11-247A12.2, RP11-38L15.3, RP11-661A12.5, RP11-218C14.8, RP11-430H10.3, RP11-977B10.2, RP5-884M6.1). We also found 19 RNA-binding proteins (FUS, CAPRIN1, DGCR8, eIF4AIII, FMRP, FUS-mutant, HuR, IGF2BP1, IGF2BP2, IGF2BP3, LIN28, LIN28A, LIN28B, PTB, TAF15, TIA1, TIAL1, U2AF65, UPF1) and 45 identified mRNAs (C22orf46, CTSL1, DDX60L, HLA-C, IFI6, IFITM1, IL15, LGALS3BP, OAS1, PARP10, PCSK7, PLSCR1, SIX5, SLC22A3, SLC4A11, TDRD7, TRANK1, C1R, C1S, C4orf33, CFH, CFI, DDX60, HLA-G, HYI, IFIT2, IFITM2, FITM3, IL12RB2, LGALS9C, MLKL, OASL, PARP12, RAC2, SERPING1, SP100, TNFSF10, TRIM21, TRIM22, TRIM5, UBE2L6, OAS2, RSAD2, TRPV4, XAF1) in the co-expression Pink module which were involved in the network (Figure 7B).

Furthermore, the Kaplan Meier-plotter analysis revealed that 8 out of the 26 identified lncRNAs in the co-expression Pink module were significantly associated with overall survival (p < 0.05), including CYTOR (p=0.0202), MIR4435-2HG (p=0.0169), RP1-137D17.2 (p=0.0340), RP11-247A12.2 (p=0.0016), RP11-646E18.4 (p=0.0286), RP11-661A12.4 (p=0.0417), RP11-661A12.5 (p=0.0134), RP11-977B10.2 (p=0.0201) (Figure 8).

RT-qPCR validation of identified molecules

qRT-PCR was used to validate the expressions of survival-associated lncRNAs and hub mRNAs resulting from WGCNA analysis, including 8 lncRNAs (CYTOR, MIR4435-2HG, RP1-137D17.2, RP11-247A12.2, RP11-646E18.4, RP11-661A12.4, RP11-661A12.5, and RP11-977B10.2) and 12 hub mRNAs (IFIT2, XAF1, UBE2L6, IFITM3, HLA-C, CTSL, ARHGDIB, LGALS3BP, IFITM1, MLKL, SERPING1, TRIM21) in cultured LC cells (Hep-2 and TU177) and one control cell (HaCaT keratinocytes) (Figure 9) [23]. The results indicated significant difference for five survival-associated lncRNAs (CYTOR, MIR4435-2HG, RP11-247A12.2, RP11-661A12.4, and RP11-661A12.5) (Figure 9A), and eight hub mRNAs (CTSL, XAF1, ARHGDIB, LGALS3BP, IFITM1, MLKL, SERPING1, and TRIM21) (Figure 9B) between LC cells and control cells.

TCGA data of LC patients

Data of LC patients was obtained from the TCGA database (http://cancergenome.nih.gov/). Level 3 RNA-seq V2 (including lncRNA and mRNA expression data) and clinical data of 99 LC patients was downloaded. The basic pretreatment method of RNA-seq data is to remove genes where the missing value (expression = 0) is more than 20%. OS analysis was performed with R survival package, according to genes expression data (cutoff value = median value of gene expression) and survival time in LC patients. 11 LC clinical traits were extracted, including age at initial pathologic diagnosis (Patients were aged 38 to 83), history of smoking (from 0.9 packs to 150 packs), history of alcohol consumption (yes or no), intermediate dimension (from 0.2cm-1.4cm), lymph node count (from 0 to 104), neck lymph node dissection (yes or no), pathologic N stage (N0, N1, N2, N3, and NX), radiation therapy (yes or no), targeted molecular therapy (yes or no), tumor status (with tumor or tumor free), and success of follow-up treatment (complete remission/response, partial remission/ response, persistent disease, progressive disease, and stable disease).

Weighted correlation network analysis of lncRNAs and mRNAs

Based on the processed TCGA data in LC patients, the top 2500 lncRNAs and 2500 mRNA were selected as good genes for further WGCNA analysis. Those 5000 genes were normalized by Limma package with Voom function, then the auxiliary data was removed and expression data was transposed for further analysis. First, all samples were checked for outliers with flashClust to construct sampleTree. The outlier in the sample was then removed based on cutHeight. A sample dendrogram and trait heatmap were visualized by WGCNA package of R software (http://www.r-project.org/) to develop networks to investigate the relationship between the corresponding sample gene expression data and clinical phenotypes [15]. The adjacency matrix aij that calculated the connection strength between each pair of nodes was measured as follows: sij = |cor(xi, xj)| aij = Sijβ. Xi and Xj were vectors of expression values for genes i and j, sij represented the Pearson’s correlation coefficient of genes I and j, aij encoded the network connection strength between genes i and j. β value was the soft-threshold (power value). Moreover, the Scale-Free Topology Fit Index (SFTFI) (scale free R²) ranging from 0 to 1 was used to determine a scale-free topology model. Choosing a set of soft-thresholding powers (range: 1 to 20) can help calculate the scale free topology model fit [16]. In this study, we defined the adjacency matrix using soft thresholding with beta=5, and the corresponding scale free R2 value was 0.87 to obtain a good scale-free topology model. In the co-expression network, genes with highly absolute correlations were clustered into the same module to generate a cluster dendrogram (The parameters are described below: TOMType = unsigned method, minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25). The dendrogram can be displayed together with the color assignment to form the network heatmap plot using the adjacency matrix algorithm method. The average linkage hierarchical clustering was conducted according to Topological Overlap Matrix (TOM)-based dissimilarity measure. Heatmap tool package was plotted to analyze the strength of network interactions. The relationships between co-expression modules and 11 LC clinical traits were calculated by the Pearson correlation coefficient and plotted by heat map. In addition, the Gene Significance (GS) was defined as mediated p-value of each gene (GS = lgP) in the linear regression between gene expression and the clinical traits.

Bioinformatics analysis provided insight into pathways and potential functions

The pathway (http://ci.smu.edu.cn/genclip3/analysis.php) and Gene Ontology (GO) (http://www.cytoscape.org/) enrichment analyses within the co-expression modules were performed to identify LC-related pathway and biological functions, with an adjusted p value < 0.05 (Benjamini-Hochberg for multiple testing). Furthermore, the LC survival-related lncRNAs in the co-expression Pink module were plotted with Rstudio.

Identification of hub molecules with molecular complex detection

The mRNA-mRNA interactions were constructed by STRING and analyzed with Cytoscape software (version 3.2.0; National Resource for Network Biology) to obtain the hub genes. The hub genes were extracted from the identified clinical-related co-expression modules using Molecular Complex Detection. MCODE plugin is a reliable method to distinguish hub molecules from non-hub molecules. The criteria for hub-molecule searching was set as the molecular complex detection (MCODE) score > 6, and statistical significance of p < 0.05 [17].

The ceRNA network and identification of integrated lncRNA-RNA binding protein-mRNA signatures

StarBase provides systematical data of the RNA-RNA and protein-RNA interaction networks from 108 CLIP-Seq (PAR-CLIP, HITS-CLIP, iCLIP, CLASH) data sets generated by 37 independent studies [48]. LncRNA-miRNA-mRNA and lncRNA-RNA binding protein-mRNA interaction networks based on the co-expression modules were decoded from the large-scale CLIP-Seq data by starBase v2.0 (http://starbase.sysu.edu.cn/starbase2/index.php). Network visualization was performed with Cytoscape 3.4.0 (http://www.cytoscape.org/).

Cell lines and cell culture

LC cell lines Hep-2, TU177 cells and normal control cell line HaCaT keratinocytes cells were purchased from Keibai Academy of Science (Nanjing, China). Hep-2 cells were cultured in EMEM medium, TU177 cells were cultured in 1640 medium, and HaCaT cells were in DMEM medium (Corning, NY, USA) supplemented with 10% fetal bovine serum (FBS, Gibco). All these cells were maintained with 5% CO2 atmosphere at 37°C [23].

RNA extraction and qRT-PCR

Total RNAs were extracted from cell lines with TRizol Reagent (Invitrogen) according to the manufacturer’s instructions. Total RNAs were reversely transcribed into cDNAs and then used to perform quantitative real-time PCR (qRT-PCR) with SYBR Premix ExTaq (TaKaRa). Beta-actin was used as an internal control for gene quantification. The RNA molecules that were assessed on the cell lines and their corresponding primers were listed in Supplementary Table 9.

Statistical analysis

All data was analyzed by R software 3.4.3 (https://www.r-project.org/). In all cases, p < 0.05 was considered statistically significance. For KEGG pathway and GO analysis, Pearson correlation coefficient was calculated, and Benjamin-Hochberg was used for multiple testing and calculated to adjust the p-value.