Prognostic role of long non-coding RNA USP30-AS1 in ovarian cancer: insights into immune cell infiltration in the tumor microenvironment

Data acquisition and preprocessing

We retrieved processed whole transcriptome data for 34 cancer cohorts, including both normal and tumor tissues, from The Cancer Genome Atlas (TCGA) database through the UCSC XENA website (https://xena.ucsc.edu/). Clinical data and gene mutation data for Ovarian Cancer (TCGA-OV) were also obtained from the UCSC XENA platform. In addition, transcriptome data for ovarian cancer were sourced from the Genotype-Tissue Expression (GTEx) database (https://www.gtexportal.org/home/index.html). These databases are open-access, allowing unrestricted access and analysis, thus obviating the need for ethical approval. We ensured that data retrieval and analysis procedures adhered to the ethical requirements of the respective databases and journals.

Transcriptomic data for tumor and normal tissue specimens were acquired from gene expression data records, following TCGA's unique specimen naming conventions. Gene expression data was standardized using Log₂(FPKM+1). Information on immune therapy cohorts was obtained from the “IMvigor210” dataset. The list of genes related to TLS (Tertiary Lymphoid Structures) was derived from literature searches, and the detailed gene list is provided in Supplementary Table 1.

Differential TLS gene selection and further analysis

Initially, we integrated ovarian cancer whole transcriptome data from GTEx and TCGA (Normal=88, Tumor=379) and employed the “limma” package to analyze the differential expression of TLS genes between tumor and normal groups. The selection criteria were set at false discovery rate (FDR) < 0.05 and | logFC | > 1. A volcano plot was generated to visually represent the differentially expressed genes. Subsequently, we utilized the protein-protein interaction visualization feature of the STRING website (https://www.string-db.org/) to illustrate the interaction relationships among differential TLS genes. Univariate Cox analysis was employed to assess the impact of differential TLS genes on prognosis, considering genes with p < 0.01 as those significantly affecting prognosis. Subsequently, a Venn diagram was employed to identify the intersection between differential TLS genes and prognostically impactful TLS genes. The resulting genes were then visualized in a heatmap to observe their expression patterns in tumor and normal tissues. Based on the median expression values of the identified genes, tumor patients were stratified into high and low expression groups. Kaplan-Meier (KM) curves were applied to analyze the overall survival differences between the two expression groups for each gene. Finally, leveraging IMvigor immune therapy cohort data, we analyzed the expression differences of each gene between effective and ineffective immune therapy groups, presenting the results through violin plots.

Non-negative matrix factorization (NMF) analysis and subsequent investigations

Building upon the TLS genes identified through the intersection of the Venn diagram in Section 2.2, we conducted NMF clustering analysis on ovarian cancer patients from TCGA. The final number of clusters was determined by selecting the value of 'k' that exhibited a “high cohesion, low coupling” pattern in the consensus matrix plot. Subsequently, Kaplan-Meier survival curves were generated for the two clusters. Following this, we conducted a comprehensive analysis of the tumor immune microenvironment (TIME) for these two clusters, encompassing the evaluation of Immune Score, Stromal Score, Estimate Score, and tumor purity. We further assessed the expression profiles of eight immune molecules (IgG, HCK, MHC-I I, LCK, STAT1, Interferon, B7-CD28, TNF) in these two clusters, visualizing differences in survival status, cancer stage, TIME, and the expression of these eight immune molecules using heatmaps. Single-sample Gene Set Enrichment Analysis (ssGSEA) was employed to investigate differences in the infiltration of 29 immune cell types between the two clusters. Additionally, a differential gene expression analysis was conducted on these two clusters, with the selection criteria set at | logFC | > 1 and P-adj < 0.05. The resultant differentially expressed genes were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathway enrichment analyses [25–27].

Construction and validation of a prognostic model based on TLS gene-associated lncRNA

Utilizing the transcriptome data of TCGA-OV tumor patients, we conducted an analysis of lncRNAs co-expressed with TLS genes, as identified through the intersection of the Venn diagram in Section 2.2. LncRNAs were selected based on the criterion of a correlation coefficient greater than 0.35 and a p-value less than 0.05 for each gene. Subsequently, Spearman correlation analysis was performed for each TLS gene with its four strongest correlated lncRNAs, and the results were visualized using scatterplots. Next, the identified lncRNAs were incorporated into The Least Absolute Shrinkage and Selection Operator (LASSO) machine learning algorithm. Ten-fold cross-validation was employed for further refinement of the selected lncRNAs, ultimately leading to the construction of a risk model. Our random seed for this analysis was set to “1997.” The risk score was defined by multiplying the expression value of each gene (i) by its corresponding coefficient (βi), as follows:

$$\begin{array}{c}\text{Risk}\text{score}=\\ {\sum }_{\text{i}=1}^{\text{n}}[expressionvalueofgen{e}_i\ast {\beta }_i]\end{array}$$ (1)

Where β_i represents the coefficient corresponding to the model gene [28–30]. Subsequently, we performed univariate Cox analysis on the model lncRNAs to assess their impact on prognosis and conducted Kaplan-Meier survival curve analysis to evaluate the influence of each lncRNA's expression on prognosis. Receiver operating characteristic curve (ROC) analysis was employed to assess the model's performance in predicting 2-year, 3-year, and 5-year outcomes. An area under the curve (AUC) greater than 0.6 indicated good predictive performance. Patients were then categorized into high-risk and low-risk groups based on the median risk score, and Kaplan-Meier analysis was performed separately for overall survival (OS) and progression-free survival (PFS) in these two groups. Furthermore, a risk factor plot was generated to illustrate how patient survival time, survival status, and the expression of model lncRNAs change with increasing risk scores. Additionally, we utilized Cytoscape software to visualize the expression correlations between TLS genes and model lncRNAs, which were presented using chord diagrams to enhance interpretability and understanding.

Establishment and validation of the prognostic nomogram

We integrated age, stage, grade, and risk score to comprehensively construct a nomogram model for predicting patient survival at 2, 3, and 5 years. To assess the accuracy of the model, we plotted calibration curves at 2, 3, and 5 years to observe the extent of deviation from the ideal values. To validate the clinical utility of the model, we performed a 5-year Decision Curve Analysis (DCA) for age, stage, grade, risk score, and the nomogram. Finally, we conducted three-dimensional Principal Component Analysis (PCA) separately for the entire transcriptome genes, TLS genes, TLS-related lncRNA genes, and model lncRNA genes. This allowed us to examine whether there were significant differences in the three-dimensional data space dimensions between the high-risk and low-risk sample groups.

Multiple analytical approaches reveal the intrinsic mechanisms of TLS-related lncRNAs

In our study, we conducted an immune subtype analysis on the TCGA-OV dataset, eliminating subtypes with limited sample sizes. We then proceeded to analyze the differences in lncRNA model risk scores among the remaining subtypes. Utilizing Sankey diagrams, we illustrated the associations between NMF subtypes, risk score risk groups, and immune subtypes. Subsequently, we performed differential analysis of risk scores between groups categorized by treatment outcomes and tumor recurrence status. We employed the “oncopredict” tool to investigate the sensitivity differences of 200 drugs between the two risk groups. Violin plots were created to visualize the sensitivity differences of three specific drugs: LCL161, Ribociclib, and Topotecan. We conducted differential analysis between the two risk groups and subjected the obtained genes to Gene Ontology (GO) analysis. The analyses were carried out using packages such as “limma”, “scales”, “ggplot2”, “ggtext”, “reshape2”, “tidyverse”, and “ggpubr”.

Furthermore, we conducted differential analysis of 16 immune cell infiltrations between the two risk groups and visualized the results using box plots. Subsequently, we employed Weighted Correlation Network Analysis (WGCNA) to identify lncRNAs co-expressed with the ssGESA scores of various immune cells. The relevant parameters included a cutline of 600 and a power value of 4. We selected the immune cell with the strongest statistical significance from the previous two risk group analyses as the target phenotype for further analysis. We identified the WGCNA module with the strongest correlation to this immune cell phenotype and intersected the lncRNAs within this module with the model lncRNAs. We conducted an analysis of immune checkpoints based on the obtained lncRNAs. Additionally, we performed differential expression analysis of the obtained lncRNAs between different risk groups, followed by Kaplan-Meier survival curve analysis for progression-free survival (PFS) and TIDE analysis for Exclusion scores in different risk groups. Based on these three analyses, we filtered out lncRNAs without statistical significance and compared the remaining lncRNAs for their expression differences between different grades and between tumor and normal tissues. We selected USP30-AS1 for further analysis.

To stratify OC patients into two groups based on the median expression value of USP30-AS1, we analyzed the most significantly mutated genes between the two groups and visualized them using waterfall plots. Subsequently, we performed multiple Gene Set Enrichment Analyses (GSEA) to observe enriched pathways in these two groups. We selected the most significantly mutated genes from the waterfall plot and retrieved their immunohistochemistry images in normal ovarian tissues and cancer tissues from The Human Protein Atlas (HPA) website (https://www.proteinatlas.org/) [31]. Finally, we generated scatter plots to illustrate the correlation between USP30-AS1 and six immune cell types in the TCGA-OV dataset.

Pan-cancer analysis of USP30-AS1

To understand the expression pattern of USP30-AS1 in other cancers, we downloaded processed transcriptomic data of 34 cancer cohorts, including both normal and tumor tissues, from the UCSC XENA website (https://xena.ucsc.edu/). Subsequently, we conducted differential analysis of USP30-AS1 expression between tumor and normal tissues in each cancer cohort and visualized the results using violin plots. We employed the “xCELL” algorithm from the R package “IOBR” (version 0.99.9) to calculate scores for 67 immune cell infiltration in each cancer cohort [32]. Next, we computed the correlation between immune cell infiltration scores and USP30-AS1 expression levels in each cancer cohort using the “Pearson” correlation method. A p-value less than 0.05 was considered statistically significant for the correlation between immune cell infiltration scores and USP30-AS1 expression levels. We visualized the correlation between each immune cell score and USP30-AS1 expression using heatmaps for each cancer cohort. Finally, we selected the cancer cohort with the highest absolute correlation coefficient between Th1 cell score and USP30-AS1 expression. We illustrated the correlation between Th1 cell score and USP30-AS1 expression using scatter plots to visualize the relationship.

Cell culture, total RNA extraction and RT-qPCR

We purchased three human ovarian cancer cell lines, SK-OV-3, OVCAR-3, CAOV-3, and one human normal ovarian epithelial cell line, IOSE-29, from the Cell Bank of the Chinese Academy of Sciences. SK-OV-3 cells were cultured in McCoy's 5A medium (HyClone, USA), OVCAR-3 cells were cultured in Roswell Park Memorial Institute (RPMI-1640) medium (HyClone, USA), CAOV-3 cells were cultured in Dulbecco's Modified Eagle Medium (DMEM, HyClone, USA), and IOSE-29 cells were cultured in a 1:1 mixture of MCDB 105 and M-199 media (HyClone, USA). All media were supplemented with 10% fetal bovine serum (FBS, BI, USA), 100 U/ml penicillin (HyClone, USA), and 100 μg/ml streptomycin (HyClone, USA). All cell lines were maintained in a cell culture incubator at 37° C with 5% CO₂.

For total RNA extraction, we employed trypsin (BI, USA) to digest the cells grown in cell culture flasks, and collected them in EP tubes. Subsequently, 900 μl of Trizol (Takara, Japan) was added to lyse the cells, followed by a 5-minute incubation on ice. Then, 200 μl of chloroform (SINOPHARM, China) was added, and vigorous shaking was performed, followed by centrifugation at low temperature for 15 minutes. The clear aqueous phase was carefully transferred to a fresh tube, and an equal volume of isopropanol and absolute ethanol (SINOPHARM, China) was added successively. After each addition, the mixture was incubated on ice for 5 minutes, followed by centrifugation at low temperature. The organic solvents were discarded, and the RNA pellet obtained was air-dried, resembling a feather-like appearance. The RNA pellet was dissolved in an appropriate amount of DEPC-treated water, and the concentration was measured. Subsequently, following the manufacturer's instructions, cDNA synthesis was performed using the PrimeScript RT Kit (TaKaRa, Japan). For subsequent real-time fluorescence quantitative PCR, we utilized the SYBR GreenER Supermix (TaKaRa, Japan) kit. The PCR reactions were conducted on a 7500 Real-Time PCR System (Thermo Fisher Scientific, USA). The PCR reaction conditions comprised an initial denaturation step at 95° C for 1 minute, followed by 35 cycles of 95° C for 90 seconds, 60° C for 30 seconds, and 72° C for 30 seconds, and a final extension step at 72° C for 10 minutes. The relative expression of USP30-AS1 was analyzed using the 2^-ΔΔCt method, with β-actin as the reference gene for normalization. The primer sequences used were as follows: β-actin: Forward: 5’- CCTGGCACCCAGCACAAT-3’, Reverse: 5’- GGGCCGGACTCGTCATAC-3’; USP30-AS1: Forward: 5rward: 5’- CCTGGCACCCA -3rward: 5rward: TGAAAACCAAGCAGCCCCA -3AA.

Statistical analysis

For statistical analysis, R (version 4.1.0) was utilized. The two groups were contrasted using the Student’s t-test or Wilcoxon test. Comparisons between several groups were evaluated using the Kruskal-Wallis test. Kaplan-Meier plots, which show survival curves, were used to compare them to the log-rank test. Results with a p-value < 0.05 were considered significant from a statistical perspective.

Differential selection and further analysis of TLS genes

The workflow of this study is illustrated in Figure 1. A total of 14 TLS genes exhibited differential expression between the tumor and normal groups. With the exception of CCL21, all other genes were upregulated in the tumor group (Figure 2A). Protein interaction networks revealed that several genes, including CCL21, CCL2, CXCL10, CCL5, CXCL13, and CCL8, exhibited close interactions with other genes (Figure 2B). Univariate Cox analysis demonstrated that CXCL9, CXCL10, CXCL11, and CXCL13 were protective factors in ovarian cancer (HR < 1, p < 0.01, Figure 2C). Subsequently, we intersected the differentially expressed TLS genes with these four genes, indicating that these four genes were both differentially expressed and prognostically protective (Figure 2D). Heatmaps displayed high expression of these four genes in the tumor group (Figure 2E). Kaplan-Meier curves indicated that patients with high expression had a more favorable prognosis (p < 0.001, Figure 2F). In the Imvigor cohort, the effective treatment group showed elevated expression of CXCL9, CXCL10, CXCL11, and CXCL13 compared to the treatment-resistant group.

Non-negative matrix factorization (NMF) analysis and further analysis

In accordance with the principle of “high cohesion, low coupling,” we categorized TCGA-OV patients into 2 clusters (Figure 3A). Kaplan-Meier (KM) curves indicated that patients in cluster 1 had a better prognosis (p = 0.037, Figure 3B). TIME (Tumor Immune Microenvironment) analysis revealed that cluster 1 had higher Immune Score, Stromal Score, and Estimate Score, whereas cluster 2 exhibited higher tumor purity (Figure 3C). The heatmap illustrated higher survival rates among patients in cluster 1. IgG, HCK, MHC-I I, LCK, B7-CD28, and TNF were expressed at higher levels in cluster 1, whereas STAT1 and Interferon were expressed at higher levels in cluster 2 (Figure 3D). Immune analysis demonstrated significant differences in the infiltration of various immune cells between the 2 clusters (p < 0.05). The immune cell infiltration scores in cluster 1 were notably higher than those in cluster 2, consistent with the TIME results (Figure 3E). Furthermore, we conducted differential gene expression analysis for the 2 clusters (Supplementary Table 2). KEGG pathway analysis of the differentially expressed genes indicated upregulation in pathways related to protein polysaccharides, chemokine signaling, ECM-receptor interaction, cytokine-cytokine receptor interaction, and PI3K-Akt signaling in cancer (Figure 3F). GO analysis revealed upregulation of differentially expressed genes in various extracellular matrix components and cellular constituents (Figure 3G).

Construction and validation of a prognostic model based on TLS gene-associated lncRNA

Based on the transcriptome data of TCGA-OV tumor patients, we conducted an analysis of long non-coding RNAs (lncRNAs) co-expressed with TLS genes, as obtained from the Venn diagram intersection in Figure 4A. Spearman analysis revealed a strong positive correlation (p < 0.001, R > 0.5, Supplementary Figure 1) between CXCL9, CXCL10, CXCL11, CXCL13, and their respective top 4 correlated lncRNAs. By employing LASSO analysis, we derived a 15-gene risk model (MICB-DT, AL078582.1, LINC01943, AC012181.1, AL365361.1, AC012236.1, USP30-AS1, AL353699.1, AP002954.1, LINC01094, TRBV11-2, AC002511.2, DTNB-AS1, LINC01857, PSMB8-AS1) (Figure 4B, 4C, and Supplementary Table 3). Univariate Cox analysis demonstrated that 13 lncRNAs significantly impacted patient prognosis, with the majority serving as protective factors. This observation was further corroborated by KM survival curves (Supplementary Figure 2). ROC curves indicated excellent diagnostic performance of the model at 2, 3, and 5 years (AUC > 0.6, Figure 4D). Patients in the high-risk group exhibited worse overall survival (OS) and progression-free survival (PFS) (p < 0.01, Figure 4E). Risk factor analysis illustrated an increasing number of deceased patients and a gradual reduction in the expression of model lncRNAs with rising risk scores (Figure 4F). There was a wide-ranging association between CXCL9, CXCL10, CXCL11, CXCL13, and the model lncRNAs (Figure 4G). Additionally, strong positive correlations were observed among the model lncRNAs themselves (Figure 4H).

Establishment and validation of the nomogram prognostic model

In our study, we incorporated age, stage, grade, and risk score to comprehensively construct a nomogram model for predicting patients' 2-year, 3-year, and 5-year survival outcomes (Figure 5A, and Supplementary Table 4). Calibration curves demonstrated the excellent predictive accuracy of the nomogram at all three time points (Figure 5B). Decision curve analysis (DCA) curves further indicated that the nomogram model yielded superior decision-making benefits (Figure 5C). Moreover, the results of a three-dimensional principal component analysis (3D PCA) revealed that under the context of lncRNA in the model, there was better differentiation between high-risk and low-risk groups of samples (Figure 5D).

Multiple analytical modes reveal the underlying mechanisms of TLS-related lncRNAs

The risk scores among TCGA immune subtypes exhibit significant differences (p < 0.0001, Figure 6A). The Sankey plot illustrates that a majority of patients in cluster 2 occupy both risk groups, while cluster 1 patients are in the minority. Patients in the low-risk group are primarily distributed within immune subtype 2, whereas those in the high-risk group are scattered across various immune subtypes (Figure 6B). The ineffective treatment group shows higher risk scores (p = 7.5e-05), and the recurrence group exhibits even higher risk scores (p = 0.0092, Figure 6C). Sensitivity analysis of 200 drugs (Supplementary Table 5) reveals that patients in the high-risk group are more sensitive to three drugs, LCL161, Ribociclib, and Topotecan, compared to the low-risk group (p < 0.0001, Figure 6D). These drugs are currently under investigation, and our results suggest their potential in ovarian cancer treatment. Gene Ontology (GO) analysis of differentially expressed genes in the two risk groups indicates that neutrophil and lymphocyte functions are prominent in biological processes (BP). Plasma cell function ranks highest in cellular components (CC), while chemotaxis tops molecular functions (MF). This underscores the significance of tumor immunity as a crucial factor distinguishing high and low-risk groups (Figure 6E). Analysis of eight immune cell infiltrations reveals significant differences in immune cell scores between high and low-risk groups, implying a strong association between the lncRNA risk formula based on TLS and various immune cells (Figure 6F).

In the low-risk group, most immune cell infiltrations are higher than those in the high-risk group, with the most significant difference observed in Th1 cells (p = 3.5e-12, Figure 7A). We conducted a correlation analysis between the six immune cells with the smallest p-values and various modules of WGCNA. Darkgreen module exhibits the highest correlation with Th1 (cor = 0.79, p = 9e-89, Figure 7B, 7C). We selected 16 lncRNAs with correlation coefficients greater than 0.4 from this module and identified six genes that intersect with model genes, demonstrating strong prognostic prediction capabilities related to Th1 (Figure 7D). The co-expression relationship between these six intersected genes and immune checkpoints is significant (Figure 7E).

Subsequent analysis focused on the six intersected genes. We employed TIDE scoring, PFS survival curves, and TIDE Exclusion scoring to screen out USP30-AS1 and Al365261.1 (Figure 8A–8C). USP30-AS1 exhibited differences between Grade 2 subgroups, while Al365261.1 did not. Consequently, we decided to proceed with the analysis of USP30-AS1 (Figure 8D). We categorized the expression of USP30-AS1 into high and low expression groups based on its median expression value. We analyzed the mutation differential genes between the two groups and presented the top 15 significant genes in a waterfall plot. LRP1B exhibited the highest mutation level (Figure 9A). Multiple GSEA analyses of the two expression groups indicated the involvement of RYR1 and CACNA1E in the calcium channel pathway: CALCIUM SIGNALING PATHWAY. This suggests that USP30-AS1 may interfere with the tumor immune response of immune cells such as Th1 by affecting the physiological opening of calcium channels, thereby promoting/inhibiting tumor immune escape (Figure 9B). Immunohistochemical analysis demonstrated higher expression of LRP1B in tumor tissues (Figure 9C). The infiltration of multiple immune cells is positively correlated with the expression of USP30-AS1 (R > 0.2, p < 0.001, Figure 9D).

Pan-cancer analysis of USP30-AS1

The results of the differential expression analysis indicate that USP30-AS1 exhibits differential expression across multiple cancer cohorts, except for uterine carcinosarcoma (UCS), pheochromocytoma and paraganglioma (PCPG), adrenocortical carcinoma (ACC), and kidney chromophobe (KICH). In most cancer types, USP30-AS1 is upregulated (P < 0.05, Figure 9A). From a macroscopic perspective, immune cell infiltration scores are inversely correlated with USP30-AS1 expression in most cancers. However, in the case of immune cells such as aDC, CD8⁺ T cells, Macrophages M1, and USP30-AS1 expression, a positive correlation is observed in most cancer cohorts (Figure 9B). Given that the differential analysis of immune cells between high and low-risk groups previously revealed the most significant statistical difference in Th1 cells between the two risk groups, we have selected the four cancer cohorts with the highest absolute correlation coefficients between Th1 cell infiltration scores and USP30-AS1 expression for scatter plot representation of their correlation. In uterine corpus endometrial carcinoma (UCEC), head and neck squamous cell carcinoma (HNSC), and Cervical squamous cell carcinoma (CESC), Th1 cell infiltration positively correlates with USP30-AS1 expression (P < 0.0001, correlation coefficient >0.3). Conversely, in glioblastomas (GBM), Th1 cell infiltration exhibits a negative correlation with USP30-AS1 expression (P < 0.0001, correlation coefficient = -0.36, Figure 9C).

USP30-AS1 is highly expressed in ovarian cancer cell lines

We conducted RT-qPCR analysis on three ovarian cancer cell lines and one normal ovarian epithelial cell line. The results revealed a significantly higher expression of USP30-AS1 in the ovarian cancer cell lines compared to the normal ovarian epithelial cell line (p < 0.01). Among them, the OVCAR-3 cell line exhibited the highest expression level (Figure 10). Hence, we conclude that USP30-AS1 is highly expressed in ovarian cancer cell lines.