InPAS is an open-source R/Bioconductor package for identifying novel proximal and distal CPSs, as well as quantifying APA from RNA-seq data of various experimental designs, using existing 3′ UTR annotations as guide. It is built on top of more than 20 other R/Bioconductor packages, including cleanUpdTSeq [72] for adjusting CPSs via a naïve Bayes classifier based on sequence features flanking CPSs, limma [73] for linear model-based differential APA usage testing, future and batchtools for parallel computing, and RSQLite for managing the metadata and intermediate files.

The InPAS workflow comprises four key steps: (1) extraction of 3′ UTR annotation, (2) manipulation and quality control of coverage data, (3) CPS search, and (4) differential APA analysis, as well as two auxiliary steps: (5) visualization of APA events, and (6) preparation of files for gene set enrichment analysis (GSEA). The key steps are shown in Additional file 2: Supplementary Fig. 1, and the entire workflow is depicted in Additional file 2: Supplementary Fig. 2, while the details on the InPAS algorithms are described in Additional file 3. Extraction of 3′ UTR annotation involves a series of steps including extraction of 3′ UTR annotation from a Gene Transfer Format (GTF) file for a reference genome, removal of any overlaps between 3′ UTRs of transcripts on opposite strands, collapse of 3′ UTRs of the same genes by their start coordinates, and extraction of downstream regions of the resulting 3′ UTRs up to a user-defined length or the nearest boundary of a downstream gene for later de novo CPS search, whichever is satisfied first. The third step, CPS search, is the most computation-intensive step, which is typically performed through four stages: (1) searching for distal CPSs (dCPSs), (2) searching for proximal CPSs (pCPSs), (3) adjusting dCPSs, and (4) adjusting pCPSs. To search for dCPSs, InPAS assesses the change pattern of the normalized merged coverage along each representative 3′ UTR and the associated extended region using a sliding window approach adopted from APAtrap [66]. Subsequently, pCPS search is performed using an optimal segmentation algorithm [74]. If the predicted dCPSs and pCPSs co-localize with annotated CPSs, they are replaced by the annotated CPSs. Otherwise, the predicted dCPSs and pCPSs can be further adjusted using a naïve Bayes classifier such as the one from the cleanUpdTSeqpackage [72]. Using a sliding window method, the classifier scores each base within a user-defined distance from the predicted CPSs based on a set of sequence features, including the occurrence of the core and auxiliary cis-elements. The initially predicted CPSs are replaced by sites of the highest scores as determined by the classifier if the highest score is greater than a user-defined threshold. In practice, users can utilize the classifier from our previously published CleanUpdTSeq package (available at https://www.bioconductor.org/packages/cleanUpdTSeq) for APA analysis of vertebrate RNA-seq data. Alternatively, users have the options to create a species- or clade-specific classifier using CleanUpdTSeq if sufficient polyadenylation data from the species or closely related species is available. Instructions for this are provided in the vignette (accessible at https://www.bioconductor.org/packages/devel/bioc/vignettes/cleanUpdTSeq/inst/doc/cleanUpdTSeq.html). If annotated polyadenylation data for a particular organism is lacking, users can either generate the necessary annotations to train a model or opt to use a classifier from a closely related species. Our published work on CleanUpdTSeq has shown that a classifier trained in zebrafish can be successfully applied to other vertebrates [72].

For each entry of 3′ UTR annotation with predicted pCPSs and dCPSs, mean coverage of the short and the long 3′ UTRs per sample are calculated. Differential APA usage between two conditions is quantified as a change in Percentage of Distal polyA site Usage Index ($\mathrm{\Delta }$PDUI) as in DaPars [64]. For differential APA analysis, InPAS is a versatile tool that can be applied to both single and multiple samples from a variety of groups across different experimental designs. Notably, for replicated experiments, InPAS employs the linear model and empirical Bayes statistical method from the limma package, which are designed to increase the statistical power for experiments with a limited number of replicates [73].

We used polyester (v1.33.0) [75] to generate 100-bp paired-end RNA-seq data with known CPSs and differential APA events as follows. Annotation for a subset of genes was extracted from the GENCODE GTF file (v34) for the human primary genome assembly (GRCh38) with the following criteria: (1) none of the selected genes overlapped; (2) each transcript had a 3′ UTR that was at least 1 kb in length and contained by a single terminal exon; (3) for multi-isoform genes, at most two isoforms were selected. If two isoforms were chosen for a gene, the two alternative 3′ UTRs should share the same start coordinate and differ in length by at least 1 kb. This process resulted in an abridged GTF file containing annotation for 5689 and 962 genes with one and two isoforms, respectively. Based on the abridged GTF file, RNA-seq data with a targeted sequencing depth of 20$\times$, 30$\times$, 40$\times$, 60$\times$, 80$\times$, 100$\times$, 120$\times$, and 140$\times$ were generated for two conditions without biological replicates. For bi-isoform genes, the relative expression proportion of one isoform ($p_1$) for each condition, was randomly sampled from a uniform distribution U(0, 1), while the relative expression proportion of the other isoform ($p_2$) was set to $1-p_1$. For single-isoform genes, the relative expression proportion was set to 1 regardless of conditions. The number of reads per transcript (c) was calculated using the formula $c=fdp/l$, where f is the length of a read pair in base pair, d is the targeted sequencing depth, p is the relative expression proportion of a given transcript, and l is the length in base of the given transcript. The relative expression proportions $p_1$and $p_2$ of a given transcript were kept the same across all sequencing depth. Read files in the FASTA format were generated by calling the simulate_experiment_countmat function. Theoretically, 566 of the 962 bi-isoform genes had significant differential APA between the two conditions ($|$$\mathrm{\Delta }$PDUI$|$ $\ge$0.2 and $|$log₂(fold change of PDUI)$|$ $\ge$0.59), which was used as ground truth to benchmark the performance of InPAS (v2.8.0), DaPars (v1.0.0) [64], DaPars2 (v2.1) [65], APAtrap (v1.0) [66], and TAPAS [67].

The simulated RNA-seq reads were mapped to the human primary genome assembly GRCh38 with the abridged GTF file as gene annotation using STAR (v2.5.3a) [76] in the basic two-pass mode. The resulting alignment files in the SAM format were converted to the BAM format and sorted by coordinates using SAMtools (v1.9) [77]. Genome coverage was calculated based on unique alignments, using the bedtools genomcov module (v2.29.0) [78] with the following settings: -bg -split. The resulting bedGraph file was used for APA analysis.

We generated two TxDb databases (TxDB_1 and TxDB_2) using the makeTxDbFromGFF function from the GenomicFeatures package. TxDB_1 was created with the abridged GTF file containing the annotation for the 5689 single-isoform genes and the short transcript isoform of the 962 bi-isoform genes, while TxDb_2 was created based on the annotation for the long transcript isoform of the 962 bi-isoform genes. We performed InPAS-based APA analysis on the bedGraph files of different coverage depths separately using each of the two databases, with the single-isoform genes in the TxDb_1 as the negative controls. For CPS searches, only extracted 3′ UTRs whose first 100-base mean coverage was at least 20$\times$ were considered. To assess the effectiveness of refining pCPSs, the naïve Bayes classifier from the cleanUpdTSeq package was used to adjust pCPS in one round of analyses but not in the other. Differential APA analysis was only conducted for transcripts with both computed pCPSs and dCPSs. CPSs predicted by InPAS were considered correct if the predicted CPSs were within 50-bp distance of the annotated CPSs. Differential APA events detected by Fisher’s exact test with false discovery rate (FDR) $\le$0.05, $|$$\mathrm{\Delta }$PDUI$|$ $\ge$0.2, and $|$log₂(fold change of PDUI)$|$ $\ge$0.59 were considered statistically significant.

The abridged GTF files containing the annotation for the 5689 single-isoform genes and that of the 962 bi-isoform genes was converted into a 12-column BED file using the ea-utils gtf2bed module (https://github.com/ExpressionAnalysis/ea-utils). The BED file was used to generate 3′ UTR annotation files using the DaPars_Extract_Anno.py script. CPS prediction was performed using the DaPars_main.py scripts with the same bedGraph file specified twice to provide coverage data for a two-group analysis in a configuration file. Threshold parameters in the configuration file were set as follows: Num_least_in_group1 = 1, Num_least_in_group2 = 1, Coverage_cutoff = 20, FDR_cutoff = 0.05, PDUI_cutoff = 0.2, Fold_change_cutoff = 0.59. APA events with FDR $\le$0.05, $|$$\mathrm{\Delta }$PDUI$|$ $\ge$0.2, and $|$log₂(fold change of PDUI)$|$ $\ge$0.59 were considered statistically significant. Differential APA analysis was performed using the DaPars_main.py scripts as described for CPS prediction, except that a pair of bedGraph files from two simulated conditions at a specific sequencing depth were specified to provide coverage data for a two-group analysis in the configuration file.

As described for DaPars-based analysis, 3′ UTR annotation files were generated except that the DaPars_Extract_Anno.py script from the DaPar2 software was used. The number of mapped reads for each sample was calculated using SAMtools (v1.9) [77]. CPS prediction was performed using the DaPars2_Multi_Sample_Multi_Chr.py script with a configuration file as required.

For individual sample-based CPS detection, distal 3′ UTRs of each sample were first predicted using the identifyDistal3UTR command, with the same BED file for DaPars-based analysis as gene annotation and individual bedgraph file provided as coverage data. CPS detection was performed using the predictAPA command with the following option settings: -g 1 -n 1.

For differential APA analysis, distal 3′ UTRs of each sample were first predicted using the identifyDistal3UTR command, with the same BED files for DaPars-based analysis as gene annotation and the two bedgraph files from the two simulated conditions at a given sequencing depth provided as coverage data. CPS detection was performed using the predictAPA command with the following option settings: -g 2 -n 1 1. Differential APA analysis was performed using the deAPA function in the deAPA package bundled with APAtrap.

TAPAS-specific annotation files were generated from the same BED files for DaPars-based analysis using a custom Python script. CPS detection was performed using the APA_sites_detection command with read length and individual bedgraph file specified. Differential APA analysis was performed using Diff_APA_site_analysis command with the predicted CPSs for each sample by the APA_sites_detection command as input.

We obtained three sets of RNA-seq data and matched 3′-end RNA-seq data from the GEO and Sequence Read Archive (SRA) databases: GSE42420 [23] and GSE94950 [50], GSE151721 and GSE151724 [79], GSE56010 [80] and SRP065825 [81]. These data were from human or mouse cell lines treated with control siRNA or siRNA against CFIm25 [23, 50], HNRNPC [80, 81], Srsf3 or Srsf7 [79]. Additionally, we downloaded RNA-seq data (E-MTAB-2836) for human cerebral cortex and liver tissues from ArrayExpress [82], as well as RNA-seq data (GSE87410) for 50 pairs of tumor and normal tissues of nine cancer types [83] from the GEO database. Aligned RNA-seq data in the BAM format for both primary tumors and matched normal tissues for the TCGA BRCA (n = 113 patients) [84] and TCGA KIRC (n = 72 patients) [85] projects were downloaded from the GDC data portal (https://portal.gdc.cancer.gov/) [86]. Additional information about these datasets is provided in Additional file 1: Supplementary Table 2.

The quality of the raw sequencing data was assessed using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) before preprocessing. For PAC-seq data (GSE94950) [50], Cutadapt [87] was used to remove TruSeq adaptor sequences, the first six nucleotides derived from the ‘Click-Adaptor’, and poly(G) tails with explicit option settings: -a NNNNAGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -m 35 -j 10 –nextseq-trim 20. Next, custom Perl scripts were used to further trim the resulting reads to remove poly(A) tails of length $\ge$15 bases (allowing one mismatch) and discard reads shorter than 25 bases after trimming or reads with less than 15 trailing A’s prior to trimming. The MACE-seq (GSE151722) and A-seq2 (SRP065825) data were preprocessed as described [6, 79, 81]. The resulting reads were mapped to the reference genome primary assembly from GENCODE (human: GRCh38; mouse: GRCm39) using STAR (v2.5.3a) [76] in the basic two-pass mode, with a GTF annotation file (human: GENCODE v34; mouse: GENCODE m33) for the corresponding primary assembly as gene annotation. Potential CPSs with at least two supporting reads were extracted using custom Perl scripts. The potential CPSs were scored using DeepPASTA, a deep learning method for CPS analysis considering both RNA sequence and secondary structures of sequence flanking potential CPSs [88]. Sites with scores $\ge$0.5 were considered as true CPSs. CPSs per library were clustered as described [6] using custom Python scripts.

The quality of the raw RNA-seq data (GSE42420, GSE151721, GSE56010, E-MTAB-2836, and GSE87410) was assessed using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc). Subsequently, Illumina sequencing adapter sequences and low-quality bases were trimmed using fastp (v0.21.0) [89] with the following explicit option settings: –cut_tail –cut_tail_window_size 4 –cut_tail_mean_quality 15 –cut_front –cut_front_window_size 1 –cut_front_mean_quality 3 –adapter_fasta adapter.fa –length_required 25. Outputs from fastp was further processed to remove trailing poly(A) or leading poly(T) of length $\ge$10 bases and only reads of lengths $\ge$25 bases after processing were retained, using custom Perl scripts. This extra trimming step aimed to increase the mapping rate of reads with trailing poly(A) or leading poly(T), which were likely derived from 3′ ends of polyadenylated transcripts. The resulting reads were paired up using the BBTools (v38.20) repair.sh scripts (https://jgi.doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/repair-guide/). The paired-end reads were mapped as described above for the 3′-end RNA-seq data. BAM files were then coordinate-sorted using SAMtools (v1.9) [77], and genome coverage was calculated using the bedtools (v2.29.0) genomcov module [78] as described above. For the TCGA data, BAM files for technical replicates of three tumor samples were first merged, and then genome coverage was calculated using the TCGA coordinate-sorted BAM files.

Due to the unavailability of fastq files for the TCGA data, we regenerated fastq files from BAM files using the SamToFastq module of the Picard Tools (v2.26.10) (https://broadinstitute.github.io/picard/). The paired-end reads from GSE42420, GSE151721, GSE56010, E-MTAB-2836, and GSE87410, and the TCGA data were mapped to the human reference transcriptome (GENCODE v34) with the whole genome as decoy, using Salmon (v1.9) [90] with explicit options setting: –seqBias –gcBias –posBias –softclip –softclipOverhangs –biasSpeedSamp 5 -l IU –validateMappings. Next, transcript-by-sample count tables for each dataset were generated using the tximport package [91]. Mean depth of coverage per transcript was calculated using the formula 2lc/L, where l is the read length in base, c is the estimated read counts for a given transcript, and L is the computed effective length in base of the given transcript. The mean depth of coverage per transcript was summed up across all relevant samples for a given differential APA analysis to obtain the total depth of coverage of a given transcript. For each pair of compared conditions, we filtered out transcripts with a total depth of coverage lower than 10, noncoding transcripts of protein-coding genes, transcripts encoding small nucleolar RNAs (snoRNAs), miRNAs, ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), rRNA pseudogenes, small activating RNAs (saRNAs), small conditional RNAs (scRNAs), small nuclear RNAs (snRNAs), and vault RNAs, as well as other transcripts shorter than 200 bases from the corresponding count table. Annotation for the remaining transcripts was extracted from the original GTF file, resulting in an abridged GTF file for the given comparison, which was used for further APA analyses.

The abridged GTF files for each comparison were converted into 12-column BED files as described above. CPS prediction and differential APA analysis of the experimental RNA-seq data were performed as described for the simulated RNA-seq data. To evaluate the accuracy of predicted CPSs by these tools, we compared the predicted CPSs from the RNA-seq data with those identified from the 3′-end RNA-seq data as references. CPS predicted from RNA-seq data (GSE42420, GSE151721, GSE56010) were considered correct if the predicted CPSs were within 50-bp distance of the CPSs detected by 3′-end RNA-seq data (GSE94950, GSE151724, and SRP065825). An APA event with an occurrence rate $>$10% across patients within each dataset (GSE87410 and the TCGA dataset) was considered recurrent.

To benchmark InPAS along with other leading tools (DaPars, DaPars2, APAtrap, and TAPAS) for APA analysis, we first simulated a series of RNA-seq datasets of different depths of coverage (20$\times$ to 140$\times$), where the expression level of each transcript of 962 bi-isoform genes and 5,689 single-isoform genes was predefined and the CPSs of all transcripts were known. Using these simulated data, we showed that the naïve Bayes model-based adjustment of CPSs improved the precision (Fig. 1a, Ref. [64, 66, 67]) and sensitivity (Fig. 1b) of InPAS in detecting proximal CPSs (pCPSs), although the effect on the sensitivity was not as pronounced. Subsequently, we evaluated the InPAS’s performance in comparison to other leading tools using the same simulated data. InPAS consistently exhibited significantly higher precision in the detection of both pCPSs and dCPSs across all tested coverage depths, outperforming all the other four tools (Fig. 1c). The overall sensitivity of InPAS was only slightly lower than the nearly identical sensitivities of DaPars and DaPars2 (Fig. 1d), though DaPars and DaPars2 naively consider the 3′-extremety of each transcript in a GTF file as an authentic dCPS without examining the actual read coverage under a given condition. The overall sensitivity of InPAS was similar to that of APAtrap at lower coverage depths (20$\times$ and 30$\times$), but much higher than that of APAtrap at higher coverage depths, and much higher than those of TAPAS at all coverage depths (Fig. 1d). When the short 3′ UTRs of the 962 bi-isoform genes were provided as 3′ UTR annotation, InPAS demonstrated higher precision and sensitivity than APAtrap in de novo predicting distal CPSs (dCPSs) at higher coverage depths, whereas none of DaPars, DaPars2, or TAPAS can detect dCPSs de novo due to their implementation limits (Fig. 1h,i). Therefore, we compared the performance of those tools in detecting pCPSs. Again, InPAS showed much higher precision and sensitivity in detecting pCPSs than all the other tools (Fig. 1e,f). With the 3′ UTRs of 5689 single-isoform genes as true negatives of APA, InPAS demonstrated high specificity in detecting pCPSs, with a specificity of greater than 91%, which was markedly higher than those of the other leading tools at all depths of coverage tested (Fig. 1g). InPAS showed higher precision and F1 scores than DaPars and APAtrap in detecting differential APA events (Fig. 1j,l), while DaPars2 lacks the functionality to perform differential APA analysis and TAPAS requires at least two replicates for each condition to perform differential APA analysis. InPAS showed lower sensitivity in detecting differential APA than APAtrap at all tested coverage depths, slightly higher sensitivity than DaPars at low coverage depths (20$\times$ to 40$\times$) but lower sensitivity than DaPars at higher coverage depths ($>$60$\times$) (Fig. 1k). Overall, our benchmarking using the simulated RNA-seq dataset showed that InPAS outperformed the other leading tools in terms of precision in predicting all CPSs as well as recovering differential APA events. Furthermore, InPAS exhibited significantly higher sensitivity, precision, and specificity in predicting pCPSs than the other tools.

Fig. 1.

Comparison of InPAS with other leading tools in identifying and quantifying APA from simulated RNA-seq data. (a,b) Precision and sensitivity of InPAS in predicting pCPSs is slightly improved by adopting the naïve Bayes classifier from the cleanUpdTSeq package to adjust the rudimentarily predicted pCPSs at different depths of coverage. (c,d) Precision and sensitivity of InPAS and other tools in predicting both dCPS and pCPSs. (e–g) Precision, sensitivity, and specificity of InPAS and other tools in predicting pCPSs. Subpanels (d–g) shared the same figure legend as (c). (h,i) Precision and sensitivity of InPAS and other tools in predicting dCPSs when annotation for the short 3′ UTRs is provided as 3′ UTR annotation. None of TAPAS, DaPars and DaPars2 can predict the dCPSs. (j–l) Precision, sensitivity and F1 score of InPAS, DaPars, APAtrap and TAPAS in recovering known differential APA events. Fisher’s exact test was used by DaPars [64] and InPAS, linear trend test by APAtrap [66], or negative binomial regression by TAPAS [67]. DaPars2 cannot perform differential APA analysis, so it was not included in the evaluation. CPS, cleavage and polyadenylation site; pCPSs, proximal CPSs; dCPS, distal CPSs; UTR, untranslated region; InPAS, Identification of Novel alternative PolyAdenylation Sites; DaPars, Dynamic analysis of Alternative PolyAdenylation from RNA-seq; APAtrap, Alternative PolyAdenylation trap; TAPAS, Tool for Alternative Polyadenylation site AnalysiS; RNA-seq, RNA sequencing.

Next, we evaluated the performance of InPAS and the other leading tools mentioned above using experimental datasets: three sets of published RNA-seq data with matched 3′-end RNA-seq data [23, 50, 79, 80, 81] (Additional file 1: Supplementary Table 2). These data were from human or mouse cell lines treated with control siRNA or siRNA against CFIm25 [23, 50], HNRNPC [80, 81], Srsf3 or Srsf7 [79]. Previous studies have shown that CFIm25 [23], HNRNPC [92, 93], SRSF3 and SRSF7 [79] can all modulate 3′ UTR length. We first determined the CPSs of each sample using the 3′-end RNA-seq data, which were used as imperfect truth for evaluating the performances of the tools for APA analysis. More than 75% of the CPSs detected in each sample overlap with the entities in the PolyASite 2.0 database [6], though the number of detected CPSs across samples under the same conditions varies drastically (Additional file 2: Supplementary Fig. 3). InPAS showed higher precision and F1 scores for most samples in predicting CPSs than other tools (Fig. 2a,c), while its sensitivity in predicting CPSs slightly lower than those of DaPars and DaPars2, which displayed the best sensitivity (Fig. 2b). Given that none of DaPars, DaPars2, and TAPAS can predict dCPSs de novo, we compared the precision of these tools in predicting pCPSs. InPAS showed the best precision in predicting pCPSs (Fig. 2d).

Fig. 2.

Comparison of InPAS with other leading tools in identifying CPSs from experimental RNA-seq data. (a,b) Precision and sensitivity of InPAS and other tools predicting both pCPSs and dCPSs. (c) F1 scores of InPAS and other tools predicting both pCPSs and dCPSs. (d) Precision of InPAS and other tools predicting pCPSs.

Besides precision, sensitivity, and specificity, we compared the functionalities, runtime, and maximal memory consumption of these tools for APA analyses using the simulated and experimental RNA-seq. Table 1 shows the differences of these tools in functionalities [64, 65, 66, 67]. APAtrap behaved as the most time-consuming tool for both predicting CPSs and differential APA analysis (Fig. 3a,b). InPAS ran only slightly slower than DaPars and DaPars2 in predicting CPSs, where TAPAS was the fastest. For differential APA analysis, APAtrap and TAPAS performed as the first and second most time-consuming tool, respectively, followed by DaPars, while InPAS was the fastest tool (Fig. 3b). In terms of maximal memory consumption, TAPAS consumed the largest amount of memory during differential APA analysis, followed by APAtrap, InPAS, and DaPars (Fig. 3c). The maximal memory required by DaPars and InPAS were small, and comparable (Fig. 3c).

Fig. 3.

CPU time and maximal memory consumption of InPAS and other tools. (a) CPU time used by each tool to predict CPSs of each sample. The y axis is on a log₁₀ scale. Each dot represents a sample. HEK293, HeLa, and P19 are names of cell lines, from which the RNA was extracted for RNA-seq. Sim, the simulated RNA-seq data using polyester. (b) CPU time used by each tool to detect differential APA events between replicates treated by a control siRNA or an siRNA against HNRNPC (HEK293), CFIm25 (HeLa), Srsf3 (P19), or Srsf7 (P19). (c) Maximal memory required by each tool to get differential APA events between replicates treated by a control siRNA or an siRNA against HNRNPC (HEK293), CFIm25 (HeLa), Srsf3 (P19), or Srsf7 (P19).

To showcase the ability of InPAS to identify differential APA in physiological states, we analyzed RNA-seq data from normal human brain and liver tissues [82]. Our analysis revealed that for 165 genes, the dCPS was used more favorably in the brain than in the liver, whereas for 22 genes, the pCPS was used more favorably in the brain than in the liver (Additional file 1: Supplementary Table 3), which is consistent with the well-documented phenomenon of widespread lengthening of 3′ UTRs in neuronal tissues via APA [94, 95]. Additional file 2: Supplementary Fig. 4 highlights several of the genes (LRP3, SRSF7, YY1, and ZNF706) that exhibited differential APA as identified by InPAS. Additionally, we observed an overrepresentation of genes related to pathways of neurodegeneration–multiple diseases (KEGG pathway hsa05022) among the genes showing differential APA in the brain compared to the liver (gene ratio = 0.20, adjusted p-value = 0.0038).

Next, we tested the ability of InPAS to identify differential APA events in heterogeneous tumor samples compared to matched normal tissue samples by analyzing two datasets, one from 50 patients with nine different cancer types [83] and the other from 113 patients with breast invasive carcinoma (BRCA) [84] and 72 patients with kidney renal clear cell carcinoma (KIRC) [85]. Although the number of differential APA events detected for each patient varied in both datasets, we identified 1162 recurrent differential events among the 50 patients with diverse cancers, and 941 recurrent differential events among the 185 patients with BRCA or KIRC, with a large fraction of differential APA events resulting in preferential use of pCPS (3′ UTR shortening) in tumor samples (Fig. 4a, Additional file 1: Supplementary Tables 4 and 5, and Additional file 2: Supplementary Fig. 5). Our findings are consistent with a previous study [64] showing widespread preference for utilizing pCPSs in cancers. The genes exhibiting differential APA in the nine cancer types include 74 genes from the COSMIC-curated Cancer Gene Census list (v96) (http://cancer.sanger.ac.uk/census), 27 tumor suppressor genes, 24 oncogenes, seven “double agents” (genes with both oncogenic and tumor-suppressor functions) [96], and 14 fusion genes (Additional file 1: Supplementary Table 4). Similarly, genes involved in differential APA in the TCGA data include 31 genes from the COSMIC-curated Cancer Gene Census list (v96) (http://cancer.sanger.ac.uk/census), nine tumor suppressor genes, 11 oncogenes, seven “double agents” and four fusion genes (Additional file 1: Supplementary Table 5). Genes involved in the p53 signaling pathways and a few other KEGG and Reactome pathways were overrepresented among the genes with differential APA between tumor and normal samples from patients of the nine cancer types (Fig. 4b,c). We also observed an overrepresentation of genes in the renal cell carcinoma pathway (KEGG hsa05211), among the genes with differential APA in the TCGA BRCA or KIRC samples compared to the matched normal tissues (gene ratio = 0.033, adjusted p-value = 0.031).

Fig. 4.

Differential APA events identified by InPAS across nine cancer types. (a) The central heatmap displays 1162 recurrent APA events across the nine cancer types. The top stacking barplot shows the number of 3′ UTR shortening (red) and lengthening (blue) events in each tumor sample compared to the matched normal tissue. The left barplot shows the proportions of samples across the nine cancer types with 3′ UTR shortening, while the right stacking shows the number of samples across the nine cancer types with 3′ UTR shortening (red) and lengthening (blue). GAC, gastric adenocarcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; HCC, hepatocellular carcinoma; SCLC, small-cell lung carcinoma; CSCC, cervical squamous cell carcinoma; ESCC, esophageal squamous cell carcinoma; SRCC, gastric signet-ring cell carcinoma; PTC, papillary thyroid carcinoma. (b,c) KEGG and Reactome pathways overrepresented among genes showing differential APA in tumor samples compared to the matched normal tissues, respectively. (d) Heatmap showing the log₂(fold change) of expression of genes involved in mRNA 3′ end processing (GO0031124) between the tumor and the matched normal samples. The significant changes (adjusted p value 0.05) are labeled with asterisks.

To investigate the mechanisms that drive APA during tumorigenesis, we examined the genes showing recurrent APA in the two cancer datasets for overlap with a collection of 127 significantly mutated genes in cancers, curated by the TCGA PanCan12 project [97]. However, only a small fraction of genes with recurrent differential APA in the two datasets—17 (1.46%) of genes in the BRCA dataset and 21 (2.23%) of genes in the KIRC dataset—were PanCan12 significantly mutated genes (Additional file 1: Supplementary Tables 4 and 5). These results suggest that the vast majority of APA events in the two cancer datasets may not be caused by known mutations of the cis-elements in 3′ UTRs. On the other hand, we found that many genes involved in 3′-end processing of pre-RNA transcripts showed aberrant expression in most tumor samples compared to their matched normal tissues (Fig. 4d, and Additional file 2: Supplementary Figs. 6,7, and 8), which is consistent with a previous study [64] and suggest dysregulation of these 3′-end processing factors might contribute to many APA events. Furthermore, our analysis suggests that the genes and the directions of their aberrant expression in tumors relative to matched normal tissues may be specific to tumor type or even patient.