In biomedbigdata/ASimulatoR: Alternative Splicing Simulator

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

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

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:

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

If you still encounter problems with package versions you can use the renv package and simply restore our lockfile.

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

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. 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" :

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.

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

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

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

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

# 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

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

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:

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

biomedbigdata/ASimulatoR documentation built on Sept. 6, 2022, 7:55 p.m.