Development and Verification of an Autophagy-Related lncRNA Signature to Predict Clinical Outcomes and Therapeutic Responses in Ovarian Cancer

Objective: Long noncoding RNAs (lncRNAs) are key regulators during ovarian cancer initiation and progression and are involved in mediating autophagy. In this study, we aimed to develop a prognostic autophagy-related lncRNA signature for ovarian cancer. Methods: Autophagy-related abnormally expressed lncRNAs were screened in ovarian cancer with the criteria values of |correlation coefficient| > 0.4 and p < 0.001. Based on them, a prognostic lncRNA signature was established. The Kaplan–Meier overall survival analysis was conducted in high- and low-risk samples in the training, verification, and entire sets, followed by receiver operating characteristics (ROCs) of 7-year survival. Multivariate Cox regression analysis was used for assessing the predictive independency of this signature after adjusting other clinical features. The associations between the risk scores and immune cell infiltration, PD-L1 expression, and sensitivity of chemotherapy drugs were assessed in ovarian cancer. Results: A total of 66 autophagy-related abnormally expressed lncRNAs were identified in ovarian cancer. An autophagy-related lncRNA signature was constructed for ovarian cancer. High-risk scores were indicative of poorer prognosis compared with the low-risk scores in the training, verification, and entire sets. ROCs of 7-year survival confirmed the well-predictive efficacy of this model. Following multivariate Cox regression analysis, this model was an independent prognostic factor. There were distinct differences in infiltrations of immune cells, PD-L1 expression, and sensitivity of chemotherapy drugs between high- and low-risk samples. Conclusions: This study constructed an autophagy-related lncRNA signature that was capable of predicting clinical outcomes and also therapeutic responses for ovarian cancer.


INTRODUCTION
Ovarian cancer represents the major cause of death among gynecological malignancies (1). Surgery followed by chemotherapy (such as platinum and taxane) remained the first-line therapeutic strategy (2). Approximately 80% of the subjects originally respond to this therapy. Nevertheless, the majority of the subjects in late stages usually experienced recurrence following chemotherapy, thereby leading to an undesirable prognosis (5-year survival < 50%). Hence, it is significant to probe into the pathogenesis of ovarian cancer and also predictive indicators for prognostic stratification.
Alterations in gene expression profiling have become fundamental laboratory tools for improving tumor diagnoses, survival outcomes, and also treatment responses, which overcome weaknesses of typical clinical and imaging features due to heterogeneity at the genetic and molecular levels (3). Dysregulated long non-coding RNAs (lncRNAs) such as LINC00189, CACNA1G-AS1, and CHRM3-AS2 have been implicated in ovarian cancer initiation and the progress, highlighting their promising functions as markers of precision medicine (4). Their correlations to survival outcomes and treatment responses are still indistinct in ovarian cancer. A few lncRNA-based expression signatures have been developed for predicting survival outcomes, status, and chemosensitivity in ovarian cancer. For instance, Zheng et al. developed a three-lncRNA signature (LOC101927151, LINC00861, and LEMD1-AS1) for predicting clinical outcomes of the subjects with ovarian cancer on the basis of copy number variation (5). Zhang et al. proposed a three-lncRNA signature (LINC01619, DLX6-AS1, and AC004943.2) that can predict survival and response of chemotherapy in ovarian cancer (6). In comparison to abundant lncRNAs identified by genome-wide studies, functionally lncRNAs require better characterization in ovarian cancer.
Autophagy, an evolutionarily conserved process, maintains cell homeostasis through a lysosomal degradation system that supports cell survival and also maintains homeostasis under various types of stress (7). Autophagy-based cell deaths provide molecular mechanisms and clinical implications upon ovarian cancer therapy (8). Thus, it is of significance to identify key regulators of autophagy for theoretical basis and clinical practice. Several studies have found that autophagy can be mediated by lncRNA regulators in ovarian cancer (9)(10)(11)(12). For instance, GAS8-AS1 inhibits the progress of ovarian cancer by activation of autophagy through binding with Beclin1 (9). Moreover, lncRNA TUG1 induces autophagy-related paclitaxel resistance via sponging miR-29b-3p in ovarian cancer (10). Also, silencing HOTAIR enhances the sensitivity to cisplatin in ovarian cancer via inhibition of cisplatin-induced autophagy (11). LncRNA highly upregulated in liver cancer exerts a carcinogenic effect Abbreviations: LncRNAs, long non-coding RNAs; TCGA, The Cancer Genome Atlas; GTEx, Genotype-Tissue Expression; FC, fold change; LASSO, least absolute shrinkage and selection operator; OS, overall survival; HR, hazard ratio; ROC, receiver operating characteristic; ESTIMATE, Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data. via targeting autophagy-related genes ATG7 and ITGB1 in an epithelial ovarian cancer (12). Thus, autophagy-related lncRNAs possess potential as prognostic indicators and therapeutic targets.
Ovarian cancer represents an insidious malignancy, which usually develops asymptomatically to late stages along with metastases, chemoresistance, and undesirable clinical outcomes (13). Autophagy is a key bioprocess during the initiation and progression of ovarian cancer (14). LncRNA regulators are involved in the autophagy process. There is still a lack of systematic analysis for identifying autophagy-related lncRNA signature for prediction of the survival outcomes of the patients with ovarian cancer. Herein, this study developed a prognostic autophagy-related lncRNA signature for ovarian cancer.

Data Retrieval and Pre-processing
Among all databases, only the Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) database has the RNA-seq expression profiling and corresponding clinical and prognostic information of ovarian cancer. Thus, we curated transcriptome data and clinical information containing age, survival time, recurrence, survival status, histologic grade, and pathologic stage of 379 ovarian cancer tissues from the TCGA database. Mutation annotation format (MAF) files of somatic mutation data of ovarian cancer samples were also downloaded from the TCGA database. Meanwhile, the RNA-seq expression profiles of 133 normal ovarian samples were obtained from the genotype tissue expression (GTEx; https://toil.xenahubs.net/download/ GTEX_phenotype.gz) database (15). The RNA-seq counts values from the TCGA and GTEx databases were normalized and preprocessed by the TCGAbiolinks package (16). Based on the HUGO Gene Nomenclature Committee (HGNC; http://www. gene.ucl.ac.uk/cgi-bin/nomenclature/searchgenes.pl), lncRNAs and mRNAs were annotated (17).

Acquisition of Abnormally Expressed lncRNAs
Abnormally expressed lncRNAs between ovarian cancer and normal ovarian specimens were screened utilizing edgeR package (http://bioconductor.org) based on gene expression data (18). Adjusted p ≤ 0.05 and |log 2 fold change (FC)| ≥ 1 were set as the criteria values of abnormally expressed lncRNAs.

Acquisition of Autophagy-Related lncRNAs
A total of 232 autophagy-related genes were retrieved from the Human Autophagy Database (HADb; http://www.autophagy. lu/). The correlation between abnormally expressed lncRNAs and autophagy-related genes was analyzed by psych package (https:// CRAN.R-project.org/package=psych) with Pearson's correlation analysis. Autophagy-related lncRNAs were screened with the criteria values of |correlation coefficient| > 0.4 and p < 0.001.

Construction of a LASSO Regression Model
By adopting the least absolute shrinkage and selection operator (LASSO) Cox regression analysis, this study constructed an autophagy-related lncRNA model utilizing the glmnet package (19). The association between the lncRNAs in this model and overall survival (OS) was assessed by univariate Cox regression analysis. LncRNAs with hazard ratio (HR) > 1 and p < 0.05 were risk factors, while those with HR < 1 and p < 0.05 were protective factors. The risk score of each sample was determined. The formula was as follows: risk score = coefficient (lncRNA1) × expression (lncRNA1) + coefficient (lncRNA2) × expression (lncRNA2) + . . . + coefficient (lncRNAn) × expression (lncRNAn). The distributions of risk scores, survival status, and disease progress were assessed in the samples of ovarian cancer.

Assessment of the Prognostic Model
All the subjects with ovarian cancer were randomly separated into the training set and the validation set at a ratio of 1:1. The subjects with ovarian cancer were separated into highand low-risk groups. The Kaplan-Meier survival analyses were carried out to assess the OS differences between high-and lowrisk groups with survival package in the training set, validation set, and the whole dataset. Time-dependent receiver operating characteristic (ROCs) curves were depicted for estimating the predictive efficacy of survival time through risk score and other clinical features (age, pathologic stage, and histologic grade) utilizing the survival ROC package. Furthermore, the differences in the risk scores between patients with different histologic grades or pathologic stages were analyzed by oneway ANOVA. Multivariate Cox regression analysis was presented for evaluating whether the risk score could be independent of other clinical features including pathologic stage and histologic grade.

Estimation of Stromal and Immune Cells in Malignant Tumor Tissues Using Expression Data
The fractions of stromal and immune cells were inferred in each ovarian cancer sample using the ESTIMATE algorithm (https:// sourceforge.net/projects/estimateproject/) (21). The differences in immune and stromal scores were assessed between high-and low-risk samples of ovarian cancer by Student' t-test.

Sensitivity to Chemotherapy Drugs
A total of 94 ovarian cancer-related chemotherapy drugs were collected from the Cancer Genome Project (CGP; version 2014). The pRRophetic package was employed to predict the clinical chemotherapeutic responses based on the gene expression Frontiers in Medicine | www.frontiersin.org profiles of ovarian cancer (22). The IC50 values of chemotherapy drugs were compared between high-and low-risk groups utilizing Student's t-test.

Somatic Mutation Analysis
The "Masked Somatic Mutation" data were obtained and processed with VarScan software (version 2), a method for the detection of somatic mutation and copy number variation in exome data (23). The maftools package (24) was adopted to analyze the MAF of the somatic variants. The "titv" function classified single nucleotide polymorphisms (SNPs) into transitions (Ti) and transversions (Tv). The overall distribution of the six different SNVs and the proportion of transitions in each sample were evaluated. The oncoplot of the top 30 mutated genes was depicted by the "ComplexHeatmap" function. The "somaticInteractions" function was used to detect mutually exclusive or co-occurring genomes, which were tested by Fisher's exact test.

Gene Set Enrichment Analysis
Pathways underlying the autophagy-related lncRNA signature were evaluated through GSEA package (25). The "c5.bp.v6.2.symbols.gm" gene set was curated from the Molecular Signatures Database, which was employed as a reference set. Terms with |nominal enrichment score (NES)| > 1.7 and nominal p < 0.05 were significantly enriched.

Screening Abnormally Expressed lncRNAs
Herein, we retrieved RNA expression profiles of ovarian cancer and normal ovarian tissues from TCGA and GTEx datasets. Totally, 590 lncRNAs with adjusted p ≤ 0.05 and |log 2 FC| ≥ 1 were abnormally expressed in 379 specimens of ovarian cancer in comparison with the normal ovarian cancer (Figures 1A-C). There were 373 upregulated (Supplementary Table 1) and 217 downregulated (Supplementary Table 2) lncRNAs for ovarian cancer. We separately displayed the top 20 upregulated ( Table 1) and downregulated ( Table 2) lncRNAs for ovarian cancer ( Figure 1D).

Construction of an Autophagy-Related lncRNA Prognostic Signature in Ovarian Cancer
The correlation between abnormally expressed lncRNAs and autophagy-related genes was evaluated by Pearson's correlation analysis, which was visualized into the heat map (Figure 2A; Supplementary Table 3). With the criteria values of |correlation coefficient| > 0.4 and p < 0.001, we selected 66 autophagy-related abnormally expressed lncRNAs, as shown in Table 3. Using LASSO regression analysis, we determined the regression coefficients of the 14 lncRNAs in the model (Figures 2B,C). The risk score of each specimen of ovarian cancer was calculated as follows: AC084018.    (Figure 2D). In Figure 2E, we visualized the distributions of the risk scores of the patients with ovarian cancer. Also, we found that the risk scores displayed associations with survival status ( Figure 2F) and disease progression ( Figure 2G).

Evaluation of the Autophagy-Related lncRNA Signature as a Robust Prognostic Factor in Ovarian Cancer
Herein, all the patients with ovarian cancer were randomly separated into the training set and validation set (both n = 187) in Table 4. Based on the median value of the risk scores, the patients were divided into high-and low-risk groups. Our data showed that there were distinct differences in OS time between the two groups in the training set (p = 3.08e-10), validation set (p = 1.91e-03), and the whole dataset (p = 9.9e-10; Figures 3A-C). Patients in the low-risk group had prolonged OS duration in comparison with those in the highrisk group. ROCs were conducted to validate the predictive performance of this signature. The area under the curves under 7-year survival were 0.717, 0.747, and 0.727 in the training set, validation set, and the whole dataset, respectively (Figures 3D-F).

The Autophagy-Related lncRNA Signature as an Independent Prognostic Factor for Ovarian Cancer
We further evaluated the correlations between the risk scores and other clinical features in specimens with ovarian cancer. As a result, lowered risk scores were detected in grade 3 than grade 2 (p = 8.55e-04; Figure 4A). Meanwhile, we found that the risk scores increased gradually as the pathologic stages increased (Figure 4B). The AUCs under 3-, 5-, and 7-year survival time were 0.602, 0.651, and 0.727 for the patients with ovarian cancer, indicating the well-predictive performance of this signature (Figure 4C). Furthermore, compared with the other clinical characteristics including age (AUC = 0.536), pathologic stage (AUC = 0.497), and histologic grade (AUC = 0.508), there was a higher AUC value under 7-year OS time for the risk score ( Figure 4D). The data suggested that this risk score might possess higher sensitivity and accuracy in predicting the prognosis of the patients with ovarian cancer. Multivariate Cox regression analysis demonstrated that this risk score could be independently predictive of the prognosis of the patients (HR: 2.31, 95% CI: 1.75-3.03, p = 2.64e-09; Figure 4E).

The Autophagy-Related lncRNA Signature Is Associated With Immune Cell Infiltration in Ovarian Cancer
Here, we adopted the CIBERSORT algorithm to infer the immune cell infiltration in samples with ovarian cancer. Figure 5A visualized the proportions of 22 kinds of immune cell components in ovarian cancer tissues. Also, we analyzed the correlations between distinct immune cells. In Figure 5B, there were strong correlations between B cells naïve and macrophages M2 (r = 0.93), between B cells memory and monocytes (r = 0.98), between T cells CD4 naïve and T cells CD4 memory resting (r = 0.98), between T cells follicular helper and Tregs (r = 0.97), between T cells follicular helper and T cells gamma delta (r = 0.96), between Tregs and macrophages M0 (r = 0.97), between T cells gamma delta and mast cells activated (r = 0.92), and between monocytes and mast cells activated (r = 0.95) in ovarian cancer tissues. Also, heatmap visualized the differences in the immune cell infiltrations between high-and low-risk samples of the ovarian cancer ( Figure 5C). Compared with the low-risk group, there were lowered infiltration levels of macrophages M1 (p < 0.0001), mast cells resting (p = 0.007), plasma cells (p = 0.003), T-cell CD8 (p = 0.016), and T-cell follicular helper (p = 0.001) and also higher infiltration levels of macrophages M2 (p = 0.003), mast cells activated (p < 0.0001), and neutrophils (p = 0.036) in the high-risk group (Figure 5D).  The Autophagy-Related lncRNA Signature Is Associated With Immune Checkpoints in Ovarian Cancer Herein, we observed whether the autophagy-related lncRNA signature was associated with the expression of immune checkpoints in ovarian cancer tissues. Higher LSECtin expression was found in the high-risk group compared with the low-risk group ( Figure 6A). Furthermore, there was distinctly decreased PD-L1 expression in the low-risk group compared with the high-risk group (p = 2.79e-04; Figure 6B). We also assessed the correlations between the risk scores and immune or stromal scores. As a result, higher immune scores and lower stromal scores were detected in the high-risk samples compared with the low-risk samples, which were not statistically significant (Figures 6C,D).

Somatic Mutation Landscapes in Ovarian Cancer
We further analyzed the somatic mutation landscapes in specimens of ovarian cancer. Both in the high-( Figure 8A) and low-risk ( Figure 8B) samples, missense mutation was the most common mutation type. TP53 was the top-ranked mutated gene, followed by TTN. C > T mutation had the highest proportion in the two groups (Figures 8C,D). Furthermore, potential druggable categories targeted mutated genes were predicted, such as the druggable genome, clinically actionable, and kinase in the high-( Figure 8E) and low-risk (Figure 8F) samples. Approximately 98.48% of samples occurred genetic mutations in the high-risk samples ( Figure 9A) and 99.28% of the occurred mutations in the low-risk samples (Figure 9B). Many genes co-occur in cancer or show strong exclusivity in their mutation patterns. Here, we visualized the close interactions between the mutated genes in the high-and low-risk samples (Figures 9C,D).

Significant Signaling Pathways Underlying the Autophagy-Related lncRNA Signature
To explore the significant signaling pathways underlying the autophagy-related lncRNA signature, this study carried out GSEA by comparing high-and low-risk groups. As shown in
As immunogenic cancers, the spontaneous anticancer immune responses may increase survival duration, and the immune escape may cut down survival duration (32). In the recent years, several immune-based strategies such as immune checkpoint inhibition, vaccination, and antigenspecific active immunotherapy have been developed in ovarian cancer (33). The immune-suppressive networks in the tumor microenvironment have been considered for immunotherapy implementation (34). The immune-related markers may be utilized for predicting the responses to immunotherapy (35). Hence, the interactions between molecules and tumor microenvironment require in-depth exploration. Our study demonstrated that the autophagy-related lncRNA signature was in relation to infiltrations of macrophages M1, mast cells resting, plasma cells, T-cell CD8, T-cell follicular helper, macrophages M2, mast cells activated, and neutrophils, indicating the crosstalk between these autophagy-related lncRNAs and immune cells. Consistently, Deng et al. identified a novel autophagy-related lncRNA model that was distinctly related to the infiltrations of immune cells in pancreatic cancer (36). Ovarian cancer represents a low immune-reactive malignancy with restricted immune cell infiltrations and extensive immunosuppressive T cell infiltrations (37). Tumors positive for PD-L1 usually exhibit a higher response to immune checkpoint inhibition therapy, and highly expressed PD-L1 is a predictor of undesirable clinical outcomes (37). Here, our data showed that the low-risk ovarian cancer patients had higher PD-L1 expression in comparison to those with high risk.
Chemotherapy (platinum and taxanes) plays a fundamental role in adjuvant therapy against ovarian cancer (38). Despite the initial response to this therapy, most of the patients diagnosed with ovarian cancer developed chemotherapy resistance (39). This resistance may be driven by a range of mechanisms. Hence, it is of importance to develop individual markers for the prediction of the sensitivity to chemotherapy. Our data showed that the risk scores were in relation to 29 chemotherapy drugs such as cisplatin, indicating that the autophagy-related lncRNAs were involved in chemotherapy sensitivity, as a previous study (11).

CONCLUSIONS
Collectively, this study provided a knowledge base of novel autophagy-related lncRNAs in ovarian cancer, improving the understanding of the functions of lncRNA on the regulation of autophagy in ovarian cancer. We established an autophagyrelated lncRNA signature as a robust prognostic marker for the prediction of survival outcomes, immunotherapy response, and chemotherapy sensitivity. Our findings may assist to precisely guide therapeutic strategies for individual patients with ovarian cancer in clinical practice.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.