Contents

1 Introduction

This vignette outlines a workflow of detecting retrotransposed transcripts (RTs) from Variant Call Format (VCF) using the svaRetro package.

2 Installation

The svaRetro package can be installed from Bioconductor as follows:

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

3 Using GRanges for structural variants: a breakend-centric data structure

This package uses a breakend-centric event notation adopted from the StructuralVariantAnnotation package. In short, breakends are stored in a GRanges object with strand used to indicate breakpoint orientation, where breakpoints are represented using a partner field containing the name of the breakend at the other side of the breakend. This notation was chosen as it simplifies the annotations of RTs which are detected at breakend-level.

4 Workflow

4.1 Loading data from VCF

VCF data is parsed into a VCF object using the readVCF function from the Bioconductor package VariantAnnotation. Simple filters could be applied to a VCF object to remove unwanted calls. The VCF object is then converted to a GRanges object with breakend-centric notations using StructuralVariantAnnotation. More information about VCF objects and breakend-centric GRanges object can be found by consulting the vignettes in the corresponding packages with browseVignettes("VariantAnnotation") and browseVignettes("StructuralVariantAnnotation").

library(StructuralVariantAnnotation)
library(VariantAnnotation)
library(svaRetro)

RT_vcf <- readVcf(system.file("extdata", "diploidSV.vcf", package = "svaRetro"))
RT_gr <- StructuralVariantAnnotation::breakpointRanges(RT_vcf, 
                                                       nominalPosition=TRUE)
head(RT_gr)

Note that StructuralVariantAnnotation requires the GRanges object to be composed entirely of valid breakpoints. Please consult the vignette of the StructuralVariantAnnotation package for ensuring breakpoint consistency.

4.2 Identifying Retrotransposed Transcripts

The package provides rtDetect to identify RTs using the provided SV calls. This is achieved by detecting intronic deletions, which are breakpoints at exon-intron (and intron-exon) boundaries of a transcript. Fusions consisting of an exon boundary and a second genomic location are reported as potential insertion sites. Due to the complexity of RT events, insertion sites can be discovered on both left and right sides, only one side, or none at all.

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(dplyr)
hg19.genes <- TxDb.Hsapiens.UCSC.hg19.knownGene

RT <- rtDetect(RT_gr, hg19.genes, maxgap=10, minscore=0.8)

The output is a list of GRanges object consisting of two sets of GRanges calls, insSite and junctions, containing candidate insertion sites and exon-exon junctions respectively. Candidate insertion sites are annotated by the source transcripts and whether exon-exon junctions are detected for the source transcripts. RT junction breakends are annotated by the UCSC exon IDs, corresponding transcripts, and NCBI gene symbols.

RT$SKA3
#> $junctions
#> GRanges object with 14 ranges and 17 metadata columns:
#>                                  seqnames    ranges strand | paramRangeID                    REF         ALT      QUAL      FILTER               sourceId                partner      svtype     svLen      insSeq    insLen       event    HOMLEN      exon                   txs  exons gene_symbol
#>                                     <Rle> <IRanges>  <Rle> |     <factor>            <character> <character> <numeric> <character>            <character>            <character> <character> <numeric> <character> <numeric> <character> <numeric> <integer>                <list> <list>      <list>
#>    MantaDEL:245251:6:6:0:0:0_bp2       13  21729832      - |           NA TCTGCAACAGATACAAATAA..           T       999        PASS MantaDEL:245251:6:6:.. MantaDEL:245251:6:6:..         DEL      -542                     0        <NA>         1    176912 uc001unt.3,uc001unv.3 176912        SKA3
#>    MantaDEL:245251:5:8:0:0:0_bp2       13  21732061      - |           NA                      G       <DEL>       999        PASS MantaDEL:245251:5:8:.. MantaDEL:245251:5:8:..         DEL     -2110        <NA>         0        <NA>         2    176913 uc001unt.3,uc001unv.3 176913        SKA3
#>    MantaDEL:245251:5:9:0:0:0_bp2       13  21734038      - |           NA                      A       <DEL>       525        PASS MantaDEL:245251:5:9:.. MantaDEL:245251:5:9:..         DEL     -1776        <NA>         0        <NA>         4    176914 uc001unt.3,uc001unv.3 176914        SKA3
#>   MantaDEL:245251:7:10:0:0:0_bp2       13  21735929      - |           NA                      T       <DEL>       539        PASS MantaDEL:245251:7:10.. MantaDEL:245251:7:10..         DEL     -1802        <NA>         0        <NA>         1    176915 uc001unt.3,uc001unv.3 176915        SKA3
#>   MantaDEL:245251:4:11:0:0:0_bp2       13  21742127      - |           NA                      A       <DEL>       999        PASS MantaDEL:245251:4:11.. MantaDEL:245251:4:11..         DEL     -6112        <NA>         0        <NA>         2    176916 uc001unt.3,uc001unv.3 176916        SKA3
#>                              ...      ...       ...    ... .          ...                    ...         ...       ...         ...                    ...                    ...         ...       ...         ...       ...         ...       ...       ...                   ...    ...         ...
#>    MantaDEL:245251:5:9:0:0:0_bp1       13  21732261      + |           NA                      A       <DEL>       525        PASS MantaDEL:245251:5:9:.. MantaDEL:245251:5:9:..         DEL     -1776        <NA>         0        <NA>         4    176913 uc001unt.3,uc001unv.3 176913        SKA3
#>   MantaDEL:245251:7:10:0:0:0_bp1       13  21734126      + |           NA                      T       <DEL>       539        PASS MantaDEL:245251:7:10.. MantaDEL:245251:7:10..         DEL     -1802        <NA>         0        <NA>         1    176914 uc001unt.3,uc001unv.3 176914        SKA3
#>   MantaDEL:245251:4:11:0:0:0_bp1       13  21736014      + |           NA                      A       <DEL>       999        PASS MantaDEL:245251:4:11.. MantaDEL:245251:4:11..         DEL     -6112        <NA>         0        <NA>         2    176915 uc001unt.3,uc001unv.3 176915        SKA3
#>    MantaDEL:245251:3:4:0:0:0_bp1       13  21742538      + |           NA                      A       <DEL>       999        PASS MantaDEL:245251:3:4:.. MantaDEL:245251:3:4:..         DEL     -3939        <NA>         0        <NA>         2    176916 uc001unt.3,uc001unv.3 176916        SKA3
#>    MantaDEL:245251:2:3:0:0:0_bp1       13  21746642      + |           NA                      T       <DEL>       999        PASS MantaDEL:245251:2:3:.. MantaDEL:245251:2:3:..         DEL     -3870        <NA>         0        <NA>         2    176917 uc001unt.3,uc001unv.3 176917        SKA3
#>   -------
#>   seqinfo: 25 sequences from an unspecified genome
#> 
#> $insSite
#> GRanges object with 2 ranges and 18 metadata columns:
#>                               seqnames    ranges strand | paramRangeID         REF             ALT      QUAL      FILTER               sourceId                partner      svtype     svLen      insSeq    insLen       event    HOMLEN  exons                   txs    rtFound rtFoundSum gene_symbol
#>                                  <Rle> <IRanges>  <Rle> |     <factor> <character>     <character> <numeric> <character>            <character>            <character> <character> <numeric> <character> <numeric> <character> <numeric> <list>                <list>     <list>  <logical>      <list>
#>   MantaBND:245251:0:3:0:0:0:0       13  21746762      + |           NA           T T[11:108585702[        49        PASS MantaBND:245251:0:3:.. MantaBND:245251:0:3:..         BND        NA                     0        <NA>         0 176918 uc001unt.3,uc001unu.3 TRUE,FALSE       TRUE        SKA3
#>   MantaBND:245251:0:3:0:0:0:1       11 108585702      - |           NA           T  ]13:21746762]T        49        PASS MantaBND:245251:0:3:.. MantaBND:245251:0:3:..         BND        NA                     0        <NA>         0   <NA>                  <NA>       <NA>       <NA>        <NA>
#>   -------
#>   seqinfo: 25 sequences from an unspecified genome

5 Visualising breakpoint pairs via circos plots

One way of visualising RT is by circos plots. Here we use the package circlize to demonstrate the visualisation of insertion site and exon-exon junctions.

To generate a simple circos plot of RT event with SKA3 transcript:

library(circlize)
rt_chr_prefix <- c(RT$SKA3$junctions, RT$SKA3$insSite)
seqlevelsStyle(rt_chr_prefix) <- "UCSC"
pairs <- breakpointgr2pairs(rt_chr_prefix)
pairs

To see supporting breakpoints clearly, we generate the circos plot according to the loci of event.

circos.initializeWithIdeogram(
    data.frame(V1=c("chr13", "chr11"),
               V2=c(21720000,108585000),
               V3=c(21755000,108586000),
               V4=c("q12.11","q24.3"),
               V5=c("gneg","gpos50")))
circos.genomicLink(as.data.frame(S4Vectors::first(pairs)), 
                   as.data.frame(S4Vectors::second(pairs)))

circos.clear()

6 SessionInfo

sessionInfo()
#> R version 4.3.0 RC (2023-04-13 r84269)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.17-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               LC_TIME=en_GB              LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] circlize_0.4.15                         dplyr_1.1.2                             TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.52.0                  AnnotationDbi_1.62.0                    svaRetro_1.6.0                          StructuralVariantAnnotation_1.16.0      VariantAnnotation_1.46.0                Rsamtools_2.16.0                        Biostrings_2.68.0                       XVector_0.40.0                          SummarizedExperiment_1.30.0             Biobase_2.60.0                          MatrixGenerics_1.12.0                   matrixStats_0.63.0                      rtracklayer_1.60.0                      GenomicRanges_1.52.0                    GenomeInfoDb_1.36.0                     IRanges_2.34.0                          S4Vectors_0.38.0                        BiocGenerics_0.46.0                     BiocStyle_2.28.0                       
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.0         blob_1.2.4               filelock_1.0.2           bitops_1.0-7             fastmap_1.1.1            RCurl_1.98-1.12          BiocFileCache_2.8.0      GenomicAlignments_1.36.0 XML_3.99-0.14            digest_0.6.31            lifecycle_1.0.3          KEGGREST_1.40.0          RSQLite_2.3.1            magrittr_2.0.3           compiler_4.3.0           rlang_1.1.0              sass_0.4.5               progress_1.2.2           tools_4.3.0              utf8_1.2.3               yaml_2.3.7               knitr_1.42               prettyunits_1.1.1        bit_4.0.5                curl_5.0.0               DelayedArray_0.26.0      xml2_1.3.3               BiocParallel_1.34.0      grid_4.3.0               fansi_1.0.4              colorspace_2.1-0         biomaRt_2.56.0           cli_3.6.1                rmarkdown_2.21           crayon_1.5.2             generics_0.1.3           httr_1.4.5               rjson_0.2.21             DBI_1.1.3               
#> [40] cachem_1.0.7             stringr_1.5.0            zlibbioc_1.46.0          assertthat_0.2.1         parallel_4.3.0           BiocManager_1.30.20      restfulr_0.0.15          vctrs_0.6.2              Matrix_1.5-4             jsonlite_1.8.4           bookdown_0.33            hms_1.1.3                bit64_4.0.5              magick_2.7.4             jquerylib_0.1.4          glue_1.6.2               codetools_0.2-19         shape_1.4.6              stringi_1.7.12           BiocIO_1.10.0            tibble_3.2.1             pillar_1.9.0             rappdirs_0.3.3           htmltools_0.5.5          GenomeInfoDbData_1.2.10  BSgenome_1.68.0          R6_2.5.1                 dbplyr_2.3.2             evaluate_0.20            lattice_0.21-8           highr_0.10               png_0.1-8                memoise_2.0.1            bslib_0.4.2              Rcpp_1.0.10              xfun_0.39                GlobalOptions_0.1.2      pkgconfig_2.0.3