Example Workflow for Bulk RNA-Seq Analysis

Limma-voom workflow

Prepare data

To download and process GEO datasets, use the prepare_geo() function. This will generate a list containing count data, sample information, and gene data.

library(easybio)
# x <- prepare_geo('gseid')

For preparing TCGA RNA-Seq data using the TCGAbiolinks`` package, use theprepare_tcga()` function. This will return a list containing count data for all samples and unstranded FPKM data for tumor samples, along with sample and feature information.

Example Workflow: TCGA CHOL Project

Three functions are designed for this workflow:

library(TCGAbiolinks)
library(SummarizedExperiment)

query <- GDCquery(
  project = "TCGA-CHOL",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification"
)
GDCdownload(query = query)
data <- GDCprepare(query = query)

lt <- prepare_tcga(data)
lt$all$sampleInfo[["group"]] <- fifelse(lt$all$sampleInfo$sample_type %ilike% "Tumor", "Tumor", "Normal")

# limma workflow
x <- dgeList(lt$all$exprCount, lt$all$sampleInfo, lt$all$featuresInfo)
x <- dprocess_dgeList(x, "group", 10)
efit <- limmaFit(x, "group")

CHOL.DEGs <- limma::topTable(fit = efit, coef = 1, number = Inf)

For a detailed overview of the Limma workflow, refer to the article: RNA-seq analysis is as easy as 1-2-3 with limma, Glimma and edgeR.

Next, visualize the differentially expressed genes using the plotVolcano() function.

# View the help page for additional arguments.
plotVolcano(data = CHOL.DEGs, x = logFC, y = -log10(adj.P.Val))
?plotVolcano

Pathway Enrichment Analysis

To perform pathway enrichment analysis such as GO and KEGG, download the r4msigdb package to access the MSigDB gene set.

devtools::install("person-c/r4msigdb")

To know more details about this package, please see r4msigdb. To get GO pathways in MSigDB:

pathwayGO <- r4msigdb::query(species = "Hs", pathway = "^GO(MF)|(BP)|(CC)_")

Gene Set Enrichment Analysis (GSEA)

The core function is adapted from the fgsea package with minor visual enhancements.

library(fgsea)
data(examplePathways)
data(exampleRanks)

Run the GSEA analysis:

fgseaRes <- fgsea(pathways = examplePathways, 
                  stats    = exampleRanks,
                  minSize  = 15,
                  maxSize  = 500)
plotGSEA(
  fgseaRes, 
  pathways = examplePathways, 
  pwayname = "5991130_Programmed_Cell_Death", 
  stats = exampleRanks, 
  save = FALSE
)
#> Warning in fsort(stats, TRUE): New parallel sort has not been implemented for
#> decreasing=TRUE so far. Using one thread.

Over-Representation Analysis

Conduct Over-Representation Analysis (ORA) using the following approach:

foraRes <- fora(examplePathways, genes=tail(names(exampleRanks), 200), universe=names(exampleRanks))

Display the results:

# Adjust the pathway position on the y-axis based on the adjusted p-value (padj)
foraRes[, pathway := factor(pathway, levels = rev(pathway))]

plotORA(data = foraRes[1:8], x = -log10(padj), y = pathway, size = overlap, fill = 'constant')