README.Rmd

---
output: github_document
---

```{r, include = FALSE}
options(width = 400)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  results = "hold"
)
```

# ASimulatoR

The goal of ASimulatoR is to simulate RNA-seq reads with alternative
splicing events. The alternative splicing events are well documented and
the true origin of each read is used for exon and juntion coverage via a
modified version of the bioconductor polyester package.

Please note that we use a custom version of Polyester, that is available
at <https://github.com/biomedbigdata/polyester>.

## Installation

You can install the current version using the
[remotes](https://cran.r-project.org/web/packages/remotes/index.html)
package:

``` r
# install.packages("remotes")
remotes::install_github("biomedbigdata/ASimulatoR")

# for a specific release e.g. v1.0.0 use
remotes::install_github("biomedbigdata/ASimulatoR@v1.0.0")
```

If you want the latest development version, just install from the dev
branch:

``` r
remotes::install_github("biomedbigdata/ASimulatoR", ref = "dev")
```

This will automatically install all imported R packages. 

There may be system requirements that you have to install manually
e.g. on Ubuntu/Debian:

``` bash
apt-get install -y libxml2-dev libcurl4-openssl-dev
```

If you still encounter problems with package versions you can use the
[renv](https://rstudio.github.io/renv/articles/packages.html) package
and simply `restore` our [lockfile](renv.lock).


## Usage

The input directory should contain a gtf/gff3 genome annotation file and one fasta file per chromosome (e.g. derived from the [Ensembl ftp
server](https://www.ensembl.org/info/data/ftp/index.html)). 

ASimulatoR accepts any gtf/gff3 genome annotation and chromosome fasta files independent of the organism.

### Quick start

This repository contains a documented example Rscript
[runASimulatoR.R](runASimulatoR.R). After installation, scripts like
this can be run from the command-line with the command:

`Rscript runASimulatoR.R /path/to/input_folder/ /path/to/output_folder/`.

### Step by step

For usage in an interactive R-session and to investigate the full
functionality of this package check the following example. 
The example is limited to the chromosome 21 for the runtime purpose. 
For even more flexibility refer to the section ["The main function: `simulate_alternative_splicing`" ](## The main function: `simulate_alternative_splicing`):

#### Creating exon supersets

Firstly, we create exon supersets by joining all exons of a gene from a
gtf/gff3 file. These supersets are then used to create splice variants.

Note: Since all exons from one gene are used to create the exon superset, you
may find that the term exon superset is used analogously to gene.

```{r createSuperset}
suppressMessages(library(ASimulatoR))

# create exon superset for genes on chromosome 21 of ensembl release 99
gtf_file = system.file('extdata', 'Homo_sapiens.GRCh38.99.21.gtf', package = 'ASimulatoR')

# by default the produced superset will be saved as .rda file into the same directory
exon_superset = get_exon_supersets(gtf_file)
exon_superset[[1]][1:5, ]
```

#### Simulating Alternative Splicing

You can find more information about the main function of this package at
the [end of the page]((## The main function: `simulate_alternative_splicing`)).

The simulator supports eight different AS events:

| es           | mes                    | ir               | a3                                  | a5                               | mee                      | afe                    | ale                   |
| ------------ | ---------------------- | ---------------- | ----------------------------------- | -------------------------------- | ------------------------ | ---------------------- | --------------------- |
| exon skiping | multiple exon skipping | intron retention | alternative 3’/acceptor splice site | alternative 5’/donor splice site | mutually exclusive exons | alternative first exon | alternative last exon |


```{r parameters}
# define your input_dir, where the annotation gtf (or the exon supersets if you have already created them) and the genome fasta files (one per chromosome) are located
# here we will use the example data
input_dir = system.file('extdata', package = 'ASimulatoR')

# define, how many groups and samples per group you analyze. Here we create a small experiment with two groups with one sample per group:
num_reps = c(1,1)

# define your outdir with NO slash
outdir = 'simulation'

# define the number of genes you want to work with. If you want all exons, do not specify this parameter or set it to NULL
# here we create splice variants from 9 exon supersets:
max_genes = 9
```

The user could define the distribution of the events by probability or
relative frequency.

  - Probability: For each superset we create an event with the
    probability mentioned in `event_prob`.
  - Frequency: Set `probs_as_freq = T`. The exon supersets are
    partitioned corresponding to the `event_prob` parameter.


```{r simulateAS}
# in this example we use relative frequencies
# here we produce eight variants with one of each AS events as well as one variant containing every event
# if probs_as_freq was FALSE, a random number would be drawn for each event-superset combination and only if it was smaller than 1/9 the AS event would be created
probs_as_freq = T
event_freq = 
    setNames(rep(1/9, 9),
             c('es', 'mes', 'ir', 'a3', 'a5', 'afe', 'ale', 'mee', 'es,ir,mes,a3,a5,afe,ale,mee'))



# we use the previously created superset to simulate splice variants from, since it is saved in the same directory as the gtf
# if no superset is found, a new one will be created
simulate_alternative_splicing(input_dir = input_dir,
                              outdir = outdir, 
                              event_probs = event_freq,
                              probs_as_freq = probs_as_freq, 
                              max_genes = max_genes,
                              num_reps = num_reps,
                              verbose = TRUE)
```

The user can also use predefined parameters for common scenarios. More information is in `?presets`:

```{r simulatePreset}
outdir_preset = 'simulation_preset'

# We will use the preset 'experiment_bias' which includes typical rna-seq biases such as positional bias arising in protocols that use cDNA fragmentation
# we can still use other arguments such as max_genes to set additional/override existing parameters
simulate_alternative_splicing(input_dir = input_dir,
                              outdir = outdir_preset,
                              preset = 'experiment_bias',
                              max_genes = max_genes,
                              verbose = FALSE)
```

#### Visualize Splice Variants

```{r visualization, dpi=300, results = "hold"}
# to visualize the splice variants we will use ggbio
suppressMessages(library(ggbio))

# firstly, we load the newly created gtf file 
gtf = rtracklayer::import('simulation/splicing_variants.gtf')

# the gene id of the variant with all events
gene_id = gtf$gene_id[grep('es,ir,mes,a3,a5,afe,ale,mee', gtf$transcript_id, fixed = T)[1]]
exons = gtf[gtf$type == 'exon' & gtf$gene_id == gene_id]
suppressWarnings(ggbio::autoplot(split(exons, exons$transcript_id)))

# have a look at the event annotation
event_anno = read.csv('simulation/event_annotation.tsv', sep = '\t')
event_anno[grepl(gene_id, event_anno$template) | grepl(gene_id, event_anno$variant), ]
```

## The main function: `simulate_alternative_splicing`

### Description

Here we describe in more detail the parameters that can be passed to this function for a better control over the simulation.

### Usage

``` r
simulate_alternative_splicing(input_dir, outdir, event_probs, ncores = 1L, ...)
```

### Arguments

| Argument      | Description                                                                                                                                                                                                                                                                                                                                                  |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
| `input_dir`   | Character path to directory containing the gtf/gff file from which splice variants are created and genome fasta files with one file per chromosome i.e. <chr_name>.fa passed to polyester                                                                                                                                                                    |
| `outdir`      | character, path to folder where simulated reads and all annotations should be written, with *no* slash at the end. By default, reads are written to current working directory.                                                                                                                                                                               |
| `event_probs` | Named list/vector containing numerics corresponding to the probabilites to create the event (combination). If `probs_as_freq` is `TRUE` `event_probs` correspond to the relative frequency of occurences for the event(combination) and in this case the sum of all frequencies has to be \<=1. No default, must not be `NULL`, except if `preset` is given. |
| `preset`      | if you want to use preset parameters one of ‘event\_partition’, ‘experiment\_bias’, ‘event\_combination\_2’. Check `?presets` for more information                                                                                                                                                                                                           |
| `ncores`      | the number of cores to be utilized for parallel generation of splice variant creation and read simulation. **This will spawn one process per core! Be aware that a lot of memory might be required for many processes.**                                                                                                                                                                                                                                                  |
| `...`         | any of several other arguments that can be used to add nuance to the simulation and splice variant creation. See section Details.                                                                                                                                                                                                                                    |

### Details

Several optional parameters can be defined to adjust the
simulation. These parameters are further used by the `simulate_experiment` function from
our [custom polyester R
package](https://github.com/biomedbigdata/polyester). 
 
 The following parameters are specific for the ASimulatoR package:
 
  - `novel_variants` : Numeric value between 0 and 1 indicating the
    percentage of splicing variants that will be suppressed in an
    additional gtf file splicing\_variants\_novel.gtf.
    
  - `write_gff` : Additionally to the gtf file containing the splice
    variants, a gff3 file with the same content will be printed to the
    outdir. Default `TRUE`

  - `max_genes` : The maximum number of genes/exon supersets to be
    included in the process of splice variant creation. Default `NULL`
    which means that all available exon supersets will be used. **This
    is a computation heavy default and you might want to adjust it\!**

  - `exon_junction_coverage` : Should the real coverage of exons,
    junctions and retained introns be written into a additional file.
    Default `TRUE`

  - `multi_events_per_exon` : Should it be possible to have more than
    one AS event at the same exon if multiple variants are created for
    the same exon superset? \!If this option is set to `TRUE` , there
    may occur unforeseen AS events that are not documented in the
    event\_annotation file\!. Default `FALSE`

  - `probs_as_freq` : Should `event_probs` be treated as relative
    frequencies instead of probabilities? Default `FALSE`

  - `save_exon_superset` : Should the exon supersets be saved to .rda
    file? Default `TRUE`

These parameters are used by the polyester function
`simulate_experiment` and have different defaults assigned in
ASimulatoR:
 
  - `fold_changes` : Currently, ASimulatoR introduces random isoform
    switches. Those can be retraced in the sim\_tx\_info.txt file
    written by polyester. We plan on improving this in the future.

  - `strand_specific` : Strand-specific simulation (1st read forward
    strand, 2nd read reverse strand with respect to transcript
    sequence). Default `TRUE` .

  - `meanmodel` : `reads_per_transcripts` as a function of transcript
    length. Always `TRUE` in ASimulatoR.

  - `frag_GC_bias`: A sample-specific GC content bias on the fragment
    level. Currently not supported in ASimulatoR: always ‘none’.

  - `verbose` : Should progress messages be printed during the
    sequencing process? Default `TRUE` .

  - `exon_junction_coverage` : Should the coverage of exons, junctions
    and retained introns be determined? Default `TRUE` .

  - `exon_junction_table` : If `exon_junction_coverage=TRUE` a
    `data.table` produced by `create_splicing_variants_and_annotation`
    to determine exon and intron coverage.

These parameters are used by the `simulate_experiment` from the polyester
R package to introduce technical biases.

  - `pcr_rate`: Fraction of fragments that will be duplicated. Reads
    from these fragments will have PCR\_DUP in the name.
  - `pcr_lambda`: If `!is.null(pcr_rate)` lambda for the poisson
    distribution to draw the number of duplicates.
  - `adapter_contamination`: If the fragment is smaller than the
    readlength, should we sequence into the `adapter_sequence`?
  - `adapter_sequence`: If `adapter_contamination`: adapter sequence
  - `bias`: Positional bias model to use when fragmenting transcripts.
    By default, all fragments from a transcript are equally likely
    (`'none'`). Other choices are `'rnaf'` and `'cdnaf'`, which mimic
    positional bias arising from different fragmentation protocols. See
    `?generate_fragments` and the polyester manuscript (Frazee et al, 2014) for details.

Note:  The user can further adjust the RNA-Seq reads simulation. Try `?simulate_experiment` to check all available parameters.

## Manuscript Use Case:

If you would like to reproduce the use case mentioned in the manuscript
you can use the corresponding preset:

``` r
# this preset uses a sequencing depth of 200 million reads
simulate_alternative_splicing('some_input_dir', 'some_output_dir', preset = 'manuscript')

# if you want to adjust the sequencing depth you can easily override it to create a lower coverage
simulate_alternative_splicing('some_input_dir', 'some_output_dir', preset = 'manuscript', seq_depth = 5e07)
```

## Value

No return, but simulated reads, a simulation info file, an alternative
splicing event annotation and exon and junction coverages are written to
`outdir` .

## References

 Alyssa C. Frazee, Andrew E. Jaffe, Ben Langmead, Jeffrey T. Leek,
Polyester: simulating RNA-seq datasets with differential transcript
expression, Bioinformatics, Volume 31, Issue 17, 1 September 2015, Pages
2778–2784, <https://doi.org/10.1093/bioinformatics/btv272>

## License

ASimulatoR Copyright (C) 2020 Manz, Quirin

This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
Public License for more details.

You should have received a copy of the GNU General Public License along
with this program. If not, see <https://www.gnu.org/licenses/>.