Working with transcripts
Introduction
Welcome to the introduction to working with transcripts in ORFik.
This vignette will walk you through how to find and manipulate the
sections of transcript you want. And how to calculate coverage over
these regions. ORFik is an R package containing various
functions for analysis of RiboSeq, RNASeq, RCP-seq, TCP-seq, Chip-seq
and Cage data, we advice you to read Importing data vignette, before
starting this one.
Motivation
Working with all the specific regions of transcripts can be hard. The standard Bioconductor packages GenomicRanges and GenomicFeatures does not solve most of these problems, so ORFik extends this ecosystem with what we call TranscriptRanges and TranscriptFeatures. We will here learn how to use ORFik effectively getting the transcript regions we want and compute coverage over those regions etc.
Getting data
Use the bundled synthetic reference from the importing data vignette.
This keeps the examples reproducible without downloading a genome. For a
real analysis, use a matching genome and annotation for your organism;
see getGenomeAndAnnotation() in the importing data
vignette.
library(ORFik)
reference_dir <- system.file("extdata/references/homo_sapiens", package = "ORFik")
output_dir <- tempfile("orfik-transcripts-")
dir.create(output_dir)
reference_files <- c("Homo_sapiens_dummy.fasta", "Homo_sapiens_dummy.gtf")
file.copy(file.path(reference_dir, reference_files), output_dir)## [1] TRUE TRUEgenome <- file.path(output_dir, reference_files[1])
gtf <- file.path(output_dir, reference_files[2])
indexFa(genome)## [1] "/tmp/Rtmp0yielu/orfik-transcripts-19fd72828651/Homo_sapiens_dummy.fasta.fai"## Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
## is not available for this TxDb objectLoading Transcript region (single point locations)
Let us get these regions:
- Transcription start sites (TSS, start of leaders)
- Translation start sites (TIS, start of CDSs)
- Translation termination sites (TTS, end of CDSs)
- Transcription end sites (TES, end of mRNA)
## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 306 +
## [2] chr2 856 -
## [3] chr3 306 +
## [4] chr4 306 +
## [5] chr5 856 -
## [6] chr6 306 +
## -------
## seqinfo: 6 sequences from an unspecified genome## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 306 +
## [2] chr2 856 -
## [3] chr3 306 +
## [4] chr4 306 +
## [5] chr5 856 -
## [6] chr6 306 +
## -------
## seqinfo: 6 sequences from an unspecified genome## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 446 +
## [2] chr2 751 -
## [3] chr3 446 +
## [4] chr4 446 +
## [5] chr5 751 -
## [6] chr6 446 +
## -------
## seqinfo: 6 sequences from an unspecified genome## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 751 +
## [2] chr2 446 -
## [3] chr3 751 +
## [4] chr4 751 +
## [5] chr5 446 -
## [6] chr6 751 +
## -------
## seqinfo: 6 sequences from an unspecified genome## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 856 +
## [2] chr2 306 -
## [3] chr3 856 +
## [4] chr4 856 +
## [5] chr5 306 -
## [6] chr6 856 +
## -------
## seqinfo: 6 sequences from an unspecified genomeLoading Transcript region (window locations)
What if we want a window around that region point? Then we also need might the mrna region information, to map to transcript coordinates and get area around that point
mrna <- loadRegion(txdb, "mrna")
## TSS
startRegion(loadRegion(txdb, "mrna"), upstream = 30, downstream = 30, is.sorted = TRUE)## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10001 chr1 306-336 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10002
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10002 chr2 826-856 -
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10003
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10003 chr3 306-336 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## ...
## <3 more elements>## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10001 chr1 416-445 +
## ENSTTEST10001 chr1 446-476 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10002
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10002 chr2 752-781 -
## ENSTTEST10002 chr2 721-751 -
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10003
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10003 chr3 416-445 +
## ENSTTEST10003 chr3 446-476 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## ...
## <3 more elements>## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10001 chr1 721-751 +
## ENSTTEST10001 chr1 752-781 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10002
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10002 chr2 446-476 -
## ENSTTEST10002 chr2 416-445 -
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10003
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10003 chr3 721-751 +
## ENSTTEST10003 chr3 752-781 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## ...
## <3 more elements>## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10001 chr1 826-856 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10002
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10002 chr2 306-336 -
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## $ENSTTEST10003
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## ENSTTEST10003 chr3 826-856 +
## -------
## seqinfo: 6 sequences from an unspecified genome
##
## ...
## <3 more elements>Extending transcript into genomic flanks
You might notice that since we are doing transcript coordinates, it will not actually extend 30 upstream from TSS, because there is no transcript there.
To extend into the genomic region do:
# TSS genomic extension
extendLeaders(startRegion(mrna[1], upstream = -1, downstream = 30, is.sorted = TRUE), 30)## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 1 range and 1 metadata column:
## seqnames ranges strand | names
## <Rle> <IRanges> <Rle> | <character>
## ENSTTEST10001 chr1 277-336 + | ENSTTEST10001
## -------
## seqinfo: 6 sequences from an unspecified genome# TES genomic extension
extendTrailers(stopRegion(mrna[1], upstream = 30, downstream = 0, is.sorted = TRUE), 30)## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 1 range and 1 metadata column:
## seqnames ranges strand | names
## <Rle> <IRanges> <Rle> | <character>
## ENSTTEST10001 chr1 826-886 + | ENSTTEST10001
## -------
## seqinfo: 6 sequences from an unspecified genomeCoverage over transcript regions
So how could we now get the coverage per nt in the TIS region of a cds?
Loading NGS data
We first need some NGS data with matching annotation, so let us now use the ORFik experiment included in the package:
## ORFik experiment: ORFik (Tjeldnes et al.)
## Libraries: 4 library types and 16 runs
## Organism: Homo sapiens
## libtype rep condition
## 1: CAGE 1 Mutant
## 2: CAGE 2 Mutant
## 3: CAGE 1 WT
## 4: CAGE 2 WT
## 5: PAS 1 Mutant
## 6: PAS 2 Mutant
## 7: PAS 1 WT
## 8: PAS 2 WT
## 9: RFP 1 Mutant
## 10: RFP 2 Mutant
## 11: RFP 1 WT
## 12: RFP 2 WT
## 13: RNA 1 Mutant
## 14: RNA 2 Mutant
## 15: RNA 1 WT
## 16: RNA 2 WTThere are 4 CAGE libraries, 4 PAS (Poly-A seq) libraries, 4 Ribo-seq libraries and 4 RNA-seq libraries. Let’s load the first Ribo-seq library in this artificial study
## ORFik experiment: ORFik (Tjeldnes et al.)
## Libraries: 1 library type and 1 runs
## Organism: Homo sapiens
## libtype
## 1: RFPLoad the Ribo-seq
Some NGS statistics, first Read length distribution
##
## 27 28 29
## 1001 961 984Then RFP peak count distribution
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 4.00 15.00 49.68 57.75 775.00Total number of unique (collapsed) and original (non collapsed reads):
## [1] "Number of collapsed: 2946 Number of non-collapsed: 146368"## [1] "duplication rate: 98 %"Total counts per region/window
TIS_window <- startRegion(loadRegion(df, "cds"), loadRegion(df, "mrna"),
is.sorted = TRUE, upstream = 20, downstream = 20)Now let’s get the total counts per window, we will do it in two identical ways. The second version uses the automatic 5’ end anchoring to create an identical output.
countOverlapsW(TIS_window, ribo, weight = "score") # score is the number of pshifted RFP peaks at respective position## ENSTTEST10001 ENSTTEST10002 ENSTTEST10003 ENSTTEST10004 ENSTTEST10005
## 1416 1818 1571 1662 1300
## ENSTTEST10006
## 1537# This is shorter version (You do not need TIS_window defined first) ->
startRegionCoverage(loadRegion(df, "cds"), ribo, loadRegion(df, "mrna"),
is.sorted = TRUE, upstream = 20, downstream = 20)## ENSTTEST10001 ENSTTEST10002 ENSTTEST10003 ENSTTEST10004 ENSTTEST10005
## 33 36 37 39 36
## ENSTTEST10006
## 38You see the first gene had 33 reads in the 20-20 window around TIS
Per nucleotide in region/window
A per nucleotide coverage is in Bioconductor called a ‘tiling’, if you want the coverage per nucleotide in that predefined window?
## RleList of length 3
## $ENSTTEST10001
## numeric-Rle of length 41 with 39 runs
## Lengths: 1 1 1 1 3 1 1 1 ... 1 1 1 1 1 1 1 1
## Values : 0 15 3 7 0 9 2 7 ... 68 7 120 20 0 155 86 9
##
## $ENSTTEST10002
## numeric-Rle of length 41 with 38 runs
## Lengths: 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1
## Values : 0 24 3 1 15 8 1 8 ... 190 255 31 64 10 11 46 285
##
## $ENSTTEST10003
## numeric-Rle of length 41 with 40 runs
## Lengths: 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1
## Values : 2 9 4 5 1 0 1 21 ... 0 80 130 14 11 185 194 17For plotting it is best to use the data.table return:
## count genes position
## <num> <int> <int>
## 1: 0 1 1
## 2: 15 1 2
## 3: 3 1 3
## 4: 7 1 4
## 5: 0 1 5
## 6: 0 1 6Let’s plot the coverage for fun, now we want frame information too:
# TIS region coverage
dt <- coveragePerTiling(TIS_window, ribo, as.data.table = TRUE, withFrames = TRUE)
# Rescale 0 to position 21
dt[, position := position - 21]
pSitePlot(dt)
Per nucleotide in region/window (split by read length)
To get coverage of a window per read length, we can anchor on a specific location, namely the 5’ end of CDSs (i.e. TISs) in this case.
# Define transcript with valid UTRs and lengths
txNames <- filterTranscripts(df, 25, 25,0, longestPerGene = TRUE)
# TIS region coverage
dt <- windowPerReadLength(loadRegion(df, "cds", txNames), loadRegion(df, "mrna", txNames), ribo)
# Remember to set scoring to scoring used above for dt
coverageHeatMap(dt, scoring = "transcriptNormalized")
To get coverage of a full region per read length (coverage of full transcript per read length), we do:
dt <- regionPerReadLength(loadRegion(df, "cds", txNames)[1], ribo,
exclude.zero.cov.grl = FALSE, drop.zero.dt = FALSE)
# Remember to set scoring to scoring used above for dt
coverageHeatMap(dt, scoring = NULL)
You see the arguments exclude.zero.cov.grl and drop.zero.dt are set to
false, default is true, which leads to considerable speed ups for large
datasets, by reducing number of 0 positions that are kept in the final
object.
Advanced details
Most coverage calculations in ORFik in some way relies on the GenomicRanges::coverage function, to understand how it all works and how to speed it up, look up the two functions:
- ORFik:::coverageByTranscriptW (Must convert to coverage)
- ORFik:::coverageByTranscriptC (Uses precomputed RLE coverage tracks, much faster)
The remaining logic is usually pure R loops, as very little computation is done in R code anyway. The output relies on data.table package for efficient computation on results and easy plotting.
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.81.1
## [3] rtracklayer_1.73.0 BiocIO_1.23.3
## [5] data.table_1.18.6.1 GenomicFeatures_1.65.0
## [7] AnnotationDbi_1.75.2 Biobase_2.73.2
## [9] fstcore_0.10.0 ORFik_1.33.1
## [11] GenomicAlignments_1.49.2 Rsamtools_2.29.0
## [13] Biostrings_2.81.9 XVector_0.53.0
## [15] GenomicRanges_1.65.4 Seqinfo_1.3.2
## [17] IRanges_2.47.5 S4Vectors_0.51.10
## [19] BiocGenerics_0.59.12 generics_0.1.4
## [21] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0
## [3] gridExtra_2.3.1 httr2_1.3.0
## [5] biomaRt_2.69.3 rlang_1.3.0
## [7] magrittr_2.0.5 otel_0.2.0
## [9] matrixStats_1.5.0 compiler_4.6.1
## [11] RSQLite_3.53.3 txdbmaker_1.9.0
## [13] png_0.1-9 vctrs_0.7.3
## [15] stringr_1.6.0 pkgconfig_2.0.3
## [17] crayon_1.5.3 fastmap_1.2.0
## [19] dbplyr_2.6.0 labeling_0.4.3
## [21] biomartr_1.0.7 rmarkdown_2.32
## [23] UCSC.utils_1.9.0 bit_4.6.0
## [25] xfun_0.61 cachem_1.1.0
## [27] cigarillo_1.3.1 GenomeInfoDb_1.49.1
## [29] jsonlite_2.0.0 progress_1.2.3
## [31] blob_1.3.0 DelayedArray_0.39.6
## [33] BiocParallel_1.47.0 parallel_4.6.1
## [35] prettyunits_1.2.0 R6_2.6.1
## [37] stringi_1.8.9 bslib_0.12.0
## [39] RColorBrewer_1.1-3 jquerylib_0.1.4
## [41] Rcpp_1.1.2 SummarizedExperiment_1.43.0
## [43] knitr_1.52 R.utils_2.13.0
## [45] BiocBaseUtils_1.15.1 Matrix_1.7-6
## [47] tidyselect_1.2.1 abind_1.4-8
## [49] yaml_2.3.12 stringfish_0.19.2
## [51] codetools_0.2-20 curl_8.0.0
## [53] lattice_0.23-1 tibble_3.3.1
## [55] withr_3.0.3 KEGGREST_1.53.6
## [57] S7_0.2.2 evaluate_1.0.5
## [59] RcppParallel_6.2.1 BiocFileCache_3.3.0
## [61] xml2_1.6.0 pillar_1.11.1
## [63] BiocManager_1.30.27 filelock_1.0.3
## [65] MatrixGenerics_1.25.0 RCurl_1.98-1.20
## [67] hms_1.1.4 ggplot2_4.0.3
## [69] scales_1.4.0 qs2_0.3.1
## [71] glue_1.8.1 maketools_1.3.2
## [73] tools_4.6.1 sys_3.4.3
## [75] locfit_1.5-9.12 buildtools_1.0.0
## [77] XML_3.99-0.24 cowplot_1.2.0
## [79] grid_4.6.1 GenomeInfoDbData_1.2.15
## [81] restfulr_0.0.17 cli_3.6.6
## [83] fst_0.9.8 S4Arrays_1.13.0
## [85] dplyr_1.2.1 gtable_0.3.6
## [87] R.methodsS3_1.8.2 DESeq2_1.53.3
## [89] sass_0.4.10 digest_0.6.39
## [91] SparseArray_1.13.2 rjson_0.2.23
## [93] farver_2.1.2 memoise_2.0.1
## [95] htmltools_0.5.9 R.oo_1.27.1
## [97] lifecycle_1.0.5 httr_1.4.9
## [99] bit64_4.8.6