Introduction

This vignette provides the implementation of the procedure described in point 7 of our Guidelines for RNA-Seq data analysis¹ protocol available from the Epigenesys website.

Briefly, it details the step necessary to: 1. create a non-redundant annotation

2. count reads per feature

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

Annotation parameters - AnnotParam

The AnnotParam class is meant to store the minimal set of information necessary to retrieve the annotation

The minimal information to provide is:

1. a datasource: a path to a file provided as a character string or if the type is biomaRt, the datasource you want to connect to, as retrieved using the biomaRt datasource function.
2. a type: one of “gff3”,“gtf”,“rda”, and “biomaRt”. gff3 is the default. If rda is used, it expects the corresponding file to contain a GRanges object by the name of gAnnot.

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 guidelines², we will to circumvent that issue create a set of synthetic transcripts.

Synthetic transcripts - createSyntheticTranscripts

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

gI <- createSyntheticTranscripts(
    "./Ptrichocarpa_210_v3.0_gene_exons.gff3.gz",
    verbose=FALSE)

writeGff3(gI,
          file="./Ptrichocarpa_210_v3.0_gene_exons_synthetic-transcripts.gff3.gz")

Bam Parameters - BamParam

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

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

RNA-Seq parameters - RnaSeqParam

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:

1. by read

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

Session Info

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.4 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-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_3.2          easyRNASeq_2.16.0 BiocStyle_2.8.0  
## 
## loaded via a namespace (and not attached):
##  [1] SummarizedExperiment_1.10.0 genefilter_1.62.0          
##  [3] progress_1.1.2              locfit_1.5-9.1             
##  [5] splines_3.5.0               lattice_0.20-35            
##  [7] htmltools_0.3.6             stats4_3.5.0               
##  [9] yaml_2.1.18                 blob_1.1.1                 
## [11] XML_3.98-1.11               survival_2.42-3            
## [13] DBI_0.8                     BiocParallel_1.14.0        
## [15] BiocGenerics_0.26.0         bit64_0.9-7                
## [17] RColorBrewer_1.1-2          matrixStats_0.53.1         
## [19] GenomeInfoDbData_1.1.0      stringr_1.3.0              
## [21] zlibbioc_1.26.0             Biostrings_2.48.0          
## [23] hwriter_1.3.2               evaluate_0.10.1            
## [25] memoise_1.1.0               latticeExtra_0.6-28        
## [27] Biobase_2.40.0              knitr_1.20                 
## [29] geneplotter_1.58.0          IRanges_2.14.0             
## [31] biomaRt_2.36.0              GenomeInfoDb_1.16.0        
## [33] parallel_3.5.0              AnnotationDbi_1.42.0       
## [35] Rcpp_0.12.16                edgeR_3.22.0               
## [37] xtable_1.8-2                backports_1.1.2            
## [39] limma_3.36.0                DelayedArray_0.6.0         
## [41] S4Vectors_0.18.0            annotate_1.58.0            
## [43] XVector_0.20.0              ShortRead_1.38.0           
## [45] bit_1.1-12                  Rsamtools_1.32.0           
## [47] digest_0.6.15               stringi_1.1.7              
## [49] DESeq_1.32.0                GenomicRanges_1.32.0       
## [51] grid_3.5.0                  rprojroot_1.3-2            
## [53] tools_3.5.0                 bitops_1.0-6               
## [55] magrittr_1.5                RCurl_1.95-4.10            
## [57] LSD_4.0-0                   RSQLite_2.1.0              
## [59] Matrix_1.2-14               prettyunits_1.0.2          
## [61] httr_1.3.1                  assertthat_0.2.0           
## [63] rmarkdown_1.9               R6_2.2.2                   
## [65] intervals_0.15.1            genomeIntervals_1.36.0     
## [67] GenomicAlignments_1.16.0    compiler_3.5.0

## [1] TRUE TRUE TRUE

Appendix