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 about 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 named "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/
    2018-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 package (source code is available at our GitHub repository). Use the bcbioRNASeq() constructor function (see example below) to create a structured S4 object that contains all of the necessary information for downstream analysis. The only required argument when creating this object is uploadDir, which specifies the path to the bcbio final upload directory. Note that bcbioRNASeq will transform all sample metadata column names to lowerCamelCase format without spaces, dashes, periods or underscores; therefore the interestingGroups argument should be specified in the same format. The values that can be used for interestingGroups are the same as the column names found in the CSV file used to bcbio_nextgen.py. Additionally, specifying the organism of the dataset is strongly recommended, and must use the full Latin name (e.g. Mus musculus); this enables automatic downloading of gene annotations from Ensembl. Once the S4 object is assigned, use the saveData() function to write the dataset to disk as an R Data file. Note that the following code block does not need to be run to reproduce the figures in this paper; its purpose is to describe how to load the data from a bcbio analysis into R using this package, facilitating use of the downstream functions.

# Non-working example demonstrating how to load a bcbio run.
# Load the pre-saved object instead (See use case below).
library(bcbioRNASeq)
bcb <- bcbioRNASeq(
    uploadDir = file.path("gse65267", "final"),
    organism = "Mus musculus",
    interestingGroups = "day"
)
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 analyses, 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 (Figure 1) that is an extension of RangedSummarizedExperiment¹⁵. The assays() slot contains Salmon quasi-alignment data imported with tximport¹³, and automatically generated DESeq2¹⁶ count transformations that provide support for quality control plots. To see all available assays in the bcbioRNASeq object listed by name, use assayNames(bcb). These matrices are described in more detail below:

Figure 1. bcbioRNASeq S4 object structure.

The bcbioRNASeq object is an S4 class extension of the RangedSummarizedExperiment container. The assays slot contains several matrices with sample IDs as column names and gene IDs as row names, including raw counts and various derivations of the raw counts, and can be accessed via the assays() function. Gene IDs tie these assays to additional information about the genes stored as a GRanges object, and can be accessed via the rowRanges() function. Similarly, sample IDs tie the assays to further sample information such as run metrics from bcbio and experimental factors also stored as a DataFrame, and can be accessed via the colData() function. Non-gene or sample specific metadata about the entire run such as tool versions and genome builds is stored as a list and can be accessed via the metadata() function.

Exporting quantified data

The package contains multiple convenience functions to extract the expression abundances described above. 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 18; 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, normalized = FALSE)
normalized <- counts(bcb, normalized = TRUE)
tpm <- counts(bcb, normalized = "tpm")
rlog <- counts(bcb, normalized = "rlog")
vst <- counts(bcb, normalized = "vst")
saveData(raw, normalized, tpm, rlog, vst)
writeCounts(raw, normalized, tpm, rlog, vst)

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 runs as part of the bcbio pipeline and 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 dashed line is used as a cutoff threshold, helping to identify samples that require additional inspection. These default suggested cutoff values are suitable for human/mouse RNA-seq data sets (based on our experience), but should be manually modified for other species. The package does not focus on automatically determining threshold values, but instead provides a mechanism by which to visually communicate the quality of the data. The default cutoffs for these thresholds can be changed using function arguments. The metrics() function allows the user to extract a data.frame containing all metrics information used by the QC report functions. 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 utility. By default, normalized counts generated with the DESeq2 varianceStabilizingTransformation() function is used to plot gene expression. This can be changed using the option named normalized in each function used to plot expression values.

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 2A–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. Some figures show lines indicating the optimal minimum number for the quality metric.

plotTotalReads(bcb)
plotMappingRate(bcb)

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

Sample classes, as defined by the interestingGroups argument, are represented by the different colors defined in the legend of each plot. Vertical dashed black lines indicate suggested cutoff values. The total reads plot (A) indicates the total number of reads sequenced per sample. There is some variability, but in all cases the coverage is higher than 20 million reads, which we consider as sufficient to give good quality gene quantification. (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 at least one mapped read. Optimal minimum 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. For this specific data, it shows no saturation yet for gene detection, meaning that more genes could be detected if coverage were increased for the samples with lower number of reads.

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 rRNA contamination rates below 10% (not shown) (Figure 2C–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 2D–E). 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. For this figure, normalized counts comes from the trimmed mean of M-values, calculated by edgeR. This is output is from the function cpm(normalized.lib.sizes = TRUE) after applying calcNormFactors to a DGEList object containing the raw counts. We would expect similar count distributions for all genes across the samples unless the library sizes or total RNA expression are different (Figure 3). The plotCountsPerGene() and plotCountDensity() functions provide two ways to visualize this comparison. Zero counts are removed in these figures and the values are log2 scale.

plotCountsPerGene(bcb)
plotCountDensity(bcb)

Figure 3. Gene count distributions.

Normalized count distributions are displayed as box plots (A) and density plots (B). The log2 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 log2 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. This plot shows how all the samples have very similar gene expression profile what makes them comparable for the differential expression analysis.

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 for downstream visualization. The plotMeanSD() function wraps the output of different variance stabilizing methods (including the varianceStabilizingTransformation() and rlog() transformations from the DESeq2 package) and plots them with the vsn package’s meanSdplot() function²⁰ (Figure 4). For the example data, the vst and rlog transformations work well with respect to reducing the noise of low expression genes, illustrated by a flatter distribution observed in these plots. The user may choose the desired transformation for their figures using the normalized argument.

plotMeanSD(bcb)

Figure 4. Variance stabilizing transformations.

Plots showing the standard deviation of normalized counts using log2() (top left), rlog() (top right), VST (varianceStabilizingTransformation()) (bottom left), and TMM (bottom right) transformations by rank(mean). The red line denotes the running median estimator. As these panels show, the rlog and VST transformations greatly reduce the standard deviation and 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 package (Figure 5).

plotDispEsts(bcb)

Figure 5. 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 (vst) or 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 expressed genes, which can consequently dominate the calculations in the similarity analysis. An rlog transformation will shrink the values of low counts towards the genes’ average across samples, without affecting the high expression genes.

Inter-correlation analysis 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 6), although this is not always the case. We can also identify potential outlier samples using a correlation heatmap, if there are samples that show low correlation with all other samples in the dataset. For more control over the graphic parameters of the correlation heatmap, other packages can be used to generate these plots by using the normalized data, accessed with counts(bcb, normalized = "vst"), as input. One example is the pheatmap() function from the pheatmap package²¹, which underlies this plot.

plotCorrelationHeatmap(bcb)

Figure 6. 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 valuable way to explore both inter- and intra-group clustering (Figure 7). As with the correlation heatmap plots, other packages can be used to generate these kinds of plots from the normalized data, which can be accessed with counts(bcb, normalized = "vst").

plotPCA(bcb, label = FALSE)

Figure 7. 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 8). 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 8). For the example data, no covariate shows a strong correlation with the PCs with FDR < 0.05, indicating that there is no need for correction for any covariates during the differential expression analysis.

plotPCACovariates(bcb, fdr = 0.1)

Figure 8. Principal component covariate correlations.

PCA is performed on transformed counts 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, non of variables show a significant correlation with the any principal component of the data. Asterisks indicates p-value < 0.05.

QC discussion

With the QC figures, outlier samples can be identified and may be removed from the set of samples. For instance, in Figure 2B, if one samples shows a very low mapping rate (< 50%) compared to others, this could indicate an issue with the RNA material or the library preparation. The same reasoning applies to Figure 2C, that shows exonic mapping rates. In this case, low mapping rates may indicate DNA contamination if the rate is very low, and should be removed if only a minority of samples are affected. Moreover, Figure 8 can reveal clustering of samples correlated to some covariates, providing guidance on which variables to add to the formula during differential expression analysis.

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 <- as(bcb, "DESeqDataSet")
design(dds) <- ∼day
dds <- DESeq(dds)
vst <- varianceStabilizingTransformation(dds)
saveData(dds, vst)

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(
    object = dds,
    contrast = c(
        factor = "day",
        numerator = "7",
        denominator = "0"
    ),
    alpha = c(0.1, 0.05)
)

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(
    object = dds,
    name = "day_7_vs_0",
    alpha = 0.05
)
saveData(res)

Mean average (MA) plot

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

plotMeanAverage(res)

Figure 9. 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 10). This function requires an Internet connection to map Ensembl gene IDs to gene names (symbols).

plotVolcano(res)

Figure 10. 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 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 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 11).

plotDEGHeatmap(results = res, counts = vst)

Figure 11. 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²⁷. Results are clustered by both row and column, with sample classes (as set with the interestingGroups parameter) highlighted across the top of the heatmap. Heatmaps are drawn internally 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 12). 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 = "vst",
    log = TRUE,
    ann = c("geneID", "geneName"),
    xs = "day"
)

Figure 12. 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 13). 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).

dds_lrt <- DESeq(dds, test = "LRT", reduced = ˜1)
res_lrt <- results(dds_lrt)

ma <- counts(bcb, "vst")[significants(res, fc = 2), ]	
res_patterns <- degPatterns(
    ma = ma,
    metadata = colData(bcb),
    time = "day",
    minc = 60
)
saveData(dds_lrt, res_lrt, res_patterns)
res_patterns[["plot"]]

Figure 13. 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. Box plots 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 (DEGreport_1.14.1 used for this plot.)

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

subset(res_patterns[["df"]], cluster == 3, 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. This function requires an Internet connection to map Ensembl gene IDs to gene names (symbols).

res_tbl <- resultsTables(res, lfc = 1)

Top tables

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

topTables(res_tbl, n = 5)

File outputs

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

If desired, variance stabilized counts (e.g. rlog, vst) can also be included 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²⁹. The functional analysis includes over-representation analysis to identify significantly enriched biological processes among a list of DE genes, and gene set enrichment analysis (GSEA) to identify pathways significantly enriched for genes with larger fold changes. The functional analysis results suggest genes/pathways that may be involved with the condition of interest; however, the results are not conclusive and all identified processes/pathways require experimental validation.

For over-representation analysis, clusterProfiler requires a vector of background genes tested for differential expression, a vector of significant DE genes, and a named vector of fold changes associated with each DE gene.

# Generate the gene IDs for the DE list of genes and the background list of genes
library(clusterProfiler)
all_genes <- rownames(res) %>%
    .[!is.na(res[["padj"]])] %>%
    as.character()
sig_genes <- significants(
    object = res,
    fc = 1,
    padj = 0.05
)

# Generate fold change values for significant results
sig_results <- as.data.frame(res)[sig_genes, ]
fold_changes <- sig_results$log2FoldChange
names(fold_changes) <- rownames(sig_results)

The enrichGO() function takes the significant gene list, the background gene list, and the ontology to test as input to perform statistical enrichment analysis of gene ontology (GO) terms using hypergeometric testing.

# Run GO enrichment analysis
ego <- enrichGO(
    gene = sig_genes,
    OrgDb = "org.Mm.eg.db",
    keyType = "ENSEMBL",
    ont = "BP"
    universe = all_genes,
    qvalueCutoff = 0.05,
    readable = TRUE
)

The GO enrichment analysis output is an S4 class object named enrichResult. The results can be extracted from this object using the slot() accessor function:

ego_summary <- slot(ego, "result") %>%
    as_tibble() %>%
    camel()

The results table contains the following columns:

There are a variety of plots that can be used to explore the enriched processes using the enrichResult object, including the dot plot, enrichment GO plot, and category netplot.

The dot plot shows the number of DE genes associated with the top 25 GO terms (size) and the P-adjusted values for these terms. The order of GO terms is based on the gene ratio (number of significant genes associated with the GO term / total number of significant genes) (Figure 14).

# Dotplot of top 25
dotplot(ego, showCategory = 25)

Figure 14. Enriched GO terms: dot plot.

The 25 GO processes with the largest gene ratios are plotted in order of gene ratio. The size of the dots represent the number of genes in the significant DE gene list associated with the GO term and the color of the dots represent the P-adjusted values (BH).

The enrichment GO plot shows the relationship between the top 25 most significantly enriched GO terms, by grouping similar terms together (Figure 15).

# Enrichment plot of top 25
enrichMap(ego, n = 25, vertex.label.cex = 0.5)

Figure 15. Enriched GO terms: enrichment GO plot.

The top 25 most significantly enriched GO terms (by P-adjusted values) are shown with similar terms grouped together. The color represents the P values relative to the other displayed terms (brighter red is more significant) and the size of the terms represents the number of genes that are significant from the significant genes list.

The category netplot shows the relationships between the genes associated with the top five most significant GO terms and the fold changes of the significant genes associated with these terms (color) (Figure 16). The size of the GO terms reflects the P values of the terms, with the more significant terms being larger.

# Cnet plot for top 5 most significant GO processes
cnetplot(
    x = ego,
    showCategory = 5,
    foldChange = fold_changes,
    vertex.label.cex = 0.5
) 
     

Figure 16. Enriched GO terms: category netplot.

The top five most significant GO terms (P-adjusted values) are plotted and connected with lines to associated DE genes. The fold changes of the significant genes associated with these terms are represented by color. The size of the GO terms reflects the P values of the terms, with the more significant terms being larger.

The Functional Analysis template also includes code for performing over-representation analysis with GO terms and code to perform GSEA analysis using KEGG and GO terms. The final step in the template is the visualization of the genes and corresponding fold changes associated with the GSEA enriched pathways using the pathview package.