The Identification of Alternative Polyadenylation in Stomach Adenocarcinomas Using the Genotype-Tissue Expression Project and the Cancer Genome Atlas– Stomach Adenocarcinoma Profiles

Introduction

As the primary subtype of stomach cancer, stomach adenocarcinoma (STAD) is the fifth most frequently diagnosed malignant tumor and the third leading cause of cancer-related deaths in the world.¹ Approximately 90% ~ 95% of all stomach cancer cases are STADs.² The incidence of STAD in China is the highest in the world, accounting for 49.9% of global cases.³ Although the survival rate of patients with STADs has greatly improved in the past 20 years due to the development of treatment methods such as surgery, chemotherapy, and targeted therapy, its prognosis remains unsatisfactory.⁴ The five-year survival rate of patients with a STAD in China is 30.2% ~ 35.9%; in European countries, the five-year survival rate is only 10% ~ 30%.⁵

Genome-wide gene expression profiles have been established including The Cancer Genome Atlas (TCGA) to better understand the impact of clinical outcomes on genetic changes in tumors. Extensive sequencing data analyses have identified mutant genes (eg, TP53, ARID1A, FAT4, KMT2D, and PIK3CA) and dysregulated functional pathways (such as neuroactive ligand–receptor interaction, pancreatic secretion, and the calcium signaling pathway) related to STAD cases. However, the mechanism of gene expression dysregulation in patients with a STAD is still not fully understood.

Alternative polyadenylation (APA) is one of the key regulatory mechanisms of gene expression in tumors. As a post-transcriptional processing mechanism, APA can produce a different three prime untranslated region (3ʹ-UTR) on messenger ribonucleic acid (mRNA) transcripts via RNA polymerase II.^6,7 The variability of the 3ʹ-UTR can change the micro-ribonucleic acid (miRNA) or RNA-binding protein regulatory sites of mRNA transcripts, thereby regulating the stability, function, expression, or subcellular localization of mRNA.^8,9 In humans, 51% ~ 79% of genes express alternative 3ʹ-UTR, indicating that the transcripts corresponding to most genes have multiple polyadenylation sites (PAS).¹⁰ Alternative polyadenylation participates in many cells’ physiological processes, such as cell replication, proliferation, migration and invasion, chromatin signaling, pluripotent cell fate, and tissue aging.^11–15 In addition, APA disorder is considered a driver of tumorigenesis.^16,17 A large number of APA dysregulation events have been reported in many types of cancer, and 3ʹ-UTR-shortened events (61% ~ 98% of total APA events) are the principal consideration in this regard.^18,19 These 3ʹ-UTR-shortened APA events regulate the expression of cancer-related genes by losing miRNA regulatory sites.^15,20 The shortened APA events also regulate the expression of tumor suppressors by reducing competitive endogenous RNA interference.^17,21 Due to a lack of control samples, existing studies on APA events in cases of STAD have primarily focused on pan-cancer APA analysis; the single-cancer APA study of STADs still requires further characterization.

We analyzed 289 STAD tumor data from the STAD study in TCGA (TCGA-STAD)²² to identify the distribution of APA in STAD. Since there was less data on normal tissue compared with cancer in TCGA related to STAD, we used data for 207 normal stomach tissue cases in the Genotype-Tissue Expression (GTEx) project as controls.²³ We found 271 genes with highly recurrent dynamic APA events associated with STAD, demonstrating the additional prognostic capabilities of APA beyond common demographic and clinical variables, thereby expanding our knowledge of the mechanisms and consequences of APA regulation during STAD tumorigenesis.

Methods

Data Collection

The RNA-Seq data for 289 patients with STAD were downloaded from the TCGA-STAD database (https://portal.gdc.cancer.gov). The RNA-Seq data for 207 normal stomach tissue cases were downloaded from the GTEx database (dbGaP Study Accession: phs000424.v7.p2. https://gtexportal.org). These RNA-Seq data were both paired-end, non-strand-specific sequencing data. The data titles and the depth and length of the reads are listed in Table S1. Clinical data, such as age, TNM stage, gender, survival time, and status were also downloaded from the TCGA database.

Bioinformatics Analysis and Statistical Methods

After Sequence Read Archive data from GTEx was converted to FASTQ data, the same Pipeline process as in TCGA (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/) was used to generate BAM files to reduce the difference between the two cohorts.

Dynamic Analysis of Alternative Polyadenylation from RNA Sequencing Analysis

The aligned BAM files (converted to bedGraph files) were generated using BEDTools (v.2.29.0). The bedGraph files were provided as input to the dynamic analysis of alternative polyadenylation from RNA sequencing (DaPars) algorithm.¹⁸ The DaPars process (developed by Xia et al) is a regression-based algorithm that performs the de novo identification of dynamic APAs from standard RNA sequencing data.¹⁸ In this study, DaPars (Python 2.7.17) was used to process the bedGraph coverage files. The dynamic poly (A) sites in the RNA-Seq data were used to de novo identify and analyze differences in APA between the 3ʹ-UTR of tumor and normal samples. The output of DaPars is the PDUI score matrix for any given gene (row) in each sample (column).

The Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis

For the specific target gene, the R package clusterProfiler²⁴ was used for the Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, as well as dot plot drawing using R-3.6.1.

The Micro-Ribonucleic Acid Binding Site

The highly conserved miRNA binding sites and their genomic location was downloaded from TargetScanHuman 7.2 (http://www.targetscan.org). Combining this list and DaPars’ prediction of the genetic coordinates of the lost 3ʹUTR, the missing highly conserved miRNA binding sites were counted.

Heterogeneity Analysis

We calculated the variance of each unique gene’s PDUI score in the tumor samples (Var [Tumor]) and normal samples (Var [Normal]) to measure their heterogeneity in the use of PAS. We plotted the variance difference between the two datasets (Var [Normal]–Var [Tumor]) using R-3.6.1. Homogeneity of variance testing was used to assess the heterogeneity of each gene.

Survival Analysis

We selected 21 genes with observably altered 3ʹ-UTRs and univariate prognostic value based on the P-value order and combined them with clinical covariates (gender, race, age, tumor stage, and surgical outcome) to conduct the study’s survival analysis. First, according to the PDUI score corresponding to the 21 genes, we selected the k-means grouping scheme with K = 3 to group patients with STADs. We used the standard Cox proportional hazard model implemented by the “survival” package in R to perform survival analysis, as well as Kaplan–Meier plots for each group of patients.²⁵ We used the Log rank test to calculate the P-value.

Results

The mRNA-Seq data corresponding to 289 tumor samples was randomly selected from the TCGA-STAD study to analyze differences in the APA spectrum for comparing tumor and normal samples. Concurrently, considering the lack of data sources for normal comparable tissue in TCGA-STAD, 207 cases of normal gastric tissue data from the GTEx project were used as controls. The corresponding sample list is given in Table S1.

Since the data obtained from the TCGA-STAD project had been aligned to the hg38 genome, the data obtained from the GTEx project were processed using the TCGA data processing scheme, and the aligned BAM file was converted into a genomic coverage file for identifying 3ʹ-UTR differences. Finally, the dynamic analysis of an alternative DaPars algorithm using RNA-Seq data was employed to perform APA analysis on the coverage data to generate a percentage of distal PDUI scores for all of the given genes in each sample.¹⁸ The PDUI scores were obtained by quantifying the lengthened and shortened 3ʹ-UTR coverage levels. The gene corresponding to the lengthened 3ʹ-UTR (distal PAS) corresponded to a PDUI score close to 1, and the gene corresponding to the shortened 3ʹ-UTR (proximal PAS) corresponded to a PDUI score close to 0. The change in APA between the tumor and normal samples was characterized by ΔPDUI, which indicated the difference between the average PDUI scores of tumors and normal samples (ΔPDUI = MeanPDUI Tumor/MeanPDUI Normal). The ΔPDUI value was also used as a measure for correlating between the 3ʹ-UTR lengthening and shortening of the tumor samples. The final output was a PDUI matrix that included 1638 non-redundant unique genes, and each column included a tumor/normal sample (289 tumors + 207 normal = 496 columns; see Table S2).

Integrative Analysis of the Cancer Genome Atlas and the Cancer Genome Atlas–Stomach Adenocarcinoma Ribonucleic Acid Sequencing Data to Identify Dynamic Polyadenylation Events in Stomach Adenocarcinomas

We used a strict ΔPDUI threshold of = ±0.2 to further screen for significant APA events between the tumor and the normal samples (Figure 1A and B) to determine the extent of APA-mediated shortened and lengthened 3′-UTR in STADs. Although the ΔPDUI that corresponded to most genes did not exceed the threshold, the number of significantly shortened APA events was higher (the red dots, N = 149) compared with the significantly lengthened APA events (the blue dots, N = 122); all of the adjusted P-values corresponding to the significant APA event were less than 0.05. According to statistics, 55% of significant APA events were shortened 3ʹ-UTR events, which was consistent with the patterns observed in multiple pan-cancer and single-cancer analyses.^17,18,26 The alternative lengths of STAD-related APA events were mainly distributed in 200–300 bp (Figure 1C). In addition, the number of shortened APA events with an alternative length greater than 500 bp was significantly larger compared with lengthened APA events (Figure 1C). We found that two genes in the Rho GTPases family of genes, CDC42 and RAC1, were significantly shortened APA events. The Rho GTPase family of enzymes is essential for multiple biological processes (including cell morphology, polarity, and migration) in many cancer types.^27–31 In particular, the Rho GTPases family is also involved in the mechanism of tumorigenesis caused by the infection of gastric epithelial cells by Helicobacter pylori³² bacteria, microorganisms that are also the primary risk factor for stomach cancer.³³ We visualized a 3ʹ-UTR coverage map of CDC42 and RAC1 between STAD and normal samples to confirm the shortening of the 3ʹ-UTR (Figure 1D). However, APA was not found in the most frequently mutated STAD genes (ie, TP53, ARID1A, FAT4, KMT2D, and PIK3CA; see Figure S1).

Figure 1 Comprehensive analysis of RNA-seq data identified APA related to STAD. (A) PDUI score graph of each gene in human tumor and normal samples. The dotted line indicates the 0.2 cutoff. The blue dots represent the 3ʹ-UTR lengthened genes, and the red dots represent the 3ʹ-UTR shortened genes. (B) Volcano graphs showing 3ʹ-UTR shortened (red) and lengthened (blue) genes, with |ΔPDUI| ≥ 0.2 (vertical dotted line). The horizontal dotted line indicates that adjusted p-value is less than 0.05. (C) Histogram shows the number of base pairs lost/gained by 3ʹ-UTR altered genes. (D) The IGV genome browser diagram depicts the 3ʹ-UTR RNA-seq density map of two 3ʹ-UTR shortened genes (CDC42 and RAC1) to highlight the difference in coverage between tumor patients (red) and normal samples (blue). (E) KEGG pathway with significant enrichment (p <0.01) associated with 3ʹ-UTR altered genes.

The KEGG pathway analysis of 271 significant 3ʹ-UTR-altered genes showed that these genes were mainly enriched in mRNA splicing and processing, as well as in oxidative phosphorylation pathways. Pan-cancer APA analysis identified similar pathways that were consistent with the presence of recurrent APA events in multiple cancer types.^17,18 However, we observed that STAD-related pathways were also enriched, including non-alcoholic fatty liver disease, and epithelial cell signaling in Helicobacter pylori and salmonella infections (Figure 1E). Therefore, the alternative 3ʹ-UTR may regulate the activity of the STAD development pathway.

Alternative Polyadenylation Drives Gene Expression Disorders in Stomach Adenocarcinomas

The shortening of the 3ʹ-UTR during tumorigenesis can escape miRNA inhibition and up-regulate the expression of its parental gene.¹⁸ We calculated the number of highly conserved miRNA binding sites in the 3ʹ-UTR-shortened region (Figure 2A) to assess a pattern for the APA-mediated loss of miRNA binding sites in patients with STADs. Accordingly, we determined that 50% of the 3′-UTR-shortened genes lost at least one highly conserved miRNA binding site, suggesting that changes in the miRNA binding site pool were common patterns of APA-mediated regulation. In addition, by deriving the number of lost highly conserved miRNA binding sites, the miRNA family with the most missing miRNA binding sites was found to be miR-133a-3p, miR-133b, and miR-30-5p. (Figure 2B); miR-133a-3p is a tumor suppressor in many cancers and can regulate cell proliferation and differentiation.^34–37 It has been reported that the overexpression of miR-133a-3p can prevent the activation of autophagy, destroy abnormal glutamine decomposition, and further inhibit the growth and metastasis of gastric cancer cells.³⁸ In addition, miR-133 is reported to induce a negative adjustment of the CDC42-PAK pathway in gastric cancer.³⁹ The analysis results of miRNA indicated that APA influenced the progress of STAD by changing the binding site-library of 3ʹ-UTR. To examine the consequences of the loss of miRNA binding sites in the 3ʹ-UTR, we calculated the differential gene expression between tumor and normal tissue types. The results indicated that the expression levels of genes with a shortened 3ʹ-UTR in tumors tended to be overexpressed (Figure 2C and D; P < 0.01, chi-square test).

Figure 2 APA drives changes in gene expression. (A) The number of genes that lost highly conserved miRNA binding sites due to shortened 3ʹ-UTR. The percentage of genes missing at least one miRNA binding site is shown above the brackets. (B) The dot plot shows the number of miRNA binding sites of the highly conserved miRNA family lost in the 3ʹ-UTR altered genes. (C) For the 3ʹ-UTR altered gene, the dot plot shows the log-fold change (threshold=1) in gene expression relative to ΔPDUI. The overexpressed genes (red dots) and underexpressed genes (blue dots) on the left represent 3ʹ-UTR shortened genes, while those on the right represent 3ʹ-UTR lengthened genes. (D) The bar graph shows the quantification of 3ʹ-UTR altered genes that are overexpressed (red) or underexpressed (blue) in STAD tumors.

The Heterogeneity of Proximal Polyadenylation Sites’ Use in Patients with Stomach Adenocarcinomas

The genetic and epigenetic mutations in tumors can promote the migration and invasion of cancer cells and, subsequently, stimulate the development of the tumor.⁴⁰ Patients with pancreatic ductal adenocarcinoma showed significant heterogeneity in the use of proximal PAS.¹⁷ To identify whether patients with a STAD will also generate alternative 3ʹ-UTR lengths through APA heterogeneity and thus promote the progress of STADs, we compared the differences in the use of proximal PAS for all APA genes between STAD and normal samples. Taking CDC42 and ARF1 as examples, the PDUI scores in normal samples were closely distributed, while the distribution of PDUI scores in tumor samples was loose (Figure 3A and B). The distribution of PDUI in tumor samples had broader coverage, indicating that patients with a STAD showed greater heterogeneity in the use of proximal PAS. Although most genes did not show significant changes between normal and tumor samples (with a threshold of 0.015), the analysis of the use of proximal PAS in the case of all genes confirmed that 273 genes showed greater heterogeneity in tumor samples (the orange color), and only three genes showed significant heterogeneity in normal samples (the blue color, Figure 3C). There was no significant difference in the number of important heterogeneity genes between the 3′-UTR-shortened/lengthened genes (Figure S2). Using KEGG enrichment analysis, a subset of 273 genes was found to have primarily been enriched in pathways such as endocytosis, lysosome, and spliceosome pathways (Figure 3D).

Figure 3 STAD patients show great heterogeneity in the usages of proximal PAS. (A and B) Plotting examples of 3ʹ-UTR shortening genes (ALDOA, RAC1). There is a tightly distributed 3ʹ-UTR-shortening gene in the normal pancreas (purple); however, the degree of proximal PAS usages varies greatly between STAD patients (orange). (C) PDUI variance plot of all genes between tumor and normal samples. Purple dots indicate genes with high variance in normal samples, and orange dots indicate genes with high variance in tumor samples. The dotted line indicates the cutoff value of ± 0.015. (D) Top 10 KEGG pathways with significant enrichment associated with 273 genes with high variance in STAD patients.

Widespread Three Prime Untranslated Region Changes in Patients with Stomach Adenocarcinoma

To further analyze the distribution of 3ʹ-UTR-altered genes in patients with STAD, we selected 235 genes from the 271 genes reviewed for changes in PDUI scores; the selected genes had corresponding PDUI values but did not have NA value (ΔPDUI)-based cluster analysis. It included 135 3′-UTR-shortened genes and 100 3′-UTR-lengthened genes (Figure 4). We found that a subset of genes (n = 34, top heatmap, Figure 4) reflected 3′-UTR-shortened genes (the color red, ΔPDUI ≤ –0.2) in more than 80% of patients. Among them, the ARF1 gene showed 3′-UTR-shortened events in 99.7% of patients. A small subset of genes (n = 18, bottom heatmap, Figure 4) was found in 3′-UTR-lengthened genes to be recurrently lengthened (the blue color, ΔPDUI ≥ 0.2) in most patients with a STAD. This showed that APA events were broadly distributed among patients with STADs. Hierarchical clustering was performed and, accordingly, patients with a STAD were divided into five categories (subgroups 1–5), which corresponded to different APA gene ΔPDUI distribution patterns. Subgroup 5 (the color purple) included 116 patients with more 3′-UTR-shortened genes and fewer 3′-UTR-lengthened genes. In contrast, subgroup 1 (the color blue) included 42 patients with more 3′-UTR-lengthened and fewer 3′-UTR-shortened genes. These subgroups were not related to the most frequently mutated STAD genes (ie, TP53, ARID1A, FAT4, KMT2D, and PIK3CA) or patient stage.

Figure 4 The alternative 3ʹ-UTR was widespread in STAD patient. The heatmap showed the genes (rows) with shortened (red) or lengthened (blue) 3ʹ-UTR in each STAD patient (column) compared to the median value of each gene in normal stomach tissue. TP53, ARID1A, FAT4, KMT2D, and PIK3CA mutations and tumor stage profiles were displayed in the context of different APA distribution patient subgroups.

The Alternative Polyadenylation Events Contributed to the Assessment of the Prognosis of Patients with Stomach Adenocarcinomas

To further verify whether the APA event added prognostic information to patients with a STAD beyond the normal demographic and clinical factors (gender, race, age, tumor stage, and surgical outcome), we screened 21 genes that showed obvious 3ʹ-UTR changes and univariate prognostic value, as well as clustered patients with STADs, based on the PDUI scores of these 21 genes. Cluster analysis revealed three patient cohorts. Patients in cohort C showed an extensive degree of shortened 3ʹ-UTRs in all 21 genes; patients in cohort A showed shortening in the second group of genes, and patients in cohort B reflected primarily distal PAS (Figure 5A). There was a significant difference in the overall survival rate between cohorts A, B, and C (P = 0.04). The survival period of patients in cohort B was significantly longer than those in cohorts A (P = 0.05) and C (P = 0.03) (Figure 5B). In addition, the FPKM values of the most frequently mutated STAD genes (ie, TP53, ARID1A, FAT4, KMT2D, and PIK3CA) were compared. The results showed that there were differences in TP53, ARID1A, KMT2D, and PIK3CA among the three cohorts (Figure S3). The KMT2D and PIK3CA genes reflected relatively low expression in cohort B. Downregulation of KMT2D could potentially inhibit the proliferation and induce the apoptosis of gastric cancer,⁴¹ while the overexpression of PIK3CA may lead to a poor gastric cancer prognosis.⁴² Therefore, the APA distribution pattern can be used as an independent prognostic indicator in STAD cases, based on its potential for being a new prognostic biomarker.

Figure 5 The APA event determined the poor prognosis cohort of STAD patients. (A) Clustered patients based on PDUI score of 3ʹ-UTR altered genes and divided them into patient group A (blue)), patient group B (black) and patient group C (green). The 3ʹ-UTR altered genes were also clustered into 3 groups (1, 2 and 3). (B) Kaplan–Meier survival plots of patient A (blue), patient B (black) and patient C (green): *p ≤0.05.

Discussion

The disturbance of gene expression is the primary feature of cancer but our current understanding of its mechanism is limited.⁴³ Recently, APA was identified as an important regulatory mechanism for the dysregulation of gene expression.^17,44 Alternative polyadenylation is widely distributed in tumorigenesis. Some APA events are widely shared in different cancer types but many events are specific to tumor types.^18,44,45 At present, APA analysis and research on STADs are limited. This work comprehensively characterized APA in patients with a STAD including its distribution, the loss of miRNA loci, changes in parental gene expression, heterogeneity, patient distribution, and prognostic information.

Tumor heterogeneity is a characteristic of malignant tumors. This means that during the growth of a tumor, its daughter cells show changes in morphological and phenotypic profiles following multiple divisions and proliferation, making the tumor different in terms of characteristics such as growth rate, invasive ability, sensitivity to drugs, and prognosis.^46–48 Tumor heterogeneity can reflect different genetic backgrounds, such as differences in chromosome number and quality, different cell-case types, variable degrees of differentiation, and the diversity of cell evolution in different clinical stages.^47,49 An in-depth study of tumor heterogeneity will help to analyze the mechanism and evolution of tumor formation and explain the reasons for the different reactivity in tumor treatments.^46,48 Our data found that some of the genes of patients with STADs showed obvious heterogeneity to the extent of using proximal or distal PAS. The heterogeneity of these genes may be one of the key mechanisms causing tumor heterogeneity. Further study of the heterogeneity of the proximal PAS will support the screening of tumor biomarkers and therapeutic targets.

Endocytosis is the method by which large molecules or other cells enter a cell. The substance that enters the cell through the endocytosis pathway will participate in a series of life activities inside the cell, which can cause malignant transformations and contribute to the formation of cancer.^50,51 The endocytic circulation of cells can promote the plasticity, infiltration, and metastasis of cancer. However, there is little research evidence that genetic changes in endoproteins cause high incidences of cancer in humans. Our data confirmed that in the case of STAD, both significant APA genes and heterogeneous APA genes were enriched in the endocytosis pathway. The 3ʹ-UTR changes to these genes may also participate in the endocytosis cycle of tumor cells and promote the progression of STADs.

We observed that the distribution of APA events varied among patients with a STAD but most of these patients had more 3′-UTR-shortened genes. Studies have shown that in proliferating cells, the upregulation of genes related to cell growth was consistent with shortening of the 3ʹ-UTR,^52,53 which was consistent with the phenomenon we observed. We found that genes with a shortened 3ʹ-UTR were more likely to overexpress, possibly due to the regulation of missing miRNA binding sites. The miRNA can suppress gene expression by combining the complementary sequence in the 3ʹ-UTR of target mRNAs to degrade them, thereby preventing their translation.⁵⁴ Therefore, APA can affect the interaction of miRNA–mRNA, which is an important aspect of post-transcriptional regulation and is considered a fundamental mediator of gene expression involved in many types of cancers.⁸ In addition, our data showed that the APA gene in patients with a STAD could be used as an independent prognostic indicator of STAD results. However, since this article is mainly based on the analysis of second-generation sequencing data, more confirmatory experiments are urgently needed; the description of results in this paper may have been overestimated due to a lack of experimental data.

Overall, our research revealed that APA events regulated STAD genes that were functionally related, thereby providing a new way for better understanding the progress of STADs and identifying new drug targets.