In this vignette, we show how to prepare reference files for estimating abundances of spliced and unspliced abundances with alignment-free methods (e.g., Salmon, alevin or kallisto). Such abundances are used as input, e.g., for estimation of RNA velocity in single-cell data.
library(eisaR)The first step is to generate a GRangesList object containing the genomic
ranges for the features of interest (spliced transcripts, and either unspliced
transcripts or intron sequences).
This is done using the getFeatureRanges() function, based on a reference GTF
file.
Here, we exemplify this with a small subset of a GTF file from
Gencode release 28.
We extract genomic ranges for spliced transcript and introns, where the latter
are defined for each transcript separately (using the same terminology as in
the BUSpaRse package).
For each intron, a flanking sequence of 50 nt is added on each side to
accommodate reads mapping across an exon/intron boundary.
gtf <- system.file("extdata/gencode.v28.annotation.sub.gtf",
package = "eisaR")
grl <- getFeatureRanges(
gtf = gtf,
featureType = c("spliced", "intron"),
intronType = "separate",
flankLength = 50L,
joinOverlappingIntrons = FALSE,
verbose = TRUE
)
#> Import genomic features from the file as a GRanges object ... OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .get_cds_IDX(mcols0$type, mcols0$phase): The "phase" metadata column contains non-NA values for
#> features of type stop_codon. This information was ignored.
#> OK
#> 'select()' returned 1:1 mapping between keys and columns
#> Extracting spliced transcript features
#> Extracting introns using the separate approachThe output from getFeatureRanges() is a GRangesList object, with all
the features of interest (both spliced transcripts and introns):
grl
#> GRangesList object of length 895:
#> $ENST00000456328.2
#> GRanges object with 3 ranges and 5 metadata columns:
#> seqnames ranges strand | exon_id exon_rank
#> <Rle> <IRanges> <Rle> | <character> <integer>
#> [1] chr1 11869-12227 + | ENSE00002234944.1 1
#> [2] chr1 12613-12721 + | ENSE00003582793.1 2
#> [3] chr1 13221-14409 + | ENSE00002312635.1 3
#> transcript_id gene_id type
#> <character> <character> <character>
#> [1] ENST00000456328.2 ENSG00000223972.5 exon
#> [2] ENST00000456328.2 ENSG00000223972.5 exon
#> [3] ENST00000456328.2 ENSG00000223972.5 exon
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
#>
#> $ENST00000450305.2
#> GRanges object with 6 ranges and 5 metadata columns:
#> seqnames ranges strand | exon_id exon_rank
#> <Rle> <IRanges> <Rle> | <character> <integer>
#> [1] chr1 12010-12057 + | ENSE00001948541.1 1
#> [2] chr1 12179-12227 + | ENSE00001671638.2 2
#> [3] chr1 12613-12697 + | ENSE00001758273.2 3
#> [4] chr1 12975-13052 + | ENSE00001799933.2 4
#> [5] chr1 13221-13374 + | ENSE00001746346.2 5
#> [6] chr1 13453-13670 + | ENSE00001863096.1 6
#> transcript_id gene_id type
#> <character> <character> <character>
#> [1] ENST00000450305.2 ENSG00000223972.5 exon
#> [2] ENST00000450305.2 ENSG00000223972.5 exon
#> [3] ENST00000450305.2 ENSG00000223972.5 exon
#> [4] ENST00000450305.2 ENSG00000223972.5 exon
#> [5] ENST00000450305.2 ENSG00000223972.5 exon
#> [6] ENST00000450305.2 ENSG00000223972.5 exon
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
#>
#> $ENST00000473358.1
#> GRanges object with 3 ranges and 5 metadata columns:
#> seqnames ranges strand | exon_id exon_rank
#> <Rle> <IRanges> <Rle> | <character> <integer>
#> [1] chr1 29554-30039 + | ENSE00001947070.1 1
#> [2] chr1 30564-30667 + | ENSE00001922571.1 2
#> [3] chr1 30976-31097 + | ENSE00001827679.1 3
#> transcript_id gene_id type
#> <character> <character> <character>
#> [1] ENST00000473358.1 ENSG00000243485.5 exon
#> [2] ENST00000473358.1 ENSG00000243485.5 exon
#> [3] ENST00000473358.1 ENSG00000243485.5 exon
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
#>
#> ...
#> <892 more elements>The metadata slot of the GRangesList object contains a list of the feature
IDs for each type of feature:
lapply(S4Vectors::metadata(grl)$featurelist, head)
#> $spliced
#> [1] "ENST00000456328.2" "ENST00000450305.2" "ENST00000473358.1"
#> [4] "ENST00000469289.1" "ENST00000607096.1" "ENST00000606857.1"
#>
#> $intron
#> [1] "ENST00000456328.2-I" "ENST00000456328.2-I1"
#> [3] "ENST00000450305.2-I" "ENST00000450305.2-I1"
#> [5] "ENST00000450305.2-I2" "ENST00000450305.2-I3"As we can see, the intron IDs are identified by a -I suffix.
Each feature is further annotated to a gene ID.
For the intronic features, the corresponding gene ID also bears the -I
suffix appended to the original gene ID.
Having separate gene IDs for spliced transcripts and introns allows, for
example, joint quantification of spliced and unspliced versions of a gene
with alevin.
Adding a suffix rather than defining a completely new gene ID allows us to
easily match corresponding spliced and unspliced feature IDs.
To further simplify the latter, the metadata of the GRangesList object
returned by getFeatureRanges() contains data.frames matching the
corresponding gene IDs (as well as transcript IDs, if unspliced transcripts
are extracted):
head(S4Vectors::metadata(grl)$corrgene)
#> spliced intron
#> 1 ENSG00000223972.5 ENSG00000223972.5-I
#> 2 ENSG00000243485.5 ENSG00000243485.5-I
#> 3 ENSG00000284332.1 ENSG00000284332.1-I
#> 4 ENSG00000268020.3 ENSG00000268020.3-I
#> 5 ENSG00000240361.2 ENSG00000240361.2-I
#> 6 ENSG00000186092.6 ENSG00000186092.6-IOnce the genomic ranges of the features of interest are extracted, we can
obtain the nucleotide sequences using the extractTranscriptSeqs() function
from the GenomicFeatures package.
In addition to the ranges, this requires the genome sequence.
This can be obtained, for example, from the corresponding BSgenome package,
or from an external FASTA file.
suppressPackageStartupMessages({
library(BSgenome)
library(BSgenome.Hsapiens.UCSC.hg38)
})
seqs <- GenomicFeatures::extractTranscriptSeqs(
x = BSgenome.Hsapiens.UCSC.hg38,
transcripts = grl
)
seqs
#> DNAStringSet object of length 895:
#> width seq names
#> [1] 1657 GTTAACTTGCCGTCAGC...ACACTGTTGGTTTCTG ENST00000456328.2
#> [2] 632 GTGTCTGACTTCCAGCA...ACAGGGGAATCCCGAA ENST00000450305.2
#> [3] 712 GTGCACACGGCTCCCAT...TGAGGGATAAATGTAT ENST00000473358.1
#> [4] 535 TCATCAGTCCAAAGTCC...GTATGATTTTAAAGTC ENST00000469289.1
#> [5] 138 GGATGCCCAGCTAGTTT...AGAATTAAGCATTTTA ENST00000607096.1
#> ... ... ...
#> [891] 178 GCGCCAGCCGGACCCCA...CGCGTATTAACGAGAG ENST00000304952.1...
#> [892] 193 GGAGATGACCGTGAGAC...GTACCGCGCCGGCTTC ENST00000481869.1-I
#> [893] 193 GGAGATGACCGTGAGAC...GTACCGCGCCGGCTTC ENST00000484667.2-I
#> [894] 352 GCGCCAGCCGGACCCCA...TCCTGGAGATGACCGT ENST00000484667.2-I1
#> [895] 826 ACCTCCGCCTGCCAGGC...AGAGATTGTGGTGAGC ENST00000458555.1-IThe resulting DNAStringSet can be written to a FASTA file and used to
generate an index for alignment-free methods such as Salmon and kallisto.
In addition to extracting feature sequences, we can also export a GTF file with the full set of features. This is useful, for example, in order to generate a linked transcriptome for later import of estimated abundances with tximeta.
exportToGtf(
grl,
filepath = file.path(tempdir(), "exported.gtf")
)In the exported GTF file, each entry of grl will correspond to a “transcript”
feature, and each individual range corresponds to an “exon” feature.
In addition, each gene is represented as a “gene” feature.
Finally, we can export a data.frame (or a tab-separated test file) providing
a conversion table between “transcript” and “gene” identifiers.
This is needed to aggregate the transcript-level abundance estimates from
alignment-free methods such as Salmon and kallisto to the gene level.
df <- getTx2Gene(grl)
head(df)
#> transcript_id gene_id
#> ENST00000456328.2 ENST00000456328.2 ENSG00000223972.5
#> ENST00000450305.2 ENST00000450305.2 ENSG00000223972.5
#> ENST00000473358.1 ENST00000473358.1 ENSG00000243485.5
#> ENST00000469289.1 ENST00000469289.1 ENSG00000243485.5
#> ENST00000607096.1 ENST00000607096.1 ENSG00000284332.1
#> ENST00000606857.1 ENST00000606857.1 ENSG00000268020.3
tail(df)
#> transcript_id gene_id
#> ENST00000304952.10-I1 ENST00000304952.10-I1 ENSG00000188290.10-I
#> ENST00000304952.10-I2 ENST00000304952.10-I2 ENSG00000188290.10-I
#> ENST00000481869.1-I ENST00000481869.1-I ENSG00000188290.10-I
#> ENST00000484667.2-I ENST00000484667.2-I ENSG00000188290.10-I
#> ENST00000484667.2-I1 ENST00000484667.2-I1 ENSG00000188290.10-I
#> ENST00000458555.1-I ENST00000458555.1-I ENSG00000224969.1-IsessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 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] stats4 parallel stats graphics grDevices utils
#> [7] datasets methods base
#>
#> other attached packages:
#> [1] BSgenome.Hsapiens.UCSC.hg38_1.4.3
#> [2] BSgenome_1.58.0
#> [3] Biostrings_2.58.0
#> [4] XVector_0.30.0
#> [5] edgeR_3.32.0
#> [6] limma_3.46.0
#> [7] QuasR_1.30.0
#> [8] Rbowtie_1.30.0
#> [9] rtracklayer_1.50.0
#> [10] GenomicFeatures_1.42.0
#> [11] AnnotationDbi_1.52.0
#> [12] Biobase_2.50.0
#> [13] GenomicRanges_1.42.0
#> [14] GenomeInfoDb_1.26.0
#> [15] IRanges_2.24.0
#> [16] S4Vectors_0.28.0
#> [17] BiocGenerics_0.36.0
#> [18] eisaR_1.2.0
#> [19] BiocStyle_2.18.0
#>
#> loaded via a namespace (and not attached):
#> [1] MatrixGenerics_1.2.0 httr_1.4.2
#> [3] splines_4.0.3 bit64_4.0.5
#> [5] GenomicFiles_1.26.0 assertthat_0.2.1
#> [7] askpass_1.1 BiocManager_1.30.10
#> [9] BiocFileCache_1.14.0 latticeExtra_0.6-29
#> [11] blob_1.2.1 GenomeInfoDbData_1.2.4
#> [13] Rsamtools_2.6.0 yaml_2.2.1
#> [15] progress_1.2.2 pillar_1.4.6
#> [17] RSQLite_2.2.1 lattice_0.20-41
#> [19] glue_1.4.2 digest_0.6.27
#> [21] RColorBrewer_1.1-2 htmltools_0.5.0
#> [23] Matrix_1.2-18 XML_3.99-0.5
#> [25] pkgconfig_2.0.3 ShortRead_1.48.0
#> [27] biomaRt_2.46.0 magick_2.5.0
#> [29] bookdown_0.21 zlibbioc_1.36.0
#> [31] purrr_0.3.4 jpeg_0.1-8.1
#> [33] BiocParallel_1.24.0 tibble_3.0.4
#> [35] openssl_1.4.3 generics_0.0.2
#> [37] ellipsis_0.3.1 SummarizedExperiment_1.20.0
#> [39] magrittr_1.5 crayon_1.3.4
#> [41] memoise_1.1.0 evaluate_0.14
#> [43] hwriter_1.3.2 xml2_1.3.2
#> [45] tools_4.0.3 prettyunits_1.1.1
#> [47] hms_0.5.3 lifecycle_0.2.0
#> [49] matrixStats_0.57.0 stringr_1.4.0
#> [51] locfit_1.5-9.4 DelayedArray_0.16.0
#> [53] compiler_4.0.3 rlang_0.4.8
#> [55] grid_4.0.3 RCurl_1.98-1.2
#> [57] VariantAnnotation_1.36.0 rappdirs_0.3.1
#> [59] bitops_1.0-6 rmarkdown_2.5
#> [61] DBI_1.1.0 curl_4.3
#> [63] R6_2.4.1 GenomicAlignments_1.26.0
#> [65] knitr_1.30 dplyr_1.0.2
#> [67] bit_4.0.4 stringi_1.5.3
#> [69] Rcpp_1.0.5 png_0.1-7
#> [71] vctrs_0.3.4 dbplyr_1.4.4
#> [73] tidyselect_1.1.0 xfun_0.18