PRAM: a novel pooling approach for discovering intergenic transcripts from large-scale RNA sequencing experiments

Abstract

Publicly available RNA-seq data is routinely used for retrospective analysis to elucidate new biology. Novel transcript discovery enabled by joint analysis of large collections of RNA-seq data sets has emerged as one such analysis. Current methods for transcript discovery rely on a '2-Step' approach where the first step encompasses building transcripts from individual data sets, followed by the second step that merges predicted transcripts across data sets. To increase the power of transcript discovery from large collections of RNA-seq data sets, we developed a novel '1-Step' approach named Pooling RNA-seq and Assembling Models (PRAM) that builds transcript models from pooled RNA-seq data sets. We demonstrate in a computational benchmark that 1-Step outperforms 2-Step approaches in predicting overall transcript structures and individual splice junctions, while performing competitively in detecting exonic nucleotides. Applying PRAM to 30 human ENCODE RNA-seq data sets identified unannotated transcripts with epigenetic and RAMPAGE signatures similar to those of recently annotated transcripts. In a case study, we discovered and experimentally validated new transcripts through the application of PRAM to mouse hematopoietic RNA-seq data sets. We uncovered new transcripts that share a differential expression pattern with a neighboring gene Pik3cg implicated in human hematopoietic phenotypes, and we provided evidence for the conservation of this relationship in human. PRAM is implemented as an R/Bioconductor package.

Figure 1.

1-Step outperforms 2-Step reconstruction methods. (A) Precision and recall of five meta-assembly methods in a benchmark test on target transcripts stratified by their maximum TPMs in the 30 ENCODE RNA-seq data sets: (1) TPM < 1 (413 transcripts); (2) 1 ≤ TPM < 10 (515 transcripts); and (3) TPM ≥ 10 (328 transcripts). (B) Comparison of target transcript GCM1 and predicted models by five meta-assembly methods.

Figure 2.

PRAM as a new computational framework predicts a valid master set of transcript models in human intergenic regions. (A) PRAM's workflow of input (cyan), intermediate (yellow), and output (green) files, with format labeled in brackets. PRAM's R functions and example parameters for each step are displayed next to arrows. (B) Distribution of GENCODE and PRAM transcripts in terms of expression levels across seven ENCODE cell lines. (C) PRAM transcript with the highest TPM had multiple complementary genomic features supporting its existence. The model ‘plcf_chr5_minus.2607.1’ had an average TPM of 245 in HeLa-S3 cells. It had high DNase-seq signals around its 5′ exon, suggesting high chromatin accessibility, and had multiple H3K4me3 ChIP-seq peaks, suggesting active transcription. Moreover, it had two RNA Pol II ChIP-seq peaks in close proximity to its transcription start site. All of these external genomic data supported the existence of this highly expressed PRAM transcript. (D,E) RAMPAGE (D) and histone modification ChIP-seq (E) signals of GENCODE and PRAM transcripts stratified by their expression levels together with ‘silent genomic regions’ defined based on H3K27me3 peaks as negative controls in all of GM12878 or K562's data sets. RAMPAGE and ChIP-seq values were derived from replicate 1 in their corresponding data sets (Supplemental Tables 13, 14). Transcripts with promoter or genomic span mappability <0.8 were excluded from D or E, respectively, due to uncertainty in their RAMPAGE or epigenetic signals. RAMPAGE and ChIP-seq signals were calculated as reads per million (RPM) and reads per kilobase per million (RPKM), respectively.

Figure 3.

Genomic features and experimental validations of PRAM mouse transcripts. (A) Workflow of applying PRAM to discover transcripts from mouse hematopoiesis-related RNA-seq data sets: input (cyan), intermediate results (yellow), and output (green). (B) PRAM transcripts CUFFp.chr12.15498 and CUFFm.chr12.33668 had multiple supporting genomic features from external data sets. (C) Semi-qRT-PCR measurements of the six PRAM models in untreated (Unt) and 48-h ß-estradiol (ß-est)-treated G1E-ER-GATA1 cells. Red dots demarcate anticipated transcript sizes. Isoforms with splice junctions distant from each other were measured separately. Gene model name prefixes were removed for brevity.

Figure 4.

Expression of PRAM transcripts correlate with the neighboring gene Pik3cg in mouse and human. (A) Expression levels of PRAM transcripts and their neighboring genes in untreated (Unt) and 48-h ß-estradiol (ß-est)-treated G1E-ER-GATA1 cells. CUFFp.chr12.15498's isoform 1 was not detected by semi-qRT-PCR and thus was not measured here. Two-tailed Student's t-test; (**) P-value < 0.01, (***) P-value < <0.001. (B) Fold changes of PRAM mouse transcripts and their neighboring genes in the RNA-seq data sets of untreated and 48-h ß-estradiol-reated G1E-ER-GATA1 cells (G1E) and wild type versus deletion of Gata2 +9.5 enhancer aorta-gonad-mesonephros (AGM). (C) Counterparts of PRAM mouse transcripts in human with their supporting genomic features. (D) Semi-qRT-PCR measurement of PRAM human transcripts and their neighboring genes. Gene model name prefixes were removed for brevity. (E,F) Correlation of gene expression levels between CUFFm.chr7.6148 with PIK3CG and PRKAR2B in K562 cells (E) and TCGA-LAML patients (F).