Introduction
The 3' untranslated region (3’-UTR) of a messenger RNA (mRNA) molecule is a critical segment located downstream of the coding sequence and plays a pivotal role in gene regulation. Polyadenylation sites (pA sites) are specific locations, primarily within the 3’-UTR but also found in introns, coding sequences, and intergenic regions, where the mRNA undergoes a process known as polyadenylation, resulting in the addition of a poly(A) tail to the pA site. Alternative polyadenylation (Mitschka & Mayr, 2022) refers to the phenomenon in which a single gene possesses multiple pA sites within its 3’-UTR, leading to the production of distinct mRNA isoforms with varying 3’-UTR lengths.
These alternative polyadenylation isoforms can have profound implications for post-transcriptional regulation, as they may influence mRNA localization, stability, and interaction with regulatory factors (Chen et al., 2024) (Mitschka & Mayr, 2022), thereby affecting the ultimate functional output of the gene. Understanding the mechanisms and consequences of alternative polyadenylation is essential to unravel the complexities of gene expression and its role in diverse biological processes.
We previously developed the software package SCAPE (Zhou et al., 2022) that can estimate alternative polyadenylation events from single-cell RNA-seq (scRNA-seq) data. Because the reads in scRNA-seq are generally not long enough to cover the pA sites, we used a Bayesian mixture model (Figure 1A) to infer the pA sites using the insert length distribution. The mixture model contained k (k > 1) Gaussian components and one uniform component. Each pA site is represented by a Gaussian component, which has three parameters: mean μ, standard deviation δ and weight π. Detailed meanings of these parameters and model specification are provided in Supplmentary Note 1. The full SCAPE analysis workflow can be summarized as the following steps: (1) extract the 3’-UTR annotation from the General Transfer Format (GTF) or General Feature Format (GFF) file; (2) map the reads to the reference genome and calculate the relative positions of reads on the corresponding 3’-UTR regions; (3) perform parameter inference for each 3’-UTR region and output the result; and (4) perform downstream tasks such as differential pA site usage analysis.
Figure 1.
(A) Original SCAPE model and improvements of SCAPE-APA. In data preprocessing, we developed a new script to extract 3’-UTR annotations from GTF files and implemented read binning. In modeling, we introduced a new likelihood function for binned reads and derived an updated inference algorithm, as detailed in Supplementary Note 1. In implementation, we rewrote SCAPE using Taichi to enhance computational efficiency. For both SCAPE-APA and SCAPE, the outputs include the estimated mean and standard deviation of pA sites (represented as α and β in the bottom panel). In the SCAPE model, pA site locations are modeled as Gaussian components. (B) Schematic illustration of a spurious pA site due to splicing of 3’-UTR. The SCAPE model assumes that each 3’-UTR is a single continuous segment; however, in reality, 3’-UTRs can be multi-segmented. This results in reads that are close together on the transcript appearing far apart on the genome, violating the SCAPE model's assumption and potentially leading to spurious pA sites. In the figure, orange represents junction reads, while gray indicates normal reads. (C) A real example of a spurious pA site. PA1 is the spurious pA site and should be merged into PA2. The bottom panel displays 3’-UTR annotations from the GTF file, showing one single-segment 3’-UTR transcript and two multi-segment 3’-UTR transcripts. PA1 and PA2 are most likely derived from the first transcript, necessitating the merging of PA1 into PA2. (D) Comparison of running time across different datasets on the same computer. Each 3’-UTR in these datasets contains more than 20,000 reads. All datasets were analyzed using a single CPU core with 10 GB of memory. (E–G) Comparison of estimated mean, standard deviation (std) and weight of estimated pA sites. The curves are the density contour of data points. The means of the Gaussian components indicate the predicted pA site locations. For a detailed explanation of these parameters, please refer to Supplementary Note 1.
However, computational speed remains a challenge for SCAPE in real data analysis. Some 3’-UTRs contain a large number of reads (>200,000), significantly slowing down the inference algorithm. SCAPE often processes most 3’-UTRs efficiently but struggles with a small subset of high-read-count 3’-UTRs. Simple parallelization strategies, such as processing chromosomes independently or using multithreading, are ineffective because SCAPE’s inference algorithm is inherently serial. Additionally, SCAPE is written in Python, which is slower than just-in-time compiled languages like Taichi (Hu et al., 2019). Code maintainability is also an issue, particularly in the inference algorithm and 3’-UTR annotation extraction scripts. To address these challenges, we developed a new version of SCAPE.
Methods
Improvements
We made significant updates to the previous implementation, as shown in Figure 1A. The biggest challenge in real data analysis is that the read number can be very large (n > 200,000) for some 3’-UTRs. This makes the original inference algorithm extremely slow. As the 3’-UTR is relatively short (~1500bp), there are many redundant reads. To solve this problem, we binned reads in the vicinity, that is, reads with similar x, l, r values in the following likelihood.
Note that the bin is 3-dimensional. The nth bin was denoted by the unique combination (xn, ln, rn) and its count by cn, where the bin sizes were set to 5bp, 10bp, 10bp, respectively. The joint likelihood of the mixture model can be rewritten as:
It can be seen that we only need to calculate the likelihood once for all cn reads in the nth bin in the new likelihood, whereas this has to be repeated cn times in the original likelihood. As the likelihood changes, the inference algorithm must be modified accordingly. The detailed mathematical derivations of the inference algorithm for the new likelihood are provided in Supplementary Note 1 (see Extended data). To further accelerate the inference, we used the Taichi language (Hu et al., 2019) to re-implement the inference algorithm, where Taichi used the just-in-time (JIT) compilation to accelerate the Python code. Because the inference algorithm is serial, the computation times are similar, regardless of whether the GPU is used.
In addition to the new likelihood and inference algorithm, we re-implemented and added new bioinformatics functionalities to streamline the workflow. We re-wrote the Python script to extract 3’-UTR annotations from GTF or GFF files based on the package “gffutils” (https://github.com/daer/gffutils), which was easier to maintain than the previous ad hoc implementation. We added a post-processing script to handle spurious pA sites owing to the junction reads. As shown in Figure 1B, some 3’-UTR may contain multiple exons/fragments, while the SCAPE model assumes the 3’-UTR to be a consecutive region on the genome, that is, we take the start of the first exon/fragment and the end of the last exon/fragment as the full 3’-UTR region. Because of this multi-segment property of the given 3’-UTR in Figure 1B, two pA sites were estimated from the data, where the orange pA site was the spurious pA site. To resolve this, we merged a pA site with the next downstream pA site if there were > 40% junction reads in all reads assigned to this pA site. Figure 1C shows a spurious pA site case in the real data, where the first estimated pA site should be merged with the second. Finally, we made a python package “scape-apa” and uploaded it to the Python Package Index (PyPI) to facilitate the installation.
We compared the speed and inference results of the old and new SCAPE software by sampling 100, 500, 1000, and 2000 3’-UTRs from public scRNA-seq data (Notaras et al., 2022). Each 3’-UTR analyzed contains at least 20,000 reads to effectively illustrate the speed difference between the new and old versions of SCAPE. Figure 1D shows the running times for the analyses of these datasets on a CPU-only computer (Dell PowerEdge C4130), with the maximum memory set to 10GB. It can be seen that the speedup is 9-fold for the dataset with 1000 3’-UTRs. Figure 1E–G shows the estimated mean, std, and weights given by the old and new SCAPE packages of the dataset with 1000 3’-UTR. We fixed the number of pA sites in each UTR based on the results of the new package. It can be seen that the majority of the data points are scattered around the y = x line for the mean and weight, whereas a larger variation is observed for the standard deviation. These differences are due to the binning operation; however, the overall comparisons suggest a highly consistent result between the old and new packages. Additionally, this benchmark focuses specifically on high-read-count 3’-UTRs, as the performance of both versions remains similar for low-read-count cases.
Spurious pA site detection is a new feature of SCAPE-APA, as the original SCAPE package does not detect spurious pA sites. In our test dataset (Notaras et al., 2022), SCAPE-APA identified 191 spurious pA sites out of 1009 total pA sites across 1000 3’-UTRs, each containing at least 20,000 reads in scRNA-seq data. In summary, we revised the statistical model for binned reads and utilized Taichi to accelerate the parameter inference of SCAPE, as well as developed useful Python scripts to process GTF files and handle spurious pA sites.
Implementation
The SCAPE-APA package has five main modules:
utils.py: This module processes GFF files from the Ensembl database into 3’-UTR annotations using the dedicated APIs of gffutils. If no 3’-UTR annotation is available for a gene, the last exon is used as its 3’-UTR.
input_processor.py: This module extracts relevant information from the input BAM files for APA event estimation. This includes variables such as x, l, r in the SCAPE model, as well as cell barcode, UMI, and junction read information.
apa_core.py: This module contains the core inference algorithm of the SCAPE model. The ApaModel class implements the core functionality, holding the parameters and results for the model. Its run() function gathers all relevant components for the inference algorithm. Details of the statistical model and inference algorithm are provided in Supplementary Note 1.
taichi_core.py: This module provides equivalent implementations of relevant functions from the inference algorithm using the Taichi language. Since Taichi is a static language, all memory for inputs and outputs must be pre-allocated before calling a Taichi function. These functions are imported by the apa_core.py module to perform the actual inference.
junction_handler.py: This module merges spurious pA sites that arise due to junction reads. If more than 40% of the reads assigned to a given pA site are junction reads, the site is labeled as spurious and merged with the next pA site.
Operation
SCAPE-APA can be installed on Linux, Mac, and Windows systems with Python 3.11 or higher. A GPU is not required to run SCAPE-APA. The basic workflow consists of five steps:
1. Generate the 3’-UTR annotation file using the “scape gen_utr_annotation” command. The input file of this step is GFF3 file of suitable genomes downloaded from Ensembl database.
2. Prepare the input files using the “scape prepare_input” command. The 3’-UTR annotation file from step 1 and the BAM file from the scRNA-seq experiment are the inputs for this step.
3. Perform pA site inference using the “scape infer_pa” command.
4. Remove spurious pA sites caused by junction reads using the “scape merge_pa” command.
5. Extract the pA count matrix using the “scape ex_pa_cnt_mat” command. The rows represent cells, and the columns represent the estimated pA sites. The pA count matrix can be fed into Seurat for differential analysis.
Discussion
Since 2019, the scale of scRNA-seq data has expanded dramatically, growing from datasets with thousands of cells to those containing millions in 2024. This surge in data volume has highlighted the limitations of the original SCAPE tool, which was not equipped to handle such large datasets. In fact, the original SCAPE was extremely slow, sometimes taking weeks to analyze a single dataset with just thousands of cells. Moreover, it contained several bugs that compromised its performance, and the reliance on ad-hoc Python scripts made maintenance difficult, especially as databases like Ensembl regularly update their annotations. Using a dedicated API to extract 3’-UTR annotations is a far more robust and reliable solution for long-term sustainability.
SCAPE-APA was developed to address these critical issues. It offers substantial improvements in speed, reliability, and maintainability, making it far better suited for the analysis of large-scale scRNA-seq datasets. By incorporating modern bioinformatic tools and creating a formal Python package, SCAPE-APA enables faster and more efficient analysis while reducing maintenance challenges. This updated version will be an invaluable tool for researchers, capable of handling the increasingly large and complex datasets in single-cell RNA sequencing, ensuring that SCAPE-APA meets the growing demands of the field.
Ethics and consent
Ethical approval and consent not required
Author contribution
G. C. developed an inference algorithm for binned reads and implemented it in Taichi. T. L. developed Python scripts for parsing the GTF file, removing spurious pA sites, and ran benchmark experiments. R. Z. provided advice regarding the 3’-UTR annotation extraction. L. C. supervised the project.
Acknowledgements
GC, TL, and LC acknowledge the computational resources provided by the Aalto Science-IT project and the CSC – IT Center for Science, Finland.
Faculty Opinions recommendedReferences
- Cheng L:
SCAPE-APA Supplementary Note 1.
Zenodo.
2024. http://www.doi.org/10.5281/zenodo.14517304
- Chen H, Wang Z, Gong L, et al.:
A distinct class of pan-cancer susceptibility genes revealed by an alternative polyadenylation transcriptome-wide association study.
Nat Commun.
2024; 15(1): 1729. PubMed Abstract
| Publisher Full Text
| Free Full Text
- Hu Y, Li TM, Anderson L, et al.:
Taichi: a language for high-performance computation on spatially sparse data structures.
ACM Trans Graph.
2019; 38(6): 1–16. Publisher Full Text
- Mitschka S, Mayr C:
Context-specific regulation and function of mRNA alternative polyadenylation.
Nat Rev Mol Cell Biol.
2022; 23(12): 779–96. PubMed Abstract
| Publisher Full Text
| Free Full Text
- Notaras M, Lodhi A, Dündar F, et al.:
Schizophrenia is defined by cell-specific neuropathology and multiple neurodevelopmental mechanisms in patient-derived cerebral organoids.
Mol Psychiatry.
2022; 27(3): 1416–34. PubMed Abstract
| Publisher Full Text
| Free Full Text
- Zhou R, Xiao X, He P, et al.:
SCAPE: a mixture model revealing single-cell polyadenylation diversity and cellular dynamics during cell differentiation and reprogramming.
Nucleic Acids Res.
2022; 50(11): e66. PubMed Abstract
| Publisher Full Text
| Free Full Text
