This vignette provides the implementation of the procedure described in point 7 of our Guidelines for RNA-Seq data analysis1 protocol available from the Epigenesys website.
Briefly, it details the step necessary to: 1. create a non-redundant annotation
count reads per feature
pre-analyse the data, i.e. assess the pertinence of the samples’ charateristics in the light of their biological provenance; i.e. in other words perform a so called “biological QA” using assessment methods such as Principal Component Analysis, Multi-dimensional Scaling, Hierarcical Clustering, etc.
The aim of this vignette is to go through these steps using the easyRNASeq package, hence the rationale of the implementation will not be discussed, albeit relevant litterature will be pointed at when necessarry.
Throughout this vignette we are going to replicate the analysis conducted in Robinson, Delhomme et al.(Robinson et al. 2014), a study looking at sexual dimorphism in Eurasian aspen.
To perform the listed steps, we need to instantiate a number of objects to store the minimal set of parameters describing the conducted RNA-Seq experiment, e.g. the BAM files location, the annotation location and type, the sequencing parameters, etc.
To get started with this process, we load the package into our R session:
library(easyRNASeq)before instantiating an AnnotParam object informing on the location and type of the annotation to be used.
The AnnotParam class is meant to store the minimal set of information necessary to retrieve the annotation
The minimal information to provide is:
In this tutorial, we will reproduce the analysis performed in Robinson, Delhomme et al. (Robinson et al. 2014). For that we will start by downloading the original annotation gff3 file for P. trichocarpa, a close related species of the trees used in the study into the current directory.
library(curl)
curl_download(url=paste0("ftp://ftp.plantgenie.org/Data/PopGenIE/",
"Populus_trichocarpa/v3.0/v10.1/GFF3/",
"Ptrichocarpa_210_v3.0_gene_exons.gff3.gz"),
destfile=,"./Ptrichocarpa_210_v3.0_gene_exons.gff3.gz")Before instantiating an “AnnotParam” object.
annotParam <- AnnotParam(
datasource="./Ptrichocarpa_210_v3.0_gene_exons.gff3.gz")This annotation file however contains multiple copy of the same exons, i.e. when exons are shared by several isoforms of a gene. This might result in so-called “multiple-counting” and as described in these guidelines2, we will to circumvent that issue create a set of synthetic transcripts.
One major caveat of estimating gene expression using aligned RNA-Seq reads is that a single read, which originated from a single mRNA molecule, might sometimes align to several features (e.g. transcripts or genes) with alignments of equivalent quality.
This, for example, might happen as a result of gene duplication and the presence of repetitive or common domains. To avoid counting unique mRNA fragments multiple times, the stringent approach is to keep only uniquely mapping reads - being aware of potential consequences, see the note below.
Not only can multiple counting arise from a biological reason, but also from technical artifacts, introduced mostly by poorly formatted gff3/gtf annotation files. To avoid this, it is best practice to adopt a conservative approach by collapsing all existing transcripts of a single gene locus into a synthetic transcript containing every exon of that gene. In the case of overlapping exons, the longest genomic interval is kept, i.e. an artificial exon is created. This process results in a flattened transcript: a gene structure with a one to one relationship.
To create such a structure, we use the createSyntheticTranscripts function on the file we just downloaded, simply by passing our annotParam object as argument.
annotParam <- createSyntheticTranscripts(annotParam,verbose=FALSE)This function returns an updated annotParam object that contains the newly created, flattened transcript annotation. This object can then be saved as an rda file for later re-use or for sharing with collaborators.
save(annotParam,
file="./Ptrichocarpa_210_v3.0_gene_exons_synthetic-transcripts_annotParam.rda")Instead of updating the annotParam object, we could have created an object of class Genome_Intervals from the genomeIntervals package, using the same function but using the actual datasource of the previous annotParam object as argument rather than the object itself.
gI <- createSyntheticTranscripts(
"./Ptrichocarpa_210_v3.0_gene_exons.gff3.gz",
verbose=FALSE)This gI object can then be exported as a gff3 file.
writeGff3(gI,
file="./Ptrichocarpa_210_v3.0_gene_exons_synthetic-transcripts.gff3.gz")Note: Ignoring multi-mapping reads may introduce biases in the read counts of some genes (such as that of paralogs or of very conserved gene families), but in the context of a conservative first analysis we are of the current opinion that they are best ignored. One should of course assess how many reads are multi-mapping (check for example the STAR output) and possibly extract them from the alignment read file to visualize them using a genome browser so as to understand where they are located and how they may affect any analysis. Based on this, one may, at a later stage, decide to relax the counting parameters to accept multi-mapping reads.
High throughput sequencing RNA-Seq data comes in a multitude of flavours, i.e. even from a single provider, protocol - e.g. strand specific, paired-end - reads characteristics - e.g. read length - will vary.
The easyRNASeq simpleRNASeq method will infer these information based on excerpts sampled from the data. However, it is always best to provide these information, as 1. the inference is done on small excerpt and can fail 2. it is always good to document an analysis
By default easyRNASeq simpleRNASeq will keep the inferred parameters over the user-provided parameters if these do not agree and emit corresponding warnings. The choice to rely on inferred parameters over user-provided one is to enforce user to cross-validate their knowledge of the data characteristics, as these are crucial for an adequate processing. Remember GIGO3.
If the automatic inference does fail, please let me know, so that I optimise it. Meanwhile, you can use the override argument to enforce the use of user-passed parameters.
To reproduce the results from Robinson, Delhomme et al. (Robinson et al. 2014), we first need to download an excerpt of the data.
We first retrieve the file listing and md5 codes
library(curl)
curl_download(url=paste0("ftp://ftp.plantgenie.org/Tutorials/RnaSeqTutorial/",
"data/star/md5.txt"),
destfile="md5.txt")In this part of the vignette, we will NOT process all the data, albeit it would be possible, but for the sake of brevity, we will only retrieve the six first datasets. We get these from the sample information contained within this package.
data(RobinsonDelhomme2014)
lapply(RobinsonDelhomme2014[1:6,"Filename"],function(f){
curl_download(url=paste0("ftp://ftp.plantgenie.org/Tutorials/",
"RnaSeqTutorial/data/star/",f),
destfile=f)
})## list()These six files - as the rest of the dataset - have been sequenced on an Illumina HiSeq 2500 in paired-end mode using a non-strand specific library protocol with a read length of 100 bp. The raw data have been processed as described in the aforementioned guidelines4 and as such have been filtered for rRNA sequences, trimmed for adapters and clipped for quality. The resulting reads (of length 50-100bp) have then been aligned using STAR.
Using these information, we finally generate the BamParam object.
bamParam <- BamParam(paired = TRUE,
stranded = FALSE)A third parameter yieldSize can be set to speed up the processing on multi-CPU or multi-core computers. It splits and processed the BAM files in chunk of size yieldSize with a default of 1M reads.
The final set of parameters we need to define encapsulate the AnnotParam and BamParam and detail how the read summarization should be performed. simpleRNASeq supports A) 2 modes of counting:
by read
by bp
the latter of which, was the default counting method the easyRNASeq function. Due to the more complex implementation required, the non-evidence of increase in counting accuracy and the extended support of the read-based approach by the mainstream, standardised Bioconductor package has led the read method to be the default in simpleRNASeq. Due to lack of time for maintenance and improvement, the bp-based method is also not recommended.
over B) 4 feature types: exon, transcript, gene or any feature provided by the user. The latter may be for example used for counting reads in promoter regions.
Given a flattened transcript structure - as created in a previous section - summarizing by transcripts or genes is equivalent. Note that using a non flattened annotation with any feature type will result in multiple counting!! i.e. the product of a single mRNA fragment will be counted for every features it overlap, hence introducing a significant bias in downstream analyses.
Given a flattened transcript structure, summarizing by exon enables the use of the resulting count table for processes such as differential exon usage analyses, as implemented in the DEXSeq package.
For the Robinson, Delhomme et al. dataset, we are interested in the gene expression, hence we create our RnaSeqParam object as follows:
rnaSeqParam <- RnaSeqParam(annotParam = annotParam,
bamParam = bamParam,
countBy = "genes",
precision = "read")## R version 3.4.1 (2017-06-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.5-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.5-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] curl_2.7 easyRNASeq_2.12.1 BiocStyle_2.4.0
##
## loaded via a namespace (and not attached):
## [1] SummarizedExperiment_1.6.3 genefilter_1.58.1
## [3] locfit_1.5-9.1 splines_3.4.1
## [5] lattice_0.20-35 htmltools_0.3.6
## [7] stats4_3.4.1 yaml_2.1.14
## [9] blob_1.1.0 XML_3.98-1.9
## [11] survival_2.41-3 rlang_0.1.1
## [13] DBI_0.7 BiocParallel_1.10.1
## [15] BiocGenerics_0.22.0 bit64_0.9-7
## [17] RColorBrewer_1.1-2 matrixStats_0.52.2
## [19] GenomeInfoDbData_0.99.0 stringr_1.2.0
## [21] zlibbioc_1.22.0 Biostrings_2.44.1
## [23] hwriter_1.3.2 evaluate_0.10.1
## [25] memoise_1.1.0 latticeExtra_0.6-28
## [27] Biobase_2.36.2 knitr_1.16
## [29] biomaRt_2.32.1 geneplotter_1.54.0
## [31] IRanges_2.10.2 GenomeInfoDb_1.12.2
## [33] parallel_3.4.1 AnnotationDbi_1.38.1
## [35] Rcpp_0.12.11 edgeR_3.18.1
## [37] xtable_1.8-2 backports_1.1.0
## [39] limma_3.32.2 DelayedArray_0.2.7
## [41] S4Vectors_0.14.3 annotate_1.54.0
## [43] XVector_0.16.0 ShortRead_1.34.0
## [45] bit_1.1-12 Rsamtools_1.28.0
## [47] digest_0.6.12 stringi_1.1.5
## [49] DESeq_1.28.0 GenomicRanges_1.28.4
## [51] grid_3.4.1 rprojroot_1.2
## [53] tools_3.4.1 bitops_1.0-6
## [55] magrittr_1.5 RCurl_1.95-4.8
## [57] LSD_3.0 RSQLite_2.0
## [59] tibble_1.3.3 Matrix_1.2-10
## [61] rmarkdown_1.6 intervals_0.15.1
## [63] genomeIntervals_1.32.0 GenomicAlignments_1.12.1
## [65] compiler_3.4.1## [1] TRUE TRUE TRUERobinson, Kathryn, Nicolas Delhomme, Niklas Mähler, Bastian Schiffthaler, Jenny Önskog, Benedicte Albrectsen, Pär Ingvarsson, Torgeir Hvidsten, Stefan Jansson, and Nathaniel Street. 2014. “Populus Tremula (European Aspen) Shows No Evidence of Sexual Dimorphism.” BMC Plant Biology 14 (1): 276–76. doi:10.1186/s12870-014-0276-5.
http://www.epigenesys.eu/en/protocols/bio-informatics/1283-guidelines-for-rna-seq-data-analysis↩
http://www.epigenesys.eu/en/protocols/bio-informatics/1283-guidelines-for-rna-seq-data-analysis↩
Garbage In Garbage Out↩
http://www.epigenesys.eu/en/protocols/bio-informatics/1283-guidelines-for-rna-seq-data-analysis↩