Accurate identification of polyadenylation sites from 3' end deep sequencing using a naive Bayes classifier

Quantifying 3'UTR length from scRNA-seq data reveals changes independent of gene expression

Although more than half of all genes generate transcripts that differ in 3'UTR length, current analysis pipelines only quantify the amount but not the length of mRNA transcripts. 3'UTR length is determined by 3' end cleavage sites (CS). We map CS in more than 200 primary human and mouse cell types and increase CS annotations relative to the GENCODE database by 40%. Approximately half of all CS are used in few cell types, revealing that most genes only have one or two major 3' ends. We incorporate the CS annotations into a computational pipeline, called scUTRquant, for rapid, accurate, and simultaneous quantification of gene and 3'UTR isoform expression from single-cell RNA sequencing (scRNA-seq) data. When applying scUTRquant to data from 474 cell types and 2134 perturbations, we discover extensive 3'UTR length changes across cell types that are as widespread and coordinately regulated as gene expression changes but affect mostly different genes. Our data indicate that mRNA abundance and mRNA length are two largely independent axes of gene regulation that together determine the amount and spatial organization of protein synthesis.

Transcript shortening via alternative polyadenylation promotes gene expression during fracture healing

Maturation of the 3' end of almost all eukaryotic messenger RNAs (mRNAs) requires cleavage and polyadenylation. Most mammalian mRNAs are polyadenylated at different sites within the last exon, generating alternative polyadenylation (APA) isoforms that have the same coding region but distinct 3' untranslated regions (UTRs). The 3'UTR contains motifs that regulate mRNA metabolism; thus, changing the 3'UTR length via APA can significantly affect gene expression. Endochondral ossification is a central process in bone healing, but the impact of APA on gene expression during this process is unknown. Here, we report the widespread occurrence of APA, which impacts multiple pathways that are known to participate in bone healing. Importantly, the progression of endochondral ossification involves global 3'UTR shortening, which is coupled with an increased abundance of shortened transcripts relative to other transcripts; these results highlight the role of APA in promoting gene expression during endochondral bone formation. Our mechanistic studies of transcripts that undergo APA in the fracture callus revealed an intricate regulatory network in which APA enhances the expression of the collagen, type I, alpha 1 (Col1a1) and Col1a2 genes, which encode the 2 subunits of the abundantly expressed protein collagen 1. APA exerts this effect by shortening the 3'UTRs of the Col1a1 and Col1a2 mRNAs, thus removing the binding sites of miR-29a-3p, which would otherwise strongly promote the degradation of both transcripts. Taken together, our study is the first to characterize the crucial roles of APA in regulating the 3'UTR landscape and modulating gene expression during fracture healing.

Proposed criteria for nevoid basal cell carcinoma syndrome in children assessed using statistical optimization

Nevoid basal cell carcinoma syndrome (NBCCS) is a tumor predisposition condition, the cardinal features of which emerge in adolescence or adulthood. Using statistical optimization, this study proposes NBCCS criteria with improved sensitivity in children less than 18 years of age. Earlier detection may lead to improved surveillance and prevention of sequelae. A survey eliciting medical history was completed by, or on behalf of, individuals with NBCCS. Based on these findings, criteria for suspicion of NBCCS in children were suggested using information from a Bernoulli naïve Bayes classifier relying on the human phenotype ontology. The sensitivity and specificity of the existing and proposed diagnostic criteria were also assessed. Participants (n = 48) reported their first signs of NBCCS appeared at a median age of 8 months, but by our retrospective analysis, they did not fulfill the current diagnostic criteria until a median age of 7 years. This study delineates the early-onset features of NBCCS and proposes criteria that should prompt consideration of NBCCS. Additionally, we demonstrate a method for quantitatively assessing the utility of diagnostic criteria for genetic disorders.

Polyadenylation-related isoform switching in human evolution revealed by full-length transcript structure

Rhesus macaque is a unique nonhuman primate model for human evolutionary and translational study, but the error-prone gene models critically limit its applications. Here, we de novo defined full-length macaque gene models based on single molecule, long-read transcriptome sequencing in four macaque tissues (frontal cortex, cerebellum, heart and testis). Overall, 8 588 227 poly(A)-bearing complementary DNA reads with a mean length of 14 106 nt were generated to compile the backbone of macaque transcripts, with the fine-scale structures further refined by RNA sequencing and cap analysis gene expression sequencing data. In total, 51 605 macaque gene models were accurately defined, covering 89.7% of macaque or 75.7% of human orthologous genes. Based on the full-length gene models, we performed a human–macaque comparative analysis on polyadenylation (PA) regulation. Using macaque and mouse as outgroup species, we identified 79 distal PA events newly originated in humans and found that the strengthening of the distal PA sites, rather than the weakening of the proximal sites, predominantly contributes to the origination of these human-specific isoforms. Notably, these isoforms are selectively constrained in general and contribute to the temporospatially specific reduction of gene expression, through the tinkering of previously existed mechanisms of nuclear retention and microRNA (miRNA) regulation. Overall, the protocol and resource highlight the application of bioinformatics in integrating multilayer genomics data to provide an intact reference for model animal studies, and the isoform switching detected may constitute a hitherto underestimated regulatory layer in shaping the human-specific transcriptome and phenotypic changes.

Divergence in alternative polyadenylation contributes to gene regulatory differences between humans and chimpanzees

While comparative functional genomic studies have shown that inter-species differences in gene expression can be explained by corresponding inter-species differences in genetic and epigenetic regulatory mechanisms, co-transcriptional mechanisms, such as alternative polyadenylation (APA), have received little attention. We characterized APA in lymphoblastoid cell lines from six humans and six chimpanzees by identifying and estimating the usage for 44,432 polyadenylation sites (PAS) in 9518 genes. Although APA is largely conserved, 1705 genes showed significantly different PAS usage (FDR 0.05) between species. Genes with divergent APA also tend to be differentially expressed, are enriched among genes showing differences in protein translation, and can explain a subset of observed inter-species protein expression differences that do not differ at the transcript level. Finally, we found that genes with a dominant PAS, which is used more often than other PAS, are particularly enriched for differentially expressed genes.

Alternative polyadenylation: methods, mechanism, function, and role in cancer

Occurring in over 60% of human genes, alternative polyadenylation (APA) results in numerous transcripts with differing 3'ends, thus greatly expanding the diversity of mRNAs and of proteins derived from a single gene. As a key molecular mechanism, APA is involved in various gene regulation steps including mRNA maturation, mRNA stability, cellular RNA decay, and protein diversification. APA is frequently dysregulated in cancers leading to changes in oncogenes and tumor suppressor gene expressions. Recent studies have revealed various APA regulatory mechanisms that promote the development and progression of a number of human diseases, including cancer. Here, we provide an overview of four types of APA and their impacts on gene regulation. We focus particularly on the interaction of APA with microRNAs, RNA binding proteins and other related factors, the core pre-mRNA 3'end processing complex, and 3'UTR length change. We also describe next-generation sequencing methods and computational tools for use in poly(A) signal detection and APA repositories and databases. Finally, we summarize the current understanding of APA in cancer and provide our vision for future APA related research.

Alternative polyadenylation mediates genetic regulation of gene expression

Little is known about co-transcriptional or post-transcriptional regulatory mechanisms linking noncoding variation to variation in organismal traits. To begin addressing this gap, we used 3’ Seq to study the impact of genetic variation on alternative polyadenylation (APA) in the nuclear and total mRNA fractions of 52 HapMap Yoruba human lymphoblastoid cell lines. We mapped 602 APA quantitative trait loci (apaQTLs) at 10% FDR, of which 152 were nuclear specific. Effect sizes at intronic apaQTLs are negatively correlated with eQTL effect sizes. These observations suggest genetic variants can decrease mRNA expression levels by increasing usage of intronic PAS. We also identified 24 apaQTLs associated with protein levels, but not mRNA expression. Finally, we found that 19% of apaQTLs can be associated with disease. Thus, our work demonstrates that APA links genetic variation to variation in gene expression, protein expression, and disease risk, and reveals uncharted modes of genetic regulation.

Evolutionarily conserved pachytene piRNA loci are highly divergent among modern humans

In the fetal mouse testis, PIWI-interacting RNAs (piRNAs) guide PIWI proteins to silence transposons but, after birth, most post-pubertal pachytene piRNAs map to the genome uniquely and are thought to regulate genes required for male fertility. In the human male, the developmental classes, precise genomic origins and transcriptional regulation of postnatal piRNAs remain undefined. Here, we demarcate the genes and transcripts that produce postnatal piRNAs in human juvenile and adult testes. As in the mouse, human A-MYB drives transcription of both pachytene piRNA precursor transcripts and messenger RNAs encoding piRNA biogenesis factors. Although human piRNA genes are syntenic to those in other placental mammals, their sequences are poorly conserved. In fact, pachytene piRNA loci are rapidly diverging even among modern humans. Our findings suggest that, during mammalian evolution, pachytene piRNA genes are under few selective constraints. We speculate that pachytene piRNA diversity may provide a hitherto unrecognized driver of reproductive isolation.

Free circular introns with an unusual branchpoint in neuronal projections

The polarized structure of axons and dendrites in neuronal cells depends in part on RNA localization. Previous studies have looked at which polyadenylated RNAs are enriched in neuronal projections or at synapses, but less is known about the distribution of non-adenylated RNAs. By physically dissecting projections from cell bodies of primary rat hippocampal neurons and sequencing total RNA, we found an unexpected set of free circular introns with a non-canonical branchpoint enriched in neuronal projections. These introns appear to be tailless lariats that escape debranching. They lack ribosome occupancy, sequence conservation, and known localization signals, and their function, if any, is not known. Nonetheless, their enrichment in projections has important implications for our understanding of the mechanisms by which RNAs reach distal compartments of asymmetric cells.

Detection of Differentially Expressed Cleavage Site Intervals Within 3' Untranslated Regions Using CSI-UTR Reveals Regulated Interaction Motifs

The length of untranslated regions at the 3' end of transcripts (3'UTRs) is regulated by alternate polyadenylation (APA). 3'UTRs contain regions that harbor binding motifs for regulatory molecules. However, the mechanisms that coordinate the 3'UTR length of specific groups of transcripts are not well-understood. We therefore developed a method, CSI-UTR, that models 3'UTR structure as tandem segments between functional alternative-polyadenylation sites (termed cleavage site intervals-CSIs). This approach facilitated (1) profiling of 3'UTR isoform expression changes and (2) statistical enrichment of putative regulatory motifs. CSI-UTR analysis is UTR-annotation independent and can interrogate legacy data generated from standard RNA-Seq libraries. CSI-UTR identified a set of CSIs in human and rodent transcriptomes. Analysis of RNA-Seq datasets from neural tissue identified differential expression events within 3'UTRs not detected by standard gene-based differential expression analyses. Further, in many instances 3'UTR and CDS from the same gene were regulated differently. This modulation of motifs for RNA-interacting molecules with potential condition-dependent and tissue-specific RNA binding partners near the polyA signal and CSI junction may play a mechanistic role in the specificity of alternative polyadenylation. Source code, CSI BED files and example datasets are available at: https://github.com/UofLBioinformatics/CSI-UTR.

Alterations in mRNA 3' UTR Isoform Abundance Accompany Gene Expression Changes in Human Huntington's Disease Brains

The huntingtin gene has two mRNA isoforms that differ in their 3' UTR length. The relationship of these isoforms with Huntington's disease is not established. We provide evidence that the abundance of huntingtin 3' UTR isoforms differs between patient and control neural stem cells, fibroblasts, motor cortex, and cerebellum. Huntingtin 3' UTR isoforms, including a mid-3' UTR isoform, have different localizations, half-lives, polyA tail lengths, microRNA sites, and RNA-binding protein sites. Isoform shifts in Huntington's disease motor cortex are not limited to huntingtin; 11% of alternatively polyadenylated genes change the abundance of their 3' UTR isoforms. Altered expression of RNA-binding proteins may be associated with aberrant isoform abundance; knockdown of the RNA-binding protein CNOT6 in control fibroblasts leads to huntingtin isoform differences similar to those in disease fibroblasts. These findings demonstrate that mRNA 3' UTR isoform changes are a feature of molecular pathology in the Huntington's disease brain.

A high-resolution mRNA expression time course of embryonic development in zebrafish

We have produced an mRNA expression time course of zebrafish development across 18 time points from 1 cell to 5 days post-fertilisation sampling individual and pools of embryos. Using poly(A) pulldown stranded RNA-seq and a 3′ end transcript counting method we characterise temporal expression profiles of 23,642 genes. We identify temporal and functional transcript co-variance that associates 5024 unnamed genes with distinct developmental time points. Specifically, a class of over 100 previously uncharacterised zinc finger domain containing genes, located on the long arm of chromosome 4, is expressed in a sharp peak during zygotic genome activation. In addition, the data reveal new genes and transcripts, differential use of exons and previously unidentified 3′ ends across development, new primary microRNAs and temporal divergence of gene paralogues generated in the teleost genome duplication. To make this dataset a useful baseline reference, the data can be browsed and downloaded at Expression Atlas and Ensembl.

Defining the 5΄ and 3΄ landscape of the Drosophila transcriptome with Exo-seq and RNaseH-seq

Cells regulate biological responses in part through changes in transcription start sites (TSS) or cleavage and polyadenylation sites (PAS). To fully understand gene regulatory networks, it is therefore critical to accurately annotate cell type-specific TSS and PAS. Here we present a simple and straightforward approach for genome-wide annotation of 5΄- and 3΄-RNA ends. Our approach reliably discerns bona fide PAS from false PAS that arise due to internal poly(A) tracts, a common problem with current PAS annotation methods. We applied our methodology to study the impact of temperature on the Drosophila melanogaster head transcriptome. We found hundreds of previously unidentified TSS and PAS which revealed two interesting phenomena: first, genes with multiple PASs tend to harbor a motif near the most proximal PAS, which likely represents a new cleavage and polyadenylation signal. Second, motif analysis of promoters of genes affected by temperature suggested that boundary element association factor of 32 kDa (BEAF-32) and DREF mediates a transcriptional program at warm temperatures, a result we validated in a fly line where beaf-32 is downregulated. These results demonstrate the utility of a high-throughput platform for complete experimental and computational analysis of mRNA-ends to improve gene annotation.

FLT1 and transcriptome-wide polyadenylation site (PAS) analysis in preeclampsia

Maternal symptoms of preeclampsia (PE) are primarily driven by excess anti-angiogenic factors originating from the placenta. Chief among these are soluble Flt1 proteins (sFlt1s) produced from alternatively polyadenylated mRNA isoforms. Here we used polyadenylation site sequencing (PAS-Seq) of RNA from normal and PE human placentae to interrogate transcriptome-wide gene expression and alternative polyadenylation signatures associated with early-onset PE (EO-PE; symptom onset < 34 weeks) and late-onset PE (LO-PE; symptom onset > 34 weeks) cohorts. While we observed no general shift in alternative polyadenylation associated with PE, the EO-PE and LO-PE cohorts do exhibit gene expression profiles distinct from both each other and from normal placentae. The only two genes upregulated across all transcriptome-wide PE analyses to date (microarray, RNA-Seq and PAS-Seq) are NRIP1 (RIP140), a transcriptional co-regulator linked to metabolic syndromes associated with obesity, and Flt1. Consistent with sFlt1 overproduction being a significant driver of clinical symptoms, placental Flt1 mRNA levels strongly correlate with maternal blood pressure. For Flt1, just three mRNA isoforms account for > 94% of all transcripts, with increased transcription of the entire locus driving Flt1 upregulation in both EO-PE and LO-PE. These three isoforms thus represent potential targets for therapeutic RNA interference (RNAi) in both early and late presentations.

Review of current methods, applications, and data management for the bioinformatics analysis of whole exome sequencing

The advent of next-generation sequencing technologies has greatly promoted advances in the study of human diseases at the genomic, transcriptomic, and epigenetic levels. Exome sequencing, where the coding region of the genome is captured and sequenced at a deep level, has proven to be a cost-effective method to detect disease-causing variants and discover gene targets. In this review, we outline the general framework of whole exome sequence data analysis. We focus on established bioinformatics tools and applications that support five analytical steps: raw data quality assessment, pre-processing, alignment, post-processing, and variant analysis (detection, annotation, and prioritization). We evaluate the performance of open-source alignment programs and variant calling tools using simulated and benchmark datasets, and highlight the challenges posed by the lack of concordance among variant detection tools. Based on these results, we recommend adopting multiple tools and resources to reduce false positives and increase the sensitivity of variant calling. In addition, we briefly discuss the current status and solutions for big data management, analysis, and summarization in the field of bioinformatics.