bcbioRNASeq: R package for bcbio RNA-seq analysis

Reading data

As noted, bcbio runs a number of tools to generate QC metrics and compute gene counts from RNA-seq data. Additional information on the bcbio RNA-seq pipeline is available on readthedocs). At the end of a bcbio run, the most important files are stored in a separate directory specified by the user in the bcbio configuration YAML file under the upload: parameter; this directory is called "final" by default. Within this directory there is a dated project directory containing quality metrics, provenance information, and data derived from the analysis that have been aggregated across all samples, e.g. count files. In addition, there is a directory corresponding to each sample that contains the binary alignment map (BAM) files and Salmon count data for that sample.

final/
2017-01-01_illumina_rnaseq/
annotated_combined.counts
bcbio-nextgen-commands.log
bcbio-nextgen.log
combined.counts
combined.dexseq
combined.dexseq.ann
data_versions.csv
multiqc/
programs.txt
project-summary.yaml
tx2gene.csv
sample1/
qc/
salmon/
sample1-ready.bam
sample1-ready.bam.bai
sample1-ready.counts
sample1-transcriptome.bam

The final upload directory generated by bcbio is used as the input for bcbioRNASeq. Once the bcbio run is complete, you can open an R session and load the bcbioRNASeq library (available from our GitHub repository). Use the loadRNASeq() function to create a structured S4 object (see below) that contains all of the necessary information for downstream analysis. The only required argument when creating this object is the full path to the final directory. We also recommend that you use the interestingGroups argument to indicate variables that are present in your metadata that are of interest for the analysis. Note that bcbioRNASeq will transform all metadata column headings to lowerCamelCase format without spaces, dashes, periods or underscores; therefore interestingGroups should be specified in the same format. Once the S4 object is created, use the saveData() function to save it as an RData (.rda) file.

> library(bcbioRNASeq)
> bcb <- loadRNASeq(
+      file.path("bcbio_example_run", "final"),
+      interestingGroups = c("genotype", "treatment"))
> saveData(bcb)

This S4 object is unique to the bcbioRNASeq package, as it contains all of the necessary data from the bcbio run required for analysis. From here, you can use various functions in bcbioRNASeq to perform analysis, make figures, and generate data tables and results files as we describe in later sections. This object is also used as the input for the R Markdown templates for report generation. First, we begin by describing the object in more detail.

Object description

We have designed a new S4 class named bcbioRNASeq, which is an extension of SummarizedExperiment¹⁵. Our S4 class adds a slot to SummarizedExperiment named bcbio that facilitates the inclusion of additional objects related to the experiment that cannot be contained in a regular SummarizedExperiment. The bcbio slot allows the incorporation of three additional data structures: the Salmon quasi-alignment data for differential expression analyses from tximport¹³, an automatically generated DESeqDataSet to provide support for quality control plots, and alternative counts generated by featureCounts. To see all available slots in the bcbioRNASeq object listed by name, you can use slotNames(bcb). Each of the slots are described in more detail below:

Exporting quantified data

Functions exist to extract expression abundances in the various formats. Outlined below are the steps to save these counts external to the bcbioRNASeq object. These steps utilize functions from the DESeq2, edgeR and tximport packages both directly as well as within wrapper functions. For discussions on RNA-seq data normalization methods and count formats see 19; we typically save at least the DESeq2 normalized counts (library size adjusted) and transcripts per million counts (gene length adjusted) for further analyses.

> raw <- counts(bcb, "raw")
> normalized <- counts(dds, "normalized")
> rlog <- counts(bcb, "rlog")
> tpm <- tpm(txi)
> saveData(raw, normalized, rlog, tpm)
> writeCounts()(raw, normalized, rlog, tpm)

Note, the raw counts generated by featureCounts are different from the Salmon-derived counts slotted in assays(), as described above. This alternative counts matrix will have a different number of rows, and is accessible using bcbio(bcb, "featureCounts").

Quality control steps

A typical RNA-seq analysis requires multiple quality control assessments at the read, alignment, sample and model level. Most of the data required to make these assessments is automatically generated by bcbio; the bcbioRNASeq package makes it easier for users to access it. For instance, the Qualimap tool run as part of the bcbio pipeline generates various metrics that can be used to assess the quality of the data and consistency across samples. The output of Qualimap is stored in the bcbioRNASeq object, and the package has several functions to visualize this output in a graphical format. These plots can be used to check data quality. Visual thresholds appear in many plots to help the user to assess quality. For example, a vertical burgundy line is used as a warning threshold, indicating that any samples with values below that require close attention. A vertical black line threshold represents the optimal value. The default cutoffs for these thresholds can be easily changed using function arguments. Using metrics(bcb) allows the user to extract a data.frame with all metrics information used by the functions in the QC report. In this way, custom figures can also be easily created using the same data but with user-preferred packages.

Below, we provide several examples of recommended QC steps for RNA-seq data with a short explanation outlining their usefulness.

Read statistics

Total reads per sample and mapping rate are metrics that can help identify imbalances in sequencing depth or failures among the samples in a dataset (Figure 1A–B). Generally, we expect to see a similar sequencing depth across all samples in a dataset and mapping rates greater than 75%. Low genomic mapping rates are indicative of sample contamination, poor sequencing quality or other artifacts.

> plotTotalReads(bcb)
> plotMappingRate(bcb)

Figure 1. Read, alignment, genomic context, and gene detection statistics.

Sample classes, as defined by the interestingGroups argument, are represented by the different colors as defined in the legend of each plot. Vertical black lines indicate optimal values while vertical burgundy lines indicate cutoffs for substandard values. The total reads plot (A) indicates the total number of reads sequenced per sample and the mapping rate plot (B) shows the percentage of reads mapping to the reference genome. Here, all samples are well within recommended ranges, having well over 25 million reads and almost 90% of reads mapping. The exonic and intronic mapping rate plots (C and D) indicate the percentage of reads mapping to exons or introns, respectively. Here, all samples are within recommended ranges, with samples from day 3 and day 7 showing higher proportions of reads mapping to intronic versus exonic regions as compared to the day 1 and normal sample classes. The genes detected plot (E) indicates the total number of genes for each sample with a least one mapped read. Optimal gene detection values will vary based on an organism’s transcriptome size. The gene detection saturation plot (F) shows the relationship between the number of reads mapped and the number of genes detected. If this trend is not linear, it indicates that the sequencing may have been saturated in terms of detecting gene expression.

Genomic context

For RNA-seq, the majority of reads should map to exons and not introns. Ideally, at least 60% of total reads should map to exons. High levels of intronic mapping may indicate high proportions of nuclear RNA or DNA contamination. Samples should also have ribosomal RNA (rRNA) contamination rates below 10% (not shown) (Figure 1C–D).

> plotExonicMappingRate(bcb)
> plotIntronicMappingRate(bcb)

Number of genes detected

Determining how many genes are detected relative to the number of mapped reads is another good way to assess the sample quality (Figure 1E–F). Ideally, all samples will have similar numbers for genes detected, and samples with higher number of mapped reads will have more genes detected. Large differences in gene detection numbers between samples can introduce biases and should be monitored at later steps for potential influence on sample clustering.

> plotGenesDetected(bcb)
> plotGeneSaturation(bcb)

Counts per gene

Comparing the distribution of normalized gene counts across samples is one way to assess sample similarity within a dataset. We would expect similar count distributions for all genes across the samples unless the library sizes or total RNA expression are different (Figure 2). The plotCountsPerGene() and plotCountDensity() functions provide two ways to visualize this comparison.

> plotCountsPerGene(bcb)
> plotCountDensity(bcb, style = "line")

Figure 2. Gene count distributions

Normalized count distributions are displayed as boxplots (A) and density plots (B). The log10 TMM-normalized counts per gene normalization method equates the overall expression levels of genes between samples under the assumption that a majority of them are not differentially expressed²⁰. Therefore, by normalizing for total RNA expression by sample, we expect the spread of the log10 TMM-normalized counts per gene to be similar for every sample. Sample classes (as set with the interestingGroups argument) are highlighted in different colors. Here, there is high similarity among the samples.

Model fitting

It is important to explore the fit of the model for a given dataset before performing differential expression analysis. The normalized and transformed data can be used to assess the variance-expression level relationship in the data, to identify which method is best at stabilizing the variance across the mean. The plotMeanSD() function wraps the output of different variance stabilizing methods (including the vst() and rlog() transformations from the DESeq2 package) and plots them with the vsn package’s meanSdplot() function²¹ (Figure 3).

> plotMeanSD(bcb)

Figure 3. Variance stabilizing transformations.

Plots show the standard deviation of normalized counts (normalizedCounts) using log2() (A), rlog() (B), and variance stabilizing (vst()) (C) transformations by rank(mean). The red line shows the running median estimator. As (B) and (C) show, the transformations greatly reduce the standard deviation, and the rlog() transformation is most effective at stabilizing the variance across the mean.

Dispersion

Another plot that is important to evaluate when performing QC on RNA-seq data is the plot of dispersion versus the mean of normalized counts. For a good dataset, we expect the dispersion to decrease as the mean of normalized counts increases for each gene. The plotDispEsts() function provides easy access to model information stored in the bcbioRNASeq object, using the plotting code provided in the DESeq2 library (Figure 4).

> plotDispEsts(bcb)

Figure 4. DESeq2 dispersion plot.

Sample similarity within groups

The QC metrics assessed up to this point are performed to get a global assessment across the dataset and look for similar trends across all samples . However, often we have a dataset in which samples can be classified into groups and it is common to interrogate how similar replicates are to each other within those groups, and the relationship between groups. To this end, bcbioRNASeq provides functions to perform this level of QC with Inter-Correlation Analyses (ICA) and Principal Components Analyses (PCA) between samples. Furthermore, we can use the results of the PCA to identify covariates that correlate with principal components.

Since these analyses are based on variance measures, it is recommended that the variance stabilized rlog transformed counts be used. Using these transformed counts minimizes large differences in sequencing depth and helps normalize all samples to a similar dynamic range . Simple logarithmic transformations of normalized count values tend to generate a lot of noise for low expressor genes, which can consequently dominate the calculations in the similarity analysis. An rlog transformation will shrink the values of low counts towards the genes’ averages across samples, without affecting the high expression genes.

Inter-correlation analysis (lCA) and hierarchical clustering

Inter-Correlation analysis allows us to look at how closely samples are related to each other by first computing pair-wise correlations between expression profiles of all samples in the dataset and then clustering based on those correlation values. Samples that are similar to one another will be highly correlated and will cluster together. We expect samples from the same group to cluster together (Figure 5), although this is not always the case. We can also identify potential outlier samples using ICA, if there are samples that show low correlation with all other samples in the dataset. For more control over the graphing parameters of the ICA heatmap, other packages can be used to generate these plots by using the normalized data, accessed with counts (bcb, normalized = "rlog"), as input. One example is the pheatmap () function from the pheatmap package²², which underlies this plot.

> plotCorrelationHeatmap(bcb)

Figure 5. Sample correlations.

All pairwise sample Pearson correlations are shown. Correlations are clustered by both row and column, with sample classes (as set with the interestingGroups argument) highlighted across the top of the heatmap. Here, the sample classes cluster well, with the normal and day 1 samples showing the highest intra-group correlations.

Principal component analysis (PCA)

PCA is a multivariate technique that allows us to summarize the systematic patterns of variations in the data²³. PCA takes the expression levels for genes and transforms them in principal component space, reducing each sample to a single point. It allows us to separate samples by expression variation , and identify potential sample outliers. The PCA plot is a great way to explore both inter- and intra-group clustering (Figure 6). As with the ICA plots, other packages can be used to generate these kinds of plots from the normalized data, which can be accessed with counts(bcb, normalized = "rlog").

> plotPCA(bcb, label = FALSE)

Figure 6. Principal component analysis.

The first two principal components of the gene expression dataset are plotted here for each of the samples. Sample classes (as set with the interestingGroups argument) are highlighted in different colors. Alternatively, sample labels can be added with the label = TRUE argument (not shown) to identify individual samples, which is particularly useful for identifying outliers. Here, we see good clustering of the samples by group with no apparent outliers.

Correlation of covariates with PCs

When there are multiple factors that can influence the results of a given experiment, it is useful to assess which of them is responsible for the most variance as determined by PCA (Figure 7). The plotPCACovariates() function passes transformed count data and metadata from the bcbioRNASeq object to the degCovariates() function of the DEGreport package²⁴. This method adapts the method described by Daily et al. for which they integrated a method to correlate covariates with principal components values to determine the importance of each factor²⁵ (Figure 7).

> plotPCACovariates(bcb, fdr = 0.1)

Figure 7. Principal component covariate correlations.

PCA is performed on transformed counts (rlog by default; but vst is also supported) and correlations between principal components and the different metadata covariates are computed. Significant correlations (FDR < 0.1) are shaded from blue (anti-correlated) to orange (correlated), with non-significant correlations shaded in gray. Here, none of variables show a significant correlation with any principal component of the data. Asterisks indicates p-value < 0.05.

Differential expression analysis using DESeq2

Once the QC is complete and the dataset looks good, the next step is to identify differentially expressed genes (DEG). For this part of the workflow, we follow instructions and guidelines from the DESeq2 vignette, using Salmon-derived abundance estimates imported with tximport. As previously noted, a Differential Expression R Markdown template is available with the bcbioRNASeq package for these steps.

The first step is to define the factors to include in the statistical analysis as a design formula. We chose to study the difference between the normal group and the day 7 group in our dataset. In our example, we have only one variable of interest; however, DESeq2 is able to model additional covariates. When additional variables are included, the last variable entered in the design formula should generally be the main condition of interest. More detailed instructions and examples are available in the DESeq2 vignette.

> # DESeqDataSet
> dds <- DESeqDataSetFromTximport(
+      txi = bcbio(bcb, "tximport"),
+      colData = colData(bcb),
+      design = formula(~group)
+     ) %>%
+      DESeq()
> saveData(dds)

Alpha level (FDR) cutoffs

The results from DESeq2 include a column for the P values associated with each gene/test as well as a column containing P values that have been corrected for multiple testing (i.e. false discovery rate values). The multiple test correction method performed by default is the Benjamini Hochberg (BH) method²⁶. Since it can be difficult to arbitrarily select an adjusted P value cutoff, the alphaSummary() function is useful for summarizing results for multiple adjusted P value or FDR cutoff values (Table 1).

> alphaSummary(dds,
+     contrast = c(
+          factor = "group",
+          numerator = "day7",
+          denominator = "normal"))

Differential expression analysis

Use the results() function to generate a DESeqResults object containing the output of the differential expression analysis. The desired BH-adjusted P value cutoff value is specified here with the alpha argument (< 0.05 shown).

> res <- results(dds,
+      contrast = c(
+          factor = "group",
+          numerator = "day7",
+          denominator = "normal"),
+      alpha = 0.05)
> saveData(res)

Mean average (MA) plot

The plotMA() function plots the mean of the normalized counts versus the log2 fold changes for all genes tested (Figure 8).

> plotMA(res)

Figure 8. MA plot.

Each point represents one gene, with mean expression levels across all samples plotted on the x-axis and the log2 fold change observed in the contrast of interest on the y-axis. Significant differentially expressed genes (with an adjusted P value less than the adjusted cutoff P value we chose earlier) are colored red²⁷.

Volcano plot

The plotVolcano() function produces a volcano plot comparing significance (here the BH-adjusted P value) for each gene against the fold change (here on a log2 scale) observed in the contrast of interest²⁸ (Figure 9).

> plotVolcano(res)

Figure 9. Volcano plot.

This plot compares the amount of gene expression change (not strictly interpretable here as these fold changes were derived from an LRT test) to the significance of that change (here plotted as the -log10 transformation of the multiple test adjusted P value), with each grey point representing a single gene. Options are available to highlight top gene candidates by shading (here, genes in the green shaded areas are bounded by a minimal fold change and -log10 adjusted P value cutoffs) and by text labeling. The two (optional) marginal plots showing the distributions of the log2 fold changes and negative log10 adjusted P values are useful in assessing cutoff choices and trade-offs.

Gene Expression Heatmap

The plotDEGHeatmap() function produces a gene expression heatmap useful for visualizing the expression of differentially expressed (DE) genes across samples. The heatmap shows only DE genes on a per-sample basis, using an additional log2 fold change cutoff. By default, this plot is scaled by row and uses the ward.D2 method for clustering²⁹. The gene expression heatmap is a nice way to explore the consistency in expression across replicates or differences in expression between sample groups for each of the DE genes (Figure 10).

> # DESeqResults, DESeqTransform
> plotDEGHeatmap(res, counts = assays(bcb)[["rlog"]])

Figure 10. Differentially Expressed Gene Heatmap.

The heatmap shows only DE genes on a per-sample basis, using an additional log2 fold change cutoff. By default, this plot is scaled by row (centered and scaled) and uses the ward.D2 method for clustering²⁹, with red and blue colors denoting higher and lower expression levels respectively. Results are clustered by both row and column. Heatmaps are drawn with the pheatmap() function of the pheatmap package²².

Specific genes plot

In addition to looking at the overall results from the differential expression analysis, it is useful to plot the expression differences for a handful of the top differentially expressed genes. This helps to check the quality of the analysis, by validating the expression for genes that are identified as significant. It is also helpful to visualize trends in expression change across the various sample groups (Figure 11). As bcbioRNASeq is integrated with the DEGreport package²⁴, we can use DEGreport’s degPlot() function to view the expression of individual DE genes.

> degPlot(
+     bcb,
+      res = res,
+      n = 3,
+      slot = "normalized",
+      # Use the ’ann’ argument for gene to symbol mapping
+     ann = c("ensgene", "symbol"),
+      xs = "group")

Figure 11. Individual gene expression patterns.

Gene expression patterns are shown for the top 3 (as selected in the function options) differentially expressed genes (by BH-adjusted P value). Normalized, transformed counts are shown for each replicate from each sample group; groups are set within the function options. Interestingly, even though our analysis compared day 7 to normal, all 3 genes show their greatest increases in gene expression at day 3, with some leveling off or relative decreases in expression at day 7.

Detecting patterns

The full set of example data is from a time course experiment, as described previously. Up to this point, we have only compared gene expression between two time points (normal and day 7), but we can also analyze the whole dataset to identify genes that show any change in expression across the different time points. As recommended by DESeq2, the best approach for this type of experimental design is to perform a likelihood ratio test (LRT) to test for differences in gene expression between any of the sample groups in the context of the time course. More information about time-course experiments and LRT is available in the DESeq2 vignette. This approach will yield a list of differentially expressed genes, but will not report how the expression is changing. Visualizing patterns of expression change amongst the significant genes is helpful in identifying groups of genes that have similar trends, which in turn can help determine a biological reason for the changes we observe (Figure 12). The DEGreport package includes the degPatterns() function, which is designed to extract and plot genes that have a similar trend across the various time points. More information about this function can be found in the DEGreport package²⁴. Note that this function works only with significant genes; significants() returns the significant genes based on log2FoldChange and padj values (0 and 0.05 respectively, by default).

> ddsLRT <- DESeqDataSetFromTximport(
+      txi = bcbio(bcb, "tximport"),
+      colData = colData(bcb),
+      design = formula(~group)) %>%
+      DESeq(test = "LRT", reduced = ~1)
> resLRT <- results(ddsLRT)
> # This step is CPU intensive
> resPatterns <- degPatterns(
+      counts(bcb, "rlog")[significants(res), ],
+      metadata = colData(bcb),
+      time = "group",
+      col = NULL)
> saveData(ddsLRT, resLRT, resPatterns)
> resPatterns[["plot"]]

Figure 12. Clustered expression patterns.

Clustered and scaled expression patterns for the top 500 differentially expressed genes. Each group represents a gene expression pattern shared among different DE genes, with the number of groups determined by the expression correlation patterns of the groups. Boxplots are shown for the expression patterns of each gene within the group to give a better idea of how well the groupings fit the expression data.

The output of the degPatterns() function is the plot, as well as a list object that contains a data.frame with two columns - the "genes" and the corresponding "cluster" number. To extract genes from a specific cluster for further analysis, the base function subset() can be utilized as follows to obtain the list of genes from cluster 8:

> subset(resPatterns[["df"]], cluster == 8, select = "genes")

Summarize analysis

Finally, at the end of the analysis a results table can be extracted containing the DE genes at a specified log fold change threshold. These results can be written to files in the specified output folder. In addition to the DE genes, a detailed summary and description of results and the output files are generated along with the cutoffs used to identify the significant genes.

> resTbl <- resultsTables(res, lfc = 1, write = TRUE, dir = deDir)

Top tables

Once the results table object is ready, the top up- and down-regulated genes (sorted by log2 fold change) can now be displayed. Here, we output the top 5 DE genes in each direction (Table 2).

> topTables(resTbl, n = 5)

File outputs

Output files from this use case include the following gene counts files (output at the count normalization step):

If desired, rlog variance stabilized counts can also be output for future use in variance based plotting methods such as PCA. DEG tables containing the differential expression results from the analysis summary step are sorted by BH-adjusted P value, and contain the following columns:

Functional analysis

To gain greater biological insight into the list of DE genes, it is helpful to perform functional analysis. We provide the Functional Analysis R Markdown template, which contains code from the clusterProfiler package³⁰, to identify potentially enriched biological processes among the DE genes. We use this tool to take the significant gene list and the background gene list (all genes tested) as input to perform statistical enrichment analysis of gene ontology (GO) terms using hypergeometric testing. The template includes a table of GO terms that are significantly enriched among the DE genes and a variety of plots summarizing the significantly enriched GO processes. The plots included are: