Polyester is an R package designed to simulate an RNA sequencing experiment. Given a set of annotated transcripts, polyester will simulate the steps of an RNA-seq experiment (fragmentation, reverse-complementing, and sequencing) and produce files containing simulated RNA-seq reads. Simulated reads can be analyzed using any of several downstream analysis tools.
In particular, Polyester was designed to simulate a case/control experiment with biological replicates. Users are able to set differential transcript expression between cases and controls. This allows users to create datasets with known differential expression, which means they can the accuracy of statistical methods for differential expression detection.
Polyester was developed with several specific features in mind:
- Simulation of differential expression at the transcript level
- Ability to set differential expression signal strength
- Simulation of small datasets, since large RNA-seq datasets can require lots of time and computing resources to analyze
- Generation of raw RNA-seq reads (as opposed to read alignments or transcript-level abundance estimates)
- Transparency/open-source code
Start R and run:
source("http://bioconductor.org/biocLite.R")
biocLite("polyester")Polyester depends on the Biostrings and IRanges libraries from Bioconductor.
You will need either:
- a reference FASTA file containing names and sequences of transcripts from which reads should be simulated. Known transcripts from human chromosome 22 (hg19 build) are available in the
datasubdirectory of this package. - or a file in GTF format denoting transcript structures, along with one FASTA file of the DNA sequence for each chromosome in the GTF file. All the FASTA files should be in the same directory. DNA sequences for some organisms can be downloaded here (sequences are in the
<organism>/<source>/<build>/Sequence/Chromosomessubdirectory, e.g.,Homo_sapiens/UCSC/hg19/Sequence/Chromosomes).
Simulating an RNA-seq experiment with Polyester requires just one function call. You can choose either simulate_experiment() or simulate_experiment_countmat().
A FASTA file called chr22.fa is provided with polyester. This file contains sequences for 918 transcripts on chromosome 22, as annotated in hg19. For this example, we will only simulate from the first 20 of these transcripts, and we will set the first 2 transcripts to be overexpressed in group A and the next 2 transcripts to be overexpressed in group B, each at a fold change of 4. A small experiment like this will only take a few seconds to run, even with many reps. A larger experiment (say, with all 918 transcripts) will run also run in a reasonable amount of time (minutes, not hours), with the exact timing depending on number of reads generated and number of reps.
To simulate a two-group experiment with 10 biological replicates in each group where the first 3 transcripts are differentially expressed with a fold change of 4, you can use code like this:
library(polyester)
library(Biostrings)
fasta_file = system.file('extdata', 'chr22.fa', package='polyester')
fasta = readDNAStringSet(fasta_file)
small_fasta = fasta[1:20]
writeXStringSet(small_fasta, 'chr22_small.fa')
fold_changes = c(4, 4, 1/4, 1/4, rep(1, 16))
outdir = 'simulated_reads'
# ~20x coverage ----> reads per transcript = length/readlength * 20
# "width" is operating on a DNAStringSet (from Biostrings)
readspertx = round(20 * width(small_fasta) / 100)
simulate_experiment('chr22_small.fa', reads_per_transcript=readspertx,
num_reps=10, fold_changes=fold_changes, outdir=outdir) The simulate_experiment function draws the number of reads to simulate from each transcript from a negative binomial model. See below for details. Depending on your use case, it may be important to account for transcript length when deciding on the baseline mean number of reads to simulate from that transcript (as we did above with readspertx).
For more flexibility, you can use the simulate_experiment_countmat function. For example, we may want to simulate timecourse data. To do this, we can explicitly specify the number of reads for each transcript (rows), at each timepoint (columns). We will again only simulate from 20 transcripts.
# set up matrix:
num_timepoints = 12
countmat = matrix(readspertx, nrow=length(small_fasta), ncol=num_timepoints)
# add spikes in expression at certain timepoints to certain transcripts:
up_early = c(1,2)
up_late = c(3,4)
countmat[up_early, 2] = 3*countmat[up_early, 2]
countmat[up_early, 3] = round(1.5*countmat[up_early, 3])
countmat[up_late, 10] = 6*countmat[up_late, 10]
countmat[up_late, 11] = round(1.2*countmat[up_late, 11])
# simulate reads:
simulate_experiment_countmat('chr22_small.fa', readmat=countmat,
outdir='timecourse_reads') In this scenario, we simulated 12 total timepoints. We also added differential expression: transcripts 1 and 2 are overexpressed (compared to baseline) at timepoints 2 and 3, with a fold change of 3 at timepoint 2 and a fold change of 1.5 at timepoint 3. Similarly, transcripts 3 and 4 are overexpressed at timepoints 10 and 11, with a fold change of 6 at timepoint 10 and a fold change of 1.2 at timepoint 11.
The simulate_experiment function draws the number of reads to simulate from each transcript from a negative binomial distribution. For this function, you need to specify:
num_reps: Number of biological replicates per experimental group (default: 10; can specify different numbers of replicates in the groups)fold_changes: A fold change for each transcript. This fold change represents the multiplicative change in the mean number of reads generated from each transcript, between the two experimental groups.reads_per_transcript: The baseline mean number of reads for each transcript.- Fold changes compare the mean number of reads in group 1 to group 2. So a fold change of 0.5 means group 2's baseline mean number of reads for this transcript is twice that of group 1.
- Long transcripts usually produce more reads in RNA-seq experiments than short ones, so you may want to specify
reads_per_transcriptas a function of transcript length - Default is 300 (regardless of transcript length).
size: controls the per-transcript mean/variance relationship. In the negative binomial distribution, the mean/variance relationship is:mean = mean + (mean^2) / size. You can specify the size for each transcript. By default, size is defined as 1/3 of the transcript's mean, which (in our experience) creates a somewhat idealized, low-variance situation. Decrease the value ofsizeto introduce more variance into your simulations.
The simulate_experiment_readmat function takes a count matrix as an argunent. Each row of this matrix represents a transcript, and each column represents a sample in the experiment. Entry i,j of the matrix specifies how many reads should be sampled from transcript i for sample j, allowing you to precisely and flexibly define the (differential) transcript expression structure for the experiment.
For both simulate_experiment and simulate_experiment_countmat, you can change these parameters:
fraglen: Mean fragment length (default 250)fragsd: Standard devation of fragment lengths (default 25)readlen: Read length (default 100)error_rate: Sequencing error rate: probability that the sequencer records the wrong nucleotide at any given base (default 0.005, uniform error model assumed)paired: Whether the reads should be paired-end (default TRUE)
This review paper (Oshlack, Robinson, and Young, Genome Biology 2010, open access) provides a good overview of the RNA sequencing process, and might be particularly useful for understanding where some of these simulation parameters come into play.
If you'd like to explore specific steps in the sequencing process (fragmentation, reverse-complementing, error-adding), the functions called within simulate_experiment are also available and individually documented in Polyester.
See ?simulate_experiment and ?simulate_experiment_countmat for details on how to change these parameters.
To create a count matrix that resembles a real dataset, use the create_read_numbers function. To run this example, you will need to install the Ballgown package from Bioconductor if you do not already have it:
source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")library(ballgown)library(ballgown)
data(bg)
countmat = fpkm_to_counts(bg, threshold=0.01, mean_rps=400000)
params = get_params(countmat)
Ntranscripts = 50
Nsamples = 10
custom_readmat = create_read_numbers(mu=params$mu, fit=params$fit, p0=params$p0, m=Ntranscripts, n=Nsamples, seed=103)## Generating data from baseline model.
The Ballgown package here is optional: the mean/variance relationship for each transcript can be estimated from any matrix of counts using get_params. You can add differential expression to the output from create_read_numbers (here, custom_readmat) and pass the resulting matrix to simulate_experiment_countmat.
A call to simulate_experiment or simulate_experiment_countmat will write FASTA files to the directory specified by the outdir argument. Reads in the FASTA file will be labeled with the transcript from which they were simulated.
If paired is true, you'll get two FASTA files per biological replicate (left mates are designated by the suffix _1.fasta; right mates by _2.fasta). If single-end reads are generated (paired=FALSE) you'll get one FASTA file per replicate.
Files will be named sample_01 through sample_N where N is the total number of replicates. The first num_reps (or num_reps[1]) samples belong to the same group in the two-group experiment scenario.
In simulate_experiment, by default, a table called sim_info.txt is written to outdir, which will contain transcript IDs, fold changes, and whether or not that transcript was set to be differentially expressed. This file could be useful for downstream analysis. If the transcript names in the FASTA file cause problems down the line (e.g., a dangling single quote from a 5'-end label), you can specify your own transcript names with the transcriptid argument. You will need to keep track of this information separately if you use simulate_experiment_countmat.
Report bugs as issues on our GitHub repository.
sessionInfo()## R version 3.1.0 (2014-04-10)
## Platform: x86_64-apple-darwin10.8.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ballgown_0.99.5 Biostrings_2.33.14 XVector_0.5.8
## [4] IRanges_1.99.28 S4Vectors_0.2.4 BiocGenerics_0.11.5
## [7] polyester_0.99.1 knitr_1.6 devtools_1.5
## [10] BiocInstaller_1.15.5
##
## loaded via a namespace (and not attached):
## [1] annotate_1.43.5 AnnotationDbi_1.27.10
## [3] BatchJobs_1.3 BBmisc_1.7
## [5] Biobase_2.25.0 BiocParallel_0.99.19
## [7] bitops_1.0-6 brew_1.0-6
## [9] checkmate_1.4 codetools_0.2-9
## [11] DBI_0.3.0 digest_0.6.4
## [13] evaluate_0.5.5 fail_1.2
## [15] foreach_1.4.2 formatR_1.0
## [17] genefilter_1.47.6 GenomeInfoDb_1.1.19
## [19] GenomicAlignments_1.1.29 GenomicRanges_1.17.40
## [21] grid_3.1.0 httr_0.5
## [23] iterators_1.0.7 lattice_0.20-29
## [25] limma_3.21.15 Matrix_1.1-4
## [27] memoise_0.2.1 mgcv_1.8-3
## [29] nlme_3.1-117 RColorBrewer_1.0-5
## [31] RCurl_1.95-4.3 Rsamtools_1.17.33
## [33] RSQLite_0.11.4 rtracklayer_1.25.16
## [35] sendmailR_1.1-2 splines_3.1.0
## [37] stringr_0.6.2 survival_2.37-7
## [39] sva_3.11.4 tools_3.1.0
## [41] whisker_0.3-2 XML_3.98-1.1
## [43] xtable_1.7-3 zlibbioc_1.11.1