Identification of 10 differentially expressed genes involved in the tumorigenesis of cervical cancer via next-generation sequencing

Uterine cervical specimens

Uterine cervical specimens were separated immediately from the cervix of 18 cervical cancer patients after radical hysterectomy. The distance between cancer sample location and normal sample location was at least 2 cm. The patients were selected following the criteria below. Inclusion criteria: (1) Those who were eligible for this study and signed the informed consent; (2) Women aged 18–75 years; (3) Preoperative stage was IB1-IIA2, and surgery was the best treatment choice for them; (4) Those who had macroscopic cervical cancer tissues and adjacent normal tissues; (5) Those who were able to cooperate according to the research protocol. Exclusion criteria: (1) Pregnant or lactating women; (2) Preoperative stage was IIB or severe, and the preferred treatment was radiotherapy or chemotherapy; (3) Those whose cervical cancer tissues and adjacent normal tissues from the isolated uterine could not be distinguished; (4) Patients with severe circulatory, immune, or respiratory diseases; (5) Those who had other uncontrolled medical conditions that the researchers deem unsuitable for enrollment; (6) Any situation that the investigators considered may increase the risk to the subject or interfere with the results. Clinical information of the three patients whose specimens were used for RNA sequencing is shown in Table 1. Table 2 exhibited the key medical records of other 15 cases whose specimens were used for confirmatory experiments. The disease status was assigned according to the 2018 FIGO classification. Paired adjacent normal tissues and cancer tissues were stored long-term at −80 °C without any additive for subsequent experiments. The study was conducted in accordance with the principles of the Declaration of Helsinki. All the specimens were obtained after approval from the Ethics Committee of Seventh Medical Center of Chinese PLA General Hospital (S2023-025-01) and written informed consent was obtained from patients in the Department of Obstetrics and Gynecology of the Seventh Medical Center of PLA General Hospital.

RNA sequencing dataset establishment

Total RNA of paired adjacent normal tissues and cancer tissues from three patients was extracted by Trizol (Invitrogen, Waltham, MA, USA) according to the manufacturer’s protocol. Therefore, we had three biological replicates. Three technical replicates were completed simultaneously. The extracted RNA of each sample was sequenced three times. Total amounts and integrity of RNA were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Then, total RNA was used as input material for the library preparation. Messenger RNA was purified from total RNA by using poly-T oligo-attached magnetic beads. After the construction of the library, the library was initially quantified by Qubit2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, MA, USA), then diluted to 1.5 ng/μl, and the insert size of the library was detected by Agilent 2100 bioanalyzer. After insert size met the expectation, RT-qPCR was used to accurately quantify the effective concentration of the library (the effective concentration of the library is higher than that of 2 nM) to ensure the quality of the library. After the library was qualified, the different libraries were pooled according to the effective concentration and the target amount of data off the machine, then sequenced by the Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA) (Fig. 1A). The end reading of 150 bp pairing is generated. Illumina Casava1.8 software (Illumina, San Diego, CA, USA) was used for basecalling. The basic principle of sequencing is to synthesize and sequence at the same time (Sequencing by Synthesis). After quality control for the reads, index of the reference genome (hg38) was built using Hisat2 (v2.0.5) and paired-end clean reads were aligned to the reference genome using Hisat2 (v2.0.5). FeatureCounts (v1.5.0-p3) was used to count the reads numbers mapped to each gene. The process of building RNA sequencing datasets was similar to what we described in our previous work (Zhao et al., 2018). Sequence data of this study have been deposited in the Gene Expression Omnibus (GEO) repository with the primary accession record GSE253690. The sequencing depth is 35.099995, and the statistical power of this experimental design calculated in R environment is 0.140838.

Identification of DEGs and functional enrichment analysis

Differential expression analysis of two groups (paired adjacent normal cervical tissues and cervical cancer tissues) of three patients was performed respectively using the edgeR package (Version 3.22.5) in the R environment (Version 4.2.1) as described in Robinson, McCarthy & Smyth (2010). The integrated DEGs union between normal group and cancer group of all three patients were analyzed using the DESeq2 package (Version 1.20.0) following the method in Love, Huber & Anders (2014) and shown in Volcano diagram. The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. We defined padj ≤ 0.05 and |log₂(foldchange)| ≥ 1 as the thresholds for significantly differential expression. We selected the top 500 DEGs of each patient through edgeR package. There were 47 intersecting genes in Venn diagram, and we selected the top 10 genes for next-step public data and experimental analysis. The heatmap of 47 genes and 10 genes were plotted using pheatmap package and ggplot2 package (Version 3.3.6).

The integrated DEGs union of the 3 patients were used for enrichment analysis. Gene Ontology (GO) enrichment analysis was implemented by the clusterProfiler R package (Version 3.8.1), in which gene length bias was corrected. GO terms with corrected p values less than 0.05 were considered significantly enriched with DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological systems, especially large-scale molecular datasets generated by genome sequencing and other high throughput experimental technologies. We used clusterProfiler R package (Version 3.8.1) to test the statistical enrichment of DEGs in KEGG pathways. The Disease Gene Network (DisGeNET) database integrates human disease-related genes. DisGeNET pathways with corrected p values less than 0.05 were considered significantly enriched by DEGs. The Reactome database brings together the various reactions and biological pathways of human model species. Reactome pathways with corrected p values less than 0.05 were considered significantly enriched by DEGs. We also used clusterProfiler R package (Version 3.8.1) to test the statistical enrichment of DEGs in the Reactome pathway and the DisGeNET pathway. The functional enrichment analysis was performed as previously described in Yu et al. (2012).

Public RNA sequencing data analysis

Public data (RNA sequencing data and clinical information) were requested from TCGA and GTEx databases provided on the UCSC Xena website (http://xenabrowser.net/datapages/). The FPKM-form RNA sequencing data were uniformly treated through Toil process (Vivian et al., 2017). The expression of selected DEGs was transformed into log₂(FPKM+1) values to compare the differences. The receiver operating characteristic (ROC) curves were analyzed by pROC package (Version 1.18.0) which was described in Robin et al. (2011). Diagnostic values were measured by area under the curve (AUC) values. Survival analysis was finished using overall survival (OS) time through survival package (Version 3.3.1) described in Groeneveld et al. (2019). Survival outcomes referred to data from the literature reported in 2018 (Liu et al., 2018). The groups were divided by the median value of the expression level of each DEG.

Quantitative real time polymerase chain reaction

Total RNA was extracted from 20–30 mg paired adjacent normal cervical tissues and cervical cancer tissues through macrodissection according to the manufacturer’s protocol of total RNA rapid extraction kit (BIOMAN, Beijing, China) within 1 week. All the equipments were purified with DEPC water. RNase free water was used for the dilution of RNA. The adsorption columns in the kit were used to remove DNA and purify RNA. The quality of RNA was assessed by the A260/A280 and A260/A230 value of each sample using Multiskan Skyhigh (Thermo Fisher Scientific, Waltham, MA, USA) and all the values were conformed to the standard. RNA (200 ng for each sample) was reverse-transcribed into cDNA by reverse transcriptase and random primer (N9) (0.1 μg/μl) in a 20 μl reaction volume with transcript one-step gDNA removal and cDNA synthesis supermix kit (TRANSGEN, Beijing, China). The reverse-transcribed reaction conditions included 25 °C for 10 min, 42 °C for 15 min and 85 °C for 5 s. cDNA was stored at −20 °C and was used for qPCR within 1 week. RT-qPCR was performed using the PerfectStart Green qPCR Supermix kit (TRANSGEN, Beijing, China) in 0.2 mL Polypropylene PCR Tube Strips (Axygen, Union City, CA, USA) via QuantStudio 3 Real-Time PCR System (ABI, Loma Linda, CA, USA). QuantStudio Design & Analysis Software (ABI, Waltham, MA, USA) was used for manual reaction setup and analysis. The PCR conditions included a denaturing step at 94 °C for 30 s and 45 cycles of 94 °C for 5 s and 60 °C for 30 s for real time plate read. The melt curve stage included 95 °C for 15 s, 60 °C for 1 min and 95 °C for 1 s. The whole reaction volume was 20 μl, and the concentration of forward and reverse primers were 10 μM. SYBR Green I was the fluorescent dye applied in this assay. GAPDH was used for data normalization. The 2^−ΔΔCt value was regarded as the expression level of each gene. The RT-qPCR was performed in investigator’s lab and repeated three times similar to what we described in our previous work (Jin et al., 2019). The results showed well repeatability. Each melt curve has only one peak and we found little outliers. Targeted gene information and specific primer sequences were described in Table S1. The manufacturer of oligonucleotides was Sangon (Shanghai, China).

Statistical analysis

RNA sequencing and public database data were conducted using the R environment (Version 4.2.1). Padj ≤ 0.05 and |log₂(foldchange)| ≥ 1 were the thresholds for DEGs. Corrected p values < 0.05 were considered significantly enriched in function enrichment analysis. Quantitative data were analyzed using Shapiro-Wilk normality test and Levene’s test. The expression level difference was calculated by Mann–Whitney U test using stats package (Version 4.2.1) and car package (Version 3.1-0). The ROC curves were analyzed by pROC package (Version 1.18.0). Cox proportional hazard model was tested and survival curves were matched through survival package (Version 3.3.1) using Kaplan-Meier methodology. Then the difference in OS of cervical cancer patients stratified by target gene expression level was tested by the Log-rank test. Data visualization was accomplished via survminer package and ggplot2 package (Version 3.3.6). The comparison of 10 DEGs expression between paired normal and cancer samples from clinical sampling in RT-qPCR was analyzed by paired t test and plotted in Graphpad Prism 9 (GRAPHPAD SOFTWARE, CA, USA). p < 0.05 was considered statistically significant.

Identification of DEGs

To discuss the potential mechanisms of cervical cancer tumorigenesis and identify the DEGs involved in these processes, we extracted the specimens from three patients (stage IB2 to IIA2) after the radical hysterectomy. Their clinical information was shown in Table 1. Then, RNA sequences were determined by NGS. The up-regulated and down-regulated DEGs were sorted out, respectively (Fig. 1A). According to the integrated NGS data of three patients, there were 875 up-regulated and 1,482 down-regulated genes in cervical cancer tissues compared to paired adjacent normal cervical tissues (Fig. 1B). DEGs were reported only if the log-fold change was >1 and the adjusted p value was less than 0.05. We focused on the top 500 DEGs of each patient, and found that their intersection contained 47 DEGs (Fig. 1C). The heatmap demonstrated that only five DEGs expressions were increased in cancer specimens, while other 42 DEGs expressions were obviously decreased (Fig. 1D). We chose the top 10 significant DEGs in the 47 intersectional DEGs for further research (Fig. 1E) and we marked their locations in the integrated DEGs union between normal group and cancer group of three patients (Fig. 1B). These genes included APOD, MASP1, ACKR1, C1QTNF7, SFRP4, HSPB6, GSTM5, IGFBP6, F10 and DCN. All of them expressed at lower levels in cervical cancer samples than in normal cervix cells.

Functional enrichment analysis of DEGs

Functional enrichment analysis revealed the potential mechanisms of the DEGs union between normal group and cancer group of three patients in the development of cervical cancer. GO analysis is associated with cellular component (CC), molecular function (MF) and biological process (BP). It was indicated that the incidence of malignant cervical tumor could be attributed to extracellular matrix, angiogenesis and circulatory system process (Fig. 2A). DisGeNET database contains a collection of genes associated with abundant human diseases. The results showed that 2357 DEGs (including 875 up-regulated DEGs and 1482 down-regulated DEGs) were more likely to generate tumor angiogenesis, systemic scleroderma and malignant neoplasm of thyroid (Fig. 2B). KEGG database reflected the reaction network of molecules and pathways of diseases. We found that PI3K-Akt signaling pathway, calcium signaling pathway, neuroactive ligand-receptor interaction, cell cycle and focal adhesion may participate in the occurrence and progression of cervical cancer (Fig. 2C). Changes in extracellular matrix organization and cell cycle checkpoints were the two most important processes associated with cervical tumors on the basis of Reactome database (Fig. 2D).

The expression levels of selected DEGs in public database

The TCGA program is the largest public database of cancer genetic information worldwide. The database provides multi-omics hereditary information of cancer patients. The GTEx project collected the sequencing data of healthy human donors. In this study, we focused only on the transcriptomic information of 10 selected DEGs to verify their expression levels in a larger population. All the selected genes (including APOD, MASP1, ACKR1, C1QTNF7, SFRP4, HSPB6, GSTM5, IGFBP6, F10 and DCN) showed lower expression levels in cervical cancer samples, and these results were statistically significant (Figs. 3A–3J).

The diagnostic performance of selected DEGs in public database

To assess the performance of these selected DEGs as biomarkers, ROC curves were generated to measure the AUC values. All the AUC values exceeded 0.90, which meant that the selected DEGs have chances to function as excellent diagnostic biomarkers (Figs. 4A–4J). In particular, the AUC values of MASP1, GSTM5 and C1QTNF7 were 0.998, 0.998 and 0.997, respectively (Figs. 4B, 4G, 4D). Other DEGs whose AUC values were greater than 0.99 included APOD, HSPB6, F10 and DCN (Figs. 4A, 4F, 4I, 4J).

Survival analysis of selected DEGs in the TCGA program

Survival analysis was performed according to the OS time. The groups were divided by the median expression value of each target gene (Figs. 5A–5J). Only the expression levels of APOD, ACKR1 and SFRP4 contributed to the variation in OS. Patients who had higher expression of APOD, ACKR1 and SFRP4 might have a longer survival time (Figs. 5A, 5C, 5E).

Selected DEGs were down-regulated in cervical cancer tissues

Fifteen pairs of cervical cancer tissues and adjacent normal cervical tissues were collected after the radical hysterectomy, and their brief clinical information were shown in Table 2. We performed RT-qPCR to determine the expressions of 10 selected DEGs in various patients (Figs. 6A–6J). GAPDH was used as the normalization gene. Similar to the analysis results of public database, the expression levels of most DEGs were down-regulated in cervical cancer tissues compared to the paired normal cervical tissues (Figs. 6A–6J). In Cases 2, 3, 8, 10, 11, 12 and 15, the expression levels of all 10 selected DEGs decreased significantly in cervical cancer tissues (Figs. S1B, S1C, S1H, S1J, S1K, S1L, S1O). The expression of DCN was almost unchanged in Case 5 (Fig. S1E), and IGFBP6 expression remained stable in both normal and cancer tissues of Case 14 (Fig. S1N). However, the expression levels of other DEGs in Case 5 and 14 were lower in cancer tissues than that in normal tissues (Figs. S1E and S1N). All data in Fig. 6 and S1 showed that the expression levels of APOD, MASP1, ACKR1 and SFRP4 were extremely consistent in most cases, demonstrating that they have a potential to be biomarkers for cervical cancer diagnosis.