PRAM: Pooling RNA-seq and Assembling Models

Peng Liu, Colin N. Dewey, and Sündüz Keleş

2023-10-24

Introduction

Pooling RNA-seq and Assembling Models (PRAM) is an R package that utilizes multiple RNA-seq datasets to predict transcript models. The workflow of PRAM contains four steps. Figure 1 shows each step with function name and associated key parameters. In addition, we provide a function named evalModel() to evaluate prediction accuracy by comparing transcript models with true transcripts. In the later sections of this vignette, we will describe each function in details.

PRAM workflow

PRAM workflow

Installation

From GitHub

Start R and enter:

devtools::install_github('pliu55/pram')

From Bioconductor

Start R and enter:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("pram")

For different versions of R, please refer to the appropriate Bioconductor release.

Quick start

Description

PRAM provides a function named runPRAM() to let you conveniently run through the whole workflow.

For a given gene annotation and RNA-seq alignments, you can predict transcript models in intergenic genomic regions:

##
## assuming the stringtie binary is in folder /usr/local/stringtie-1.3.3/
##
pram::runPRAM(  in_gtf, in_bamv, out_gtf, method='plst',
                stringtie='/usr/loca/stringtie-1.3.3/stringtie')

Examples

PRAM has included input examples files in its extdata/demo/ folder. The table below provides a quick summary of all the example files.

runPRAM()’s input example files.
input argument file name(s)
in_gtf in.gtf
in_bamv SZP.bam, TLC.bam

You can access example files by system.file() in R, e.g. for the argument in_gtf, you can access its example file by

system.file('extdata/demo/in.gtf', package='pram')
#> [1] "/tmp/Rtmp6fws3p/Rinst31175651c7b1f8/pram/extdata/demo/in.gtf"

Below shows the usage of runPRAM() with example input files:

in_gtf = system.file('extdata/demo/in.gtf', package='pram')

in_bamv = c(system.file('extdata/demo/SZP.bam', package='pram'),
            system.file('extdata/demo/TLC.bam', package='pram') )

pred_out_gtf = tempfile(fileext='.gtf')

##
## assuming the stringtie binary is in folder /usr/local/stringtie-1.3.3/
##
pram::runPRAM(  in_gtf, in_bamv, pred_out_gtf, method='plst',
                stringtie='/usr/loca/stringtie-1.3.3/stringtie')

Define intergenic genomic ranges: defIgRanges()

Description

To predict intergenic transcripts, we must first define intergenic regions by defIgRanges(). This function requires a GTF file containing known gene annotation supplied for its in_gtf argument. This GTF file should contain an attribue of gene_id in its ninth column. We provided an example input GTF file in PRAM package: extdata/gtf/defIGRanges_in.gtf.

In addition to gene annotation, defIgRanges() also requires user to provide chromosome sizes so that it would know the maximum genomic ranges. You can provide one of the following arguments:

By default, defIgRanges() will define intergenic ranges as regions 10 kb away from any known genes. You can change it by the radius argument.

Example

pram::defIgRanges(system.file('extdata/gtf/defIgRanges_in.gtf', package='pram'),
                genome = 'hg38')
#> GRanges object with 456 ranges and 0 metadata columns:
#>                 seqnames          ranges strand
#>                    <Rle>       <IRanges>  <Rle>
#>     [1]             chr1 30001-248956422      *
#>     [2]             chr2     1-242193529      *
#>     [3]             chr3     1-198295559      *
#>     [4]             chr4     1-190214555      *
#>     [5]             chr5     1-181538259      *
#>     ...              ...             ...    ...
#>   [452] chrUn_KI270753v1         1-62944      *
#>   [453] chrUn_KI270754v1         1-40191      *
#>   [454] chrUn_KI270755v1         1-36723      *
#>   [455] chrUn_KI270756v1         1-79590      *
#>   [456] chrUn_KI270757v1         1-71251      *
#>   -------
#>   seqinfo: 455 sequences from an unspecified genome; no seqlengths

Prepare input RNA-seq alignments: prepIgBam()

Description

Once intergenic regions were defined, prepIgBam() will extract corresponding RNA-seq alignments from input BAM files. In this way, transcript models predicted at later stage will solely from intergenic regions. Also, with fewer RNA-seq alignments, model prediction will run faster.

Three input arguments are required by prepIgBam():

Example

finbam =system.file('extdata/bam/CMPRep2.sortedByCoord.raw.bam', 
                    package='pram')

iggrs = GenomicRanges::GRanges('chr10:77236000-77247000:+')

foutbam = tempfile(fileext='.bam')

pram::prepIgBam(finbam, iggrs, foutbam)
#> Extracted RNA-seq alignments are saved in /tmp/RtmpNcxstO/file311c7611bad553.bam

Build transcript models: buildModel()

Description

buildModel() predict transcript models from RNA-seq BAM file(s). This function requires two arguments:

Transcript prediction methods

buildModel() has implemented seven transcript prediction methods. You can specify it by the method argument with one of the keywords: plcf, plst, cfmg, cftc, stmg, cf, and st. The first five denote meta-assembly methods that utilize multiple RNA-seq datasets to predict a single set of transcript models. The last two represent methods that predict transcript models from a single RNA-seq dataset.

The table below compares prediction steps for these seven methods. By default, buildModel() uses plcf to predict transcript models.

(#tab:methods)Prediction steps of the seven buildModel() methods
method meta-assembly preparing RNA-seq input building transcripts
assembling transcripts
plcf yes pooling alignments Cufflinks no
plst yes pooling alignments StringTie no
cfmg yes no Cufflinks Cuffmerge
cftc yes no Cufflinks TACO
stmg yes no StringTie StringTie-merge
cf no no Cufflinks no
st no no StringTie no

Required external software

Depending on your specified prediction method, buildModel() requires external software: Cufflinks, StringTie and/or TACO, to build and/or assemble transcript models. You can either specify the software location using the cufflinks, stringtie, and taco arguments in buildModel(), or simply leave these three arugments undefined and let PRAM download them for you automatically. The table below summarized software versions buildModel() would download when required software was not specified. Please note that, for macOS, pre-compiled Cufflinks binary versions 2.2.1 and 2.2.0 appear to have an issue on processing BAM files, therefore we recommend to use version 2.1.1 instead.

(#tab:software)buildModel()-required software and recommended version
software Linux binary macOS binary required by
Cufflinks, Cuffmerge v2.2.1 v2.1.1 plcf, cfmg, cftc, and cf
StringTie, StringTie-merge v1.3.3b v1.3.3b plst, stmg, and st
TACO v0.7.0 v0.7.0 cftc

Example

fbams = c(  system.file('extdata/bam/CMPRep1.sortedByCoord.clean.bam', 
                        package='pram'),
            system.file('extdata/bam/CMPRep2.sortedByCoord.clean.bam', 
                        package='pram') )

foutgtf = tempfile(fileext='.gtf')

##
## assuming the stringtie binary is in folder /usr/local/stringtie-1.3.3/
##
pram::buildModel(fbams, foutgtf, method='plst',
                stringtie='/usr/loca/stringtie-1.3.3/stringtie')

Select transcript models: selModel()

Description

Once transcript models were built, you may want to select a subset of them by their genomic features. selModel() was developed for this purpose. It allows you to select transcript models by their total number of exons and total length of exons and introns.

selModel() requires two arguments:

By default: selModel() will select transcript models with \(\ge\) 2 exons and \(\ge\) 200 bp total length of exons and introns. You can change the default using the min_n_exon and min_tr_len arguments.

Example

fin_gtf = system.file('extdata/gtf/selModel_in.gtf', package='pram')

fout_gtf = tempfile(fileext='.gtf')

pram::selModel(fin_gtf, fout_gtf)
#> Selected transcript models are saved in /tmp/RtmpNcxstO/file311c76315422cf.gtf

Evaluate transcript models: evalModel()

Motivation

After PRAM has predicted a number of transcript models, you may wonder how accurate these models are. To answer this question, you can compare PRAM models with real transcripts (i.e., positive controls) that you know should be predicted. PRAM’s evalModel() function will help you to make such comparison. It will calculate precision and recall rates on three features of a transcript: exon nucleotides, individual splice junctions, and transcript structure (i.e., whether all splice junctions within a transcript were constructed in a model).

Input

evalModel() requires two arguments:

The two arguments can be in multiple formats:

Output

The output of evalModel() is a data.table object, where columns are evaluation results and each row is three transcript features.

evalModel() output columns
column name representation
feat transcript feature
ntp number of true positives (TP)
nfn number of false negatives (FN)
nfp number of false positives (FP)
precision precision rate: \(\frac{TP}{(TP+FP)}\)
recall recall rate: \(\frac{TP}{(TP+FN)}\)
evalModel() output rows
feature name representation
exon_nuc exon nucleotide
indi_jnc individual splice junction
tr_jnc transcript structure

Example

fmdl = system.file('extdata/benchmark/plcf.tsv', package='pram')
ftgt = system.file('extdata/benchmark/tgt.tsv',  package='pram')

mdldt = data.table::fread(fmdl, header=TRUE, sep="\t")
tgtdt = data.table::fread(ftgt, header=TRUE, sep="\t")

pram::evalModel(mdldt, tgtdt)
#>        feat     ntp    nfn  nfp precision    recall
#> 1: exon_nuc 1723581 109337 8424 0.9951363 0.9403481
#> 2: indi_jnc    2889    162  192 0.9376826 0.9469027
#> 3:   tr_jnc    1138    118  252 0.8187050 0.9060510

Session Info

Below is the output of sessionInfo() on the system on which this document was compiled.

#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> 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       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] GenomicAlignments_1.38.0 XVector_0.42.0           GenomicRanges_1.54.0    
#> [4] S4Vectors_0.40.0         IRanges_2.36.0          
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.6-1.1              jsonlite_1.8.7             
#>  [3] rjson_0.2.21                compiler_4.3.1             
#>  [5] crayon_1.5.2                SummarizedExperiment_1.32.0
#>  [7] Biobase_2.62.0              Rsamtools_2.18.0           
#>  [9] Biostrings_2.70.0           bitops_1.0-7               
#> [11] parallel_4.3.1              jquerylib_0.1.4            
#> [13] BiocParallel_1.36.0         yaml_2.3.7                 
#> [15] fastmap_1.1.1               lattice_0.22-5             
#> [17] R6_2.5.1                    S4Arrays_1.2.0             
#> [19] pram_1.18.0                 GenomeInfoDb_1.38.0        
#> [21] knitr_1.44                  BiocGenerics_0.48.0        
#> [23] XML_3.99-0.14               DelayedArray_0.28.0        
#> [25] MatrixGenerics_1.14.0       GenomeInfoDbData_1.2.11    
#> [27] bslib_0.5.1                 rlang_1.1.1                
#> [29] cachem_1.0.8                xfun_0.40                  
#> [31] sass_0.4.7                  SparseArray_1.2.0          
#> [33] cli_3.6.1                   zlibbioc_1.48.0            
#> [35] digest_0.6.33               grid_4.3.1                 
#> [37] rtracklayer_1.62.0          evaluate_0.22              
#> [39] data.table_1.14.8           codetools_0.2-19           
#> [41] abind_1.4-5                 stats4_4.3.1               
#> [43] RCurl_1.98-1.12             restfulr_0.0.15            
#> [45] rmarkdown_2.25              BiocIO_1.12.0              
#> [47] matrixStats_1.0.0           tools_4.3.1                
#> [49] htmltools_0.5.6.1