High-grade serous cancer (HGSC) is the most common subtype of ovarian cancer. HGSC is highly aggressive with poor patient outcomes, and a deeper understanding of HGSC tumorigenesis could help guide future treatment development. To systematically characterize the underlying pathologic mechanisms and intratumoral heterogeneity in human HGSC, we used an optimized single-cell multiomics sequencing technology to simultaneously analyze somatic copy-number alterations (SCNA), DNA methylation, chromatin accessibility, and transcriptome in individual cancer cells. Genes associated with interferon signaling, metallothioneins, and metabolism were commonly upregulated in ovarian cancer cells. Integrated multiomics analyses revealed that upregulation of interferon signaling and metallothioneins was influenced by both demethylation of their promoters and hypomethylation of satellites and LINE1, and potential key transcription factors regulating glycolysis using chromatin accessibility data were uncovered. In addition, gene expression and DNA methylation displayed similar patterns in matched primary and abdominal metastatic tumor cells of the same genetic lineage, suggesting that metastatic cells potentially preexist in the subclones of primary tumors. Finally, the lineages of cancer cells with higher residual DNA methylation levels and upregulated expression of CCN1 and HSP90AA1 presented greater metastatic potential. This study characterizes the critical genetic, epigenetic, and transcriptomic features and their mutual regulatory relationships in ovarian cancer, providing valuable resources for identifying new molecular mechanisms and potential therapeutic targets for HGSC.

Significance:

Integrated analysis of multiomic changes and epigenetic regulation in high-grade serous ovarian cancer provides insights into the molecular characteristics of this disease, which could help improve diagnosis and treatment.

Ovarian cancer is one of the leading causes of gynecologic cancer death (1), due in part to its high heterogeneity. High-grade serous cancer (HGSC) is the most aggressive and common subtype of ovarian cancer (comprising approximately 70% of all cases). HGSC is generally diagnosed at an advanced stage, presenting with extensive metastases and massive ascites. Despite an initial clinical response to the surgery and chemotherapy, the majority of patients will relapse with fatal outcomes. Hence, understanding the pathologic mechanisms of HGSC initiation, promotion, and progression will be of great value.

Intratumor heterogeneity (ITH) is one of the most challenging obstacles to study HGSC (2). Previous studies have revealed the heterogeneity of HGSC in single-nucleotide variations (SNV), somatic copy-number alterations (SCNA), aberrant transcriptomic programs, and tumor subtypes (3–5). However, these studies are usually based on bulk sequencing, which permits only tissue-level resolution and has limited ability to illustrate ITH. In addition, tumor development involves not only alterations in the genome and transcriptome, but also extensive epigenetic changes, which remain to be comprehensively analyzed in HGSC.

Single-cell sequencing provides a powerful tool for the study of ITH. Most studies have used single-cell RNA sequencing (scRNA-seq) to characterize different phenotypic cell states of cancer cells (6–9). However, to accurately explore the dynamic changes of tumor promotion (i.e., the changes between normal tissues and primary tumors) and progression (i.e., the changes between primary tumors and metastases), reconstruction of the genetic lineages of cancer cells is necessary. Naturally occurring genomic mutations, like SNVs and SCNAs, leave the real signatures to trace evolutionary histories of cancer cells, which can be detected by single-cell DNA sequencing. Besides tracing the genetic lineages and revealing the phenotypic changes during tumorigenesis, we also wondered the epigenetic alterations and their regulation on phenotypes during tumor promotion and progression. Thus, we improved single-cell multiomics sequencing technology (scCOOL-seq; ref. 10) to simultaneously assess the SCNAs, DNA methylome, chromatin accessibility, and transcriptome in the same individual cells. Using improved scCOOL-seq, we reconstructed the genetic lineages using SCNAs. Tracing the genetic lineages, we characterized the multiomics changes as well as the epigenetic regulation of phenotypes from normal cells to primary tumor cells and further to matched metastatic cells. Our study offers deeper insights for HGSC and provides a valuable resource for developing new therapeutic strategies.

Study approval

This study was approved by the ethics committee of the Peking University People's Hospital (License No. 2019HPB034-01) and was conducted in accordance with the Declaration of Helsinki. All patients provided the written informed consent before the surgery.

Clinical samples

No patient received any tumor-related treatment before surgery. Fallopian tubes (FT) were obtained from patients who underwent hysterectomy with bilateral salpingo-oophorectomy (BSO) for benign and malignant gynecologic diseases (patient OC20, OC21, OC25, OC26, and OC28) and risk-reducing salpingo-oophorectomy (RRSO) surgery (patient OC22). Only epithelial cells obtained from pathologically normal FT tissues and proved to have no SCNAs using improved scCOOL-seq were considered as normal control cells.

Single cell preparation

The fresh specimens were digested into single-cell suspensions using the MACS Tissue Dissociation Kit (Miltenyi Biotec, #130-095-929). Epithelial cells of tumor tissues were isolated by magnetic activated cell sorting (MACS; CD326 Microbeads, Miltenyi Biotec, #130-061-101, RRID: AB_2832928). The live normal fallopian tube epithelial (FTE) cells were sorted by FACS using the anti-human CD326 antibody (BioLegend, #324208, RRID: AB_756082) and 7AAD (BD Pharmingen, #559925).

Library construction of improved scCOOL-seq

The sorted single cells were first treated by the GpC methyltransferase (NEB, #M0227L). Next, each cell was incubated with prewashed carboxylic acid dynabeads (Invitrogen, #65011) to separate the nuclei from the RNA-containing cytoplasm according to the scTrio-seq2 protocol (11).

The supernatant containing mRNA was processed following the modified STRT-seq protocol (12). The nuclear DNA fraction was processed following a previously described method (13). Briefly, the nuclei with magnetic Dynabeads were treated with bisulfite (Zymo Research, #D5044) and were then amplified using random primers with anchor sequence. After purification by AMPure XP beads, the products were amplified by index primers for 16 cycles. The libraries were finally sequenced on the Illumina HiSeq 4000 platform.

Gene knockout, gene overexpression, and cell viability assay

The human ovarian cancer cell line SKOV3 (RRID: CVCL_0532) and ES2 (RRID: CVCL_3509) were purchased from ATCC. Cell lines were authenticated by short tandem repeat (STR) analysis and were mycoplasma negative. Cells within 20 passages were used in the study. We used the clustered regularly interspaced short palindromic repeats (CRISPR)-Cas9 to knockout CCN1. The guide RNA (gRNA) sequences are as follows: gRNA1: CGCGCACTTGGGCGCCTCCA; gRNA2: GCAGATCCCCTTCAGAGCGG. Overexpression of HSPA6 was performed by transfecting the overexpression plasmid using Lipofectamine 3000 Reagent (Invitrogen, #L3000150) as described previously (14). The HSP90AA1 inhibitor TAS-116 was purchased from Selleck Chemicals (#S7716). Cell viability was analyzed by the CellTiter-Glo 2.0 Assay (Promega, #G9241).

Wound healing assay

When cells grew to full confluence, a wound was scratched using a 200-μL micropipette tip. The cells were washed with PBS to remove the cell debris and then grown in serum-free medium. The wound healing results were recorded using a phase-contrast microscope (Leica), and the migration areas were calculated using the ImageJ software.

Transwell assays

Transwell assays were performed with transwell chambers (Coring, #3422) with or without coating with Matrigel (Corning, #356231) to measure the invasion and migration ability, respectively. Cells were harvested in serum-free medium and seeded into the upper chambers. After 24-hour incubation, the cells in the lower chambers were fixed and stained, and cells in the upper chambers were removed with a cotton swab. The photographs were taken using a microscope and the cell numbers were counted using ImageJ software.

Data availability

The raw data have been deposited in the Genome Sequence Archive (GSA) under the accession number HRA000360 and are available upon reasonable request. The processed data associated with SCNAs, DNA methylome, chromatin accessibility, and transcriptome have been submitted to the Gene Expression Omnibus (GEO) database under the accession number GSE189955.

Code availability

Customized codes for data analysis are available at https://github.com/hlxie.

Single-cell multiomics sequencing reveals tumor heterogeneity in HGSC

Here, using improved scCOOL-seq, we profiled normal fallopian tubes, primary tumors, and matched metastases to unveil the molecular characteristics associated with tumor development (Fig. 1A; Supplementary Fig. S1A). To cover ITH as much as possible, the tumor tissues were processed with a multiregional sampling strategy. To avoid variations resulting from treatment effects and disease stages, only patients who were treatment-naïve and had stage III or IV disease were enrolled in our study (Supplementary Table S1).

Figure 1.

Identification of cancer cells and normal fallopian tube epithelial cells. A, The workflow diagram illustrates the sampling and analysis strategies. B, Uniform manifold approximation and projection (UMAP) plot of all single cells we sequenced. Colors indicate cell types (left), patients (middle), and sampling regions (right). LN, lymph node; RL, round ligament. C, Global SCNA patterns inferred from RNA data of single cells. Each row of the heatmap represents a single cell. The color bar on the left represents the patient origin of each single cell. D, The global SCNA patterns profiled by DNA methylome data of single cells at 1-M resolution (top) and schematic diagram of the evolutionary histories of genetic lineages during tumorigenesis (bottom) for OC09. Each row of the top panel represents a single cell. Each column represents a chromosome. The representative focal aberrations within chromosome arms used to define genetic lineages are marked. The global SCNA patterns of other patients were introduced in Supplementary Fig. S2.

Figure 1.

Identification of cancer cells and normal fallopian tube epithelial cells. A, The workflow diagram illustrates the sampling and analysis strategies. B, Uniform manifold approximation and projection (UMAP) plot of all single cells we sequenced. Colors indicate cell types (left), patients (middle), and sampling regions (right). LN, lymph node; RL, round ligament. C, Global SCNA patterns inferred from RNA data of single cells. Each row of the heatmap represents a single cell. The color bar on the left represents the patient origin of each single cell. D, The global SCNA patterns profiled by DNA methylome data of single cells at 1-M resolution (top) and schematic diagram of the evolutionary histories of genetic lineages during tumorigenesis (bottom) for OC09. Each row of the top panel represents a single cell. Each column represents a chromosome. The representative focal aberrations within chromosome arms used to define genetic lineages are marked. The global SCNA patterns of other patients were introduced in Supplementary Fig. S2.

Close modal

We first sequenced the RNA compartment of single cells and obtained single-cell transcriptomes for a total of 5,329 cells (Supplementary Table S1). On average, 4,701 genes were detected per cell (Supplementary Fig. S1B). We distinguished cancer cells and normal FTE cells from stromal cells in three steps. First, according to clustering analysis and well-known cell type-specific marker gene expression, cells can be classified as epithelial cells, T cells, macrophages, and fibroblasts (Fig. 1B; Supplementary Fig. S1C). Second, we used the R package InferCNV (15) to infer SCNAs based on averaged expression levels across chromosomal intervals (Fig. 1C). Next, we sequenced the nuclear fractions for some of the epithelial cells from tumor tissues with inferred SCNAs and epithelial cells from FTs without inferred SCNAs. We profiled the SCNA patterns using the DNA methylome data of these selected cells (called “genomic SCNAs” for short; Fig. 1D; Supplementary Figs. S1D and S1E; Supplementary Fig. S2). Finally, we defined the epithelial cells with specific genomic SCNAs as cancer cells and epithelial FT cells without genomic SCNAs as normal FTE cells. Cancer cells display great interpatient heterogeneities as previously reported (6–9), whereas FTE cells and stromal cells exhibit low interpatient heterogeneities (Fig. 1B).

Reconstruction of the genetic lineages of cancer cells

We reconstructed genetic lineages of cancer cells in the same patient using the improved scCOOL-seq genome data. Diverse genetic alterations in cancer cells provide hallmarks for lineage tracing, such as SNVs and SCNAs. Considering that HGSC tumors frequently carry extensive SCNAs, we used SCNAs to construct the genetic lineages. SCNAs were classified into regional aberrations and focal aberrations (16). Focal aberrations frequently occur within chromosome arms and are generally considered irreversible (16, 17). Therefore, we used the breakpoints of subchromosomal focal aberrations to reconstruct the genetic lineages of cancer cells within the same patient (Fig. 1D; Supplementary Fig. S2; ref. 11). Cells containing more than three identical intra-chromosomal breakpoints were considered to belong to the same genetic lineage (11). Our data show that though HGSC carry extensive SCNAs, the composition of genetic lineages in HGSC is relatively simple, compared with colorectal cancer (11). The majority of HGSC patients (8/9) have only one or two lineages of cancer cells (Fig. 1D; Supplementary Fig. S2). Especially for OC16, more than 200 single cells collected from 10 sampling regions all belonged to the same lineage (Supplementary Fig. S2; Supplementary Table S1).

Dynamic changes of RNA expression from normal FTE cells to primary tumor cells

In recent years, the notion that HGSC originates from fallopian tube epithelium has been supported by increasing amounts of evidences (18). Thus, we used FTE cells as controls for the following analyses. According to the genetic lineages, all cancer cells in a given patient originated from the same ancestor cancer cell (Fig. 1D; Supplementary Fig. S2). Thus, we compared all primary tumor cells from the same patient with normal FTE cells to explore the phenotypic changes associated with tumor promotion.

As shown in Fig. 2A, 187 genes were downregulated in primary tumors. Gene ontology (GO) analysis shows that the pathways of eukaryotic translation process and protein ubiquitin were enriched (Fig. 2B), implying great alterations in protein levels. As an important organ to produce and response to hormones, the expression of genes related to hormone responses was also aberrant in cancer cells.

Figure 2.

Dynamic changes of RNA expression from normal FTE to primary tumors. A, Heatmap showing DEGs between normal FTE and primary tumor cells. Each row represents a DEG. The color bar at the top represents sampling regions and patient origins. B, The GO enrichment of DEGs. Red, upregulated genes in primary tumor cells; blue, downregulated genes in primary tumor cells. C and E, Scatterplots of the OXPHOS and glycolysis scores of cells in our study colored by sampling regions (C) and patients (E). The colors of patients are the same as A. D and F, 2D density plots of the OXPHOS and glycolysis scores of cells in our study (D) and patients in TCGA project (F). G, IHC determination of ISG15 expression in FTE, non-HGSC, and HGSC tumors. Scale bar, 40 μm.

Figure 2.

Dynamic changes of RNA expression from normal FTE to primary tumors. A, Heatmap showing DEGs between normal FTE and primary tumor cells. Each row represents a DEG. The color bar at the top represents sampling regions and patient origins. B, The GO enrichment of DEGs. Red, upregulated genes in primary tumor cells; blue, downregulated genes in primary tumor cells. C and E, Scatterplots of the OXPHOS and glycolysis scores of cells in our study colored by sampling regions (C) and patients (E). The colors of patients are the same as A. D and F, 2D density plots of the OXPHOS and glycolysis scores of cells in our study (D) and patients in TCGA project (F). G, IHC determination of ISG15 expression in FTE, non-HGSC, and HGSC tumors. Scale bar, 40 μm.

Close modal

On the other hand, 361 genes were upregulated in primary tumors (Fig. 2A). GO analysis revealed that these differentially expressed genes (DEG) were mainly enriched in pathways associated with VEGFA-VEGFR2 signaling (such as VEGFA, PGK1, and NRP2), metal ion homeostasis (such as MT1E, MT1F, MT1G, MT1H, MT1X, and MT2A), and regulation of viral process (such as ISG15, BST2, STAT1, and IFI27; Fig. 2B). Of note, genes associated with both of oxidative phosphorylation (OXPHOS; such as ATP5MC3, ATP5ME, and ATP6AP1) and aerobic glycolysis (such as ALDOA, ENO1, and LDHA) were upregulated. To investigate whether this phenomenon is due to the existence of two cell subtypes (one upregulated OXPHOS, and the other one upregulated glycolysis), we calculated OXPHOS and glycolysis scores for each cell (see Supplementary Materials and Methods). The FTE cells and cancer cells could be clearly distinguished, whereas cancer cells could not be further separated into subtypes and there were minute variations between patients, suggesting that cancer cells simultaneously elevated the expression of both OXPHOS and aerobic glycolytic genes compared with FTE cells (Fig. 2CE). Similarly, the patients of TCGA also cannot be divided into two groups according to the expression of OXPHOS and glycolysis (Fig. 2F).

To further clarify the typical characteristics of HGSC, we also compared HGSC with four non-HGSC tumors. We found that HGSC tumors could be clearly separated from non-HGSC tumors by principal component analysis (PCA), demonstrating that HGSC had very distinct transcriptional patterns (Supplementary Figs. S3A and S3B). Furthermore, the differential expression analysis shows that expression of metallothionein genes and interferon signaling pathway genes were also increased in HGSC when compared with non-HGSC tumors (Supplementary Figs. S3C and S3D). Moreover, we performed IHC staining of ISG15 and BST2 to verify this result at protein levels (Fig. 2G; Supplementary Figs. S3E and S3F).

The regulation of SCNAs and DNA methylation on gene expression

HGSC frequently carried strong amplification on chromosome 8 (Chr8; ref. 4). We wondered whether the 361 DEGs upregulated in primary tumors may be resulted from the amplification of Chr8. We identified that 46 genes of these 361 genes were on Chr8, which is the most significantly enriched chromosome for these DEGs (P = 2.7E−14, hypergeometric test; Fig. 3A). Interestingly, these 46 genes were most enriched in the pathways associated with electron transport chain and OXPHOS (Fig. 3B). We also found that the upregulation of these 46 genes was associated with shorter overall survival time (Fig. 3C). In addition, we also observed that MYC (located on Chr8), which regulated cell division, OXPHOS and aerobic glycolysis, was upregulated in 10 out of the 12 patients with HGSC we analyzed (Supplementary Fig. S4A). The results suggested that the amplification of Chr8 may contribute to tumorigenesis of HGSC through amplifying the master regulator MYC and increasing the dosages of these 46 genes related to OXPHOS.

Figure 3.

The regulation of SCNAs and DNA methylation on RNA expression. A, The enrichment of chromosomes using 361 DEGs upregulated in primary tumors. B and C, GO enrichment (B) and survival analysis (C) using the 46 upregulated DEGs on chromosome 8. D, DNA methylation levels (1-kb tile) of satellites and LINE1. E and F, Top, Venn plot showing the DEGs whose expression levels were correlated with the methylation levels of satellites and LINE1. E, Upregulated DEGs whose expression levels were negatively correlated with the methylation levels of satellites and LINE1. F, Downregulated DEGs whose expression levels were positively correlated with the methylation levels of satellites and LINE1. Bottom, the GO enrichment analysis of the genes that may be potentially regulated by both of satellite and LINE1 methylation. G, Examples of genes derived from E. mRNA, the gene expression of the corresponding gene; WCG, the DNA methylation levels of satellites or LINE1. H, Examples of genes derived from F. I, Survival analysis of gene sets derived from E and F. The mean expression of selected genes from gene sets were used to group patients. J, The GO enrichment analysis of upregulated DEGs (left) and downregulated DEGs (right) whose expression levels were negatively correlated with their promoter methylation. K, Genome browser view showing that a distal region (WCG7649168: chr14: 94,577,079-94,583,033) potentially regulated the expression of IFI27. The left bottom boxplot shows the expression of IFI27. The middle bottom panel shows the zoomed-in views of a known enhancer, which is overlapped with WCG7649168 and the chromatin accessibility of corresponding locations. The right bottom boxplot shows the DNA methylation levels of WCG7649168. PT, primary tumor.

Figure 3.

The regulation of SCNAs and DNA methylation on RNA expression. A, The enrichment of chromosomes using 361 DEGs upregulated in primary tumors. B and C, GO enrichment (B) and survival analysis (C) using the 46 upregulated DEGs on chromosome 8. D, DNA methylation levels (1-kb tile) of satellites and LINE1. E and F, Top, Venn plot showing the DEGs whose expression levels were correlated with the methylation levels of satellites and LINE1. E, Upregulated DEGs whose expression levels were negatively correlated with the methylation levels of satellites and LINE1. F, Downregulated DEGs whose expression levels were positively correlated with the methylation levels of satellites and LINE1. Bottom, the GO enrichment analysis of the genes that may be potentially regulated by both of satellite and LINE1 methylation. G, Examples of genes derived from E. mRNA, the gene expression of the corresponding gene; WCG, the DNA methylation levels of satellites or LINE1. H, Examples of genes derived from F. I, Survival analysis of gene sets derived from E and F. The mean expression of selected genes from gene sets were used to group patients. J, The GO enrichment analysis of upregulated DEGs (left) and downregulated DEGs (right) whose expression levels were negatively correlated with their promoter methylation. K, Genome browser view showing that a distal region (WCG7649168: chr14: 94,577,079-94,583,033) potentially regulated the expression of IFI27. The left bottom boxplot shows the expression of IFI27. The middle bottom panel shows the zoomed-in views of a known enhancer, which is overlapped with WCG7649168 and the chromatin accessibility of corresponding locations. The right bottom boxplot shows the DNA methylation levels of WCG7649168. PT, primary tumor.

Close modal

DNA methylation is one of the most important epigenomic modifications, which undergoes extensive reprogramming during tumorigenesis and may impact gene expression (19, 20). We then analyzed the potentially regulatory effects of repeat element methylation, promoter methylation, and distal region methylation on gene expression.

We calculated the DNA methylation levels of whole genome and different annotated genomic elements (Fig. 3D; Supplementary Figs. S4B and S4C). For all patients, both satellites and long interspersed nuclear element 1 (LINE1) exhibited the most dramatic DNA demethylation (Fig. 3D), even for the cancer cells with global DNA methylation levels comparable with FTE (OC04, OC08, and lineage A1 in OC05; Supplementary Fig. S4B), indicating that strong demethylation of satellites and LINE1 is a hallmark of ovarian cancer. Next, to explore which DEGs were potentially regulated by DNA methylation of these genomic elements, we calculated the correlations between RNA expression and DNA methylation of these genomic elements (Supplementary Table S2). RNA expression of 61 DEGs upregulated in cancer cells were negatively correlated with DNA methylation levels of both satellites and LINE1 (Fig. 3E). Interestingly, the higher expression of genes associated with metallothioneins bind metals and interferon signaling, the two typical upregulated pathways in HGSC, was potentially affected by demethylation of satellites and LINE1 (Fig. 3E and G). In addition, we found that the upregulation of these 61 genes was associated with poorer progression-free survival (PFS) in HGSC (Fig. 3I). As for the DEGs downregulated in primary tumor cells, we revealed that the genes associated with translation elongation may result from demethylation of these two repetitive genomic elements (Fig. 3F and H); and the downregulation these genes predicted shortened PFS (Fig. 3I).

DNA methylation of gene promoters is one of the most important regulations of gene expression. Correlation and GO analysis show that genes associated with metallothioneins bind metals and interferon signaling are not only potentially regulated by satellite and LINE1 methylation, but also by their promoter methylation (Fig. 3J). Interestingly, the metabolic pathway of OXPHOS was significantly enriched but glycolysis was not (Fig. 3J), indicating that the expression of different metabolic pathways was regulated by different mechanisms.

The DNA methylation of distal regulatory regions, such as enhancers, may also be involved in gene expression regulation. Thus, we used MICMIC (21) to identify methylation regulation networks of distal regulatory regions. We totally detected 16,941 regulatory region-target pairs (Supplementary Table S3). Of these pairs, 53.9% were anticorrelated with the expression of target genes. Most target genes were potentially regulated by DNA methylation of less than five regulatory regions. These results demonstrate that some genes are not only potentially regulated by the methylation of their promoters, but also by methylation of distal regulatory regions. For example, the upregulated gene IFI27 (Fig. 3K), whose expression was negatively correlated with its promoter DNA methylation (Supplementary Table S3), was also negatively correlated with methylation of a distal region [∼235 kb away from IFI27 transcriptional start site (TSS); correlation, −0.26; P value, 4.8E−13]. In addition, this distal region is overlapped with a known ovarian enhancer (Fig. 3K) and its flanking regions exhibit higher chromatin accessibility in primary tumor cells compared with FTE cells, highlighting its putative regulatory roles of IFI27 expression.

Chromatin accessibility and the key transcription factors in HGSC

Binding of transcription factors (TF) to the cis-regulatory elements of a gene is a fundamental mechanism of transcription regulation, and TF access is often regulated by chromatin states. We found that the global chromatin accessibility of cancer cells tended to be more open than that of FTE cells (Fig. 4A), whereas the accessibility of regions adjacent to TSSs was less open in cancer cells than that of FTE cells (Fig. 4B). In addition, cancer cells had less and shorter nucleosome-depleted regions (NDR) than FTE cells (Fig. 4C), indicating that the cancer cells had less open chromatin states of focal open chromatin regions. The results of different genomic elements show that the accessibility of CpG islands (CGI), promoters, and microsatellite regions evidently decreased in cancer cells (Fig. 4D).

Figure 4.

Characteristics of chromatin accessibility in HGSC. A, Chromatin accessibility of the whole genome in cancer cells and FTE cells. B, Normalized chromatin accessibility around TSS (±7 kb) in different patients. Each line is colored by different patient, which is the same as the colors in C. C, The number and length of NDRs in FTE and cancer cells, colored by patients. D, The normalized GCH levels of different genomic elements. E, Heatmap representing chromVAR bias-corrected deviations of the TFs whose binding motif presented different accessibility across FTE cells and primary tumor cells. F, Scatterplot showing the TFs whose binding motifs exhibit different deviation scores predicted by chromVAR and may regulate the DEGs inferred by Lisa. Top, TFs were more active in primary tumor cells and its target genes were upregulated in primary tumor cells. Bottom, TFs were more inert in primary tumor cells and its target genes were downregulated in primary tumor cells. The colors of the dots represent the fold changes of the expression in cancer cells and FTE cells.

Figure 4.

Characteristics of chromatin accessibility in HGSC. A, Chromatin accessibility of the whole genome in cancer cells and FTE cells. B, Normalized chromatin accessibility around TSS (±7 kb) in different patients. Each line is colored by different patient, which is the same as the colors in C. C, The number and length of NDRs in FTE and cancer cells, colored by patients. D, The normalized GCH levels of different genomic elements. E, Heatmap representing chromVAR bias-corrected deviations of the TFs whose binding motif presented different accessibility across FTE cells and primary tumor cells. F, Scatterplot showing the TFs whose binding motifs exhibit different deviation scores predicted by chromVAR and may regulate the DEGs inferred by Lisa. Top, TFs were more active in primary tumor cells and its target genes were upregulated in primary tumor cells. Bottom, TFs were more inert in primary tumor cells and its target genes were downregulated in primary tumor cells. The colors of the dots represent the fold changes of the expression in cancer cells and FTE cells.

Close modal

To explore potential key regulators in cancer cells, we performed TF binding motif enrichment of NDRs using chromVAR (22). In total, motifs of 107 TFs were more open in cancer cells, whereas 126 were more closed (Fig. 4E). To further explore which TFs were more likely to participate in DEG expression regulation, we used Lisa (23) to infer the transcriptional regulators of the DEGs (Fig. 4F). Of note, NFE2L2 was the most significantly enriched TF predicted using DEGs and simultaneously its binding motifs were more open in cancer cells than FTE evaluated by NDRs. Previous studies have verified that NFE2L2 can promote ovarian cancer growth (24) and it is associated with resistance to chemotherapy (25, 26). However, the expression level of NFE2L2 has no difference between cancer cells and FTE cells. Of these 107 TFs, expression levels of four TFs in cancer cells are significantly higher than that in FTE cells: FOXK1, TFAP2C, NR2F6, and DDIT3. FOXK1 is an important regulator to induce aerobic glycolysis (27, 28), which was previously proved to facilitate cell proliferation and metastasis in ovarian cancer (27, 29). Considering that the glycolysis pathway was upregulated in HGSC and the binding motifs of its transcriptional regulator FOXK1 exhibited high accessibility, we suggested that FOXK1 could be an attractive candidate therapeutic target for HGSC. For the 126 TFs with more closed chromatin states, several members of the GATA family were enriched. However, GATA1-5 were almost not expressed in cancer cells and FTE cells (Supplementary Fig. S5), whereas only GATA6 was expressed and its expression levels were different between cancer cells and FTE cells. Therefore, we considered that GATA6 played more important roles during tumorigenesis of HGSC among these GATA members. Previous studies had shown that loss of GATA6 expression led to nuclear deformation in ovarian cancer (30, 31). Here, we also illustrated that GATA6 may also regulate tumorigenesis process of HSGC through affecting the chromatin states of its target genes.

The gene expression profiles of abdominal metastases are similar to their matched primary tumors

After describing several key differences between FTE and primary tumors, we next explored the differences between primary tumors and matched metastases. On the basis of genetic lineage analyses, we compared the primary tumor cells and matched metastases within the same genetic lineage in the same patient.

Results show that very few genes were differentially expressed for the patients with abdominal metastases (Fig. 5A) and primary tumor cells exhibited close relationships to its matched metastatic cancer cells (Fig. 5B). These results propose that the microenvironment of abdominal metastases did not strongly affect the global RNA expression patterns of cancer cells, although the cancer cells have invaded and proliferated in a different organ, and imply that the metastatic capacity was probably acquired at an early stage during tumorigenesis. In contrast, previous studies based on bulk RNA profiling have reported several differences between primary tumors and metastases (32–34). Such contradiction may be caused by different technologies used in the researches. The DEGs between primary tumors and matched metastases revealed by bulk profiling may reflect the differences of different genetic lineages of cancer cells or different cell type compositions of the primary tumors and metastases, but not the actual differences of cancer cells from the same genetic lineage in the primary tumors and matched metastases.

Figure 5.

Dynamic changes of RNA expression from primary tumors to metastases. A, UpSet plot visualization of the DEGs between primary tumors and matched metastases. Left, genes upregulated in metastases; right, genes downregulated in metastases. B, Shared nearest neighbor (SNN) clustering of cancer cells of patients with metastases based on RNA expression. C and F, Western blot analysis to verify the overexpression of HSPA6 in SKOV3 cells (C) and ES2 cells (F). D and G, The migration and invasion ability of SKOV3 cells (D) and ES2 cells (G) after HSPA6 overexpression was assessed by transwell assay. Scale bar, 100 μm. E, The migration ability of SKOV3 cells after HSPA6 overexpression was assessed by wound healing assay. Scale bar, 200 μm. ***, P < 0.001; ****, P < 0.0001.

Figure 5.

Dynamic changes of RNA expression from primary tumors to metastases. A, UpSet plot visualization of the DEGs between primary tumors and matched metastases. Left, genes upregulated in metastases; right, genes downregulated in metastases. B, Shared nearest neighbor (SNN) clustering of cancer cells of patients with metastases based on RNA expression. C and F, Western blot analysis to verify the overexpression of HSPA6 in SKOV3 cells (C) and ES2 cells (F). D and G, The migration and invasion ability of SKOV3 cells (D) and ES2 cells (G) after HSPA6 overexpression was assessed by transwell assay. Scale bar, 100 μm. E, The migration ability of SKOV3 cells after HSPA6 overexpression was assessed by wound healing assay. Scale bar, 200 μm. ***, P < 0.001; ****, P < 0.0001.

Close modal

Although only a few DEGs between metastases and primary tumors were shared by different patients (Fig. 5A), these genes may provide some clues to identify critical genes promoting or suppressing metastasis. For example, FN1 exhibited higher expression in metastasis of 2 patients (Fig. 5A). Previous studies have proved that FN1 can promote epithelial–mesenchymal transition and metastasis in serous ovarian cancer (35, 36). We next tested the functions of one of the five genes (HSPA6) downregulated in metastases shared by 2 patients (Fig. 5A). The results showed that overexpression of HSPA6 inhibited migration and invasion of the cancer cells, suggesting that its downregulation promotes tumor metastasis of HGSC (Fig. 5CG).

Cancer cells maintained DNA methylation levels and patterns during abdominal metastasis

We found that global DNA methylation levels of cancer cells are stable during abdominal metastasis in most patients (Fig. 6A). Thus, we speculated that cancer cells maintained the global methylation characteristics of their genetic lineages during metastasis, so that the DNA methylation patterns can be used to describe cancer lineage histories. We then performed unsupervised hierarchical clustering of DNA methylation (see Supplementary Materials and Methods), and found that DNA methylation patterns produce similar lineage histories to those inferred by SCNAs (Fig. 6B; Supplementary Fig. S6). In contrast, different genetic lineages of cancer cells in the same patient could not be distinguished based on the clustering of the transcriptome and chromatin accessibility (Supplementary Fig. S6).

Figure 6.

Dynamic changes of DNA methylation from primary tumors to metastases. A, Single-cell global DNA methylation levels of primary tumor cells and matched metastases within the same lineage. Each point represents a single cell, colored by genetic lineages. The lines at the bottom, middle, and top represent 25 percentile, median, and 75%ile values, respectively. “Pri_tum” and “Met” represent primary tumor and metastasis, respectively. B, Clustering of cancer cells based on DNA methylation for OC09 and OC14, colored by genetic lineages (left) and sampling regions (right). Pri_tum, primary tumor; Met_LN, lymph node metastasis; Met_Ome, omental metastasis.

Figure 6.

Dynamic changes of DNA methylation from primary tumors to metastases. A, Single-cell global DNA methylation levels of primary tumor cells and matched metastases within the same lineage. Each point represents a single cell, colored by genetic lineages. The lines at the bottom, middle, and top represent 25 percentile, median, and 75%ile values, respectively. “Pri_tum” and “Met” represent primary tumor and metastasis, respectively. B, Clustering of cancer cells based on DNA methylation for OC09 and OC14, colored by genetic lineages (left) and sampling regions (right). Pri_tum, primary tumor; Met_LN, lymph node metastasis; Met_Ome, omental metastasis.

Close modal

Intratumor heterogeneities of primary tumors provide insights into the metastasis process

Reconstruction of genetic lineages could unveil ITH and provide precise routes of metastasis. Different lineages not only have different SCNAs, DNA methylation levels, and RNA expression patterns but may also have different metastatic potential. For example, both patient OC09 and OC14 harbor two genetic lineages (lineage A1 and A2) of cancer cells, and the lineage A1 had established secondary tumors, but the other lineage (lineage A2) had not (Fig. 6A). Therefore, we wondered what kind of primary tumor cells possess greater metastatic potential. To explore the common upregulated genes in lineage A1 for patient OC09 and OC14, we performed differential expression analysis and revealed that CCN1 and HSP90AA1 were upregulated in the metastasized lineage A1 for both OC09 and OC14 (Fig. 7A), which prompted us to test if inhibition of CCN1 and HSP90AA1 can suppress the metastasis of ovarian cancer cells. We knocked out CCN1 using the CRISPR-Cas9 system in ovarian cancer cell line SKOV3 and ES2 (Fig. 7B; Supplementary Fig. S7B). Results showed that knockout of CCN1 significantly restrained the migration and invasion of cancer cells by wound healing assay and transwell assays (Fig. 7C and D; Supplementary Figs. S7A and S7C). We inhibited HSP90AA1 using the HSP90α inhibitor TAS-116. We revealed that TAS-116 can dose-dependently inhibit ovarian cancer cell migration and invasion (Fig. 7EG; Supplementary Figs. S7D–S7F). In addition, TAS-116 could significantly suppress the growth of cancer cells in a time- and dose-dependent manner (Fig. 7H; Supplementary Fig. S7G). Together, our results demonstrate that CCN1 and HSP90AA1 may participate in the metastasis of HGSC and may serve as potential therapeutic targets for ovarian cancer.

Figure 7.

Intratumor heterogeneities of primary tumors reveal key genes involved in the metastasis. A, Venn plot showing the number of DEGs between lineage A1 and lineage A2 for OC09 and OC14. Left, upregulated genes in lineage A1; right, upregulated genes in lineage A2. B, Western blot analysis to verify the knockout of CCN1 in SKOV3 cells. C, The migration of SKOV3 cells knocked out of CCN1 was detected by wound healing assay, which was compared with the cells treated with the nontarget gRNA. D, The migration and invasion of SKOV3 cells after CCN1 knockout were assessed by transwell assay. Scale bar, 100 μm. E, The migration of SKOV3 cells treated with the indicated concentrations of TAS-116 was evaluated using wound healing assay. F and G, The migration and invasion of SKOV3 cells treated with the indicated concentrations of TAS-116 were assessed by transwell assay. Scale bar, 100 μm. H, The cell viability of SKOV3 cells treated with the indicated concentrations of TAS-116 for 24, 48, and 72 hours. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Figure 7.

Intratumor heterogeneities of primary tumors reveal key genes involved in the metastasis. A, Venn plot showing the number of DEGs between lineage A1 and lineage A2 for OC09 and OC14. Left, upregulated genes in lineage A1; right, upregulated genes in lineage A2. B, Western blot analysis to verify the knockout of CCN1 in SKOV3 cells. C, The migration of SKOV3 cells knocked out of CCN1 was detected by wound healing assay, which was compared with the cells treated with the nontarget gRNA. D, The migration and invasion of SKOV3 cells after CCN1 knockout were assessed by transwell assay. Scale bar, 100 μm. E, The migration of SKOV3 cells treated with the indicated concentrations of TAS-116 was evaluated using wound healing assay. F and G, The migration and invasion of SKOV3 cells treated with the indicated concentrations of TAS-116 were assessed by transwell assay. Scale bar, 100 μm. H, The cell viability of SKOV3 cells treated with the indicated concentrations of TAS-116 for 24, 48, and 72 hours. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Close modal

To further unveil the shared upregulated pathways of the metastasized lineage between OC09 and OC14, we performed gene set enrichment analysis (GSEA; ref. 37). Using hallmark gene sets, we observed hallmarks of interferon gamma and alpha responses were significantly enriched in the metastasized lineage A1 (Fig. 8AD). Consistently, results based on ontology gene sets show that metastasized lineage A1 exhibited higher expression associated with virus responses (Fig. 8E and F). In addition, lineage A1 also upregulated genes involved in blood vessel morphogenesis (Fig. 8E and F). In OC09, we further showed that the binding motifs of TFs associated with interferon signaling were more open in lineage A1 than A2 (Fig. 8G).

Figure 8.

Intratumor heterogeneities of primary tumors reveal critical pathways involved in the metastasis. A and C, Heatmaps of DEGs between primary tumor cells of lineage A1 and lineage A2 for OC14 (A) and OC09 (C). B and D, GSEA analysis for OC14 (B) and OC09 (D) showing the enriched pathways of lineage A1 using hallmark gene sets. E and F, GSEA analysis for OC14 (E) and OC09 (F) showing the enriched pathways of lineage A1 using ontology gene sets. G, Heat map showing that the binding motifs of TFs involved in interferon responses is more open in lineage A1 than lineage A2 for OC09. H, Enrichment analysis of annotated genomic elements on hypomethylated tiles in lineage A2 compared with lineage A1.

Figure 8.

Intratumor heterogeneities of primary tumors reveal critical pathways involved in the metastasis. A and C, Heatmaps of DEGs between primary tumor cells of lineage A1 and lineage A2 for OC14 (A) and OC09 (C). B and D, GSEA analysis for OC14 (B) and OC09 (D) showing the enriched pathways of lineage A1 using hallmark gene sets. E and F, GSEA analysis for OC14 (E) and OC09 (F) showing the enriched pathways of lineage A1 using ontology gene sets. G, Heat map showing that the binding motifs of TFs involved in interferon responses is more open in lineage A1 than lineage A2 for OC09. H, Enrichment analysis of annotated genomic elements on hypomethylated tiles in lineage A2 compared with lineage A1.

Close modal

With respect to DNA methylation, we identified that, for both OC09 and OC14, the lineage with lower residual DNA methylation levels (lineage A2) did not metastasized. To further investigate which genomic regions were demethylated in nonmetastasized lineage A2, we performed enrichment analysis of the differentially methylated tiles between lineage A1 and A2. Results showed that the tiles with decreased DNA methylation in lineage A2 were strongly enriched in intergenic regions, satellites, and LTR (Fig. 8H).

Besides OC09 and OC14, primary tumors of four other patients possessed at least two genetic lineages, all of which metastasized or did not. Although these lineages present comparable metastatic potential, differential gene expression analysis shows that the phenotypes of these genetic lineages were distinct (Supplementary Figs. S8A–S8H).

Here, we introduce an improved single-cell multiomics sequencing method, and provide the first systematic exploration about inter- and intratumor heterogeneities of SCNAs, DNA methylome, chromatin accessibility, and the transcriptome of primary tumors and matched metastases in HGSC at single-cell and single-base resolution.

Our single-cell multiomics sequencing technology exhibits several advantages in studying tumor biology. First, cell identities can be more accurately determined by multiomics data. Most studies distinguished malignant cells from nonmalignant cells using scRNA-seq data, mainly based on the gene expression of cell type-specific markers and RNA-inferred SCNAs (38, 39). In contrast, using the single-cell multiomics sequencing, we can not only determine the cell identities via RNA data but also profile higher resolution SCNA patterns by DNA sequencing data. Second, we can explore the dynamic changes of DNA methylome, chromatin states, and transcriptome during tumor development after reconstructing the genetic trajectories of cancer cells. Similarly, using the multiomics technology, Clark and colleagues utilized RNA data to reconstruct developmental trajectories of mouse embryonic stem cells and explored the epigenetic and transcriptomic properties of different cell states along these trajectories (40). Third, we can uncover potential inter-omic regulatory relationships, which are very difficult to identify by bulk sequencing or separate single-cell analyses of different molecular layers.

On the basis of genetic lineage analysis, we explored the dynamics of the epigenome and transcriptome during tumor promotion. We noticed that the metabolic pathways were highly reprogrammed in HGSC (Fig. 2B). A century ago, Warburg and colleagues proposed that tumors mainly metabolized glucose through glycolysis rather than OXPHOS (41). However, we showed that both OXPHOS and glycolysis were increased in ovarian cancer cells. Consistently, recent studies based on cell lines and mouse models have reported that ovarian cancer cells displayed higher rates of OXPHOS and glycolysis than normal epithelial cells (42, 43). We further revealed that these two pathways may be regulated by different mechanisms. The upregulation of genes related to OXPHOS may be due to either the Chr8 amplification to increase their dosages, or upregulate their upstream regulators (e.g., MYC), or their promoter demethylation; while the binding motif of the key TF involved in glycolysis, FOXK1, exhibited high accessibility in HGSC and may upregulate the glycolysis-related genes.

Taking advantage of single-cell multiomics sequencing technology, we revealed the high ITH in HGSC. We show that most tumors were composed of diverse lineages of cancer cells. These lineages vary in DNA methylation levels, RNA expression levels, and metastatic potential. The results of patient OC09 and OC14 show that the lineage with a higher residual DNA methylation level escaped from the primary tumor and invaded other tissues, whereas the other lineage, with a lower residual DNA methylation level, did not (Fig. 6A). Furthermore, the metastasized and probably more aggressive lineage exhibited upregulation of CCN1, HSP90AA1, and genes associated with virus response and angiogenesis (such as VEGFA). CCN1 has been proved to participate in cell proliferation, inflammation, and angiogenesis, and its overexpression was associated with poor prognosis of ovarian cancer (44, 45). Previous studies have shown that inhibition of HSP90AA1 may potentate the ability of chemotherapy (46–48). Here we revealed that CCN1 and HSP90AA1 may also be involved in ovarian cancer metastasis. While performing end point analysis to compare the primary tumors and metastases is a common means to investigate the mechanisms of metastasis, this approach cannot provide insights into the early events leading to metastasis. In our study, we used single-cell multiomics sequencing to compare primary tumor cells with different metastatic potential, which provided a new strategy for the study of the mechanisms of metastasis and help to better understand the metastasis process.

Next, we inspected the epigenetic and transcriptomic changes associated with tumor progression (Figs. 5 and 6). The results show that the global DNA methylation levels were maintained during metastasis. Unsupervised hierarchical clustering of DNA methylation shows that the DNA methylation patterns produce similar lineage histories to those deduced from SCNA patterns. As a covalent modification of genome, DNA methylation is relatively stable and involved in establishment and maintenance of cell type- and tissue-specific features during mitosis and development (49–51). Thus, despite the changes in DNA methylation at a relatively early stage of tumorigenesis, ovarian cancer cells retained the majority of their lineage-specific DNA methylation patterns during later tumor progression. The convergent pattern of DNA methylation and SCNAs was also demonstrated in prostate cancer (52). Another study has reported that phyloepigenetic patterns could recapitulate the phylogenetic histories reconstructed by SNVs (53). Therefore, the convergence of SCNAs, SNVs, and DNA methylome changes may exist in a variety of cancers.

In summary, using our single-cell multiomics sequencing, we obtained high-resolution profiles of HGSC, dissected its extensive inter- and intratumor heterogeneities, identified many potentially important variations during tumorigenesis, and explored the regulatory networks in HGSC. Our study provides highly valuable resources for understanding the molecular characteristics of HGSC, as well as for improving the diagnosis and individualized therapy of this disease.

K. Kee reports grants from Tsinghua University during the conduct of the study. No disclosures were reported by the other authors.

Y. Wang: Conceptualization, resources, formal analysis, investigation, visualization, methodology, writing–original draft. H. Xie: Resources, data curation, software, formal analysis, visualization, methodology, writing–review and editing. X. Chang: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–review and editing. W. Hu: Formal analysis, investigation, writing–review and editing. M. Li: Formal analysis, investigation, methodology, writing–review and editing. Y. Li: Resources, investigation, methodology, writing–review and editing. H. Liu: Validation, visualization, methodology. H. Cheng: Investigation, methodology. S. Wang: Methodology. L. Zhou: Methodology. D. Shen: Methodology. S. Dou: Methodology. R. Ma: Methodology. Y. Mao: Methodology. H. Zhu: Methodology. X. Zhang: Methodology. Y. Zheng: Methodology. X. Ye: Methodology. L. Wen: Conceptualization, supervision. K. Kee: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–review and editing. H. Cui: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–review and editing. F. Tang: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, methodology, project administration, writing–review and editing.

The authors thank the Beijing Advanced Innovation Centre for Genomics, the Computing Platform of the Center for Life Science, and the National Center for Protein Sciences at Peking University for supporting our work to F. Tang. K. Kee and H. Cui were supported by the National Key Research and Development Program of China 2018YFA0107703 and 2016YFA0201400, respectively. X. Chang was supported by the National Natural Science Foundation of China (81971360 and 81671431).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Torre
LA
,
Trabert
B
,
DeSantis
CE
,
Miller
KD
,
Samimi
G
,
Runowicz
CD
, et al
.
Ovarian cancer statistics, 2018
.
CA Cancer J Clin
2018
;
68
:
284
96
.
2.
Kim
S
,
Han
Y
,
Kim
SI
,
Kim
H-S
,
Kim
SJ
,
Song
YS
.
Tumor evolution and chemoresistance in ovarian cancer
.
NPJ Precis Oncol
2018
;
2
:
20
.
3.
Bashashati
A
,
Ha
G
,
Tone
A
,
Ding
J
,
Prentice
LM
,
Roth
A
, et al
.
Distinct evolutionary trajectories of primary high-grade serous ovarian cancers revealed through spatial mutational profiling
.
J Pathol
2013
;
231
:
21
34
.
4.
Network CGAR
.
Integrated genomic analyses of ovarian carcinoma
.
Nature
2011
;
474
:
609
.
5.
Verhaak
RG
,
Tamayo
P
,
Yang
J-Y
,
Hubbard
D
,
Zhang
H
,
Creighton
CJ
, et al
.
Prognostically relevant gene signatures of high-grade serous ovarian carcinoma
.
J Clin Invest
2013
;
123
:
517
25
.
6.
Hao
Q
,
Li
J
,
Zhang
Q
,
Xu
F
,
Xie
B
,
Lu
H
, et al
.
Single-cell transcriptomes reveal heterogeneity of high-grade serous ovarian carcinoma
.
Clin Transl Med
2021
;
11
:
e500
.
7.
Izar
B
,
Tirosh
I
,
Stover
EH
,
Wakiro
I
,
Cuoco
MS
,
Alter
I
, et al
.
A single-cell landscape of high-grade serous ovarian cancer
.
Nat Med
2020
;
26
:
1271
9
.
8.
Shih
AJ
,
Menzin
A
,
Whyte
J
,
Lovecchio
J
,
Liew
A
,
Khalili
H
, et al
.
Identification of grade and origin specific cell populations in serous epithelial ovarian cancer by single cell RNA-seq
.
PLoS One
2018
;
13
:
e0206785
.
9.
Winterhoff
BJ
,
Maile
M
,
Mitra
AK
,
Sebe
A
,
Bazzaro
M
,
Geller
MA
, et al
.
Single cell sequencing reveals heterogeneity within ovarian cancer epithelium and cancer associated stromal cells
.
Gynecol Oncol
2017
;
144
:
598
606
.
10.
Guo
F
,
Li
L
,
Li
J
,
Wu
X
,
Hu
B
,
Zhu
P
, et al
.
Single-cell multi-omics sequencing of mouse early embryos and embryonic stem cells
.
Cell Res
2017
;
27
:
967
88
.
11.
Bian
S
,
Hou
Y
,
Zhou
X
,
Li
X
,
Yong
J
,
Wang
Y
, et al
.
Single-cell multiomics sequencing and analyses of human colorectal cancer
.
Science
2018
;
362
:
1060
3
.
12.
Li
L
,
Dong
J
,
Yan
L
,
Yong
J
,
Liu
X
,
Hu
Y
, et al
.
Single-cell RNA-seq analysis maps development of human germline cells and gonadal niche interactions
.
Cell Stem Cell
2017
;
20
:
858
73
.
13.
Fan
X
,
Lu
P
,
Wang
H
,
Bian
S
,
Wu
X
,
Zhang
Y
, et al
.
Integrated single-cell multiomics analysis reveals novel candidate markers for prognosis in human pancreatic ductal adenocarcinoma
.
Cell Discov
2022
;
8
:
13
.
14.
Chen
Y
,
Li
H
,
Cheng
H-y
Rui-Qiong
M
,
Ye
X
,
Cui
H
, et al
.
Fibrinogen alpha chain is up-regulated and affects the pathogenesis of endometriosis
.
Reprod Biomed Online
2019
;
39
:
893
904
.
15.
Tickle
T
,
Tirosh
I
,
Georgescu
C
,
Brown
M
,
Haas
B
.
inferCNV of the Trinity CTAT Project
.
Cambridge, MA
:
Broad Institute of MIT, Harvard
;
2019
.
16.
Zack
TI
,
Schumacher
SE
,
Carter
SL
,
Cherniack
AD
,
Saksena
G
,
Tabak
B
, et al
.
Pan-cancer patterns of somatic copy number alteration
.
Nat Genet
2013
;
45
:
1134
40
.
17.
Beroukhim
R
,
Mermel
CH
,
Porter
D
,
Wei
G
,
Raychaudhuri
S
,
Donovan
J
, et al
.
The landscape of somatic copy-number alteration across human cancers
.
Nature
2010
;
463
:
899
905
.
18.
Labidi-Galy
SI
,
Papp
E
,
Hallberg
D
,
Niknafs
N
,
Adleff
V
,
Noe
M
, et al
.
High grade serous ovarian carcinomas originate in the fallopian tube
.
Nat Commun
2017
;
8
:
1093
.
19.
Baylin
SB
,
Esteller
M
,
Rountree
MR
,
Bachman
KE
,
Schuebel
K
,
Herman
JG
.
Aberrant patterns of DNA methylation, chromatin formation and gene expression in cancer
.
Hum Mol Genet
2001
;
10
:
687
92
.
20.
Sharma
S
,
Kelly
TK
,
Jones
PA
.
Epigenetics in cancer
.
Carcinogenesis
2010
;
31
:
27
36
.
21.
Tong
Y
,
Sun
J
,
Wong
CF
,
Kang
Q
,
Ru
B
,
Wong
CN
, et al
.
MICMIC: identification of DNA methylation of distal regulatory regions with causal effects on tumorigenesis
.
Genome Biol
2018
;
19
:
73
.
22.
Schep
AN
,
Wu
B
,
Buenrostro
JD
,
Greenleaf
WJ
.
chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data
.
Nat Methods
2017
;
14
:
975
8
.
23.
Qin
Q
,
Fan
J
,
Zheng
R
,
Wan
C
,
Mei
S
,
Wu
Q
, et al
.
Lisa: inferring transcriptional regulators through integrative modeling of public chromatin accessibility and ChIP-seq data
.
Genome Biol
2020
;
21
:
32
.
24.
Manandhar
S
,
Choi
B-h
,
Jung
K-A
,
Ryoo
I-G
,
Song
M
,
Kang
SJ
, et al
.
NRF2 inhibition represses ErbB2 signaling in ovarian carcinoma cells: implications for tumor growth retardation and docetaxel sensitivity
.
Free Radical Biol Med
2012
;
52
:
1773
85
.
25.
Liu
P
,
Dodson
M
,
Fang
D
,
Chapman
E
,
Zhang
DD
.
NRF2 negatively regulates primary ciliogenesis and hedgehog signaling
.
PLoS Biol
2020
;
18
:
e3000620
.
26.
Xia
M
,
Yu
H
,
Gu
S
,
Xu
Y
,
Su
J
,
Li
H
, et al
.
p62/SQSTM1 is involved in cisplatin resistance in human ovarian cancer cells via the Keap1-Nrf2-ARE system
.
Int J Oncol
2014
;
45
:
2341
8
.
27.
Sun
H
,
Wang
H
,
Wang
X
,
Aoki
Y
,
Wang
X
,
Yang
Y
, et al
.
Aurora-A/SOX8/FOXK1 signaling axis promotes chemoresistance via suppression of cell senescence and induction of glucose metabolism in ovarian cancer organoids and cells
.
Theranostics
2020
;
10
:
6928
.
28.
Sukonina
V
,
Ma
H
,
Zhang
W
,
Bartesaghi
S
,
Subhash
S
,
Heglind
M
, et al
.
FOXK1 and FOXK2 regulate aerobic glycolysis
.
Nature
2019
;
566
:
279
83
.
29.
Li
L
,
Gong
M
,
Zhao
Y
,
Zhao
X
,
Li
Q
.
FOXK1 facilitates cell proliferation through regulating the expression of p21, and promotes metastasis in ovarian cancer
.
Oncotarget
2017
;
8
:
70441
.
30.
Capo-chichi
CD
,
Cai
KQ
,
Testa
JR
,
Godwin
AK
,
Xu
X-X
.
Loss of GATA6 leads to nuclear deformation and aneuploidy in ovarian cancer
.
Mol Cell Biol
2009
;
29
:
4766
77
.
31.
Cai
KQ
,
Caslini
C
,
Capo-Chichi
CD
,
Slater
C
,
Smith
ER
,
Wu
H
, et al
.
Loss of GATA4 and GATA6 expression specifies ovarian cancer histological subtypes and precedes neoplastic transformation of ovarian surface epithelia
.
PLoS One
2009
;
4
:
e6454
.
32.
Sallinen
H
,
Janhonen
S
,
Plnen
P
,
Niskanen
H
,
Kaikkonen
MU
.
Comparative transcriptome analysis of matched primary and distant metastatic ovarian carcinoma
.
BMC Cancer
2019
;
19
:
1121
.
33.
Mitra
S
,
Tiwari
K
,
Podicheti
R
,
Pandhiri
T
,
Mitra
AK
.
Transcriptome profiling reveals matrisome alteration as a key feature of ovarian cancer progression
.
Cancers
2019
;
11
:
1513
.
34.
Brodsky
AS
,
Andrew
F
,
Miller
DH
,
Souriya
V
,
Shannon
ML
,
Hsin-Ta
W
, et al
.
Expression profiling of primary and metastatic ovarian tumors reveals differences indicative of aggressive disease
.
PLoS One
2014
;
9
:
e94476
.
35.
Liang
H
,
Yu
M
,
Yang
R
,
Zhang
L
,
Zhang
L
,
Zhu
D
, et al
.
A PTAL-miR-101-FN1 axis promotes EMT and invasion-metastasis in serous ovarian cancer
.
Mol Ther Oncolytics
2020
;
16
:
53
62
.
36.
Lou
X
,
Han
X
,
Jin
C
,
Tian
W
,
Yu
W
,
Ding
D
, et al
.
SOX2 targets fibronectin 1 to promote cell migration and invasion in ovarian cancer: new molecular leads for therapeutic intervention
.
OMICS
2013
;
17
:
510
8
.
37.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
.
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
Proc Natl Acad Sci U S A
2005
;
102
:
15545
50
.
38.
Puram
SV
,
Tirosh
I
,
Parikh
AS
,
Patel
AP
,
Yizhak
K
,
Gillespie
S
, et al
.
Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer
.
Cell
2017
;
171
:
1611
24
.
39.
Neftel
C
,
Laffy
J
,
Filbin
MG
,
Hara
T
,
Shore
ME
,
Rahme
GJ
, et al
.
An integrative model of cellular states, plasticity, and genetics for glioblastoma
.
Cell
2019
;
178
:
835
49
.
40.
Clark
SJ
,
Argelaguet
R
,
Kapourani
C-A
,
Stubbs
TM
,
Lee
HJ
,
Alda-Catalinas
C
, et al
.
scNMT-seq enables joint profiling of chromatin accessibility DNA methylation and transcription in single cells
.
Nat Commun
2018
;
9
:
1
9
.
41.
Warburg
O
,
Wind
F
,
Negelein
E
.
The metabolism of tumors in the body
.
J Gen Physiol
1927
;
8
:
519
30
.
42.
Anderson
AS
,
Roberts
PC
,
Frisard
MI
,
McMillan
RP
,
Brown
TJ
,
Lawless
MH
, et al
.
Metabolic changes during ovarian cancer progression as targets for sphingosine treatment
.
Exp Cell Res
2013
;
319
:
1431
42
.
43.
Dier
U
,
Shin
D-H
,
Hemachandra
LMP
,
Uusitalo
LM
,
Hempel
N
.
Bioenergetic analysis of ovarian cancer cell lines: profiling of histological subtypes and identification of a mitochondria-defective cell line
.
PLoS One
2014
;
9
:
e98479
.
44.
Shen
H
,
Cai
M
,
Zhao
S
,
Wang
H
,
Li
M
,
Yao
S
, et al
.
CYR61 overexpression associated with the development and poor prognosis of ovarian carcinoma
.
Med Oncol
2014
;
31
:
117
.
45.
Gery
S
,
Xie
D
,
Yin
D
,
Gabra
H
,
Miller
C
,
Wang
H
, et al
.
Ovarian carcinomas: CCN genes are aberrantly expressed and CCN1 promotes proliferation of these cells
.
Clin Cancer Res
2005
;
11
:
7243
54
.
46.
Banerji
U
,
Sain
N
,
Sharp
SY
,
Valenti
M
,
Asad
Y
,
Ruddle
R
, et al
.
An in vitro and in vivo study of the combination of the heat shock protein inhibitor 17-allylamino-17-demethoxygeldanamycin and carboplatin in human ovarian cancer models
.
Cancer Chemother Pharmacol
2008
;
62
:
769
78
.
47.
Maloney
A
,
Clarke
PA
,
Naaby-Hansen
S
,
Stein
R
,
Koopman
J-O
,
Akpan
A
, et al
.
Gene and protein expression profiling of human ovarian cancer cells treated with the heat shock protein 90 inhibitor 17-allylamino-17-demethoxygeldanamycin
.
Cancer Res
2007
;
67
:
3239
53
.
48.
Sain
N
,
Krishnan
B
,
Ormerod
MG
,
De Rienzo
A
,
Liu
WM
,
Kaye
SB
, et al
.
Potentiation of paclitaxel activity by the HSP90 inhibitor 17-allylamino-17-demethoxygeldanamycin in human ovarian carcinoma cell lines with high levels of activated AKT
.
Mol Cancer Ther
2006
;
5
:
1197
208
.
49.
Lokk
K
,
Modhukur
V
,
Rajashekar
B
,
Märtens
K
,
Mägi
R
,
Kolde
R
, et al
.
DNA methylome profiling of human tissues identifies global and tissue-specific methylation patterns
.
Genome Biol
2014
;
15
:
3248
.
50.
Schultz
MD
,
He
Y
,
Whitaker
JW
,
Hariharan
M
,
Mukamel
EA
,
Leung
D
, et al
.
Human body epigenome maps reveal noncanonical DNA methylation variation
.
Nature
2015
;
523
:
212
6
.
51.
Wan
J
,
Oliver
VF
,
Wang
G
,
Zhu
H
,
Zack
DJ
,
Merbs
SL
, et al
.
Characterization of tissue-specific differential DNA methylation suggests distinct modes of positive and negative gene expression regulation
.
BMC Genomics
2015
;
16
:
49
.
52.
Brocks
D
,
Assenov
Y
,
Minner
S
,
Bogatyrova
O
,
Simon
R
,
Koop
C
, et al
.
Intratumor DNA methylation heterogeneity reflects clonal evolution in aggressive prostate cancer
.
Cell Rep
2014
;
8
:
798
806
.
53.
Mazor
T
,
Pankov
A
,
Johnson
BE
,
Hong
C
,
Hamilton
EG
,
Bell
RJ
, et al
.
DNA methylation and somatic mutations converge on the cell cycle and define similar evolutionary histories in brain tumors
.
Cancer Cell
2015
;
28
:
307
17
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.

Supplementary data