Identifying Active and Alternative Promoters from RNA-Seq data with proActiv

Joseph Lee, Deniz Demircioğlu, Jonathan Göke

Summary

Most human genes have multiple promoters that control the expression of distinct isoforms. The use of these alternative promoters enables the regulation of isoform expression pre-transcriptionally. Alternative promoters have been found to be important in a wide number of cell types and diseases.

proActiv is a method that enables the analysis of promoters from RNA-Seq data. proActiv uses aligned reads as input, and then generates counts and normalized promoter activity estimates for each annotated promoter. These estimates can then be used to identify which promoter is active, which promoter is inactive, and which promoters change their activity across conditions.

Here we present a quick start guide to using proActiv, and a detailed workflow for identifying active promoters and alternative promoters across 2 conditions.

If you use proActiv in your research, please cite:

Demircioğlu, Deniz, et al. “A Pan-cancer Transcriptome Analysis Reveals Pervasive Regulation through Alternative Promoters.” Cell 178.6 (2019): 1465-1477.

Contents

Quick Start: Quantifying promoter activity with proActiv

proActiv estimates promoter activity from RNA-Seq data. Promoter activity is defined as the total amount of transcription initiated at each promoter. proActiv takes as input either BAM files or junction files (TopHat2 or STAR), and a promoter annotation object of the relevant genome. An optional argument condition can be supplied, describing the condition corresponding to each input file. Here we demonstrate proActiv with STAR junction files (Human genome GRCh38 Gencode v34) as input. Due to size constraints, analysis is restricted to a subset of chr1 (10,000,000-30,000,000).

library(proActiv)

## List of STAR junction files as input
files <- list.files(system.file('extdata/vignette/junctions', 
                                package = 'proActiv'), full.names = TRUE)
## Vector describing experimental condition
condition <- rep(c('A549','HepG2'), each=3)
## Promoter annotation for human genome GENCODE v34 restricted to a subset of chr1
promoterAnnotation <- promoterAnnotation.gencode.v34.subset

result <- proActiv(files = files, 
                   promoterAnnotation = promoterAnnotation,
                   condition = condition)

result is a summarizedExperiment object which can be accessed as follows:

A complete workflow to identify alternative promoter usage

Here we present a complete step-by-step workflow for analyzing promoter activity with proActiv and for identifying alternative promoter usage across samples from different conditions. We will compare samples from 2 different cell lines (A549 and HepG2) to identify alternative promoters.

Preparing input data

proActiv uses RNA-Seq data to quantify promoter activity. Users have the option of using as input either BAM files, or tab-delimited junction files that are generated when TopHat2 or STAR is used for read alignment.

Below, we demonstrate running proActiv with input STAR junction files. This data is taken from the SGNEx project, and restricted to the chr1:10,000,000-30,000,000 region. The reference genome used for alignment is Gencode v34 (GRCh38). These files can be found in extdata/vignette/junctions:

Preparing promoter annotations

In order to quantify promoter activity, proActiv uses a set of promoters based on genome annotations. proActiv allows the creation of a promoter annotation object for any genome from a TxDb object or from a GTF file with the preparePromoterAnnotation function. Users have the option to either pass the file path of the GTF/GFF or TxDb to be used, or use the TxDb object directly as input. proActiv includes pre-calculated promoter annotations for the human genome (GENCODE v34). However, due to size constraints, the annotation is restricted to the chr1:10,000,000-30,000,000 region. Users can build full annotations by downloading GTF files from GENCODE page and following the steps below.

We demonstrate creating the restricted promoter annotation for the Human genome (GENCODE v34) with both GTF and TxDb:

The PromoterAnnotation object has 3 slots:

When creating promoter annotations for other species, one can call genomeStyles to identify the species argument to be used:

Running proActiv

Once promoters in the genome are identified, proActiv estimates promoter activity at each annotated promoter. Here, we load pre-calculated promoter annotation for GENCODE Release 34. We also supply the experimental condition to proActiv. This information allows proActiv to summarize results across conditions.

result is a SummarizedExperiment object with assays as raw/normalized promoter counts, absolute/relative promoter activity and gene expression:

The rowData slot stores a promoter-gene ID mapping and promoter position (5’ to 3’) for each promoter by gene. Mean absolute promoter activity and gene expression for each condition is also summarized here. Promoters are also categorized into three classes. Promoters with activity < 0.25 are classified as inactive, while the most active promoters of each gene are classified as major promoters. Promoters active at lower levels are classified as minor promoters.

promoterId geneId seqnames start strand internalPromoter promoterPosition txId A549.mean A549.gene.mean HepG2.mean HepG2.gene.mean A549.class HepG2.class
1 ENSG00000000938.13 chr1 27624062 - TRUE 4 ENST0000…. 0.000000 0.00000 0.000000 0.00000 NA NA
2 ENSG00000000938.13 chr1 27626240 - FALSE 3 ENST0000…. 0.000000 0.00000 0.000000 0.00000 Inactive Inactive
3 ENSG00000000938.13 chr1 27626569 - FALSE 2 ENST0000…. 0.000000 0.00000 0.000000 0.00000 Inactive Inactive
4 ENSG00000000938.13 chr1 27635185 - FALSE 1 ENST0000…. 0.000000 0.00000 0.000000 0.00000 Inactive Inactive
5 ENSG00000001460.18 chr1 24379823 - TRUE 7 ENST0000…. 4.721803 20.91493 4.319706 18.81489 NA NA
6 ENSG00000001460.18 chr1 24391679 - TRUE 6 ENST0000…. 4.207194 20.91493 3.999170 18.81489 NA NA

For cleaner downstream analysis, one can remove single-exon transcripts for which promoter activity is not quantified. result can be filtered as such:

Identifying alternative promoters

To identify alternative promoter usage, proActiv implements a linear model, regressing both absolute promoter activity and relative promoter activity against sample conditions. A change in absolute promoter activity indicates absolute differential usage of the promoter. However, this does not imply alternative usage or a switch in major promoter, as the relative usage of promoters between conditions may still be the same. For this reason, information from relative promoter activity must also be taken into account. Alternative promoter usage is thus characterized by a significant change in absolute and relative promoter activity across conditions.

Candidates for alternative promoter usage can be identified by calling getAlternativePromoters. This function takes in the result object returned by proActiv, and referenceCondition, the reference condition to be compared to all other samples.

Detection of alternative promoters with the following parameters:

Analysis and visualization of alternative promoter usage

Here we offer several visualizations of the data returned by the workflow above.

Alternative Promoter Usage

To visualize genes with alternative promoter usage across conditions, the function boxplotPromoters can be used. This function outputs boxplots of gene expression, absolute and relative promoter activity for a gene of interest. boxplotPromoters takes the following arguments:

Below, we call boxplotPromoters with RAP1GAP:

The plots object stores the plot for absolute promoter activity, relative promoter activity and gene expression in the first, second and third slot respectively. The plot generated by boxplotPromoters reflect the alternative usage of promoters 137 and 141 as called by getAlternativePromoters.

Promoter category proportion

Here, we visualize the categorization of annotated promoters in the two cell lines. The proportions between the categories are similar across the two cell lines, with majority of the promoters being inactive.

Major/minor promoters by position

Analysis of major:minor promoter proportions against promoter position. The analysis is restricted to multi-promoter genes with at least one active promoter. Below, we generate the plot for cell line HepG2. In general, the major:minor promoter proportion decreases with increasing promoter position.

Major promoter activity and gene expression

Comparison of major promoter activity and gene expression, calculated by summing over all promoters. Single promoter genes lie on the diagonal. Multi-promoter genes lie to the right of the diagonal. Below, we generate the plot for cell line HepG2. This plot suggests that a single major promoter does not often fully explain gene expression, with minor promoters also contributing to gene expression.

t-SNE

We generate a t-SNE plot with all active promoters. Expectedly, replicates from each cell line cluster together.

Getting help

Questions and issues can be raised at the Bioconductor support site: https://support.bioconductor.org. Ensure your posts are tagged with proActiv.

Alternatively, issues can be raised at the proActiv Github repository: https://github.com/GoekeLab/proActiv.

Citing proActiv

If you use proActiv, please cite:

Demircioğlu, Deniz, et al. “A Pan-cancer Transcriptome Analysis Reveals Pervasive Regulation through Alternative Promoters.” Cell 178.6 (2019): 1465-1477.

Session information

#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] Rtsne_0.16     ggplot2_3.3.6  gridExtra_2.3  proActiv_1.8.0
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7                matrixStats_0.62.0         
#>   [3] bit64_4.0.5                 filelock_1.0.2             
#>   [5] RColorBrewer_1.1-3          progress_1.2.2             
#>   [7] httr_1.4.4                  GenomeInfoDb_1.34.0        
#>   [9] tools_4.2.1                 bslib_0.4.0                
#>  [11] utf8_1.2.2                  R6_2.5.1                   
#>  [13] KernSmooth_2.23-20          DBI_1.1.3                  
#>  [15] BiocGenerics_0.44.0         colorspace_2.0-3           
#>  [17] withr_2.5.0                 tidyselect_1.2.0           
#>  [19] prettyunits_1.1.1           DESeq2_1.38.0              
#>  [21] bit_4.0.4                   curl_4.3.3                 
#>  [23] compiler_4.2.1              cli_3.4.1                  
#>  [25] Biobase_2.58.0              xml2_1.3.3                 
#>  [27] DelayedArray_0.24.0         labeling_0.4.2             
#>  [29] rtracklayer_1.58.0          sass_0.4.2                 
#>  [31] caTools_1.18.2              scales_1.2.1               
#>  [33] genefilter_1.80.0           rappdirs_0.3.3             
#>  [35] stringr_1.4.1               digest_0.6.30              
#>  [37] Rsamtools_2.14.0            rmarkdown_2.17             
#>  [39] XVector_0.38.0              pkgconfig_2.0.3            
#>  [41] htmltools_0.5.3             MatrixGenerics_1.10.0      
#>  [43] highr_0.9                   dbplyr_2.2.1               
#>  [45] fastmap_1.1.0               rlang_1.0.6                
#>  [47] RSQLite_2.2.18              farver_2.1.1               
#>  [49] BiocIO_1.8.0                jquerylib_0.1.4            
#>  [51] generics_0.1.3              jsonlite_1.8.3             
#>  [53] gtools_3.9.3                BiocParallel_1.32.0        
#>  [55] dplyr_1.0.10                RCurl_1.98-1.9             
#>  [57] magrittr_2.0.3              GenomeInfoDbData_1.2.9     
#>  [59] Matrix_1.5-1                Rcpp_1.0.9                 
#>  [61] munsell_0.5.0               S4Vectors_0.36.0           
#>  [63] fansi_1.0.3                 lifecycle_1.0.3            
#>  [65] stringi_1.7.8               yaml_2.3.6                 
#>  [67] SummarizedExperiment_1.28.0 zlibbioc_1.44.0            
#>  [69] gplots_3.1.3                BiocFileCache_2.6.0        
#>  [71] grid_4.2.1                  blob_1.2.3                 
#>  [73] parallel_4.2.1              crayon_1.5.2               
#>  [75] lattice_0.20-45             Biostrings_2.66.0          
#>  [77] splines_4.2.1               GenomicFeatures_1.50.0     
#>  [79] annotate_1.76.0             hms_1.1.2                  
#>  [81] KEGGREST_1.38.0             locfit_1.5-9.6             
#>  [83] knitr_1.40                  pillar_1.8.1               
#>  [85] GenomicRanges_1.50.0        rjson_0.2.21               
#>  [87] geneplotter_1.76.0          codetools_0.2-18           
#>  [89] biomaRt_2.54.0              stats4_4.2.1               
#>  [91] XML_3.99-0.12               glue_1.6.2                 
#>  [93] evaluate_0.17               data.table_1.14.4          
#>  [95] png_0.1-7                   vctrs_0.5.0                
#>  [97] gtable_0.3.1                assertthat_0.2.1           
#>  [99] cachem_1.0.6                xfun_0.34                  
#> [101] xtable_1.8-4                restfulr_0.0.15            
#> [103] survival_3.4-0              tibble_3.1.8               
#> [105] GenomicAlignments_1.34.0    AnnotationDbi_1.60.0       
#> [107] memoise_2.0.1               IRanges_2.32.0             
#> [109] ellipsis_0.3.2