Abstract

RNA-seq analysis involves multiple steps from processing raw sequencing data to identifying, organizing, annotating, and reporting differentially expressed genes. bcbio is an open source, community-maintained framework providing automated and scalable RNA-seq methods for identifying gene abundance counts. We have developed bcbioRNASeq, a Bioconductor package that provides ready-to-render templates and wrapper functions to post-process bcbio output data. bcbioRNASeq automates the generation of high-level RNA-seq reports, including identification of differentially expressed genes, functional enrichment analysis and quality control analysis.

Introduction

For a high-level overview of our bcbio RNA-seq analysis pipeline, including detailed explanation of the bcbioRNASeq S4 class definition, first consult our workflow paper published in F1000 Research (Steinbaugh et al. 2018). This vignette is focused on more advanced usage and edge cases that a user may encounter when attempting to load a bcbio dataset and perform downstream quality control analysis.

Note: if you use bcbioRNASeq in published research, please include this citation:

citation("bcbioRNASeq")
## 
## To cite bcbioRNASeq in publications use:
## 
##   Steinbaugh MJ, Pantano L, Kirchner RD, Barrera V, Chapman BA,
##   Piper ME, Mistry M, Khetani RS, Rutherford KD, Hoffman O,
##   Hutchinson JN, Ho Sui SJ. (2018). bcbioRNASeq: R package for
##   bcbio RNA-seq analysis. F1000Research, 6:1976. URL
##   https://f1000research.com/articles/6-1976/v2. DOI
##   10.12688/f1000research.12093.2.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {{bcbioRNASeq}: {R} package for bcbio {RNA-seq} analysis},
##     author = {Michael J. Steinbaugh and Lorena Pantano and Rory D. Kirchner and Victor Barrera and Brad A. Chapman and Mary E. Piper and Meeta Mistry and Radhika S. Khetani and Kayleigh D. Rutherford and Oliver Hofmann and John N. Hutchinson and Shannan J. Ho Sui},
##     journal = {F1000Research},
##     volume = {6},
##     number = {1976},
##     year = {2018},
##     url = {https://f1000research.com/articles/6-1976/v2},
##     doi = {10.12688/f1000research.12093.2},
##   }
library(bcbioRNASeq)

Loading bcbio data

The bcbioRNASeq constructor function is the main interface connecting bcbio] output data to interactive use in R. It is highly customizable and supports a number of options for advanced use cases. Consult the documentation for additional details.

help(topic = "bcbioRNASeq", package = "bcbioRNASeq")

Upload directory

We have designed the constructor to work as simply as possible by default. The only required argument is uploadDir, the path to the bcbio final upload directory specified with upload: in the YAML configuration. Refer to the bcbio configuration documentation for detailed information on how to set up a bcbio run, which is outside the scope of this vignette.

For example, let’s load up the example bcbio dataset stored internally in the package.

uploadDir <- system.file("extdata", "bcbio", package = "bcbioRNASeq")

bcbio outputs RNA-seq data in a standardized directory structure, which is described in detail in our workflow paper.

list.files(path = uploadDir, full.names = FALSE, recursive = TRUE)
##  [1] "2018-03-18_GSE65267-merged/bcbio-nextgen-commands.log"
##  [2] "2018-03-18_GSE65267-merged/bcbio-nextgen.log"         
##  [3] "2018-03-18_GSE65267-merged/combined.counts"           
##  [4] "2018-03-18_GSE65267-merged/programs.txt"              
##  [5] "2018-03-18_GSE65267-merged/project-summary.yaml"      
##  [6] "2018-03-18_GSE65267-merged/tx2gene.csv"               
##  [7] "control_rep1/salmon/quant.sf"                         
##  [8] "control_rep2/salmon/quant.sf"                         
##  [9] "control_rep3/salmon/quant.sf"                         
## [10] "fa_day7_rep1/salmon/quant.sf"                         
## [11] "fa_day7_rep2/salmon/quant.sf"                         
## [12] "fa_day7_rep3/salmon/quant.sf"                         
## [13] "sample_metadata.csv"

Counts level

By default, bcbioRNASeq imports counts at gene level, which are required for standard differential expression analysis (level = "genes"). For pseudo-aligned counts (e.g. Salmon, Kallisto, Sailfish) (Bray et al. 2016; Patro, Mount, and Kingsford 2014; Patro et al. 2017), tximport (Soneson, Love, and Robinson 2016) is used internally to aggregate transcript-level counts to gene-level counts, and generates length-scaled transcripts per million (TPM) values. For aligned counts processed with featureCounts (Liao, Smyth, and Shi 2014) (e.g. STAR, HISAT2) (Dobin et al. 2013; Dobin and Gingeras 2016; Kim, Langmead, and Salzberg 2015), these values are already returned at gene level, and therefore not handled by tximport. Once the gene-level counts are imported during the bcbioRNASeq call, the DESeq2 package (Love, Huber, and Anders 2014) is then used to generate an internal DESeqDataSet from which we derive normalized and variance-stabilized counts.

object <- bcbioRNASeq(uploadDir = uploadDir, level = "genes")
print(object)
## bcbioRNASeq 0.3.28
## uploadDir: /tmp/Rtmpal4Y6i/temp_libpath819e49ffa1bb/bcbioRNASeq/extdata/bcbio
## dates(2): [bcbio] 2018-03-18; [R] 2019-09-17
## level: genes
## caller: salmon
## organism: Mus musculus
## interestingGroups: sampleName
## class: RangedSummarizedExperiment 
## dim: 100 6 
## metadata(27): allSamples bcbioCommandsLog ... wd yaml
## assays(6): counts aligned ... normalized vst
## rownames(100): ENSMUSG00000000001 ENSMUSG00000000003 ...
##   ENSMUSG00000062661 ENSMUSG00000074340
## rowData names(0):
## colnames(6): control_rep1 control_rep2 ... fa_day7_rep2
##   fa_day7_rep3
## colData names(26): averageInsertSize averageReadLength ...
##   treatment x5x3Bias

Alternatively, if you want to perform transcript-aware analysis, such as differential exon usage or splicing analysis, transcript-level counts can be obtained using level = "transcripts". Note that when counts are loaded at transcript level, TPMs are generated with tximport internally, but no additional normalizations or transformations normally calculated for gene-level counts with DESeq2 are generated.

object <- bcbioRNASeq(uploadDir = uploadDir, level = "transcripts", fast = TRUE)
print(object)
## bcbioRNASeq 0.3.28
## uploadDir: /tmp/Rtmpal4Y6i/temp_libpath819e49ffa1bb/bcbioRNASeq/extdata/bcbio
## dates(2): [bcbio] 2018-03-18; [R] 2019-09-17
## level: transcripts
## caller: salmon
## organism: Mus musculus
## interestingGroups: sampleName
## class: RangedSummarizedExperiment 
## dim: 100 6 
## metadata(27): allSamples bcbioCommandsLog ... wd yaml
## assays(3): counts avgTxLength tpm
## rownames(100): ENSMUST00000000001 ENSMUST00000000003 ...
##   ENSMUST00000000674 ENSMUST00000000687
## rowData names(0):
## colnames(6): control_rep1 control_rep2 ... fa_day7_rep2
##   fa_day7_rep3
## colData names(26): averageInsertSize averageReadLength ...
##   treatment x5x3Bias

Expression callers

Since bcbio is flexible and supports a number of expression callers, we have provided advanced options in the bcbioRNASeq constructor to support a variety of workflows using the caller argument. Salmon, Kallisto, and Sailfish counts are supported at either gene or transcript level. Internally, these are loaded using tximport. STAR and HISAT2 aligned counts processed with featureCounts are also supported, but only at gene level.

object <- bcbioRNASeq(uploadDir = uploadDir, caller = "star", fast = TRUE)
print(object)
## bcbioRNASeq 0.3.28
## uploadDir: /tmp/Rtmpal4Y6i/temp_libpath819e49ffa1bb/bcbioRNASeq/extdata/bcbio
## dates(2): [bcbio] 2018-03-18; [R] 2019-09-17
## level: genes
## caller: star
## organism: Mus musculus
## interestingGroups: sampleName
## class: RangedSummarizedExperiment 
## dim: 100 6 
## metadata(27): allSamples bcbioCommandsLog ... wd yaml
## assays(1): counts
## rownames(100): ENSMUSG00000000001 ENSMUSG00000000003 ...
##   ENSMUSG00000000563 ENSMUSG00000000567
## rowData names(0):
## colnames(6): control_rep1 control_rep2 ... fa_day7_rep2
##   fa_day7_rep3
## colData names(26): averageInsertSize averageReadLength ...
##   treatment x5x3Bias

Sample selection and metadata

If you’d like to load up only a subset of samples, this can be done easily using the samples argument. Note that the character vector declared here must match the description column specified in the sample metadata. Conversely, if you’re working with a large dataset and you simply want to drop a few samples, this can be accomplished with the censorSamples argument.

When working with a bcbio run that has incorrect or outdated metadata, the simplest way to fix this issue is to pass in new metadata from an external spreadsheet (CSV or Excel) using the sampleMetadataFile argument. Note that this can also be used to subset the bcbio dataset, similar to the samples argument (see above), based on the rows that are included in the spreadsheet.

sampleMetadataFile <- file.path(uploadDir, "sample_metadata.csv")
stopifnot(file.exists(sampleMetadataFile))
import(file = sampleMetadataFile)
##                   fileName  description  treatment day
## 1 control_rep1_R1.fastq.gz control_rep1    control   0
## 2 control_rep2_R1.fastq.gz control_rep2    control   0
## 3 control_rep3_R1.fastq.gz control_rep3    control   0
## 4 fa_day7_rep1_R1.fastq.gz fa_day7_rep1 folic acid   7
## 5 fa_day7_rep2_R1.fastq.gz fa_day7_rep2 folic acid   7
## 6 fa_day7_rep3_R1.fastq.gz fa_day7_rep3 folic acid   7
object <- bcbioRNASeq(
    uploadDir = uploadDir,
    sampleMetadataFile = sampleMetadataFile,
    fast = TRUE
)
sampleData(object)
## DataFrame with 6 rows and 5 columns
##                   day                 fileName   sampleName  treatment
##              <factor>                 <factor>     <factor>   <factor>
## control_rep1        0 control_rep1_R1.fastq.gz control_rep1    control
## control_rep2        0 control_rep2_R1.fastq.gz control_rep2    control
## control_rep3        0 control_rep3_R1.fastq.gz control_rep3    control
## fa_day7_rep1        7 fa_day7_rep1_R1.fastq.gz fa_day7_rep1 folic acid
## fa_day7_rep2        7 fa_day7_rep2_R1.fastq.gz fa_day7_rep2 folic acid
## fa_day7_rep3        7 fa_day7_rep3_R1.fastq.gz fa_day7_rep3 folic acid
##              interestingGroups
##                       <factor>
## control_rep1      control_rep1
## control_rep2      control_rep2
## control_rep3      control_rep3
## fa_day7_rep1      fa_day7_rep1
## fa_day7_rep2      fa_day7_rep2
## fa_day7_rep3      fa_day7_rep3

Genome annotations

When analyzing a dataset against a well-annotated genome, we recommend importing the corresponding metadata using AnnotationHub and ensembldb. This functionality is natively supported in the bcbioRNASeq constructor with using the organism, ensemblRelease, and genomeBuild arguments. For example, with our internal bcbio dataset, we’re analyzing counts generated against the Ensembl Mus musculus GRCm38 genome build (release 87). These parameters can be defined in the object load call to ensure that the annotations match up exactly with the genome used.

object <- bcbioRNASeq(
    uploadDir = uploadDir,
    level = "genes",
    organism = "Mus musculus",
    genomeBuild = "GRCm38",
    ensemblRelease = 87
)

This will return a GRanges object using the GenomicRanges package (Lawrence et al. 2013), which contains coordinates and rich metadata for each gene or transcript. These annotations are accessible with the rowRanges and rowData functions defined in the SummarizedExperiment package (Huber et al. 2015).

summary(rowRanges(object))
## [1] "GRanges object with 100 ranges and 7 metadata columns"
summary(rowData(object))
##    Length     Class      Mode 
##         7 DataFrame        S4

Alternatively, transcript-level annotations can also be obtained automatically using this method.

When working with a dataset generated against a poorly-annotated or non-standard genome, we provide a fallback method for loading gene annotations from a general feature format (GFF) file with the gffFile argument. If possible, we recommend providing a general transfer format (GTF) file, which is identical to GFF version 2. GFFv3 is more complicated and non-standard, but Ensembl GFFv3 files are also supported.

If your dataset contains transgenes (e.g. EGFP, TDTOMATO) or spike-ins (e.g. ERCCs), these features can be defined with the transgeneNames and spikeNames arguments, which will automatically populate the rowRanges slot with placeholder metadata.

We recommend loading up data per genome in its own bcbioRNASeq object when possible, so that rich metadata can be imported easily. In the edge case where you need to look at multiple genomes simultaneously, set organism = NULL, and bcbioRNASeq will skip the gene annotation acquisition step.

Refer to the the GenomicRanges and SummarizedExperiment package documentation for more details on working with the genome annotations defined in the rowRanges slot of the object. Here are some useful examples:

seqnames(object)
## factor-Rle of length 100 with 64 runs
##   Lengths:  1  1  1  1  1  1  1  1  5  1 ...  1  1  1  2  1  1  1  1  2  1
##   Values :  3  X 16 11  6 13  4  9 11 17 ...  6  5  9 11  8 11  7 15  2  3
## Levels(86): 3 X 16 ... CHR_MG3836_PATCH CHR_MG3496_PATCH CHR_MG3172_PATCH
ranges(object)
## IRanges object with 100 ranges and 0 metadata columns:
##                          start       end     width
##                      <integer> <integer> <integer>
##   ENSMUSG00000000001 108107280 108146146     38867
##   ENSMUSG00000000003  77837901  77853623     15723
##   ENSMUSG00000000028  18780447  18811987     31541
##   ENSMUSG00000000049 108343354 108414396     71043
##   ENSMUSG00000000058  17281185  17289115      7931
##                  ...       ...       ...       ...
##   ENSMUSG00000048583 142650766 142666816     16051
##   ENSMUSG00000055022  92051165  92341967    290803
##   ENSMUSG00000061689 156613705 156764363    150659
##   ENSMUSG00000062661  31245823  31295989     50167
##   ENSMUSG00000074340 105973711 105987423     13713
strand(object)
## factor-Rle of length 100 with 42 runs
##   Lengths: 3 7 1 3 2 2 1 2 3 1 2 6 2 3 2 ... 3 1 4 1 1 2 2 2 4 2 2 1 4 2 4
##   Values : - + - + - + - + - + - + - + - ... + - + - + - + - + - + - + - +
## Levels(3): + - *

Variance stabilization

During the bcbioRNASeq constructor call, log2 variance stabilizaton of gene-level counts can be calculated automatically, and is recommended. This is performed internally by the DESeq2 package, using the varianceStabilizingTransformation and/or rlog functions. These transformations will be slotted into assays.

summary(counts(object, normalized = "vst"))
##   control_rep1     control_rep2     control_rep3     fa_day7_rep1   
##  Min.   : 4.919   Min.   : 4.919   Min.   : 4.919   Min.   : 4.919  
##  1st Qu.: 4.919   1st Qu.: 4.919   1st Qu.: 4.919   1st Qu.: 4.919  
##  Median : 6.051   Median : 5.658   Median : 6.009   Median : 6.237  
##  Mean   : 6.713   Mean   : 6.706   Mean   : 6.715   Mean   : 6.762  
##  3rd Qu.: 7.990   3rd Qu.: 8.239   3rd Qu.: 7.815   3rd Qu.: 8.239  
##  Max.   :12.215   Max.   :12.972   Max.   :11.818   Max.   :11.764  
##   fa_day7_rep2     fa_day7_rep3   
##  Min.   : 4.919   Min.   : 4.919  
##  1st Qu.: 5.101   1st Qu.: 4.919  
##  Median : 6.216   Median : 6.185  
##  Mean   : 6.754   Mean   : 6.789  
##  3rd Qu.: 8.033   3rd Qu.: 7.930  
##  Max.   :12.086   Max.   :11.997

Use case dataset

To demonstrate the functionality and configuration of the package, we have taken an experiment from the Gene Expression Omnibus (GEO) public repository of expression data to use as an example use case. The RNA-seq data is from a study of acute kidney injury in a mouse model (GSE65267) (Craciun et al. 2016). The study aims to identify differentially expressed genes in progressive kidney fibrosis and contains samples from mouse kidneys at several time points (n = 3, per time point) after folic acid treatment. From this dataset, we are using a subset of the samples for our use case: before folic acid treatment, and 1, 3, 7 days after treatment.

For the vignette, we are loading a pre-computed version of the example bcbioRNASeq object used in our workflow paper.

object <- import(file = pasteURL(
    "github.com", "hbc", "bcbioRNASeq",
    "raw", "f1000v2", "data", "bcb.rda",
    protocol = "https"
))
object <- updateObject(object)
print(object)
## bcbioRNASeq 0.3.28
## uploadDir: /n/data1/cores/bcbio/bcbioRNASeq/F1000v2/GSE65267-merged/final
## dates(2): [bcbio] 2018-03-18; [R] 2018-06-01
## level: genes
## caller: salmon
## organism: Mus musculus
## interestingGroups: day
## class: RangedSummarizedExperiment 
## dim: 51652 12 
## metadata(28): version level ... subset previousVersion
## assays(6): counts tpm ... rlog vst
## rownames(51652): ENSMUSG00000000001 ENSMUSG00000000003 ...
##   ENSMUSG00000114967 ENSMUSG00000114968
## rowData names(7): geneID geneName ... entrezID broadClass
## colnames(12): control_rep1 control_rep2 ... fa_day7_rep2
##   fa_day7_rep3
## colData names(25): sampleName day ... x5x3Bias xGC

Sample metadata

For reference, let’s take a look at the sample metadata. By comparison, using colData will return all sample-level metadata, including our quality control metrics generated by bcbio. We recommend instead using sampleData with the clean = TRUE argument in reports, which only returns factor columns of interest.

sampleData(object)
## DataFrame with 12 rows and 7 columns
##                sampleName      day replicate   strain   tissue  treatment
##                  <factor> <factor>  <factor> <factor> <factor>   <factor>
## control_rep1 control_rep1        0         1   Balb/c   kidney    control
## control_rep2 control_rep2        0         2   Balb/c   kidney    control
## control_rep3 control_rep3        0         3   Balb/c   kidney    control
## fa_day1_rep1 fa_day1_rep1        1         1   Balb/c   kidney folic acid
## fa_day1_rep2 fa_day1_rep2        1         2   Balb/c   kidney folic acid
## ...                   ...      ...       ...      ...      ...        ...
## fa_day3_rep2 fa_day3_rep2        3         2   Balb/c   kidney folic acid
## fa_day3_rep3 fa_day3_rep3        3         3   Balb/c   kidney folic acid
## fa_day7_rep1 fa_day7_rep1        7         1   Balb/c   kidney folic acid
## fa_day7_rep2 fa_day7_rep2        7         2   Balb/c   kidney folic acid
## fa_day7_rep3 fa_day7_rep3        7         3   Balb/c   kidney folic acid
##              interestingGroups
##                       <factor>
## control_rep1                 0
## control_rep2                 0
## control_rep3                 0
## fa_day1_rep1                 1
## fa_day1_rep2                 1
## ...                        ...
## fa_day3_rep2                 3
## fa_day3_rep3                 3
## fa_day7_rep1                 7
## fa_day7_rep2                 7
## fa_day7_rep3                 7

Interesting groups

Groups of interest to be used for coloring and sample grouping in the quality control plots can be defined in the bcbioRNASeq object using the interestingGroups argument in the bcbioRNASeq constructor call. Assignment method support with the interestingGroups function is also provided, which can modify the groups of interest after the object has been created.

The interestingGroups definition defaults to sampleName, which is automatically generated from the bcbio description metadata for demultiplexed bulk RNA-seq samples. In this case, the samples will be grouped and colored uniquely.

Interesting groups must be defined using a character vector and refer to column names defined in the colData slot of the object. Note that the bcbioRNASeq package uses lower camel case formatting for column names (e.g. “sampleName”), so the interesting groups should be defined using camel case, not snake (e.g. “sample_name”) or dotted case (e.g. “sample.name”).

This approach was inspired by the DESeq2 package, which uses the argument intgroup in some functions, such as plotPCA for labeling groups of interest. We took this idea and ran with it for a number of our quality control functions, which are described in detail below.

interestingGroups(object) <- c("treatment", "day")

Quality control

Our workflow paper describes the quality control process in detail. In addition, our Quality Control R Markdown template contains notes detailing each metric. Here were are going into more technical detail regarding how to customize the appearance of the plots. Note that all quality control plotting functions inherit the interestingGroups defined inside the bcbioRNASeq object. You can also change this dynamially for each function call using the interestingGroups argument.

Read counts

plotTotalReads(object)

plotMappedReads(object)

Mapping rates

Note that the overall mapping rate is relatively low per sample, while the exonic mapping rate is acceptable. High quality samples should have low intronic mapping rate, and high values are indicative of sample degradation and/or contamination.

plotMappingRate(object)

plotExonicMappingRate(object)

plotIntronicMappingRate(object)

rRNA mapping rate

Note that the samples here have a high rRNA mapping rate. This can be indicative of the polyA enrichment or ribo depletion protocol not having removed all ribosomal RNA (rRNA) transcripts. This will reduce the number of biologically meaningful reads in the experiment and is best avoided.

plotRRNAMappingRate(object)

5’->3’ bias

RNA-seq data can have specific biases at either the 5’ or 3’ end of sequenced fragments. It is common to see a small amount of bias, especially if polyA enrichment was performed, or if there is any sample degradation. If a large amount of bias is observed here, be sure to analyze the samples with a Bioanalyzer and check the RIN scores.

plot5Prime3PrimeBias(object)

Gene distributions

plotFeaturesDetected(object)

plotGeneSaturation(object)

plotCountsPerFeature(object, geom = "boxplot")

plotCountsPerFeature(object, geom = "density")

Sample similarity

We can visualize sample similarity with principal component analysis (PCA) and hierarchical clustering of sample correlations. These functions support multiple normalization methods with the normalized argument.

plotPCA(object)

plotCorrelationHeatmap(object)

Export

The counts accessor function returns raw counts by default. Multiple normalization methods are accessible with the normalized argument.

raw_counts <- counts(object, normalized = FALSE)
normalized_counts <- counts(object, normalized = TRUE)
tpm <- counts(object, normalized = "tpm")
log2_vst_counts <- counts(object, normalized = "vst")

We are providing bcbioRNASeq method support for the export function, which automatically exports matrices defined in assays and corresponding metadata (e.g. rowData, colData) to disk in a single call.

files <- export(object)
print(files)
## $assays
## $assays$counts
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/counts.csv"
## 
## $assays$tpm
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/tpm.csv"
## 
## $assays$avgTxLength
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/avgTxLength.csv"
## 
## $assays$normalized
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/normalized.csv"
## 
## $assays$rlog
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/rlog.csv"
## 
## $assays$vst
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/assays/vst.csv"
## 
## 
## $colData
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/colData.csv"
## 
## $rowData
## [1] "/data00/draco/acidbase/packages/bcbioRNASeq/vignettes/object/rowData.csv"

Differential expression

bcbioRNASeq integrates with DESeq2 and edgeR for differential expression.

DESeq2

To prepare our dataset for analysis, we need to coerce the bcbioRNASeq object to a DESeqDataSet. We have defined an S4 coercion method that uses the as function in the package.

dds <- as(object, "DESeqDataSet")
print(dds)
## class: DESeqDataSet 
## dim: 51652 12 
## metadata(11): version caller ... date sessionInfo
## assays(1): counts
## rownames(51652): ENSMUSG00000000001 ENSMUSG00000000003 ...
##   ENSMUSG00000114967 ENSMUSG00000114968
## rowData names(7): geneID geneName ... entrezID broadClass
## colnames(12): control_rep1 control_rep2 ... fa_day7_rep2
##   fa_day7_rep3
## colData names(25): sampleName day ... x5x3Bias xGC

Since both bcbioRNASeq and DESeqDataSet classes extend RangedSummarizedExperiment, internally we coerce the original bcbioRNASeq object to a RangedSummarizedExperiment and then use the DESeqDataSet constructor, which requires a SummarizedExperiment with integer counts.

The source code for our bcbioRNASeq to DESeqDataSet coercion is accessible with getMethod.

getMethod(
    f = "coerce",
    signature = signature(from = "bcbioRNASeq", to = "DESeqDataSet")
)

We also provide a parameterized starter R Markdown template for standard DESeq2 differential expression, supporting analysis and visualization of multiple contrasts inside a single report. In our workflow paper, we also describe an example LRT analysis.

Now you can perform differential expression analysis with DESeq2. The authors of that package have provided a number of detailed, well documented references online:

edgeR

To prepare our dataset for analysis, we need to coerce the bcbioRNASeq object to a DGEList. Similar to our DESeq2 approach, we have defined an S4 coercion method that hands off to edgeR.

dge <- as(object, "DGEList")
summary(dge)
##         Length Class      Mode   
## counts  619824 -none-     numeric
## samples      3 data.frame list   
## offset  619824 -none-     numeric

The source code for our bcbioRNASeq to DGEList coercion is accessible with getMethod.

getMethod(
    f = "coerce",
    signature = signature(from = "bcbioRNASeq", to = "DGEList")
)

Functional analysis

We provide a starter R Markdown template for gene set enrichment analysis (GSEA) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (Kanehisa and Goto 2000; Subramanian et al. 2005), which leverages the functionality of the clusterProfiler package (Yu et al. 2012). This workflow is described in detail in our workflow paper, and is outside the scope of this advanced use vignette. For more information on functional analysis, consult the clusterProfiler vignette or Stephen Turner’s excellent DESeq2 to fgsea workflow.

R session information

utils::sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.3.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] bcbioRNASeq_0.3.28          basejump_0.11.15           
##  [3] SummarizedExperiment_1.14.1 DelayedArray_0.10.0        
##  [5] BiocParallel_1.18.1         matrixStats_0.55.0         
##  [7] Biobase_2.44.0              GenomicRanges_1.36.1       
##  [9] GenomeInfoDb_1.20.0         IRanges_2.18.2             
## [11] S4Vectors_0.22.1            BiocGenerics_0.30.0        
## [13] BiocStyle_2.12.0           
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.4               Hmisc_4.2-0                  
##   [3] AnnotationHub_2.16.1          BiocFileCache_1.8.0          
##   [5] plyr_1.8.4                    lazyeval_0.2.2               
##   [7] splines_3.6.1                 DESeqAnalysis_0.2.8          
##   [9] ggplot2_3.2.1                 digest_0.6.20                
##  [11] ensembldb_2.8.0               htmltools_0.3.6              
##  [13] magrittr_1.5                  checkmate_1.9.4              
##  [15] memoise_1.1.0                 cluster_2.1.0                
##  [17] syntactic_0.2.6               limma_3.40.6                 
##  [19] transformer_0.2.7             readr_1.3.1                  
##  [21] annotate_1.62.0               Biostrings_2.52.0            
##  [23] R.utils_2.9.0                 pkgdown_1.4.1                
##  [25] prettyunits_1.0.2             colorspace_1.4-1             
##  [27] blob_1.2.0                    rappdirs_0.3.1               
##  [29] ggrepel_0.8.1                 xfun_0.9                     
##  [31] dplyr_0.8.3                   jsonlite_1.6                 
##  [33] tximport_1.12.3               crayon_1.3.4                 
##  [35] RCurl_1.95-4.12               genefilter_1.66.0            
##  [37] zeallot_0.1.0                 survival_2.44-1.1            
##  [39] glue_1.3.1                    gtable_0.3.0                 
##  [41] zlibbioc_1.30.0               XVector_0.24.0               
##  [43] UpSetR_1.4.0                  SingleCellExperiment_1.6.0   
##  [45] scales_1.0.0                  pheatmap_1.0.12              
##  [47] edgeR_3.26.8                  DBI_1.0.0                    
##  [49] Rcpp_1.0.2                    xtable_1.8-4                 
##  [51] progress_1.2.2                htmlTable_1.13.1             
##  [53] foreign_0.8-72                bit_1.1-14                   
##  [55] Formula_1.2-3                 htmlwidgets_1.3              
##  [57] httr_1.4.1                    RColorBrewer_1.1-2           
##  [59] acepack_1.4.1                 bioverbs_0.2.10              
##  [61] pkgconfig_2.0.2               XML_3.98-1.20                
##  [63] R.methodsS3_1.7.1             nnet_7.3-12                  
##  [65] dbplyr_1.4.2                  locfit_1.5-9.1               
##  [67] labeling_0.3                  tidyselect_0.2.5             
##  [69] rlang_0.4.0                   later_0.8.0                  
##  [71] AnnotationDbi_1.46.1          munsell_0.5.0                
##  [73] tools_3.6.1                   cli_1.1.0                    
##  [75] RSQLite_2.1.2                 ggridges_0.5.1               
##  [77] evaluate_0.14                 stringr_1.4.0                
##  [79] yaml_2.2.0                    knitr_1.24                   
##  [81] bit64_0.9-7                   fs_1.3.1                     
##  [83] acidplots_0.2.16              bcbioBase_0.6.10             
##  [85] purrr_0.3.2                   AnnotationFilter_1.8.0       
##  [87] mime_0.7                      R.oo_1.22.0                  
##  [89] grr_0.9.5                     biomaRt_2.40.4               
##  [91] brio_0.3.8                    compiler_3.6.1               
##  [93] rstudioapi_0.10               curl_4.0                     
##  [95] interactiveDisplayBase_1.22.0 geneplotter_1.62.0           
##  [97] tibble_2.1.3                  stringi_1.4.3                
##  [99] GenomicFeatures_1.36.4        desc_1.2.0                   
## [101] lattice_0.20-38               ProtGenerics_1.16.0          
## [103] Matrix_1.2-17                 vctrs_0.2.0                  
## [105] pillar_1.4.2                  BiocManager_1.30.4           
## [107] data.table_1.12.2             cowplot_1.0.0                
## [109] bitops_1.0-6                  Matrix.utils_0.9.7           
## [111] httpuv_1.5.2                  rtracklayer_1.44.4           
## [113] R6_2.4.0                      latticeExtra_0.6-28          
## [115] bookdown_0.13                 promises_1.0.1               
## [117] gridExtra_2.3                 sessioninfo_1.1.1            
## [119] MASS_7.3-51.4                 assertthat_0.2.1             
## [121] freerange_0.2.5               DESeq2_1.24.0                
## [123] rprojroot_1.3-2               withr_2.1.2                  
## [125] GenomicAlignments_1.20.1      Rsamtools_2.0.0              
## [127] GenomeInfoDbData_1.2.1        hms_0.5.1                    
## [129] grid_3.6.1                    rpart_4.1-15                 
## [131] rmarkdown_1.15                DelayedMatrixStats_1.6.1     
## [133] goalie_0.3.8                  shiny_1.3.2                  
## [135] base64enc_0.1-3

References

The papers and software cited in our workflows are available as a shared library on Paperpile.

Bray, Nicolas L, Harold Pimentel, Páll Melsted, and Lior Pachter. 2016. “Near-Optimal Probabilistic RNA-seq Quantification.” Nat. Biotechnol. 34 (5) (May): 525–527. doi:10.1038/nbt.3519. http://dx.doi.org/10.1038/nbt.3519.

Craciun, Florin L, Vanesa Bijol, Amrendra K Ajay, Poornima Rao, Ramya K Kumar, John Hutchinson, Oliver Hofmann, et al. 2016. “RNA Sequencing Identifies Novel Translational Biomarkers of Kidney Fibrosis.” J. Am. Soc. Nephrol. 27 (6) (June): 1702–1713. doi:10.1681/ASN.2015020225. http://dx.doi.org/10.1681/ASN.2015020225.

Dobin, Alexander, and Thomas R Gingeras. 2016. “Optimizing RNA-Seq Mapping with STAR.” Methods Mol. Biol. 1415: 245–262. doi:10.1007/978-1-4939-3572-7\_13. http://dx.doi.org/10.1007/978-1-4939-3572-7_13.

Dobin, Alexander, Carrie A Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R Gingeras. 2013. “STAR: ultrafast Universal RNA-seq Aligner.” Bioinformatics 29 (1) (January): 15–21. doi:10.1093/bioinformatics/bts635. http://dx.doi.org/10.1093/bioinformatics/bts635.

Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nat. Methods 12 (2) (February): 115–121. doi:10.1038/nmeth.3252. http://dx.doi.org/10.1038/nmeth.3252.

Kanehisa, M, and S Goto. 2000. “KEGG: kyoto Encyclopedia of Genes and Genomes.” Nucleic Acids Res. 28 (1) (January): 27–30. https://www.ncbi.nlm.nih.gov/pubmed/10592173.

Kim, Daehwan, Ben Langmead, and Steven L Salzberg. 2015. “HISAT: a Fast Spliced Aligner with Low Memory Requirements.” Nat. Methods 12 (4) (April): 357–360. doi:10.1038/nmeth.3317. http://dx.doi.org/10.1038/nmeth.3317.

Lawrence, Michael, Wolfgang Huber, Hervé Pagès, Patrick Aboyoun, Marc Carlson, Robert Gentleman, Martin T Morgan, and Vincent J Carey. 2013. “Software for Computing and Annotating Genomic Ranges.” PLoS Comput. Biol. 9 (8) (August): e1003118. doi:10.1371/journal.pcbi.1003118. http://dx.doi.org/10.1371/journal.pcbi.1003118.

Liao, Yang, Gordon K Smyth, and Wei Shi. 2014. “featureCounts: an Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7) (April): 923–930. doi:10.1093/bioinformatics/btt656. http://dx.doi.org/10.1093/bioinformatics/btt656.

Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biol. 15 (12): 550. doi:10.1186/s13059-014-0550-8. http://dx.doi.org/10.1186/s13059-014-0550-8.

Patro, Rob, Geet Duggal, Michael I Love, Rafael A Irizarry, and Carl Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nat. Methods 14 (4) (April): 417–419. doi:10.1038/nmeth.4197. http://dx.doi.org/10.1038/nmeth.4197.

Patro, Rob, Stephen M Mount, and Carl Kingsford. 2014. “Sailfish Enables Alignment-Free Isoform Quantification from RNA-seq Reads Using Lightweight Algorithms.” Nat. Biotechnol. 32 (5) (May): 462–464. doi:10.1038/nbt.2862. http://dx.doi.org/10.1038/nbt.2862.

Soneson, Charlotte, Michael I Love, and Mark D Robinson. 2016. “Differential Analyses for RNA-seq: transcript-Level Estimates Improve Gene-Level Inferences.” F1000Res. 4 (December). doi:10.12688/f1000research.7563.1. https://f1000research.com/articles/4-1521/v1/pdf.

Steinbaugh, Michael J, Lorena Pantano, Rory D Kirchner, Victor Barrera, Brad A Chapman, Mary E Piper, Meeta Mistry, et al. 2018. “bcbioRNASeq: R Package for Bcbio RNA-seq Analysis.” F1000Res. 6 (June). doi:10.12688/f1000research.12093.2. https://f1000research.com/articles/6-1976/v2/pdf.

Subramanian, Aravind, Pablo Tamayo, Vamsi K Mootha, Sayan Mukherjee, Benjamin L Ebert, Michael A Gillette, Amanda Paulovich, et al. 2005. “Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc. Natl. Acad. Sci. U. S. A. 102 (43) (October): 15545–15550. doi:10.1073/pnas.0506580102. http://dx.doi.org/10.1073/pnas.0506580102.

Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters.” OMICS 16 (5) (May): 284–287. doi:10.1089/omi.2011.0118. http://dx.doi.org/10.1089/omi.2011.0118.