Discovery of alternative polyadenylation dynamics from single cell types.

Alternative polyadenylation (APA) occurs in the process of mRNA maturation by adding a poly(A) tail at different locations, resulting increased diversity of mRNA isoforms and contributing to the complexity of gene regulatory network. Benefit from the development of high-throughput sequencing technologies, we could now delineate APA profiles of transcriptomes at an unprecedented pace. Especially the single cell RNA sequencing (scRNA-seq) technologies provide us opportunities to interrogate biological details of diverse and rare cell types. Despite increasing evidence showing that APA is involved in the cell type-specific regulation and function, efficient and specific laboratory methods for capturing poly(A) sites at single cell resolution are underdeveloped to date. In this review, we summarize existing experimental and computational methods for the identification of APA dynamics from diverse single cell types. A future perspective is also provided.

Introduction

In eukaryotic cells, mature mRNA is generated from precursors that undergo multiple molecular processing, including 5′ end capping, splicing, and 3′ end cleavage and polyadenylation [1]. The location of the cleavage on the nascent RNA, followed by the synthesis of a poly(A) tail, is termed polyadenylation site or poly(A) site. It is well known that most protein-coding genes are transcribed into diverse mRNA isoforms through differential usage of different poly(A) sites. This phenomenon is known as alternative polyadenylation (APA) [2]. More than 50% of genes in animals and plants contain ≥2 poly(A) sites [3], [4], [5]. Most APAs occur in 3′ untranslated regions (UTR-APA), resulting in an inclusion or exclusion of specific sequences in 3′ UTRs. The 3′ UTRs tend to harbor many cis-elements, including poly(A) signals, mRNA stability regulatory elements, RNA export and localization signals, as well as miRNA targets [6], [7]. Thus, the dynamic usages of UTR-APA could influence mRNA stability, translational efficiency, nuclear export and cytoplasmic localization [6], [7]. In some cases, 3′ UTRs are involved in differentiation of protein functionalities through UTR-APA-mediated protein–protein interactions [8], [9], [10]. Moreover, some APAs are also found in 5′ UTR, intron and protein coding regions (CR-APA), resulting in truncated proteins or proteins with altered functions [2], [6], [11]. Inhibition of non-canonical poly(A) sites could protect mRNA integrity and globally regulate APA profiles [12], [13], and inactivation of tumor-suppressor-genes were widely observed in tumors through intronic polyadenylation [14].

APA is regulated by the polyadenylation machinery and associated proteins [15], [16], [17], [18], and is increasingly found to be tissue- and cell type-specific [19], [20], [21], [22], [23], [24], [25]. Through single-molecule sequencing, abundant tissue- and species-specific poly(A) sites were observed in maize and sorghum [19]. While APA in 3′ UTR of Pax3 was proved to control the fate of muscle stem cell and muscle function under homeostatic conditions [20], some differential APA events were found having less impact on muscle-specific expression profiles comparing to noncoding RNAs [21]. Differential APA responses to Cadmium toxicity were monitored across different root cell types in Arabidopsis [22]. The mouse brain cells showed substantial 3′ UTR dynamics across diverse cell types, especially, neurons globally preferred using distal, while microglia exhibited preferences of proximal, poly(A) sites [5], [23], [24]. Cell type specificities of APA profiles discovered from multiple tumor types show promising prognostic potential in certain cancers and cardiac diseases [26], [27], [28]. However, our understandings of APA regulation and function in different cell types is far from clear [2], [7]. Interrogation of APA profiles from different cells is the first step towards revealing such a mechanism.

To profile the cell type-specific APA or differentially expressed APA among diverse cell types, several kinds of wet lab approaches have been developed, and could be classified into three categories: (1) Sort and purify cells into different cell types before library construction and followed by a normal 3′ end sequencing [22], [29]; (2) Extract poly(A) profile of specific cell types from intact tissues based on co-immunoprecipitation (co-IP) of mRNA 3′ ends [23], [24]; (3) Integration of conventional bulk 3′ end sequencing and single cell RNA-seq protocols to interrogate the cell types and poly(A) profile simultaneously [30]. Meanwhile, in order to take full advantage of massive volume of existing single cell RNA-seq (scRNA-seq) data, most of which adopted a 3′ end-enriched sequencing strategy and were mainly used to delineate the cell heterogeneity based on gene expression, several computational tools (scDAPA [31], scAPA [32], and Sierra [33]) were recently developed to decipher the cell heterogeneity at the APA level and demonstrated promising efficacy.

In this review, we summarize these experimental methods detecting APAs from individual cell types, and computational programs delineating APA sites and dynamics from 3′ end-enriched scRNA-seq data.

Experimental methods for cell type-specific APA identification

To comprehensively and precisely obtain APA profiles of different cell types, the best way is designing appropriately experimental approaches to differentiate cell populations and capture mRNA fragments with poly(A) tails directly, followed by library construction and high-throughput sequencing. The currently available experimental methods devoting in cell type APA identification could be classified into the following three categories as shown in Fig. 1.

Fig. 1

Experimental protocols for profiling APA in cell type-specific manner. (A) A cell sorting method using PAT-seq; (B) The cTag-PAPERCLIP using crosslinking immunoprecipitation and GFP tagging; (C) The BAT-seq using cellular and molecular barcodes. See more details in the text.

Cell sorting methods

These methods first dissociate complex tissues into individual cells which is then sorted into pure cell populations before performing 3′ end sequencing [22], [29]. Tissues mincing could be finished by mechanical disaggregation (scalpels or tweezers) with or without enzymatic dissociation by collagenase (Col II, Col IV, Col V, or Col XI) plus DNase [34]. For plants, it requires chopping of tissues into small pieces, followed by a cellulase cocktail digestion [22]. After that, individual cells could be sorted into pure cell populations using fluorescence-activated cell sorting (FACS) [35], [36], micromanipulation (single cells are individually aspirated from a cell population under a microscope using glass micropipettes) [37], optical tweezers (integration of imaging based cell selection and laser beam to manipulate single cells) [38], or microfluidics [39] etc. Total mRNAs of individual cell types are then extracted for library construction and sequencing using a NGS platform, e.g. Illumina or Ion Torrent.

Based on the adopted library construction strategies, 3′ end sequencing protocols could be classified into two categories: oligo(dT) primer-based protocols and RNA manipulation-based protocols. Oligo(dT) primer-based protocols mainly utilize oligo(dT) primer to directly capture poly(A) fragments of mRNA molecules, followed by reverse transcription to produce cDNA for sequencing, such as 3′-seq [29], A-seq [40], PAS-seq [41], and PAT-seq [42] (Fig. 1A). These protocols are easy to implement and time saving, but typically accompany with internal priming (annealing to internal adenine-rich sequences), sequencing desynchronization, and bias arising from enzymatic digestions. To overcome the internal priming problem, several RNA manipulation-based protocols were developed by incorporating multiple RNA manipulation steps. For example, 3P-seq adopts a splint adaptor to ligate complete mRNAs, and RNase T1/H digestions to obtain target fragments for sequencing [43]; 3′READS employs chimeric U5 and T45 and a stringent primer washing in library construction [44]. However, they are generally time consuming, involving multiple steps of RNA manipulation, and prone to bias of the cleavage efficiency of the enzymes.

Limitations of cell sorting methods include: (1) potential introduction of unexpected cellular stress during tissue dissociation, cell type purification and laser illumination, thus potentially alter the gene expression and polyadenylation profiles; (2) the need of a high volume sample; (3) not applicable to tissues hard to be dissociated or containing rare cells that are not efficiently collected during sorting.

Crosslinking immunoprecipitation and GFP tagging methods

To address the impact of manipulation and/or cell sorting on gene expression, Hwang et al., proposed cTag-PAPERCLIP based on crosslinking immunoprecipitation (CLIP) and green fluorescent protein (GFP) tagging to extract poly(A) profile of pure cell populations from intact tissues [23], [24]. Taking advantage that poly(A)-binding protein (PABP) has high affinity to poly(A) tails (not internal adenine-rich sequences) [45], PABP was used to pull down the poly(A) tail contained mRNA fragments before performing reverse transcription reactions with oligo(dT) [46] (Fig. 1B). This approach also eliminated potential internal priming issues resulting from oligo(dT) primers annealing to internal adenine-rich sequences in mRNAs. In order to profile specific cell type in mixed tissues, mouse models conditionally and selectively expressing GFP-tagged poly(A)-binding proteins were generated using the Cre/Lox technology [47] (Fig. 1B).

The technical advantages of cTag-PAPERCLIP include less extracellular stress in tissue dissociation compared to manual manipulation and FACS, and the elimination of potential internal priming events without in silico filtering. Limitations include that cTag-PAPERCLIP cannot handle with unknown or rare cell populations, and the organism must express tagged-PABP in a cell-type specific manner. Since PABP plays multiple key roles in regulating mRNA localization, turnover and translation, it is unclear how GFP-tag may negatively impact the folding structures and many essential functions of PABP, thus altering cellular physiology and state [48]. Additionally, some RNA-binding proteins are aggregation-prone (like PABPN1), overexpression may result in aggregation and depletion of themselves (hence bias for APA), and the aggregation is also cell-type specific [17]. Moreover, it is time consuming for plants and/or those organisms without an efficient transformation protocol. Similar to cell sorting methods, it is also not applicable to cases where sample amounts are limited.

Cellular and molecular barcoding methods

The rapid growth of single cell RNA-seq (scRNA-seq) technology provides us unprecedented opportunity to delineate the cell state and diversity, and changes our understanding of organisms and organs by tracing the widespread cell heterogeneity [49], [50]. It’s conceivable to incorporate the scRNA-seq protocol into conventional 3′ end sequencing protocols. Velten et al. proposed a method BAT-seq [30] by combining a bulk 3′ end sequencing protocol TIF-seq [51] and a cellular and molecular barcode-based scRNA-seq protocol [52]. In these protocols, cellular barcodes were introduced to differentiate cell identities of mRNA molecules, and unique molecular identifiers (UMI) were added for absolute quantification of gene expression [52], [53]. Cellular and UMI barcoding could be accomplished by reverse transcription through plate-based (STRT, CEL-seq2) or droplet-based protocols (Drop-seq, MARS-seq) [54], [55]. Reverse transcription took place inside each droplet/well with a single cell, after which cells were pooled in bulk for initial amplification of cDNAs. Fragments of mRNA with poly(A) tails were captured using a biotinylated tag, followed by 3′ end library construction and Illumina paired-end sequencing (Fig. 1C).

Undoubtedly, the marriage of scRNA-seq and 3′ end sequencing is the ideal way to delineate APA dynamics from diverse cell types. Unknown and rare cells or sub-type cell populations could be detected based on gene expression profiles and cell-type level gene markers. Furthermore, RNA molecules and usages of distinct poly(A) sites could be precisely quantified by the embedded cellular and UMI barcodes, eliminating potential biases introduced by unbalanced PCR amplification efficiencies and/or diverse mRNA fragment lengths [52], [53]. However, the single cell-based method BAT-seq encompasses several limitations, including a low sensitivity (limited detection of ~5% of RNA molecules), complicated manipulation and computation steps [6].

Computational methods for cell type-specific APA identification

Although a serial of experimental methods were developed to directly profile the APA of multiple species and tissues in the past decade, the availability of 3′ end sequencing data and their coverage of different cell types are still limited and incomparable to the conventional bulk RNA-seq data [6], [56]. To remedy the lag of APA studies and take advantage of the massive volume of existing RNA-seq data, a number of computational methods have been developed to identify and quantify APA from conventional bulk RNA-seq data (reviewed in [56], [57]). In terms of sequencing strategies, reads from conventional bulk RNA-seq data would be distributed evenly across individual transcripts, fluctuations of read coverages resulted from diverse isoforms (including APA) could be detected by specific computational approaches, such as DaPars [27], APAtrap[58], ChangePoint [59], IntMAP [60], 3USS [61], and TECtool [62].

Similar to the bulk RNA-seq, massive scRNA-seq data are continuously produced in an unprecedented speed and coverage. Since efficacious single cell-based 3′ end sequencing protocol is still lacking, in order to explore the APA dynamics between different cell types or biological conditions, several computational methods were proposed using scRNA-seq data. According to their adopted quantification strategies, these methods could be clustered into two categories: peak calling-based methods and density distribution-based methods (Table 1).

Table 1

Characteristics of currently available computational methods for detecting APA dynamics from scRNA-seq data.

	Strategy	Quantificationmethod	Statistical method	scRNA-seqdata type
scAPA	Peak calling-based	log_2C_1+1〈C+1〉	Chi-squared test	3′ end
Sierra	Peak calling-based	log_2E×p_11⋯p_n1g_1+1+1	Wilcoxon rank-sum test	3′ end
scDAPA	Density distribution-based	∑n=1NpnA-pnB2	Wilcoxon rank-sum test	3′ end,full-length

Peak calling-based methods

Unlike conventional bulk RNA-seq, most of scRNA-seq methods introduces cellular and molecular barcodes into library constructions, reads from a specific cell and molecule will be uniquely barcoded and be unbiasedly quantified [54], [55]. Thus, typically only the 3′ end or 5′ end most fragments will be sequenced in scRNA-seq protocols, yielding sequenced reads that are enriched in locations around the transcription start or polyadenylation sites of individual isoforms in high probability. Taking it into account, scAPA [32] and Sierra [33] were proposed recently to quantify the cell type-specific APA regulation from 3′ enriched scRNA-seq data by a peak calling strategy (Table 1, Fig. 2).

Fig. 2

Illustration of the peak calling-based and density distribution-based methods in APA dynamics identification using 3′ enriched scRNA-seq data. Low coverage peaks are missed by peak calling-based methods, and overlapping peaks resulted from usage of adjacent poly(A) sites cannot be separated (top panel); usage of adjacent poly(A) sites are divided into separate bins by density distribution-based methods, while concrete number of poly(A) sites cannot be determined (bottom panel).

These methods generally take the aligned-BAM file and the cell cluster annotation file as inputs. Aligned BAM file and cell cluster annotation are obtained through routine scRNA-seq analysis pipelines and tools (e.g. 10× Cell Ranger [53], STAR [63], Seurat [64], and SCran [65]). Key steps of these methods include:

Divide the mixed BAM files into individual BAM files of reads from a single cell type. Assignment and deduplication of reads could be done according to their attached cellular and UMI barcodes [66], and cell type classification results obtained from analyzing the homogeneity of gene expression profiles of cells [64], [65].

Identify alignment peaks from aligned-BAM files. Since the result of read alignments of scRNA-seq data is similar to that of ChIP-seq data, the enriched and accumulated alignments are shown as lots of individual peaks (Fig. 2). Generally, an individual peak presents a preference of a specific poly(A) site of the gene, multiple peaks of a gene indicates the occurrence of APA in the corresponding gene. To conduct the peak identification, scAPA utilizes Homer fndPeaks [67], a popular peak calling tool of ChIP-seq data analysis, to locate the enriched regions. On the other hand, Sierra employs a splice-aware peak caller to overcome the potential spliced-peaks caused by introns [33].

After peak identification, scAPA utilizes the BEDTools merge [68] to combine overlapping peaks, followed by a Gaussian finite mixture modeling [69] to separate potential adjacent peaks caused by the fuzzy distribution of aligned-reads. In contrast to scAPA, Sierra takes a similarity score, considering both the distance and lengths of two peaks, to determine whether the two peaks should be merged or not.

Annotate and quantify individual peaks. All peaks called by scAPA are annotated with a predefined 3′ UTR extracted from human and mouse genomes. Numbers of unique reads located in annotated peaks of each cell type are then counted by featureCounts of the Rsubread package [70]. While Sierra takes the GTF annotation file and genomic coordinates of peaks as inputs to annotate the genomic features of identified peaks, and quantify the UMI number of each peak using the GenomicAlignments package [71].

Detect APA dynamics from cell type-specific APA profiles. Several distinct approaches were proposed in detecting the differential usage of APA in genes between samples in cell type-specific manner. In scAPA, the chi-squared test is used to evaluate the significance of difference, and a proximal poly(A) site usage index (proximal PUI, Table 1) was proposed to quantify the relative usage of the most proximal poly(A) site within a gene. Sierra employs the Wilcoxon rank-sum test to detect the 3′ UTR lengthening or shortening events, and calculates a relative expression level (R, Table 1) for each peak in each cell.

Peak calling method is not sensitive to overlapping peaks resulting from reads of isoforms with adjacent poly(A) sites (Fig. 2) [72], [73], only ~5% of poly(A) sites adjacent to each other within 200 bp could be differentiated by scAPA. For distance ranges in 200–300 bp, the recognition rate is ~30% [32]. Reads spreading in a wide range and consequently a low coverage will be missed by peak calling methods in high probability. Moreover, the high dropout rate and low sequencing throughput in single cells of currently available scRNA-seq technologies also hamper the discovery precision and depth in identification of APA dynamics. Thus, not all valid reads could be recovered or accurately split by the peak calling methods, and not all APA events could be discovered from these scRNA-seq data.

Density distribution-based methods

Through the reads generated by 3′ enriched scRNA-seq enrich in the upstream regions of poly(A) sites, their distributions generally cover a rather wide region and in a fuzzy mode [33], [74]. To overcome the fuzzy distribution, scDAPA utilizes a non-parametric density distribution method to quantify the difference between two compared groups, rather than calling the precise locations and boundaries of peaks [31]. The scDAPA also takes aligned-BAM/SAM file and cell cluster annotation file as inputs, and consists of three core modules:

Extract and divide reads by cell type classifications. The mapped reads should be in BAM/SAM format, and reads are well tagged with cell barcodes and UMIs. The scDAPA employs an automatic shell script to deduplicate and extract valid reads based on the 10× Cell Ranger BAM Tag system (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam) and the default alignment Tag system [75]. Valid reads will be separated into individual SAM files according to the cell type classification result including cell barcodes and their corresponding cell types/groups.

Extract and annotate 3′ ends from aforementioned reads. The coordinates of the 3′ most ends of these retrieved reads are extracted using BEDTools [68] and annotated by genes with a provided GTF file. All the reads overlapping with a genic region (including both exonic and intronic regions) will be reserved, others located in intergenic regions will be discarded.

Identify APA dynamics between two different cell types or the same cell type in different conditions. Instead of assigning 3′ enriched reads into separate peaks, scDAPA divides the genic region into equal-size bins (histograms), calculate the density distribution of each bins by gene and by cell type, and directly quantify the difference using a site distribution difference (SDD, Table 1) index (Fig. 2). The significance of difference is measured by the Wilcoxon rank-sum test, and adjusted by the Benjamini Hochberg method if multiple statistical tests were performed [76].

Since the performance of histogram quantification is sensitive to the bin size [77], scDAPA chooses the 100 bp as default bin-size through evaluating a sequential of bin-sizes and inspecting the distances between 3′ ends of scRNA-seq data and authentic poly(A) site annotations [31], [74]. Furthermore, the scDAPA is also applicable to scRNA-seq data sequencing full-length of mRNA molecules, in terms of only 3′ most ends are used in density quantification [54]. The limitation of density distribution-based method is that it quantifies relative difference between groups, rather than quantifies the usage of individual poly(A) sites. Besides, it also faces the innate drawback of currently available scRNA-seq technologies, a low gene capturing ratio, and many APA events cannot be detected. Thus, the identified APA events cannot fully represent the APA landscape of the cell. Before better technologies become available, other means to improve data processing (e.g. [79]) could be implemented to reach this goal.

Summary and outlook

APA is increasingly recognized as an important regulator in many processes, including mRNA translation, stability, nuclear export and localization [2], [7]. Tissue- and cell type- specific regulation of APA were also widely observed and involved in cell activation, proliferation, development and oncogenesis [23], [24], [56], [78], implying a complex role of APA. To decipher and understand the regulatory mechanism of APA, it’s necessary to investigate the APA atlas at a cellular resolution.

Here, we summarize the recent development of experimental methods in discovery of APA dynamics in cell type-specific manner, as well as computation method utilizing the scRNA-seq data. The available experimental protocols have largely advanced the cell type-specific study of APA, however, these methods still contain many inherent drawbacks as aforementioned. The cell sorting method may introduce unexpected stress in sample preparation, and affect the cellular states and APA profiles consequently; the crosslinking immunoprecipitation and GFP tagging method requires designing specify transgene model expressing GFP-tagged PABP in specific cell types, that limits its wide application in various cell populations; the cellular and molecular barcoding method could quantify gene expression and APA profiles from single cells by employing a scRNA-seq protocol, allowing discerning rare and unknown cell populations possible. Collectively, the single cell-based 3′ end sequencing method is of the most promising application prospect. However, technical ameliorations of this type of combinatory protocol, such as consideration of a more powerful and stable single cell platform 10× Chromium, to generate an easy-to-implement protocol of high accuracy and robustness are pressing needs.

In the past decade, a massive volume of scRNA-seq data, covering a wide variety of cell types and biological conditions, were generated. Computational approaches mining 3′ end enriched scRNA-seq data expand their application scope, and largely fill the gap in study of cell type-specific APA. Theses computational methods have been successfully applied on several 3′ enriched scRNA-seq data sets, e.g., a global 3′ UTR shortening were observed in activated T cells compared to naïve T cells using scAPA [32], and aberrant APA dynamics were observed in lung cancer and acute myeloid leukemia by scAPA and scDAPA [32], [74], respectively. However, the low resolution of 3′ enriched scRNA-seq in depicting poly(A) sites hampers the peak calling-based methods in accurate separation of peaks generated by close poly(A) sites; density distribution-based method cannot quantify usages of individual poly(A) site. The inherent low gene capturing ratio of scRNA-seq protocols also limits the observation depth of these computational methods. To improve the discovery accuracy and reduce the false discovery rate, potential ways include integrating the reference annotation of poly(A) sites, such as PolyA_DB 3 [3] and PlantAPAdb [4], to guide the identification and quantification of poly(A) site usages from scRNA-seq data, and improving the sequencing throughput of single cells and/or the number of cell inputs. All in all, more innovatively and efficaciously computational approaches for estimating cell type-specific APA usage from scRNA-seq data sets are still in high demand.

Declaration of competing interest

The authors report no conflict of interest.

Author contribution statement

Congting Ye: Conceptualization, Visualization, Writing - original draft, Funding acquisition. Juncheng Lin: Writing - review & editing. Qingshun Quinn Li: Conceptualization, Funding acquisition, Writing - review & editing.