Integrated analysis of lymphocyte infiltration-associated lncRNA for ovarian cancer via TCGA, GTEx and GEO datasets

Data collection and preprocessing

The RNA-seq and miRNA expression profiles of ovarian cancer were downloaded from The Cancer Genome Atlas (TCGA), while the RNA-seq datasets of normal ovarian tissue were searched from the Genotype-Tissue Expression (GTEx) (Carithers et al., 2015). The RNA-seq datasets of counts value from TCGA and GTEx were normalized and processed with TCGAbiolinks package (Mounir et al., 2019) (Table S1). It contains RNA expression profile of 363 high grade serous ovarian cancer and 108 normal tissues. The miRNA datasets from TCGA contain microRNA expression profiles of 479 HGSOC and 8 normal tissues (Table S2). On the basis of the information recorded in the HUGO Gene Nomenclature Committee (HGNC; http://www.genenames.org/) (Eyre et al., 2006), the lncRNAs and protein-coding genes in RNA expression profiles were annotated. Figure S1 illustrates the workflow for bioinformatics based analysis.

Differential RNA expression analysis

The differentially expressed protein-coding genes/lncRNAs (DEGs/DELs) from HGSOC and normal ovarian tissues were selected using the edgeR (Version 3.22.5; http://www.bioconductor.org/pack-ages/release/bioc/html/edgeR.html), limma (Version 3.30.0; http://www.bioconductor.org/packages/release/bioc/html/limma.html) and DESeq2 (Version 1.20.0; http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) package of R software (Nikolayeva & Robinson, 2014; Carvalho et al., 2007; Anders & Huber, 2010). Subsequently, we combined the differentially expressed genes, acquired from the three packages, to get the convergence gene set. On the other hand, differentially expressed miRNAs (DEMs) were identified using the limma package. A false discovery rate (FDR) lt 0.05 and —logFC (fold change)— > 1 were set as the criteria value for DEGs, DEMs and DELs. The expression intensity and direction of DEGs, DELs, and DEMs were represented by heatmaps using the pheatmap R package (Version: 1.0.12; https://cran.r-project.org/web/packages/pheatmap) based on Euclidean distance.

Validation of DEL, DEM, and DEG expressions

We chose four ovarian cancer datasets from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Three of them (GSE18520: 53 vs 10, GSE27651: 22 vs 6, GSE54388: 16 vs 6) were based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) while another one (GSE47841: 12 vs 9) was from GPL14613 platform (Affymetrix Multispecies miRNA-2 Array). In order to improve the number of samples (91 HGSOC samples vs 22 normal samples), we integrated all samples of three datasets (GSE18520, GSE27651, GSE54388), and then removed the batch effect so as to avoid generating less reliable results using sva package (Leek et al., 2012) (Table S3). Principal Component Analysis (PCA) was performed on the gene expression profiles obtained from GEO, to reduce the dimensionality and evaluate the independence of normal and tumor samples. The expressions of DELs, DEMs and DEGs were validated in microarray datasets from GEO that detected expression profiles containing mRNA and lncRNA (GSE18520, GSE27651, GSE54388) as well as miRNA (GSE47841) between the HGSOC and normal tissues. The difference of expression was tested by t-test. A value of p <  0.05 was considered statistically significant.

Function enrichment analysis

Biological process of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs was conducted by using the clusterProfiler R package (Version 6.8; http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (Yu et al., 2012) to predict their underlying functions. A p < 0.05 was considered to be statistically significant.

Survival analysis

The expressions of DEG/DEL/DEM were extracted and merged with the survival information in each sample (Table S4). Univariate Cox regression analysis was implemented to perform prognostic analysis using the R survival package (Version 2.43.3; https://cran.r-project.org/package=survival). The survival result was expressed as Kaplan–Meier (K-M) curve, with higher expression of a specific gene versus lower expression of this gene in a patient, the cutoff value determined by survminer (Version 0.4.3; https://cran.r-project.org/package=survminer). Again a p < 0.05 was considered as related to prognostic significance.

lncRNA–miRNA–mRNA ceRNA regulatory network construction

The starBase database (http://starbase.sysu.edu.cn/) (Li et al., 2014) was used to study the lncRNA-miRNA interactions. The target genes of miRNAs in the DELs–DEMs interaction network were predicted using the miRTarBase database. The negative interaction pairs between DEMs and DEGs/DELs were integrated to construct the DELs–DEMs–DEGs ceRNA network using the Cytoscape software (Version 3.6.1; http://www.cytoscape.org/) (Kohl, Wiese & Warscheid, 2011). Furthermore, we analyzed functional enrichment of the genes in the ceRNA network with clusterProfiler R package.

Immune infiltrates analysis

We analyzed the correlation between lncRNAs in ceRNA network and the abundance of immune cell infiltrates; including B cells, CD⁴⁺ T cells, CD⁸⁺ T cells, neutrophils, macrophages, and dendritic cells using Tumor IMmune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/) (Li et al., 2017). And the tumor purity was calculated using the estimate R package (https://r-forge.r-project.org/projects/estimate/), In addition, correlations between the lncRNAs and marker genes of tumor-infiltrating lymphocytes were explored using Spearman’s correlation. The gene markers of tumor-infiltrating lymphocytes included markers of CD⁸⁺ T cells, B cells, general T cells, tumor-associated macrophages (TAMs), M1 macrophages, M2 macrophages, monocytes, neutrophils, natural killer (NK) cells, dendritic cells (DC), T helper 1 (Th1) cells, T helper 2 (Th2) cells, T helper 17 (Th17) cells, follicular helper T (Tfh) cells, regulatory T cell (Tregs) and exhausted T cells, which were referenced in prior studies (Sousa & Maatta, 2016). The gene expression level is displayed with log2 transformation.

Differential expression analysis

A total of 629 DELs (343 downregulated and 286 upregulated), 7, 514 DEGs (3, 319 downregulated and 4, 195 upregulated) and 336 DEMs (165 downregulated and 171 upregulated) were identified between HGSOC and normal ovarian tissues based on the threshold of FDR<0.05 and —logFC—>1. The expression changes of DELs (Fig. S2A), DEMs (Fig. S2B), and DEGs (Fig. S2C) were displayed by heatmap analysis.

Validation of DEL, DEM, and DEG expressions

To validate the expressions of detected differently expressed genes above, we utilized datasets from GEO. The results of PCA indicated that normal samples vs. tumor samples in the merged datasets had a significant difference (Fig. S3). After combining gene expressions in GEO, 111 DELs (77 downregulated and 34 upregulated), 2,620 DEGs (1,535 downregulated and 1,085 upregulated) while 186 DEMs (101 downregulated and 85 upregulated) were identified for subsequent analysis.

Functional annotation for DEGs

GO and KEGG pathway analysis were performed to explore the potential biological functions of DEGs. The results, through KEGG pathway enrichment analysis, indicated that upregulated genes were mainly involved in cell cycle, oxidative phosphorylation, homologous recombination, cytokine-cytokine receptor interaction and IL-17 signaling pathway, while downregulated genes mainly enriched in cGMP-PKG signaling pathway, complement and coagulation cascades, TGF-β signaling pathway and Ras signaling pathway, as shown in Fig. 1A and Table S5. Furthermore, the GO biological process terms were performed. In line with the KEGG enrichment results, regulation of cell cycle phase transition and positive regulation of cell cycle were obtained for the upregulated genes, while transmembrane receptor protein serine/threonine kinase signaling pathway and regulation of membrane lipid distribution was enriched in the downregulated genes (Fig. 1B and Table S6).

Prognosis prediction for DELs, DEMs, and DEGs

We further performed the survival analysis for DELs, DEMs, and DEGs base on TCGA data and clinical information to excavate the prognostic genes. The results showed that 33 DELs (25 upregulated, 8 downregulated), 134 DEMs (76 upregulated, 58 downregulated), and 1612 DEGs (949 upregulated, 663 downregulated) were significantly associated with overall survival (OS) (p < 0.05) (Table S7).

ceRNA network construction

We then constructed ceRNA networks to explore the interactions between DELs, DEMs and DEGs. The 33 OS-related DELs were mapped into the starBase database to predict their target microRNAs. Comparing with the OS-related DEMs, lncRNA/miRNA interacted networks were constructed. Subsequently, the target genes of selected DEMs were searched using miRTarBase database for identified DEMs–DEGs interaction. After taking the change of expression and prognosis information into consideration, several DELs-DEMs-DEGs interactions including 4 DELs (3 upregulated, 1 downregulated), 8 DEMs (3 upregulated, 5 downregulated) and 44 DEGs (15 upregulated, 29 downregulated) were predicted and resultantly a DELs–DEMs–DEGs ceRNA network was established, as shown in Fig. 2.

Function enrichment analysis showed the genes in the ceRNA network participate in several KEGG pathways; including T cell receptor signaling pathway, MicroRNAs in cancer progression pathway, FoxO signaling pathway, Autophagy-animal and mTOR signaling pathways (Table S8). The DEL–DEM–DEG interactions associated with signaling pathways are shown in Fig. 3, which involved 3 DELs (2 upregulated and 1 downregulated), 4 DEMs (2 upregulated and 2 downregulated), and 7 DEGs (4 upregulated and 3 downregulated).

Further K-M curves of genes in the ceRNA network were plotted. The relationships between their respective expressions and the prognostic values were in accordance with our hypothesis. The results showed that patients with a high expression of DELs (LINC00858 and FTX) and DEGs (PAK4, ZEB1, STK11 and SETD7) indicated a poor survival prognosis, meanwhile the high expressions of DEMs (has-miR-324-5p, has-miR-150-5p and has-miR-342-3p) in the ceRNA axes presented and good survival prognosis. Similarly, the high expressions of LINC00665 and the related protein coding genes (VAV3, E2F3, RRAGD) played their roles as suppressor genes, while the high expressions of has-miR-449-5p acted as an oncogene (Fig. 4 and Fig. S4).

lncRNA expression is correlated with lymphocyte infiltration level in HGSOC

Since the functional analysis indicated that the mRNAs in ceRNA network enriched in immune-related pathways, such as T-cell receptor signaling pathway, Autophagy –animal, and mTOR signaling pathway, whether the lnRNAs played potential roles in lymphocyte infiltration? Tumor-Infiltrating Lymphocyte Grade has been proven to be an independent predictor of sentinel lymph node status, as well as survival in cancers (Ohtani, 2007; Azimi et al., 2012). However, genomic approaches that are commonly used for the analysis of lymphocyte infiltration in clinical tumor samples can be severely influenced by tumor purity (Yoshihara et al., 2013). We put the expression profiles into TIMER website after eliminating the influence of tumor purity.

The results showed that expression levels of three lncRNAs all have significantly positive correlations with infiltrating levels of tumor purity (Fig. S5, first column). FTX exhibited significant negative correlation with infiltrating levels of CD⁸⁺ T cells (r =  − 0.13, P = 0.012), neutrophils (r =  − 0.12, P = 0.018), and macrophages (r =  − 0.17, P = 0.00091) as shown in Fig. 5A. LINC00855 also has negative correlations with neutrophils (r =  − 0.12, P = 0.021) and DCs (r =  − 0.11, P = 0.043) (Fig. 5B). On the other hand, LINC00665 demonstrated significant positive correlation with infiltrating levels of CD⁴⁺ T cells (r = 0.11, P = 0.044), negative correlations with CD⁸⁺ T cells (r =  − 0.16, P = 0.0026), neutrophils (r =  − 0.27, P = 1.8^e−07), macrophages (r =  − 0.11, P = 0.029) and DCs (r =  − 0.26, P = 3.7e−07) (Fig. 5C). These findings suggest that three lncRNAs, LINC00855, FTX and LINC00665, play a specific role in lymphocyte infiltration in HGSOC, especially LINC00665.

Correlation analysis between lncRNA expression and immune marker sets

We further investigated the correlations between lncRNAs and marker sets of different immune cells, including CD⁸⁺T cells, T cells (general), B cells, monocytes, TAMs, M1 and M2 macrophages, neutrophils, NK cells, and DCs. Different type of T cells were also analyzed, which include Th1 cells, Th2 cells, Tfh cells, Th17 cells, Tregs, as well as exhausted T cells (Table 1 and Table S9).

After adjusting correlation by purity, the results revealed that most immune markers of various immune cells and T cells are significantly correlated with the FTX and LINC00665 expression level. However, only 22 gene markers were found to be correlated with the expression level of LINC00885.

It is also found that the expression levels of most markers of monocytes, TAMs, M2 macrophages have strong correlations with LINC00665 expression.

As is shown in Fig. 5, high LINC00665 expression correlates with low infiltration level of DCs. This result can be explained by the following findings that DC markers such as HLA-DPB1, HLA-DQB1, HLA-DRA, HLA-DPA1, CD1C (BDCA-1) and ITGAX (CD11c) show significant negative correlations with LINC00665 expression. We can further conclude that there is a significant negative correlation between LINC00665 and DC infiltration. This finding is rather interesting because tumor metastasis can be aggravated by DC reducing CD⁸⁺ T cell cytotoxicity and increasing Treg cells (Sawant et al., 2012).

It is also found that there are significant correlations between FTX, LINC00665 and marker genes of Treg and T cell exhaustion, such as FOXP3, TGFB1 (TGFβ), IL10, TBX21 (T-bet), PDCD1 (PD-1), CTLA4 and HAVCR2 (TIM-3). FOXP3 is a critical marker of Treg cells which suppress cytotoxic T cells attacking tumor cells (Facciabene, Motz & Coukos, 2012). We found that FTX and LINC00665 expressions have strong negative correlations with TIM-3, which is known as a critical factor that regulates T cell exhaustion. This finding suggests that these two lncRNA expressions are of great importance in TIM-3 mediating T cell exhaustion (Table 1). Accordingly, it is further confirmed that FTX and LINC00665 are specifically correlated with immune infiltrating cells, suggesting that these two lncRNAs play vital roles in inhibiting immune escape in the ovarian microenvironment.