---
title: "Example Workflow for Bulk RNA-Seq Analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{example-bulk-rna-seq-workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


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

```{r}
library(easybio)
# x <- prepare_geo('gseid')
```

For preparing TCGA RNA-Seq data using the `TCGAbiolinks`` package, use the `prepare_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:

```r
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](https://doi.org/10.12688/f1000research.9005.3).

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

```r
# 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.

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

To know more details about this package, please see [r4msigdb](https://github.com/person-c/r4msigdb). To get GO pathways in MSigDB:

```r
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.

```{r}
library(fgsea)
data(examplePathways)
data(exampleRanks)
```

Run the GSEA analysis:

```{r,fig.width=4, fig.height=2.7}
fgseaRes <- fgsea(pathways = examplePathways, 
                  stats    = exampleRanks,
                  minSize  = 15,
                  maxSize  = 500)
plotGSEA(
  fgseaRes, 
  pathways = examplePathways, 
  pwayname = "5991130_Programmed_Cell_Death", 
  stats = exampleRanks, 
  save = FALSE
)
```

#### Over-Representation Analysis

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

```{r}
foraRes <- fora(examplePathways, genes=tail(names(exampleRanks), 200), universe=names(exampleRanks))
```

Display the results:

```{r, fig.width=8, fig.height=2}
# 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')
```
