# The sRNADiff package

#### Matthias Zytnicki

#### 4 January 2019

sRNADiff aims at finding differentially expressed short RNAs with or without annotation. To do so, sRNADiff uses up to 4 different strategies:

- annotation-based,
- naive,
- HMM,
- clustering.

## Example

Raw data have been downloaded from the GEO data set GSE62830, provided in Viollet et al., 2015 (Viollet et al. 2015). Adapters were removed with fastx_clipper and mapped with bowtie2 (Salzberg and Langmead 2012) on the human genome version GRCh38.

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

The file `data.csv` contains three columns for each BAM file:

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

    dir         <- system.file("extdata", package="srnadiff", mustWork = TRUE)
    data        <- read.csv(file.path(dir, "data.csv"))
    gtfFile     <- file.path(dir, "Homo_sapiens.GRCh38.76.gtf.gz")
    annotation  <- readWholeGenomeAnnotation(gtfFile)
    bamFiles    <- file.path(dir, data$FileName)
    replicates  <- data$SampleName
    conditions  <- factor(data$Condition)
    exp         <- sRNADiffExp(annotation, bamFiles, replicates, conditions)
    diffRegions <- runAll(exp)

For your conveniance, the sample can be loaded with an only command `exp <- sRNADiffExample()`, so the script boils down to:

    exp         <- sRNADiffExample()
    diffRegions <- runAll(exp)

## Rationale and scope of the package

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.

## Explanation of the different strategies

Each strategy segments the genome its own way and finds potential differentially expressed regions. Every strategy makes assumption of the profiles of the differentially expressed regions. These assumptions differ between strategies, and thus it is useful to combined them.

### Annotation

This is the standard pipe-line, and the method simply provides the list of miRNA genes extracted from the given GTF file.

If you do not provide a GTF file, this step is skipped.

Unfortunately, there is no clear standard for an annotation format. We provide here four ways to parse the annotation file.

#### 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(dir, "Homo_sapiens.GRCh38.76.gtf.gz")
    annotation  <- readWholeGenomeAnnotation(gtfFile)

#### miRBase annotation

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 function parses a GFF3 miRBase file, and extracts the precursor miRNAs.

    gffFile    <- file.path(dir, "mirbase21_GRCh38.gff3")
    annotation <- readMiRBasePreAnnotation(gffFile)

As a result, 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:

    annotation <- readMiRBaseMatureAnnotation(gffFile)

#### Other format

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

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

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`).

### Naive

This strategy is the simplest, most commonly used, to find transcription units. It proceeds in two steps.

The first step merges all the BAM files. Regions with a minimum number coverage (summing all the reads in all input BAM files) are kept.

The second steps merges regions that are within a given distance (also provided by the user). The default distance can be change using this function:

    exp <- setMergeDistance(exp, 100)

This method is efficient when a gene is differentially expressed on the whole transcript.

### HMM

This strategy first computes the coverage of every input BAM file of every nucleotide of the genome. If two genome nucleotides have exactly the same coverage profile, we merge them. We give these coverage profiles to DESeq2 (Love, Huber, and Anders 2014), and compute the p-values. We thus have a p-value for each position of the genome. We then use a very crude hidden Markov model (HMM) on each chromosome of the genome. The HMM has two states: “not differentially expressed”, and “differentially expressed”. Each state uses a binomial distribution, the probability to emit a p-value < \(t\), and the probability to emit a p-value >= \(t\) (where \(t\) is a threshold that can be tuned, initially set to 0.1). We did not try to optimize the parameters of the HMM, because it takes times, and do not improve the results. Instead, the HMM is used a crude (but effective) way to segment the genome.

The parameters of the HMM can be changed using these functions:

    exp <- setTransitionProbabilities(exp, noDiffToDiff = 0.0001, diffToNoDiff = 0.01)
    exp <- setEmissionProbabilities(exp, probability = 0.9)
    exp <- setEmissionThreshold(exp, threshold = 0.1)

where:

- `noDiffToDiff` is the probability to move from the “not differentially expressed” state to the “differentially expressed” state,
- `diffToNoDiff` is the probability to move from the “differentially expressed” state to the “not differentially expressed” state,
- `probability` is the probability to emit a p-value < \(t\) in the “differentially expressed” state, and a p-value >= \(t\) in the “not differentially expressed” state, where \(t\) is the next threshold,
- `threshold` is the threshold that limits each state.

This method is efficient when a region is highly differentially expressed, even if it is concealed in a longer transcript.

### Slicing

This stretegy proceeds in three steps:

- It normalizes counts using a CPM strategy developped in EdgeR.
- The mean distribution of the coverage of the replicates is computed for each condition.
- For each condition, we “slice” the mean coverage distribution: we add horizontal lines at height h (where h is the minimum depth parameter defined in the general parameter section), \(h+1\), \(h+2\), …, \(+\infty\), and the intervals where the distribution and the horizontal lines meet is a putative region.
- These putative regions are the output of the strategy.

Obviously, many putative regions may have very similar start and end points. If two regions differ only by a few base pairs, the smallest regions will be discarded.

The default values can be changed using these functions:

    exp <- setMinDifferences(exp, differences = 50)

This method is especially efficient when a differentially expressed region is large, even though the fold change is moderate.

### 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:

    exp <- setMinDepth(exp, depth = 10)
    exp <- setSizes(exp, minSize = 15, maxSize = 1000)

### Combination of the strategies

#### 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 function

    exp <- setStrategies(exp, annotation = TRUE, naive = FALSE, hmm = TRUE, slice = TRUE)

where `annotation`, `naive`, `hmm` and `slice` are logical variables, set to `FALSE` if you want to skip that step.

#### 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:

    exp <- setMinOverlap(exp, minOverlap = 20)

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

### Extracting regions

The regions, with corresponding information provided by DESeq2 (mean expression, fold-change, p-value, adjusted p-value, etc.), can be extracted with this command:

    exp     <- sRNADiffExample()
    exp     <- runAll(exp)
    regions <- regions(exp, pvalue = 0.5)

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

## Visualization

You can get the read coverage profile of your selected regions using the `plotRegion` method. For instance, if `exp` is your `srnadiff` object, you can type:

    plotRegion(exp, regions[1])

Each line provides the read density on an input sample. The dark grey rectangle emphasizes the selected region.

## Misc

### Naming the conditions

The conditions can be sorted (`WT` before `mutant`, `untreated` before `treated`) this way:

    conditions  <- factor(c("wt", "wt", "ko", "ko"), levels=c("wt", "ko"))

### Using several cores

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

    exp <- setNThreads(exp, nThreads = 4)

### 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.

# Session information

    devtools::session_info()

    ## - Session info ----------------------------------------------------------
    ##  setting  value
    ##  version  R version 3.5.2 (2018-12-20)
    ##  os       Windows Server 2012 R2 x64
    ##  system   x86_64, mingw32
    ##  ui       RTerm
    ##  language (EN)
    ##  collate  C
    ##  ctype    English_United States.1252
    ##  tz       America/New_York
    ##  date     2019-01-04
    ##
    ## - Packages --------------------------------------------------------------
    ##  package              * version   date       lib source
    ##  acepack                1.4.1     2016-10-29 [2] CRAN (R 3.5.1)
    ##  annotate               1.60.0    2019-01-04 [2] Bioconductor
    ##  AnnotationDbi          1.44.0    2019-01-04 [2] Bioconductor
    ##  assertthat             0.2.0     2017-04-11 [2] CRAN (R 3.5.2)
    ##  backports              1.1.3     2018-12-14 [2] CRAN (R 3.5.1)
    ##  base64enc              0.1-3     2015-07-28 [2] CRAN (R 3.5.0)
    ##  bindr                  0.1.1     2018-03-13 [2] CRAN (R 3.5.1)
    ##  bindrcpp               0.2.2     2018-03-29 [2] CRAN (R 3.5.1)
    ##  Biobase                2.42.0    2019-01-04 [2] Bioconductor
    ##  BiocGenerics           0.28.0    2019-01-04 [2] Bioconductor
    ##  BiocManager            1.30.4    2018-11-13 [2] CRAN (R 3.5.1)
    ##  BiocParallel           1.16.5    2019-01-05 [2] Bioconductor
    ##  BiocStyle            * 2.10.0    2019-01-05 [2] Bioconductor
    ##  biomaRt                2.38.0    2019-01-05 [2] Bioconductor
    ##  Biostrings             2.50.2    2019-01-05 [2] Bioconductor
    ##  bit                    1.1-14    2018-05-29 [2] CRAN (R 3.5.2)
    ##  bit64                  0.9-7     2017-05-08 [2] CRAN (R 3.5.2)
    ##  bitops                 1.0-6     2013-08-17 [2] CRAN (R 3.5.2)
    ##  blob                   1.1.1     2018-03-25 [2] CRAN (R 3.5.2)
    ##  callr                  3.1.1     2018-12-21 [2] CRAN (R 3.5.2)
    ##  checkmate              1.8.5     2017-10-24 [2] CRAN (R 3.5.1)
    ##  cli                    1.0.1     2018-09-25 [2] CRAN (R 3.5.1)
    ##  cluster                2.0.7-1   2018-04-13 [2] CRAN (R 3.5.2)
    ##  colorspace             1.3-2     2016-12-14 [2] CRAN (R 3.5.1)
    ##  crayon                 1.3.4     2017-09-16 [2] CRAN (R 3.5.2)
    ##  data.table             1.11.8    2018-09-30 [2] CRAN (R 3.5.1)
    ##  DBI                    1.0.0     2018-05-02 [2] CRAN (R 3.5.2)
    ##  DelayedArray           0.8.0     2019-01-05 [2] Bioconductor
    ##  desc                   1.2.0     2018-05-01 [2] CRAN (R 3.5.1)
    ##  DESeq2                 1.22.2    2019-01-05 [2] Bioconductor
    ##  devtools               2.0.1     2018-10-26 [2] CRAN (R 3.5.1)
    ##  digest                 0.6.18    2018-10-10 [2] CRAN (R 3.5.2)
    ##  dplyr                  0.7.8     2018-11-10 [2] CRAN (R 3.5.1)
    ##  evaluate               0.12      2018-10-09 [2] CRAN (R 3.5.2)
    ##  foreign                0.8-71    2018-07-20 [2] CRAN (R 3.5.2)
    ##  Formula                1.2-3     2018-05-03 [2] CRAN (R 3.5.0)
    ##  fs                     1.2.6     2018-08-23 [2] CRAN (R 3.5.1)
    ##  genefilter             1.64.0    2019-01-04 [2] Bioconductor
    ##  geneplotter            1.60.0    2019-01-04 [2] Bioconductor
    ##  GenomeInfoDb           1.18.1    2019-01-05 [2] Bioconductor
    ##  GenomeInfoDbData       1.2.0     2018-12-21 [2] Bioconductor
    ##  GenomicAlignments      1.18.1    2019-01-05 [2] Bioconductor
    ##  GenomicFeatures        1.34.1    2019-01-05 [2] Bioconductor
    ##  GenomicRanges          1.34.0    2019-01-05 [2] Bioconductor
    ##  ggplot2                3.1.0     2018-10-25 [2] CRAN (R 3.5.1)
    ##  glue                   1.3.0     2018-07-17 [2] CRAN (R 3.5.2)
    ##  gridExtra              2.3       2017-09-09 [2] CRAN (R 3.5.1)
    ##  gtable                 0.2.0     2016-02-26 [2] CRAN (R 3.5.1)
    ##  Hmisc                  4.1-1     2018-01-03 [2] CRAN (R 3.5.1)
    ##  hms                    0.4.2     2018-03-10 [2] CRAN (R 3.5.2)
    ##  htmlTable              1.13      2019-01-02 [2] CRAN (R 3.5.2)
    ##  htmltools              0.3.6     2017-04-28 [2] CRAN (R 3.5.1)
    ##  htmlwidgets            1.3       2018-09-30 [2] CRAN (R 3.5.1)
    ##  httr                   1.4.0     2018-12-11 [2] CRAN (R 3.5.2)
    ##  IRanges                2.16.0    2019-01-04 [2] Bioconductor
    ##  knitr                * 1.21      2018-12-10 [2] CRAN (R 3.5.2)
    ##  labeling               0.3       2014-08-23 [2] CRAN (R 3.5.0)
    ##  lattice                0.20-38   2018-11-04 [2] CRAN (R 3.5.2)
    ##  latticeExtra           0.6-28    2016-02-09 [2] CRAN (R 3.5.1)
    ##  lazyeval               0.2.1     2017-10-29 [2] CRAN (R 3.5.1)
    ##  locfit                 1.5-9.1   2013-04-20 [2] CRAN (R 3.5.1)
    ##  magrittr               1.5       2014-11-22 [2] CRAN (R 3.5.2)
    ##  Matrix                 1.2-15    2018-11-01 [2] CRAN (R 3.5.2)
    ##  matrixStats            0.54.0    2018-07-23 [2] CRAN (R 3.5.2)
    ##  memoise                1.1.0     2017-04-21 [2] CRAN (R 3.5.2)
    ##  munsell                0.5.0     2018-06-12 [2] CRAN (R 3.5.1)
    ##  nnet                   7.3-12    2016-02-02 [2] CRAN (R 3.5.2)
    ##  pillar                 1.3.1     2018-12-15 [2] CRAN (R 3.5.1)
    ##  pkgbuild               1.0.2     2018-10-16 [2] CRAN (R 3.5.1)
    ##  pkgconfig              2.0.2     2018-08-16 [2] CRAN (R 3.5.2)
    ##  pkgload                1.0.2     2018-10-29 [2] CRAN (R 3.5.1)
    ##  plyr                   1.8.4     2016-06-08 [2] CRAN (R 3.5.1)
    ##  prettyunits            1.0.2     2015-07-13 [2] CRAN (R 3.5.2)
    ##  processx               3.2.1     2018-12-05 [2] CRAN (R 3.5.1)
    ##  progress               1.2.0     2018-06-14 [2] CRAN (R 3.5.2)
    ##  ps                     1.3.0     2018-12-21 [2] CRAN (R 3.5.2)
    ##  purrr                  0.2.5     2018-05-29 [2] CRAN (R 3.5.1)
    ##  R6                     2.3.0     2018-10-04 [2] CRAN (R 3.5.2)
    ##  RColorBrewer           1.1-2     2014-12-07 [2] CRAN (R 3.5.0)
    ##  Rcpp                   1.0.0     2018-11-07 [2] CRAN (R 3.5.2)
    ##  RCurl                  1.95-4.11 2018-07-15 [2] CRAN (R 3.5.1)
    ##  remotes                2.0.2     2018-10-30 [2] CRAN (R 3.5.1)
    ##  reshape2               1.4.3     2017-12-11 [2] CRAN (R 3.5.1)
    ##  rlang                  0.3.0.1   2018-10-25 [2] CRAN (R 3.5.2)
    ##  rmarkdown            * 1.11      2018-12-08 [2] CRAN (R 3.5.1)
    ##  rpart                  4.1-13    2018-02-23 [2] CRAN (R 3.5.2)
    ##  rprojroot              1.3-2     2018-01-03 [2] CRAN (R 3.5.1)
    ##  Rsamtools              1.34.0    2019-01-05 [2] Bioconductor
    ##  RSQLite                2.1.1     2018-05-06 [2] CRAN (R 3.5.2)
    ##  rstudioapi             0.8       2018-10-02 [2] CRAN (R 3.5.1)
    ##  rtracklayer            1.42.1    2019-01-05 [2] Bioconductor
    ##  S4Vectors              0.20.1    2019-01-04 [2] Bioconductor
    ##  scales                 1.0.0     2018-08-09 [2] CRAN (R 3.5.1)
    ##  sessioninfo            1.1.1     2018-11-05 [2] CRAN (R 3.5.1)
    ##  snow                   0.4-3     2018-09-14 [2] CRAN (R 3.5.2)
    ##  srnadiff             * 1.2.1     2019-01-04 [1] Bioconductor
    ##  stringi                1.2.4     2018-07-20 [2] CRAN (R 3.5.2)
    ##  stringr                1.3.1     2018-05-10 [2] CRAN (R 3.5.2)
    ##  SummarizedExperiment   1.12.0    2019-01-05 [2] Bioconductor
    ##  survival               2.43-3    2018-11-26 [2] CRAN (R 3.5.2)
    ##  testthat               2.0.1     2018-10-13 [2] CRAN (R 3.5.1)
    ##  tibble                 2.0.0     2019-01-04 [2] CRAN (R 3.5.2)
    ##  tidyselect             0.2.5     2018-10-11 [2] CRAN (R 3.5.1)
    ##  usethis                1.4.0     2018-08-14 [2] CRAN (R 3.5.1)
    ##  withr                  2.1.2     2018-03-15 [2] CRAN (R 3.5.1)
    ##  xfun                   0.4       2018-10-23 [2] CRAN (R 3.5.2)
    ##  XML                    3.98-1.16 2018-08-19 [2] CRAN (R 3.5.1)
    ##  xtable                 1.8-3     2018-08-29 [2] CRAN (R 3.5.1)
    ##  XVector                0.22.0    2019-01-05 [2] Bioconductor
    ##  yaml                   2.2.0     2018-07-25 [2] CRAN (R 3.5.2)
    ##  zlibbioc               1.28.0    2019-01-05 [2] Bioconductor
    ##
    ## [1] C:/Users/biocbuild/bbs-3.8-bioc/tmpdir/RtmpglU5o9/Rinst156067a0135f
    ## [2] C:/Users/biocbuild/bbs-3.8-bioc/R/library

# References

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

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15 (12):550.

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

Viollet, Coralie, David A. Davis, Martin Reczko, Joseph M. Ziegelbauer, Francesco Pezzella, Jiannis Ragoussis, and Robert Yarchoan. 2015. “Next-Generation Sequencing Analysis Reveals Differential Expression Profiles of Mirna-mRNA Target Pairs in Kshv-Infected Cells.” PLOS ONE 10:1–23.