“Finding differentially expressed unannotated genomic regions from RNA-seq data with srnadiff”

Matthias Zytnicki and Ignacio González

Package

srnadiff 1.14.0

1 Version info

• R version: R version 4.1.1 (2021-08-10)
• Bioconductor version: 3.14
• Package version: 1.14.0

2 Abstract

srnadiff is an R package that finds differently expressed regions from RNA-seq data at base-resolution level without relying on existing annotation. To do so, the package implements the identify-then-annotate methodology that builds on the idea of combining two pipelines approach: differential expressed regions detection and differential expression quantification.

3 Introduction

There is no real method for finding differentially expressed short RNAs. The most used method focusses on miRNAs, and only uses a standard RNA-Seq pipe-line on these genes.

However, annotated tRF, siRNAs, piRNA, etc. are thus out of the scope of these analyses. Several ad hoc method have been used, and this package implements a unifying method, finding differentially expressed genes or regions of any kind.

The srnadiff package implements two major methods to produce potential differentially expressed regions: the HMM and IR method. Briefly, these methods identify contiguous base-pairs in the genome that present differential expression signal, then these are regrouped into genomic intervals called differentially expressed regions (DERs).

Once DERs are detected, the second step in a sRNA-diff approach is to quantify the signification of these. To do so, reads (including fractions of reads) that overlap each expressed region are counted to arrive at a count matrix with one row per region and one column per sample. Then, this count matrix is analyzed using the standard workflow of DESeq2 for differential expression of RNA-seq data, assigning a p-value to each candidate DER.

The main functions for finds differently expressed regions are srnadiffExp and srnadiff. The first one creates an S4 class providing the infrastructure (slots) to store the input data, methods parameters, intermediate calculations and results of an sRNA-diff approach. The second one implement four methods to find candidate differentially expressed regions and quantify the statistic signification of the finded regions.

This vignette explains the basics of using srnadiff by showing an example, including advanced material for fine tuning some options. The vignette also includes description of the methods behind the package.

3.1 Citing srnadiff

We hope that srnadiff will be useful for your research. Please use the following information to cite srnadiff and the overall approach when you publish results obtained using this package, as such citation is the main means by which the authors receive credit for their work. Thank you!

Zytnicki, M., and I. González. (2021). “Finding differentially expressed sRNA-Seq regions with srnadiff.” .

3.2 How to get help for srnadiff

Most questions about individual functions will hopefully be answered by the documentation. To get more information on any specific named function, for example MIMFA, you can bring up the documentation by typing at the R.

help("srnadiff")

or

?srnadiff

The authors of srnadiff always appreciate receiving reports of bugs in the package functions or in the documentation. The same goes for well-considered suggestions for improvements. If you’ve run into a question that isn’t addressed by the documentation, or you’ve found a conflict between the documentation and what the software does, then there is an active community that can offer help. Send your questions or problems concerning srnadiff to the Bioconductor support site at .

Please send requests for general assistance and advice to the support site, rather than to the individual authors. It is particularly critical that you provide a small reproducible example and your session information so package developers can track down the source of the error. Users posting to the support site for the first time will find it helpful to read the posting guide at the Bioconductor help page.

3.3 Quick start

A typical sRNA-diff session can be divided into three steps:

1. Data preparation: In this first step, a convenient R object of class srnadiffExp is created containing all the information required for the two remaining steps. The user needs to provide a vector with the full paths to the BAM files (which should be coordinate sorted), a data.frame with sample and experimental design information and optionally annotated regions as a GRanges object.
2. Performing srnadiff: Using the object created in the first step the user can perform srnadiff to find potential DERs and quantify the statistic signification of these.
3. Visualization of the results: The DERs obtained in the second step are visualized by plotting the coverage information surrounding genomic regions.

A typical srnadiff session might look like the following. Here we assume that bamFiles is a vector with the full paths to the BAM files and the sample and experimental design information are stored in a data frame sampleInfo.

4 Using srnadiff

4.1 Installation

We assume that the user has the R program (see the R project) already installed.

The srnadiff package is available from the Bioconductor repository. To be able to install the package one needs first to install the core Bioconductor packages. If you have already installed Bioconductor packages on your system then you can skip the two lines below.

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

Once the core Bioconductor packages are installed, you can install the srnadiff package by

BiocManager::install("srnadiff")

Load the srnadiff package in your R session:

library(srnadiff)

A list of all accessible vignettes and methods is available with the following command:

help.search("srnadiff")

4.2 Data overview

To help demonstrate the functionality of srnadiff, the package includes datasets of published by (Viollet et al. 2015).

Briefly, these data consist of three replicates of sRNA-Seq of SLK (human) cell lines, and three replicates of SLK cell lines infected with Kaposi’s sarcoma associated herpesvirus. The analysis shows that several loci are repressed in the infected cell lines, including the 14q32 miRNA cluster.

Raw data have been downloaded from the GEO data set GSE62830. Adapters were removed with fastx_clipper and mapped with (Salzberg and Langmead 2012) on the human genome version GRCh38.

This data is restricted to a small locus on chr14. It uses the whole genome annotation (with coding genes, etc.) and extracts miRNAs.

The file dataInfo.csv contains three columns for each BAM file:

• the file name (FileName)
• a human readable sample name (SampleName)
• a condition, e.g. WT (Condition)

4.3 Data preparation: the srnadiffExp object

The first step in an sRNA-diff approach is to create a srnadiffExp object. srnadiffExp is an S4 class providing the infrastructure to store the input data, methods parameters, intermediate calculations and results of a sRNA-diff approach. This object will be also the input of the visualization function.

4.3.1 Content

The object srnadiffExp will usually be represented in the code here as srnaExp. To build such an object the user needs the following:

1. paths to the BAM files: a vector with the full paths to the sample BAM files, which should be coordinate sorted. Please use the sortBam function of the Rsamtools package to sort the reads.
2. sample information: a data.frame with three columns labelled FileName, SampleName and Condition. The first column is the BAM file name (without extension), the second column the sample name, and the third column the condition to which sample belongs. Each row describes one sample.
3. annotation: optionally annotation information. A GRanges object containing annotated regions.

Here, we demonstrate how to construct a srnadiffExp object from (Viollet et al. 2015).

## Specifiy path to data files
basedir <- system.file("extdata", package="srnadiff", mustWork=TRUE)

## Read sample information file, and create a data frame
sampleInfo <- read.csv(file.path(basedir, "dataInfo.csv"))

## Vector with the full paths to the BAM files to use
bamFiles <- file.path(basedir, sampleInfo$FileName)

## Creates an srnadiffExp object
srnaExp <- srnadiffExp(bamFiles, sampleInfo)

4.3.2 Adding an annotation

Optionally, if annotation information is available as a GRanges object, annotReg say, then a srnadiffExp object can be created by

srnaExp <- srnadiffExp(bamFiles, sampleInfo, annotReg)

or by

srnaExp <- srnadiffExp(bamFiles, sampleInfo)
annotReg(srnaExp) <- annotReg

4.3.3 Directionality

srnadiff chooses a reference condition, which is the control, or wild type condition. The results reflect the comparison of the other condition versus the reference condition. By default, the reference condition is the first condition, when the conditions are sorted by the alphabetical order. It is thus advisable to specify the reference condition.

As mentionned in DESeq2, you can select the reference condition using factors:

sampleInfo$Condition <- factor(sampleInfo$Condition,
                               levels = c("control", "infected"))

Internally, srnadiff (through DESeq2) computes the log-fold-change of the comparison. The log-fold-change is positive when the reference condition is less expressed than the other condition.

4.3.4 The srnaExp object

A summary of the srnaExp object can be seen by typing the object name at the R prompt

srnaExp

## Object of class srnadiffExp.
##  Sample information
##         FileName SampleName Condition
## 1 SRR1634756.bam      ctr_1   control
## 2 SRR1634757.bam      ctr_2   control
## 3 SRR1634758.bam      ctr_3   control
## 4 SRR1634759.bam      inf_1  infected
## 5 SRR1634760.bam      inf_2  infected
## 6 SRR1634761.bam      inf_3  infected

For your conveniance and illustrative purposes, an example of an srnadiffExp object can be loaded with an only command, so the script boils down to:

srnaExp <- srnadiffExample()

The srnadiffExp object in this example was constructed by:

basedir    <- system.file("extdata", package="srnadiff", mustWork=TRUE)
sampleInfo <- read.csv(file.path(basedir, "dataInfo.csv"))
gtfFile    <- file.path(basedir, "Homo_sapiens.GRCh38.76.gtf.gz")
annotReg   <- readAnnotation(gtfFile, feature="gene", source="miRNA")
bamFiles   <- file.path(basedir, sampleInfo$FileName)
srnaExp    <- srnadiffExp(bamFiles, sampleInfo, annotReg)

4.4 Read annotation

srnadiff offers the readAnnotation function related to loading annotation data. This accepts two annotation format files: GTF and GFF formats. Specification of GTF/GFF format can be found at the UCSC dedicated page.

readAnnotation reads and parses content of GTF/GFF files and stores annotated genomic features (regions) in a GRanges object. This has three main arguments: the first argument indicates the path, URL or connection to the GTF/GFF annotation file. The second and third argument, feature and source respectively, are of type character string, these specify the feature and attribute type used to select rows in the GTF/GFF annotation which will be imported. feature and source can be NULL, in this case, no selection is performed and all content into the file is imported.

4.4.1 Extraction of putative regions using an GTF annotation file

This method simply provides the genomic regions corresponding to the annotation file that is optionally given by the user. It can be a set of known miRNAs, siRNAs, piRNAs, genes, or a combination thereof.

4.4.2 Whole genome file annotation

This GTF file can be found in the central repositories (NCBI, Ensembl) and contains all the annotation found in an organism (coding genes, tranposable element, etc.). The following function reads the annotation file and extracts the miRNAs. Annotation files may have different formats, but this command has been tested on several model organisms (including human) from Ensembl.

gtfFile  <- file.path(basedir, "Homo_sapiens.GRCh38.76.gtf.gz")
annotReg <- readAnnotation(gtfFile, feature="gene", source="miRNA")

4.4.3 Extraction of precursor miRNAs using a miRBase-formatted GFF file

miRBase (Kozomara and Griffiths-Jones 2014) is the central repository for miRNAs. If your organism is available, you can download their miRNA annotation in GFF3 format (check the “Browse” tab). The following code parses a GFF3 miRBase file, and extracts the precursor miRNAs.

gffFile  <- file.path(basedir, "mirbase21_GRCh38.gff3")
annotReg <- readAnnotation(gffFile, feature="miRNA_primary_transcript")

4.4.4 Extraction of mature miRNAs using a miRBase-formatted GFF file

In the previous example, the reads will be counted per pre-miRNA, and the 5’ and 3’ arms, the miRNA and the miRNA* will be merged in the same feature. If you want to separate the two, use:

gffFile  <- file.path(basedir, "mirbase21_GRCh38.gff3")
annotReg <- readAnnotation(gffFile, feature="miRNA")

4.4.5 Other format

When the previous functions do not work, you can use your own parser with:

annotation <- readAnnotation(gtfFile, source="miRNA", feature="gene")

The source parameter keeps all the lines such that the second field matches the given parameter (e.g. miRNA). The feature parameter keeps all the lines such that the third field matches the given parameter (e.g. gene). The name of the feature will be given by the tag name (e.g. gene_name). source, feature and name can be NULL. In this case, no selection is performed on source or feature. If name is null, then a systematic name is given (annotation_N).

4.5 Performing sRNA-diff

The main function for performing an sRNA-diff analysis is srnadiff, this the wrapper for running several key functions from this package. srnadiff implement four methods to produce potential DERs: the annotation, naive, hmm and IR method (see bellow). Once potential DERs are detected, the second step in srnadiff is to quantify the statistic signification of these.

srnadiff has three main arguments. The first argument is an instance of class srnadiffExp. The second argument is of type character vector, it specify the segmentation methods to use, one of annotation, naive, hmm, IR or combinations thereof. The default all, all methods are used. The third arguments is of type list, it contain named components for the methods parameters to use. If missing, default parameter values are supplied. Details about the methods parameters are further described in the manual page of the parameters function and in Methods to produce differentially expressed regions section.

We then performs an sRNA-diff analysis on the input data contained in srnaExp by

srnaExp <- srnadiff(srnaExp)

srnadiff returns an object of class srnadiffExp again containing additional slots for:

• regions
• parameters
• countMatrix

5 Working with the srnadiffExp object

Once the srnadiffExp object is created the user can use the methods defined for this class to access the information encapsulated in the object.

By example, the sample information is accessed by

sampleInfo(srnaExp)

##         FileName SampleName Condition
## 1 SRR1634756.bam      ctr_1   control
## 2 SRR1634757.bam      ctr_2   control
## 3 SRR1634758.bam      ctr_3   control
## 4 SRR1634759.bam      inf_1  infected
## 5 SRR1634760.bam      inf_2  infected
## 6 SRR1634761.bam      inf_3  infected

For accessing the chromosomeSize slot

chromosomeSizes(srnaExp)

##        14 
## 107043718

The list of parameters can be exported by the function parameters

parameters(srnaExp)

5.1 Extracting regions

The regions, together with log2-fold-changes and adjusted p-values, can be extracted with this command:

regions <- regions(srnaExp, pvalue=0.5)

where padj is the (adjusted) p-value threshold. The output in a GenomicRanges object, and the information is accessible with the mcols() function.

You can export the regions to a BED file with the rtracklayer function export.

library(rtracklayer)
export(regions, "file.bed")

You can either export to GFF, GFF3, or GTF formats. Simply change the file suffix.

5.2 Data visualization

An insightful way of looking at the results of srnadiff is to investigate how the coverage information surrounding finded regions are distributed on the genomic coordinates.

plotRegions provides a flexible genomic visualization framework by displaying tracks in the sense of the Gviz package. Given a region (or regions), four separate tracks are represented:

1. `GenomeAxisTrack, a horizontal axis with genomic coordinate tickmarks for reference location to the displayed genomic regions;
2. GeneRegionTrack, if theannot` argument is passed, a track displaying all gene and/or sRNA annotation information in a particular region;
3. `AnnotationTrack, regions are plotted as simple boxes if no strand information is available, or as arrows to indicate their direction; and
4. `DataTrack, plot the sample coverages surrounding the genomic regions.

The sample coverages can be plotted in various different forms as well as combinations thereof. Supported plotting types are:

The sample coverages can be plotted in various different forms as well as combinations thereof. Supported plotting types are:

`p: simple dot plot;

`l: lines plot;

`b: combination of dot and lines plot;

`a: lines plot of the sample-groups average (i.e., mean) values;

`confint: confidence intervals for average values.

The default visualization for results from srnadiff is a lines plot of the sample-groups average.

plotRegions(srnaExp, regions(srnaExp)[1])

6 Methods behind srnadiff

6.1 Pre-processing data

As input, srnadiffExp expects BAM files as obtained, e.g., from RNA-Seq or another high-throughput sequencing experiment. Reading and processing of BAM files uses current Bioconductor infrastructure for processing sequencing reads: RSamtools, IRanges and GenomicRanges libraries. At this stage BAM files are summarized into base-resolution coverage and stored in a run-length encoding format in order to enhance computational performance. Run-length encoding is a compact way to store an atomic vector as a pairs of vectors (value, length). It is based on the rle function from the base R package.

As a second pre-processing step, srnadiffExp estimate the size factors (the effective library size) from the coverage data, such that count values in each sample coverage can be brought to a common scale by dividing by the corresponding size (normalization) factor. This step is also called normalization, its purpose is to render coverages (counts) from different samples, which may have been sequenced to different depths, comparable. Normalization factors are estimated using the median ratio method described by Equation 5 in (Anders and Huber 2010).

6.2 Methods to produce differentially expressed regions

6.2.1 HMM method: hmm

The first step in HMM method is quantifying the evidence for differential expression at the base-resolution level. To do this, srnadiff use the common approach in comparative analysis of transcriptomics data: test the null hypothesis that the logarithmic fold change between condition groups for a nucleotide expression is exactly zero.

The next step in the HMM approach enforces a smoothness assumption over the state of nucleotides: differential expression does not randomly switch along the chromosome, rather, continuous regions of RNA are either “differentially expressed” or “not”. This is captured with a hidden Markov model (HMM) with binary latent state corresponding to the true state of each nucleotide: differentially expressed or not differentially expressed.

The observations of the HMM are then the empirical p-values arising from the differential expression analysis corresponding to each nucleotide position. Modelling p-values directly enabled us to define the emission of each state as follows: the differentially expressed state emits a p-value \(< t\) with probability \(p\), and the not differentially expressed state emits a p-value \(\geqslant t\) with probability \(1-p\), where \(t\) is a real number between 0 and 1.

The HMM approach normally needs emission, transition, and starting probabilities values. They can be tuned by the user according to the overall p-values from differential analysis. We then run the Viterbi algorithm [ref] in order to finding the most likely sequence of states from the HMM. This essentially segments the genome into regions, where a region is defined as a set of consecutive bases showing a common expression signature. A region of bases with differentially expressed state is referred as an expressed region and is given as output of the method.

To run the HMM approach, srnadiff first form a large matrix, with rows corresponding to bases, columns corresponding to samples and entries are the coverage from a nucleotide of a particular sample. This count matrix is then analyzed as into feature-level counts using the feature-level RNA-seq differential expression analysis from DESeq2. In practice, the p-value is not computed for every nucleotide. Nucleotides for which the sum of the coverage across all samples is less than a threshold are given a p-value of 1, because these poorly expressed bases are unlikely to provide a differentially expressed sRNA.

The parameters for the HMM method are:

noDiffToDiff: Initial transition probability from “no differentially expressed” state to “differentially expressed” state.

diffToNoDiff: Initial transition probability from “differentially expressed” state to no “differentially expressed” state.

emission: Is the probability to emit a p-value \(<t\) in the “differentially expressed” state, and a p-value \(\geq t\) in the “not differentially expressed” state.

emissionThreshold: Is the threshold \(t\) that limits each state.

This parameters can be changed using using the assignment function parameters<-

parameters(srnaExp) <- list(noDiffToDiff=0.01, emissionThreshold=0.2)

6.2.2 IR method: IR

In this approach, for each base, the average from the normalized coverage is calculated across all samples into each condition. This generates a vector of (normalized) mean coverage expression per condition. These two vectors are then used to compute per-nucleotide log-ratios (in absolute value) across the genome. For the computed log-ratio expression, the method uses a sliding threshold h that run across the log-ratio levels identifying bases with log-ratio value above of h. Regions of contiguous bases passing this threshold are then analyzed using an adaptation of Aumann and Lindell algorithm for irreducibility property (Aumann and Lindell 2003).

The minimun sliding threshold, minLogFC, used in the IR method can be changed using the assignment function parameters<-

parameters(srnaExp) <- list(minLogFC=1)

6.2.3 Naive method: naive

This method is the simplest, gived a fixed threshold h, contiguous bases with log-ratio expression (in absolute value) passing this threshold are then considered as candidate differentially expressed regions.

The fixed threshold, cutoff, used in this method can be changed using the assignment function parameters<-

parameters(srnaExp) <- list(cutoff=1.5)

6.3 Quantifying DERs

The result of the ER step is a list of genomic regions which were chosen with a specific use, to quantify their expression for subsequent testing for differential expression. The selected regions are then quantified using the summarizeOverlaps function of the GenomicAlignments package. Notice that a read can overlap two different regions (e.g. extracted from the HMM and the IR methods), and thus can be counted twice for the quantification. The result is ultimately a matrix with rows corresponding to ERs and columns corresponding to samples; entries of this matrix are the number of aligned reads from a particular sample that overlap a particular region. Then, this count matrix is analyzed using the standard DESeq2 workflow for differential expression of RNA-seq data, assigning a p-value to each DER.

6.4 General parameters

The three last strategies can be tuned by specifying:

• the minimum and maximum regions sizes,
• the minimum depth of the most expressed condition.

The default values can be changed using these functions:

parameters(srnaExp) <- list(minDepth=1)
parameters(srnaExp) <- list(minSize=15, maxSize=1000)

6.5 Combination of strategies

6.5.1 Choice of the strategies

All the regions given by each strategies are then combined into a list of regions. You can choose not to use some strategies, use the parameter segMethod of the function srnadiff.

srnaExp <- srnadiffExample()
srnaExp <- srnadiff(srnaExp, segMethod=c("hmm", "IR"))

6.5.2 Quantification of the features

The selected regions are then quantified using of the summarizeOverlaps function of the GenomicAlignments package. Notice that a read can overlap two different regions (e.g. extracted from the naive and the slicing methods), and thus can be counted twice for the quantification.

You can adjust the minimum number of overlapping nucleotides between a read and a region to declare a hit, using:

parameters(srnaExp) <- list(minOverlap=1000)

DESeq2 is then used to get the adjusted p-values of these regions.

6.6 Using an other method to compute adjusted p-values

p-values are computed twice in the process: first, when preparing the HMM observed values; second, when comparing putative regions.

The defaut method uses DESeq2. You can choose instead two other methods: edgeR, or baySeq.

The method can be chosen using the parameter diffMethod of the function srnadiff:

srnaExp <- srnadiffExample()
srnaExp <- srnadiff(srnaExp, diffMethod="edgeR")

or

srnaExp <- srnadiffExample()
srnaExp <- srnadiff(srnaExp, diffMethod="baySeq")

## ......

The parameter diffMethod is case-insensitive.

7 Misc

7.1 Using several cores

The quantification and differential expression steps can be accelerated using several cores and the following command:

exp <- setNThreads(exp, nThreads=4)

7.2 Troubleshooting

While installing the package, if the compiler complains and says

#error This file requires compiler and library support for the ISO C++ 2011 standard.
This support is currently experimental, and must be enabled with the -std=c++11 or -std=gnu++11 compiler options.

Add this line

Sys.setenv("PKG_CXXFLAGS"="-std=c++11")

before installing the package.

8 Session information

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.1.1 (2021-08-10)
##  os       Ubuntu 20.04.3 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  C                           
##  ctype    en_US.UTF-8                 
##  tz       America/New_York            
##  date     2021-10-26                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version  date       lib source        
##  abind                  1.4-5    2016-07-21 [2] CRAN (R 4.1.1)
##  annotate               1.72.0   2021-10-26 [2] Bioconductor  
##  AnnotationDbi          1.56.0   2021-10-26 [2] Bioconductor  
##  AnnotationFilter       1.18.0   2021-10-26 [2] Bioconductor  
##  assertthat             0.2.1    2019-03-21 [2] CRAN (R 4.1.1)
##  backports              1.2.1    2020-12-09 [2] CRAN (R 4.1.1)
##  base64enc              0.1-3    2015-07-28 [2] CRAN (R 4.1.1)
##  baySeq                 2.28.0   2021-10-26 [2] Bioconductor  
##  Biobase                2.54.0   2021-10-26 [2] Bioconductor  
##  BiocFileCache          2.2.0    2021-10-26 [2] Bioconductor  
##  BiocGenerics         * 0.40.0   2021-10-26 [2] Bioconductor  
##  BiocIO                 1.4.0    2021-10-26 [2] Bioconductor  
##  BiocManager          * 1.30.16  2021-06-15 [2] CRAN (R 4.1.1)
##  BiocParallel           1.28.0   2021-10-26 [2] Bioconductor  
##  BiocStyle            * 2.22.0   2021-10-26 [2] Bioconductor  
##  biomaRt                2.50.0   2021-10-26 [2] Bioconductor  
##  Biostrings             2.62.0   2021-10-26 [2] Bioconductor  
##  biovizBase             1.42.0   2021-10-26 [2] Bioconductor  
##  bit                    4.0.4    2020-08-04 [2] CRAN (R 4.1.1)
##  bit64                  4.0.5    2020-08-30 [2] CRAN (R 4.1.1)
##  bitops                 1.0-7    2021-04-24 [2] CRAN (R 4.1.1)
##  blob                   1.2.2    2021-07-23 [2] CRAN (R 4.1.1)
##  bookdown               0.24     2021-09-02 [2] CRAN (R 4.1.1)
##  BSgenome               1.62.0   2021-10-26 [2] Bioconductor  
##  bslib                  0.3.1    2021-10-06 [2] CRAN (R 4.1.1)
##  cachem                 1.0.6    2021-08-19 [2] CRAN (R 4.1.1)
##  callr                  3.7.0    2021-04-20 [2] CRAN (R 4.1.1)
##  checkmate              2.0.0    2020-02-06 [2] CRAN (R 4.1.1)
##  cli                    3.0.1    2021-07-17 [2] CRAN (R 4.1.1)
##  cluster                2.1.2    2021-04-17 [2] CRAN (R 4.1.1)
##  colorspace             2.0-2    2021-06-24 [2] CRAN (R 4.1.1)
##  crayon                 1.4.1    2021-02-08 [2] CRAN (R 4.1.1)
##  curl                   4.3.2    2021-06-23 [2] CRAN (R 4.1.1)
##  data.table             1.14.2   2021-09-27 [2] CRAN (R 4.1.1)
##  DBI                    1.1.1    2021-01-15 [2] CRAN (R 4.1.1)
##  dbplyr                 2.1.1    2021-04-06 [2] CRAN (R 4.1.1)
##  DelayedArray           0.20.0   2021-10-26 [2] Bioconductor  
##  desc                   1.4.0    2021-09-28 [2] CRAN (R 4.1.1)
##  DESeq2                 1.34.0   2021-10-26 [2] Bioconductor  
##  devtools               2.4.2    2021-06-07 [2] CRAN (R 4.1.1)
##  dichromat              2.0-0    2013-01-24 [2] CRAN (R 4.1.1)
##  digest                 0.6.28   2021-09-23 [2] CRAN (R 4.1.1)
##  dplyr                  1.0.7    2021-06-18 [2] CRAN (R 4.1.1)
##  edgeR                  3.36.0   2021-10-26 [2] Bioconductor  
##  ellipsis               0.3.2    2021-04-29 [2] CRAN (R 4.1.1)
##  ensembldb              2.18.0   2021-10-26 [2] Bioconductor  
##  evaluate               0.14     2019-05-28 [2] CRAN (R 4.1.1)
##  fansi                  0.5.0    2021-05-25 [2] CRAN (R 4.1.1)
##  fastmap                1.1.0    2021-01-25 [2] CRAN (R 4.1.1)
##  filelock               1.0.2    2018-10-05 [2] CRAN (R 4.1.1)
##  foreign                0.8-81   2020-12-22 [2] CRAN (R 4.1.1)
##  Formula                1.2-4    2020-10-16 [2] CRAN (R 4.1.1)
##  fs                     1.5.0    2020-07-31 [2] CRAN (R 4.1.1)
##  genefilter             1.76.0   2021-10-26 [2] Bioconductor  
##  geneplotter            1.72.0   2021-10-26 [2] Bioconductor  
##  generics               0.1.1    2021-10-25 [2] CRAN (R 4.1.1)
##  GenomeInfoDb         * 1.30.0   2021-10-26 [2] Bioconductor  
##  GenomeInfoDbData       1.2.7    2021-09-23 [2] Bioconductor  
##  GenomicAlignments      1.30.0   2021-10-26 [2] Bioconductor  
##  GenomicFeatures        1.46.0   2021-10-26 [2] Bioconductor  
##  GenomicRanges        * 1.46.0   2021-10-26 [2] Bioconductor  
##  ggplot2                3.3.5    2021-06-25 [2] CRAN (R 4.1.1)
##  glue                   1.4.2    2020-08-27 [2] CRAN (R 4.1.1)
##  gridExtra              2.3      2017-09-09 [2] CRAN (R 4.1.1)
##  gtable                 0.3.0    2019-03-25 [2] CRAN (R 4.1.1)
##  Gviz                   1.38.0   2021-10-26 [2] Bioconductor  
##  highr                  0.9      2021-04-16 [2] CRAN (R 4.1.1)
##  Hmisc                  4.6-0    2021-10-07 [2] CRAN (R 4.1.1)
##  hms                    1.1.1    2021-09-26 [2] CRAN (R 4.1.1)
##  htmlTable              2.3.0    2021-10-12 [2] CRAN (R 4.1.1)
##  htmltools              0.5.2    2021-08-25 [2] CRAN (R 4.1.1)
##  htmlwidgets            1.5.4    2021-09-08 [2] CRAN (R 4.1.1)
##  httr                   1.4.2    2020-07-20 [2] CRAN (R 4.1.1)
##  IRanges              * 2.28.0   2021-10-26 [2] Bioconductor  
##  jpeg                   0.1-9    2021-07-24 [2] CRAN (R 4.1.1)
##  jquerylib              0.1.4    2021-04-26 [2] CRAN (R 4.1.1)
##  jsonlite               1.7.2    2020-12-09 [2] CRAN (R 4.1.1)
##  KEGGREST               1.34.0   2021-10-26 [2] Bioconductor  
##  knitr                * 1.36     2021-09-29 [2] CRAN (R 4.1.1)
##  lattice                0.20-45  2021-09-22 [2] CRAN (R 4.1.1)
##  latticeExtra           0.6-29   2019-12-19 [2] CRAN (R 4.1.1)
##  lazyeval               0.2.2    2019-03-15 [2] CRAN (R 4.1.1)
##  lifecycle              1.0.1    2021-09-24 [2] CRAN (R 4.1.1)
##  limma                  3.50.0   2021-10-26 [2] Bioconductor  
##  locfit                 1.5-9.4  2020-03-25 [2] CRAN (R 4.1.1)
##  magick                 2.7.3    2021-08-18 [2] CRAN (R 4.1.1)
##  magrittr               2.0.1    2020-11-17 [2] CRAN (R 4.1.1)
##  Matrix                 1.3-4    2021-06-01 [2] CRAN (R 4.1.1)
##  MatrixGenerics         1.6.0    2021-10-26 [2] Bioconductor  
##  matrixStats            0.61.0   2021-09-17 [2] CRAN (R 4.1.1)
##  memoise                2.0.0    2021-01-26 [2] CRAN (R 4.1.1)
##  munsell                0.5.0    2018-06-12 [2] CRAN (R 4.1.1)
##  nnet                   7.3-16   2021-05-03 [2] CRAN (R 4.1.1)
##  pillar                 1.6.4    2021-10-18 [2] CRAN (R 4.1.1)
##  pkgbuild               1.2.0    2020-12-15 [2] CRAN (R 4.1.1)
##  pkgconfig              2.0.3    2019-09-22 [2] CRAN (R 4.1.1)
##  pkgload                1.2.3    2021-10-13 [2] CRAN (R 4.1.1)
##  png                    0.1-7    2013-12-03 [2] CRAN (R 4.1.1)
##  prettyunits            1.1.1    2020-01-24 [2] CRAN (R 4.1.1)
##  processx               3.5.2    2021-04-30 [2] CRAN (R 4.1.1)
##  progress               1.2.2    2019-05-16 [2] CRAN (R 4.1.1)
##  ProtGenerics           1.26.0   2021-10-26 [2] Bioconductor  
##  ps                     1.6.0    2021-02-28 [2] CRAN (R 4.1.1)
##  purrr                  0.3.4    2020-04-17 [2] CRAN (R 4.1.1)
##  R6                     2.5.1    2021-08-19 [2] CRAN (R 4.1.1)
##  rappdirs               0.3.3    2021-01-31 [2] CRAN (R 4.1.1)
##  RColorBrewer           1.1-2    2014-12-07 [2] CRAN (R 4.1.1)
##  Rcpp                   1.0.7    2021-07-07 [2] CRAN (R 4.1.1)
##  RCurl                  1.98-1.5 2021-09-17 [2] CRAN (R 4.1.1)
##  remotes                2.4.1    2021-09-29 [2] CRAN (R 4.1.1)
##  restfulr               0.0.13   2017-08-06 [2] CRAN (R 4.1.1)
##  rjson                  0.2.20   2018-06-08 [2] CRAN (R 4.1.1)
##  rlang                  0.4.12   2021-10-18 [2] CRAN (R 4.1.1)
##  rmarkdown            * 2.11     2021-09-14 [2] CRAN (R 4.1.1)
##  rpart                  4.1-15   2019-04-12 [2] CRAN (R 4.1.1)
##  rprojroot              2.0.2    2020-11-15 [2] CRAN (R 4.1.1)
##  Rsamtools              2.10.0   2021-10-26 [2] Bioconductor  
##  RSQLite                2.2.8    2021-08-21 [2] CRAN (R 4.1.1)
##  rstudioapi             0.13     2020-11-12 [2] CRAN (R 4.1.1)
##  rtracklayer          * 1.54.0   2021-10-26 [2] Bioconductor  
##  S4Vectors            * 0.32.0   2021-10-26 [2] Bioconductor  
##  sass                   0.4.0    2021-05-12 [2] CRAN (R 4.1.1)
##  scales                 1.1.1    2020-05-11 [2] CRAN (R 4.1.1)
##  sessioninfo            1.1.1    2018-11-05 [2] CRAN (R 4.1.1)
##  srnadiff             * 1.14.0   2021-10-26 [1] Bioconductor  
##  stringi                1.7.5    2021-10-04 [2] CRAN (R 4.1.1)
##  stringr                1.4.0    2019-02-10 [2] CRAN (R 4.1.1)
##  SummarizedExperiment   1.24.0   2021-10-26 [2] Bioconductor  
##  survival               3.2-13   2021-08-24 [2] CRAN (R 4.1.1)
##  testthat               3.1.0    2021-10-04 [2] CRAN (R 4.1.1)
##  tibble                 3.1.5    2021-09-30 [2] CRAN (R 4.1.1)
##  tidyselect             1.1.1    2021-04-30 [2] CRAN (R 4.1.1)
##  usethis                2.1.2    2021-10-25 [2] CRAN (R 4.1.1)
##  utf8                   1.2.2    2021-07-24 [2] CRAN (R 4.1.1)
##  VariantAnnotation      1.40.0   2021-10-26 [2] Bioconductor  
##  vctrs                  0.3.8    2021-04-29 [2] CRAN (R 4.1.1)
##  withr                  2.4.2    2021-04-18 [2] CRAN (R 4.1.1)
##  xfun                   0.27     2021-10-18 [2] CRAN (R 4.1.1)
##  XML                    3.99-0.8 2021-09-17 [2] CRAN (R 4.1.1)
##  xml2                   1.3.2    2020-04-23 [2] CRAN (R 4.1.1)
##  xtable                 1.8-4    2019-04-21 [2] CRAN (R 4.1.1)
##  XVector                0.34.0   2021-10-26 [2] Bioconductor  
##  yaml                   2.2.1    2020-02-01 [2] CRAN (R 4.1.1)
##  zlibbioc               1.40.0   2021-10-26 [2] Bioconductor  
## 
## [1] /tmp/RtmpsqJJzl/Rinst21a2832b784e2e
## [2] /home/biocbuild/bbs-3.14-bioc/R/library

References

Anders, S., and W. Huber. 2010. “Differential Expression Analysis for Sequence Count Data.” Genome Biology 11 (10): R106.

Aumann, Y., and Y. Lindell. 2003. “A statistical theory for quantitative association rules.” Journal of Intelligent Information Systems 20: 255–83.

Kozomara, A., and S. Griffiths-Jones. 2014. “Base: annotating high confidence microRNAs using deep sequencing data.” NAR 42: D68–D73.

Salzberg, S., and B. Langmead. 2012. “Fast gapped-read alignment with Bowtie 2.” Nature Methods 9: 357–59.

Viollet, C., Davis D. A., M. Reczko, J. M. Ziegelbauer, F. Pezzella, J. Ragoussis, and R. Yarchoan. 2015. “Next-generation sequencing analysis reveals differential expression profiles of miRNA-mRNA target pairs in KSHV-infected cells.” PLOS ONE 10: 1–23.