# automatic RNA-Seq present/absent gene expression calls generation

#### Julien Wollbrett, Marc Robinson-Rechavi, Frederic Bastian

#### 2020-01-17

`BgeeCall` is a collection of functions that uses Bgee expertise to create RNA-Seq gene expression present/absent calls.

The `BgeeCall` package allows to:

- Generate calls of presence/absence of expression at the gene level. You can generate these calls for any RNA-Seq samples as long as reference intergenic sequences have been generated by the Bgee team or by the Bgee community.
- Download reference intergenic sequences for species available in Bgee.
- Generate calls of presence/absence of expression at the transcript level (beta version)

If you find a bug or have any issues with `BgeeCall` please write a bug report in our GitHub issues manager.

## How present/absent calls are generated

In Bgee present/absent gene expression calls for RNA-seq are generated using a threshold specific to each RNA-Seq library, calculated using reads mapped to reference intergenic regions. This is unlike the more usual use of an arbitrary threshold below which a gene is not considered as expressed (e.g log2(TPM) = 1).

### Bgee database

Bgee is a database to retrieve and compare gene expression patterns in multiple animal species and produced from multiple data types (RNA-Seq, Affymetrix, in situ hybridization, and EST data). It notably integrates RNA-Seq libraries for 29 species.

### Reference intergenic regions

Reference intergenic regions are defined in the Bgee RNA-Seq pipeline. Candidate intergenic regions are defined using gene annotation data. For each species, over all available libraries, reads are mapped to these intergenic regions with kallisto, as well as to genes. This “intergenic expression” is deconvoluted to distinguish reference intergenic from non annotated genes, which have higher expression. Reference intergenic regions are then defined as intergenic regions with low expression level over all RNA-Seq libraries, relative to genes. This step allows not to consider regions wrongly considered as intergenic because of potential gene annotation quality problem as intergenic. For more information please refer to the Bgee RNA-Seq pipeline.

### Threshold of present/absent

BgeeCall pipeline allows to download reference intergenic regions resulting from the expertise of the Bgee team. Moreover BgeeCall allows to use these reference intergenic regions to automatically generate gene expression calls for your own RNA-Seq libraries as long as the species is integrated to Bgee The present/absent abundance threshold is calculated for each library using the formula :

\[ \frac {proportion\ of\ reference\ intergenic\ present}{proportion\ of\ protein\ coding\ present} = 0.05 \]

## Installation

In R:

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

## How to use the BgeeCall package

BgeeCall is highly tunable. Do not hesitate to have a look at the reference manual to have a precise descripton of all slots of the 4 main S4 classes (AbundanceMetadata, KallistoMetadata, BgeeMetadata and UserMetadata) or of all available functions. BgeeCall needs kallisto to run. If you do not have kallisto installed you will find more information how to install it here

### Load the package

    library(BgeeCall)

### Quick start

With the BgeeCall package it is easy to generate present/absent gene expression calls. The most time comsuming task of this calls generation is the generation of the kallisto transcriptome index. As the time needed for this step depend on the size of the transcriptome, we choose, as an example, the smallest transcriptome file among all species available on Bgee (C. elegans). To generate these calls you will need :

- a transcriptome
- gene annotations
- your RNA-Seq reads in fastq files

For this vignette we created a toy fastq file example based on the SRX099901 library using the ShortRead R package

    library("ShortRead")
    # keep 48.000 reads
    sampler <- FastqSampler(file.path("absolute_path","/SRX099901/SRR350955.fastq.gz"), 48000)
    set.seed(1); SRR350955 <- yield(sampler)
    writeFastq(object = SRR350955,
               file =file.path( "absolute_path","SRX099901_subset", "SRR350955_subset.fastq.gz"),
               mode = "w", full = FALSE, compress = TRUE)

In this example we used the Bioconductor AnnotationHub to load transcriptome and gene annotations but you can load them from wherever you want.

    ah <- AnnotationHub::AnnotationHub()
    ah_resources <- AnnotationHub::query(ah, c("Ensembl", "Caenorhabditis elegans", "84"))
    annotation_object <- ah_resources[["AH50789"]]
    transcriptome_object <- rtracklayer::import.2bit(ah_resources[["AH50453"]])

Once you have access to transcriptome, gene annotations and your RNA-Seq library, an object of class `UserMetadata` has to be created.

    # create an object of class UserMetadata and specify the species ID
    user_BgeeCall <- new("UserMetadata", species_id = "6239")
    # import annotation and transcriptome in the user_BgeeCall object
    # it is possible to import them using an S4 object (GRanges, DNAStringSet) or a file (gtf, fasta)
    user_BgeeCall <- setAnnotationFromObject(user_BgeeCall, annotation_object, "WBcel235_84")
    user_BgeeCall <- setTranscriptomeFromObject(user_BgeeCall, transcriptome_object, "WBcel235")
    # provide path to the directory of your RNA-Seq library
    user_BgeeCall <- setRNASeqLibPath(user_BgeeCall,
                                      system.file("extdata",
                                                  "SRX099901_subset",
                                                  package = "BgeeCall"))

And that’s it… You can run the generation of your present/absent gene expression calls

    calls_output <- generate_calls_workflow(userMetadata = user_BgeeCall)

    #> Querying Bgee to get intergenic release information...
    #> Note: importing `abundance.h5` is typically faster than `abundance.tsv`
    #> reading in files with read_tsv
    #> 1
    #> summarizing abundance
    #> summarizing counts
    #> summarizing length
    #> Note: importing `abundance.h5` is typically faster than `abundance.tsv`
    #> reading in files with read_tsv
    #> 1
    #> summarizing abundance
    #> summarizing counts
    #> summarizing length
    #> Generate present/absent expression calls.
    #>
    #> TPM cutoff for which 95% of the expressed genes are coding found at TPM = 4.64724e-06

Each analyze generates 5 files and return path to each one of them.

- calls_tsv_path : path to main tsv file with TPM, count, length, biotype, type, and presence/absence of expression summarized at gene level (or at the transcript level if it was requested)

    head.DataTable(x = read.table(calls_output$calls_tsv_path, header = TRUE), n = 5)
    #>               id abundance counts  length        biotype  type    call
    #> 1 WBGene00000001   27.0993      2 1556.20 protein_coding genic present
    #> 2 WBGene00000002    0.0000      0 1761.00 protein_coding genic  absent
    #> 3 WBGene00000003    0.0000      0 1549.00 protein_coding genic  absent
    #> 4 WBGene00000004   13.8730      1 1519.92 protein_coding genic present
    #> 5 WBGene00000005    0.0000      0 1487.00 protein_coding genic  absent

- cutoff_info_file_path : path to tsv summary of the analyze containing the proportion of gene, protein coding and intergenic defined as expressed. It also contains the library ID and the present/absent TPM threshold

    read.table(calls_output$cutoff_info_file_path)
    #>                              V1                V2
    #> 1                     libraryId  SRX099901_subset
    #> 2                     cutoffTPM       4.64724e-06
    #> 3        proportionGenicPresent  25.3210509597495
    #> 4            numberGenicPresent              5580
    #> 5                   numberGenic             22037
    #> 6       proportionCodingPresent  27.1433462121583
    #> 7           numberPresentCoding              5550
    #> 8                  numberCoding             20447
    #> 9   proportionIntergenicPresent 0.490998363338789
    #> 10      numberIntergenicPresent                21
    #> 11             numberIntergenic              4277
    #> 12 ratioIntergenicCodingPresent              0.05

- abundance_tsv : path to tsv kallisto quant output file

    head.DataTable(x = read.table(calls_output$abundance_tsv, header = TRUE), n = 5)
    #>    target_id length eff_length est_counts     tpm
    #> 1 Y110A7A.10   1787    1556.20          2 27.0993
    #> 2    F27C8.1   1940    1761.00          0  0.0000
    #> 3    F07C3.7   1728    1549.00          0  0.0000
    #> 4    F52H2.2   1739    1519.92          1 13.8730
    #> 5 T13A10.10a   1734    1555.00          0  0.0000

- TPM_distribution_path : path to plot in pdf reprensting density distribution of TPM values for all sequences, protein coding sequences, and intergenic sequences. The grey line corresponds to TPM threshold used to generate present/absent calls.

    openPDF(calls_output$TPM_distribution_path)

- S4_slots_summary : path to tsv file containing a summary of values used for the most important slots of the three S4 classes (UserMetadata, KallistoMetadata, and BgeeMetadata).

    read.table(calls_output$S4_slots_summary, header = TRUE, sep = "\t")
    #>                                     Slot.name
    #> 1                 AbundanceMetadata@tool_name
    #> 2                     AbundanceMetadata@txOut
    #> 3           AbundanceMetadata@ignoreTxVersion
    #> 4                    AbundanceMetadata@cutoff
    #> 5  AbundanceMetadata@read_size_kmer_threshold
    #> 6             BgeeMetadata@intergenic_release
    #> 7                     UserMetadata@species_id
    #> 8                     UserMetadata@reads_size
    #> 9                UserMetadata@rnaseq_lib_path
    #> 10            UserMetadata@transcriptome_name
    #> 11               UserMetadata@annotation_name
    #> 12           UserMetadata@simple_arborescence
    #> 13                                 output_dir
    #>                                                                                                                           Slot.value
    #> 1                                                                                                                           kallisto
    #> 2                                                                                                                              FALSE
    #> 3                                                                                                                              FALSE
    #> 4                                                                                                                               0.05
    #> 5                                                                                                                                 50
    #> 6                                                                                                                                0.1
    #> 7                                                                                                                               6239
    #> 8                                                                                                                                 51
    #> 9                             C:/Users/biocbuild/bbs-3.11-bioc/tmpdir/Rtmp6HFiwZ/Rinst317429f2590b/BgeeCall/extdata/SRX099901_subset
    #> 10                                                                                                                          WBcel235
    #> 11                                                                                                                       WBcel235_84
    #> 12                                                                                                                              TRUE
    #> 13 C:/Users/biocbuild/bbs-3.11-bioc/tmpdir/Rtmp6HFiwZ/Rinst317429f2590b/BgeeCall/extdata/intergenic_0.1/all_results/SRX099901_subset

### Generate present/absent calls for more than one RNA-Seq library

You will potentialy be also interested to generate present/absent calls on different RNA-Seq libraries, potentially on different species, or using The main function `generate_presence_absence()` allows to generate present/absent calls from a UserMetadata object but also from a data frame or a tsv file depending on the arguments of the function you use. Please choose one of the three following arguments : - userMetadata : Allows to generate present/absent calls for one RNA-Seq library using one object of the class UserMetadata.
- userDataFrame : Provide a dataframe where each row correspond to one present/absent call generation. It allows to generate present/absent calls on different libraries, species, transcriptome, genome annotations, etc. - userFile : Similar to userDataFrame except that the information are stored in a tsv file. A template of this file called `userMetadataTemplate.tsv` is available at the root of the package.

Columns of the dataframe or the tsv file are :

- species_id : The NCBI ID of the species.
- run_ids : The runs of the RNA-Seq library you want to use for the generation of the calls. Allows to generate expression calls for a subset of the runs of one RNA-Seq library as described in Generate calls for a subset of RNA-Seq runs. If not interested by this option, leave the column empty.
- reads_size : The size of the reads of your RNA-Seq library.
- rnaseq_lib_path : Path to the directory containing all fastq files generated for this library. This directory can only contains single-end runs or paired-end runs.
- transcriptome_path : path to the transcriptome file.
- annotation_path : path to the genome annotation file. Works with GTF of GFF3 files.
- working_path : path to the working directory where results will be stored. Using the same working directory for different RNA-Seq libraries of the same species will allow to reuse previously generated data like the custom transcriptome index (generated from both transcriptome and reference intergenic sequences). By default the working path is defined by the `getwd()` function and correspond to the working directory of your R session. If not interested by this option, leave the column empty.
- output_directory : Both species results and RNA-Seq libraries results are by default stored at the same place using the value of the `working_path` column. However, this column allows you to define a different output_directory for RNA-Seq results. For instance it allows you to save calls information directly in the RNA-Seq directory. If not interested by this option, leave the column empty.

Once the file has been fill in expression calls can be generated with :

    calls_output <- generate_calls_workflow(userFile = "path_to_your_file.tsv")

### Reference intergenic sequences

#### Releases of reference intergenic sequences

Different releases of reference intergenic sequences are available. It is possible to list all these releases :

    list_intergenic_release()
    #> Downloading release information of reference intergenic sequences...
    #>     release releaseDate                             FTPURL
    #> 1       0.1  2018-12-21 ftp://ftp.bgee.org/intergenic/0.1/
    #> 2       0.2  2019-02-07 ftp://ftp.bgee.org/intergenic/0.2/
    #> 3 community  2019-07-22
    #> 4    custom  2019-07-22
    #>                                                    referenceIntergenicFastaURL
    #> 1 ftp://ftp.bgee.org/intergenic/0.1/ref_intergenic/SPECIES_ID_intergenic.fa.gz
    #> 2 ftp://ftp.bgee.org/intergenic/0.2/ref_intergenic/SPECIES_ID_intergenic.fa.gz
    #> 3
    #> 4
    #>   minimumVersionBgeeCall
    #> 1                  0.9.9
    #> 2                  0.9.9
    #> 3                  1.1.0
    #> 4                  1.1.0
    #>                                                                                                                                                                                                     description
    #> 1                                                                                                                                                                  intergenic regions used to generate Bgee 14.
    #> 2                                                                          cleaned intergenic sequences based on release 0.1 (remove blocks of Ns longer than 100 and sequences containing more than 5% of Ns).
    #> 3                                                                           Release allowing to access to all reference intergenic sequences generated by the community and not present in Bgee for the moment.
    #> 4 Release allowing to use your own FASTA reference intergenic sequences. When this release is selected BgeeCall will use the sequences at UserMetadata@custom_intergenic_path to generate present/absent calls.
    #>                                                                            messageToUsers
    #> 1
    #> 2                         be careful, this intergenic release has not been tested by Bgee
    #> 3 These reference intergenic sequences have not been generated by Bgee. Use with caution.
    #> 4                              You decided to use your own reference intergenic sequences

It is then possible to choose one specific release to create a `BgeeMetadata` object. Always use the setter method `setIntergenicRelease()` when changing the release of an already existing BgeeMetadata object.

    # create BgeeMetadata object and define one reference intergenic release
    bgee <- new("BgeeMetadata", intergenic_release = "0.1")
    #> Querying Bgee to get intergenic release information...
    # change the reference intergenic release of your BgeeMetadata object
    bgee <- setIntergenicRelease(bgee, "0.2")
    #> IMPORTANT : be careful, this intergenic release has not been tested by Bgee

By default the reference intergenic release used when a `BgeeMetadata` object is created is the last stable one created by the Bgee team.

#### Core reference intergenic from Bgee

`Core reference intergenic` releases are created by the Bgee team when a lot of new RNA-Seq libraries have been manually curated for already existing species and/or for new species. These releases are the only ones with a release number (e.g “0.1”). Each of these releases contains reference intergenic sequences for a list of species. Bgee reference intergenic sequences have been generated using Bgee team expertise. The RNA-Seq libraries were manually curated as healthy and wild type. Quality Control have been done along all steps of generation of these sequences. Reference intergenic sequences have been selected from all potential intergenic regions (see Bgee pipeline). BgeeCall allows to generate gene expression call from Bgee reference intergenic sequences for any RNA-Seq libraries as long as these sequences have been generated by the Bgee team. A tsv file containing all species available for current release of reference intergenic is available here. This file also contains a column describing the number of RNA-Seq libraries used to generated the reference intergenic sequences of each species. It is also possible to list in R all species for which Bgee reference intergenic sequences have been created :

    list_bgee_ref_intergenic_species(myBgeeMetadata = bgee)
    #>    speciesId             speciesName numberOfLibraries   genomeVersion
    #> 1       9606             Homosapiens              5026       GRCh38.p5
    #> 2      10090             Musmusculus               133       GRCm38.p4
    #> 3       9544           Macacamulatta                90         MMUL1.0
    #> 4       7955              Daniorerio                67          GRCz10
    #> 5       8364       Xenopustropicalis                66          JGI4.2
    #> 6       6239   Caenorhabditiselegans                50        WBcel235
    #> 7       9031            Gallusgallus                45         Galgal4
    #> 8      10116        Rattusnorvegicus                36        Rnor_6.0
    #> 9       9913               Bostaurus                33          UMD3.1
    #> 10     13616    Monodelphisdomestica                19         monDom5
    #> 11      9258 Ornithorhynchusanatinus                17           OANA5
    #> 12      7240      Drosophilasimulans                17 GCA_000259055.1
    #> 13      9598          Pantroglodytes                15      CHIMP2.1.4
    #> 14      7237 Drosophilapseudoobscura                14 GCA_000001765.2
    #> 15      7227  Drosophilamelanogaster                14           BDGP6
    #> 16      9593          Gorillagorilla                13       gorGor3.1
    #> 17      9597             Panpaniscus                12      CHIMP2.1.4
    #> 18      9823               Susscrofa                10     Sscrofa10.2
    #> 19     10141          Caviaporcellus                 9 Felis_catus_6.2
    #> 20      9685              Feliscatus                 9         cavPor3
    #> 21      7230    Drosophilamojavensis                 8         EquCab2
    #> 22      9796           Equuscaballus                 8 GCA_000005175.1
    #> 23      9986    Oryctolaguscuniculus                 6         eriEur1
    #> 24      9615    Canislupusfamiliaris                 6       CanFam3.1
    #> 25      9365      Erinaceuseuropaeus                 6       OryCun2.0
    #> 26      7244       Drosophilavirilis                 4 GCA_000005245.1
    #> 27     28377      Anoliscarolinensis                 4       AnoCar2.0
    #> 28      7217     Drosophilaananassae                 4 GCA_000005975.1
    #> 29      7245        Drosophilayakuba                 4 GCA_000005115.1

#### Community reference intergenic

If you want to use BgeeCall on a species for which Bgee does not provide reference intergenic sequences you have the possibility to create them by yourself and share them with the Bgee community by following all steps of this tutorial. Do not forget that the number of RNA-Seq libraries is a key point to the generation of precise reference intergenic sequences. It is possible to list in R all species for which reference intergenic sequences have been created by the community using the following code

    list_community_ref_intergenic_species()
    #>   speciesId numberOfLibraries annotationVersion genomeVersion kallistoVersion
    #> 1     13686               243            Si_gnG        Si_gnG          0.44.0
    #>                                                                                      url
    #> 1 https://zenodo.org/api/files/5492ff2f-91a3-4101-8d67-78b8f8625cc6/ref_intergenic.fa.gz

If reference intergenic sequences of the species you are interested in are available only from the community release it is then possible to use this release to generate your present/absent calls

    # create a BgeeMetadata object using the community release
    bgee <- new("BgeeMetadata", release = "community")
    calls_output <- generate_calls_workflow(bgeeMetadata = bgee, userMetadata = user_BgeeCall)

#### Your own reference intergenic

If you generated your own reference intergenic sequences follwowing this tuorial but did not share them for the moment (do not forget to do it…), it is also possible to use BgeeCall with a file containing the sequences. In this case you need to select the custom release and provide the path to the file containing reference intergenic sequences :

    bgee <- new("BgeeMetadata", release = "community", custom_intergenic_path = "path/to/custom/ref_intergenic.fa.gz")
    calls_output <- generate_calls_workflow(bgeeMetadata = bgee, userMetadata = user_BgeeCall)

### Generate present/absent calls at transcript level (beta version)

kallisto generates TPMs at the transcript level. In the Bgee pipeline we summarize this expression at the gene level to calculate our present/absent calls. In `BgeeCall` it is now possible to generate present/absent calls at the transcript level. Be careful when using this feature as it has not been tested for the moment. To generate such calls you only have to create one object of the class `KallistoMetadata` and edit the value of one attribute

    kallisto <- new("KallistoMetadata", txOut = TRUE)
    calls_output <- generate_calls_workflow(myAbundanceMetadata = kallisto, userMetadata = user_BgeeCall)

### Tune how to use kallisto

#### Download or reuse your own kallisto

By default `BgeeCall` suppose that kallisto is installed. If kallisto is not installed on your computer you can either :

- let BgeeCall automatically download the version 0.45 of kallisto. BgeeCall will use it to quantify abundance of transcripts. It will only be used by this package and will have no impact on your potential already existing version of kallisto.

    kallisto <- new("KallistoMetadata", install_kallisto = TRUE)
    calls_output <- generate_calls_workflow(myAbundanceMetadata = kallisto, userMetadata = user_BgeeCall)

- You can download and install kallisto by following the instructions on the official website : http://pachterlab.github.io/kallisto/download

#### Edit kallisto quant attributes

By default kallisto is run with the same parameters that we use in the RNA-Seq Bgee pipeline:

- single end : “-t 1 –single -l 180 -s 30 –bias”
- paired end : “-t 1 –bias”

It is possible to modify them and use your favourite kallisto parameters

    kallisto <- new("KallistoMetadata", single_end_parameters = "-t 3 --single -l 150 -s 30", pair_end_parameters = "-t 2 -b --seed 36")
    calls_output <- generate_calls_workflow(myAbundanceMetadata = kallisto, userMetadata = user_BgeeCall)

#### Choose between two kmer size

By default 2 indexes with 2 different kmer sizes can be used by `BgeeCall` The default kmer size of kallisto (31) is used for libraries with reads length equal or larger than 50 bp. A kmer size of 15 is used for libraries with reads length smaller than 50 bp. We decided not to allow to tune kmers size because the generation of the index is time consuming and index generation takes even more time with small kmers size (< 15bp). However it is possible to modify the threshold of read length allowing to choose between default and small kmer size.

    # libraries with reads smaller than 70bp will use the index with kmer size = 15
    kallisto <- new("KallistoMetadata", read_size_kmer_threshold = 70)
    calls_output <- generate_calls_workflow(myAbundanceMetadata = kallisto, userMetadata = user_BgeeCall)

### Generate calls for a subset of RNA-Seq runs

By default gene expression calls are generated using all runs of the RNA-Seq library. It is possible to select only a subset of these runs.

    # RNA-Seq run SRR350955_subsetof from the RNA-Seq library will be used to generate the calls
    user_BgeeCall <- setRunIds(user_BgeeCall, c("SRR350955_subset"))
    calls_output <- run_from_object(myUserMetadata = user_BgeeCall)

When run IDs are selected, the name output directory combine the library ID and all selected run IDs. In our example the expression calls will be stored in the directory `SRX099901_SRR350955_subset`.

### Modify present/absent threshold

By default the threshold of present/absent is calculated with the formula :

    proportion of ref intergenic present / proportion of protein coding present = 0.05

This 0.05 corresponds to the ratio used in the Bgee pipeline. However it is possible to edit this value. Be careful when editing this value as it has a big impact on your present absent.

    kallisto <- new("KallistoMetadata", cutoff = 0.1)

### Generate calls with a simple arborescence of directories

By default the arborescence of directories created by `BgeeCall` is as simple as possible. the results will be created using the path `working_path/intergenic_release/all_results/libraryId`. Generating present/absent gene expression calls for the same RNA-Seq library using different transcriptome or annotation versions using this arborescence will overwrite previous results. The `UserMetadata` class has an attribute `simple_arborescence` that is `TRUE` by default. If `FALSE`, a complexe arborescence of directories containing the name of the annotation and transcriptome files will be created. This complex arborescence will then allow to generate present/absent calls for the same library using different version of transcriptome or annotaiton.

    user_BgeeCall <- setSimpleArborescence(userObject = user_BgeeCall, simpleArborescence = FALSE)
    calls_output <- run_from_object(myUserMetadata = user_BgeeCall)

### Change directory where calls are saved

By default directories used to save present/absent calls are subdirectories of `UserMetadata@working_path`. However it is possible to select the directory where you want the calls to be generated.

    user_BgeeCall <- setOutputDir(user_BgeeCall, "path/to/calls/for/this/library/")

This output directory will only contains results generated at the RNA-Seq library level. All data generated at species level are still stored using the `UserMetadata@working_path`. They can then still be reused to generate calls from other libraries of the same species.

#Session Info

    sessionInfo()
    #> R Under development (unstable) (2019-12-14 r77569)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server 2012 R2 x64 (build 9600)
    #>
    #> 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] parallel  stats4    stats     graphics  grDevices utils     datasets
    #> [8] methods   base
    #>
    #> other attached packages:
    #> [1] rtracklayer_1.47.0   GenomicRanges_1.39.1 GenomeInfoDb_1.23.7
    #> [4] IRanges_2.21.3       S4Vectors_0.25.10    BiocGenerics_0.33.0
    #> [7] BgeeCall_1.3.1
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] Biobase_2.47.2                httr_1.4.1
    #>  [3] bit64_0.9-7                   jsonlite_1.6
    #>  [5] AnnotationHub_2.19.5          shiny_1.4.0
    #>  [7] assertthat_0.2.1              interactiveDisplayBase_1.25.0
    #>  [9] askpass_1.1                   BiocManager_1.30.10
    #> [11] BiocFileCache_1.11.4          blob_1.2.0
    #> [13] BSgenome_1.55.2               GenomeInfoDbData_1.2.2
    #> [15] Rsamtools_2.3.2               yaml_2.2.0
    #> [17] progress_1.2.2                BiocVersion_3.11.1
    #> [19] pillar_1.4.3                  RSQLite_2.2.0
    #> [21] backports_1.1.5               lattice_0.20-38
    #> [23] glue_1.3.1                    digest_0.6.23
    #> [25] promises_1.1.0                XVector_0.27.0
    #> [27] htmltools_0.4.0               httpuv_1.5.2
    #> [29] Matrix_1.2-18                 XML_3.98-1.20
    #> [31] pkgconfig_2.0.3               biomaRt_2.43.1
    #> [33] zlibbioc_1.33.0               xtable_1.8-4
    #> [35] purrr_0.3.3                   later_1.0.0
    #> [37] BiocParallel_1.21.2           tibble_2.1.3
    #> [39] openssl_1.4.1                 SummarizedExperiment_1.17.1
    #> [41] GenomicFeatures_1.39.2        mime_0.8
    #> [43] magrittr_1.5                  crayon_1.3.4
    #> [45] memoise_1.1.0                 evaluate_0.14
    #> [47] tools_4.0.0                   prettyunits_1.1.0
    #> [49] hms_0.5.3                     matrixStats_0.55.0
    #> [51] stringr_1.4.0                 Rhdf5lib_1.9.0
    #> [53] DelayedArray_0.13.2           AnnotationDbi_1.49.0
    #> [55] Biostrings_2.55.4             compiler_4.0.0
    #> [57] rlang_0.4.2                   rhdf5_2.31.3
    #> [59] grid_4.0.0                    RCurl_1.95-4.13
    #> [61] tximport_1.15.6               rappdirs_0.3.1
    #> [63] bitops_1.0-6                  rmarkdown_2.0
    #> [65] DBI_1.1.0                     curl_4.3
    #> [67] R6_2.4.1                      GenomicAlignments_1.23.1
    #> [69] knitr_1.27                    dplyr_0.8.3
    #> [71] fastmap_1.0.1                 bit_1.1-15.1
    #> [73] zeallot_0.1.0                 readr_1.3.1
    #> [75] stringi_1.4.5                 Rcpp_1.0.3
    #> [77] vctrs_0.2.1                   dbplyr_1.4.2
    #> [79] tidyselect_0.2.5              xfun_0.12