# Generating reference files for spliced and unspliced abundance estimation with alignment-free methods

Charlotte Soneson

#### 2021-05-19

# 1 Introduction

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)

# 2 Generate feature ranges

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 approach

The 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.frame`s 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-I

# 3 Extract feature sequences

Once 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-I

The resulting `DNAStringSet` can be written to a FASTA file and used to generate an index for alignment-free methods such as Salmon and kallisto.

# 4 Generate an expanded GTF file

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.

# 5 Generate a transcript-to-gene mapping

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-I

# 6 Session info

    sessionInfo()
    #> R version 4.1.0 RC (2021-05-10 r80283)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 17763)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.1252
    #> [3] LC_MONETARY=English_United States.1252
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.1252
    #>
    #> 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.60.0
    #>  [3] Biostrings_2.60.0
    #>  [4] XVector_0.32.0
    #>  [5] edgeR_3.34.0
    #>  [6] limma_3.48.0
    #>  [7] QuasR_1.32.0
    #>  [8] Rbowtie_1.32.0
    #>  [9] rtracklayer_1.52.0
    #> [10] GenomicFeatures_1.44.0
    #> [11] AnnotationDbi_1.54.0
    #> [12] Biobase_2.52.0
    #> [13] GenomicRanges_1.44.0
    #> [14] GenomeInfoDb_1.28.0
    #> [15] IRanges_2.26.0
    #> [16] S4Vectors_0.30.0
    #> [17] BiocGenerics_0.38.0
    #> [18] eisaR_1.4.0
    #> [19] BiocStyle_2.20.0
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] bitops_1.0-7                matrixStats_0.58.0
    #>  [3] bit64_4.0.5                 RColorBrewer_1.1-2
    #>  [5] filelock_1.0.2              progress_1.2.2
    #>  [7] httr_1.4.2                  GenomicFiles_1.28.0
    #>  [9] tools_4.1.0                 bslib_0.2.5.1
    #> [11] utf8_1.2.1                  R6_2.5.0
    #> [13] DBI_1.1.1                   tidyselect_1.1.1
    #> [15] prettyunits_1.1.1           bit_4.0.4
    #> [17] curl_4.3.1                  compiler_4.1.0
    #> [19] SGSeq_1.26.0                DelayedArray_0.18.0
    #> [21] bookdown_0.22               sass_0.4.0
    #> [23] rappdirs_0.3.3              stringr_1.4.0
    #> [25] digest_0.6.27               Rsamtools_2.8.0
    #> [27] rmarkdown_2.8               Rhisat2_1.8.0
    #> [29] jpeg_0.1-8.1                pkgconfig_2.0.3
    #> [31] htmltools_0.5.1.1           MatrixGenerics_1.4.0
    #> [33] highr_0.9                   dbplyr_2.1.1
    #> [35] fastmap_1.1.0               rlang_0.4.11
    #> [37] RSQLite_2.2.7               jquerylib_0.1.4
    #> [39] BiocIO_1.2.0                generics_0.1.0
    #> [41] hwriter_1.3.2               jsonlite_1.7.2
    #> [43] BiocParallel_1.26.0         dplyr_1.0.6
    #> [45] VariantAnnotation_1.38.0    RCurl_1.98-1.3
    #> [47] magrittr_2.0.1              GenomeInfoDbData_1.2.6
    #> [49] Matrix_1.3-3                Rcpp_1.0.6
    #> [51] fansi_0.4.2                 lifecycle_1.0.0
    #> [53] stringi_1.6.2               yaml_2.2.1
    #> [55] SummarizedExperiment_1.22.0 zlibbioc_1.38.0
    #> [57] BiocFileCache_2.0.0         grid_4.1.0
    #> [59] blob_1.2.1                  crayon_1.4.1
    #> [61] lattice_0.20-44             splines_4.1.0
    #> [63] hms_1.1.0                   KEGGREST_1.32.0
    #> [65] magick_2.7.2                locfit_1.5-9.4
    #> [67] knitr_1.33                  pillar_1.6.1
    #> [69] igraph_1.2.6                RUnit_0.4.32
    #> [71] rjson_0.2.20                biomaRt_2.48.0
    #> [73] XML_3.99-0.6                glue_1.4.2
    #> [75] evaluate_0.14               ShortRead_1.50.0
    #> [77] latticeExtra_0.6-29         BiocManager_1.30.15
    #> [79] png_0.1-7                   vctrs_0.3.8
    #> [81] purrr_0.3.4                 assertthat_0.2.1
    #> [83] cachem_1.0.5                xfun_0.23
    #> [85] restfulr_0.0.13             tibble_3.1.2
    #> [87] GenomicAlignments_1.28.0    memoise_2.0.0
    #> [89] ellipsis_0.3.2