11_plot_ora.Rmd

# (PART\*) Part IV: Plotting {-}

# Plot ORA {#plot-ora-1}

```{r include=FALSE}
library(knitr)
library(ggplot2)
library(igraph)
library(ggraph)
library(patchwork)
library(genekitr)
```

After selecting interested terms or pathways from `genORA` or `genGSEA` result, user could pass the data frame to `plotEnrich`, which includes many **ready-made plot types**, including barplot, dotplot, heatmap, wego-like plot, chord plot, network, wordcloud etc.

It is worth mentioning that **almost all plots are based on `r CRANpkg("ggplot2")`** and `plot_theme` function could easily change their border, legend, label etc., which helps user make own plot.

## Get ORA result

> For more details, please refer to [chapter8](#ora-basic-usage)

```{r}
# 1st step: prepare input IDs
data(geneList, package = "genekitr")
entrez_id <- names(geneList)[abs(geneList) > 2]
head(entrez_id, 5)

# 2nd step: prepare gene set
hg_gs <- geneset::getGO(org = "human",ont = "bp")

# 3rd step: ORA analysis
ego <- genORA(entrez_id, geneset = hg_gs, p_cutoff = 0.01, q_cutoff = 0.01)
# next we only show ten sample terms 
ego <- ego[1:10, ]
head(ego)
```


## Bar Plot

Bar plot is the simplest way to show enriched terms. The x-axis is enrichment metric (e.g. Gene ratio; Fold enrichment); the y-axis is selected terms. The bar color represents statistical value.

> **What's the difference between pvalue and p.adjust?**
`p.adjust` has seven adjustment methods: <https://www.rdocumentation.org/packages/stats/versions/3.5.0/topics/p.adjust>.
The adjusted p-value is always the p-value, multiplied with some factor: adj.p = f * p The actual size of this factor f depends on the strategy used to correct for multiple testing.
By the way, the "q-value" stands for "false discovery rate (FDR)" method.


The basic arguments are: 

- `term_metric`: The x-axis could be "Count", "GeneRatio", "FoldEnrich" or "RichFactor"
- `stats_metric`: Statistic value of "p.adjust", "pvalue" or "qvalue"
- `up_color`: Color for stronger statistic value (e.g. pvalue 0.01)
- `down_color`: Color for weaker statistic value (e.g. pvalue 1)
- `wrap_length`: Wrap term text if longer than this number

```{r include= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "bar", main_text_size = 15, legend_text_size = 13)
p2 <- plotEnrich(ego,
  plot_type = "bar", term_metric = "GeneRatio", stats_metric = "pvalue",
  main_text_size = 15, legend_text_size = 13
)
p3 <- plotEnrich(ego,
  plot_type = "bar", up_color = "#E69056", down_color = "#325CAC",
  main_text_size = 15, legend_text_size = 13
)
p4 <- plotEnrich(ego,
  plot_type = "bar", wrap_length = 25,
  main_text_size = 15, legend_text_size = 13
)
pp <- p1 + p2 + p3 + p4
```


```{r eval = F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "bar")
p2 <- plotEnrich(ego, plot_type = "bar", term_metric = "GeneRatio", stats_metric = "pvalue")
p3 <- plotEnrich(ego, plot_type = "bar", up_color = "#E69056", down_color = "#325CAC")
p4 <- plotEnrich(ego, plot_type = "bar", wrap_length = 25)
p1 + p2 + p3 + p4 + plot_annotation(tag_levels = "A")
```


(ref:barplotScap) Bar plot of enrichment analysis.

(ref:barplotCap) **Bar plot of enrichment analysis.** default (A), modify metrics (B),  modify color (C) and modify term length (D).

```{r barplot, fig.width=20, fig.height=10, fig.cap="(ref:barplotCap)", fig.scap="(ref:barplotScap)", echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 30))
```

## Bubble Plot

Enriched gene sets are shown as bubbles with different sizes.

The x-axis is statistic value and the y-axis is "Fold Enrichment".

The basic arguments are: 

- `stats_metric`: Statistic value of "p.adjust", "pvalue" or "qvalue"
- `scale_ratio`: Change bubble size. Default is 1.


```{r include= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "bubble", main_text_size = 15, legend_text_size = 10)
p2 <- plotEnrich(ego,
  plot_type = "bubble", scale_ratio = 0.5, stats_metric = "qvalue",
  main_text_size = 15, legend_text_size = 10
)
pp <- p1 / p2
```

```{r eval= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "bubble")
p2 <- plotEnrich(ego, plot_type = "bubble", 
                 scale_ratio = 0.5, stats_metric = "qvalue")
p1 / p2 + plot_annotation(tag_levels = "A")
```

(ref:bubbleplotScap) Bubble plot of enrichment analysis.

(ref:bubbleplotCap) **Bubble plot of enrichment analysis.** default (A), modify bubble size (B).

```{r bubbleplot, fig.height=12, fig.width=10, fig.cap="(ref:bubbleplotCap)", fig.scap="(ref:bubbleplotScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 30))
```

## Dot Plot

Similar with bar plot, dot plot is also widely used in enrichment analysis plotting. Like bubble plot, dot size represents gene number of enriched term.

```{r include= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "dot", main_text_size = 15, legend_text_size = 10)
p2 <- plotEnrich(ego,
  plot_type = "dot", scale_ratio = 1.5, stats_metric = "pvalue",
  main_text_size = 15, legend_text_size = 10,
  term_metric = "RichFactor"
)
pp <- p1 + p2
```

```{r eval= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "dot")
p2 <- plotEnrich(ego,
  plot_type = "dot",
  scale_ratio = 1.5,
  stats_metric = "pvalue",
  term_metric = "RichFactor"
)
p1 + p2 + plot_annotation(tag_levels = "A")
```

(ref:dptplotScap) Dot plot of enrichment analysis.

(ref:dotplotCap) **Dot plot of enrichment analysis.** default (A), modify bubble size (B).

```{r dotplot, fig.height=10, fig.width=20, fig.cap="(ref:dotplotCap)", fig.scap="(ref:dotplotScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 30))
```

After analyzing group enrichment analysis for [ORA](#ora-group-enrichment), we could use `dotplot` to show the result:

```{r include=FALSE}
id_100 <- c(head(names(geneList), 50), tail(names(geneList), 50))
two_groups <- list(
  exp_group = c(rep("up", 50), rep("down", 50)),
  time_group = c(rep("time1", 40), rep("time2", 60))
)
gs <- geneset::getGO(org = "human",ont = "mf")
gego2 <- genORA(id_100, geneset = gs, group_list = two_groups)
```

(ref:dptplotGroupScap) Dot plot of group enrichment analysis.

(ref:dotplotGroupCap) **Dot plot of group enrichment analysis.** number in round brackets shows total gene number in selected pathways

```{r dotplotGroup, fig.height=6, fig.width=10, fig.cap="(ref:dotplotGroupCap)", fig.scap="(ref:dotplotGroupScap)",echo = T, dpi=300}
plotEnrich(gego2,
           plot_type = 'dot',
           scale_ratio = 2, # dot size
           main_text_size = 10,
           legend_text_size =8,
           n_term = 6) # show terms
```



## Lollipop Plot

Lollipop is like combination of barplot and dotplot

```{r include= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "lollipop", main_text_size = 15, legend_text_size = 10)
p2 <- plotEnrich(ego,
  plot_type = "lollipop", scale_ratio = 1.5, stats_metric = "pvalue",
  main_text_size = 15, legend_text_size = 10,
  term_metric = "RichFactor",
  up_color = "#a32a31",
  down_color = "#f7dcca"
)
pp <- p1 + p2
```

```{r eval= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "lollipop")
p2 <- plotEnrich(ego,
  plot_type = "lollipop",
  scale_ratio = 1.2,
  stats_metric = "pvalue",
  term_metric = "RichFactor",
  up_color = "#a32a31",
  down_color = "#f7dcca"
)
p1 + p2 + plot_annotation(tag_levels = "A")
```

(ref:lolliplotScap) Lollipop plot of enrichment analysis.

(ref:lolliplotCap) **Lollipop plot of enrichment analysis.** defult (A), selected genes (B)

```{r lolliplot, fig.height=10, fig.width=20, fig.cap="(ref:lolliplotCap)", fig.scap="(ref:lolliplotScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 30))
```

## Heatmap Plot

Heatmap plot shows interactions between enriched terms and their genes.
If fold change is given, heatmap will add color for up and down-regulated genes.


```{r include= F}
library(patchwork)

logfc = geneList
p1 <- plotEnrich(ego, plot_type = "geneheat")
show_gene = c('JUN','SOX2','CD24','TLR4')
p2 <- plotEnrich(ego, plot_type = "geneheat", show_gene = show_gene)
p3 <- plotEnrich(ego, plot_type = "geneheat", show_gene = show_gene, fold_change = logfc)
pp <- p1 / p2 / p3
```

```{r eval= F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "geneheat")
show_gene = c('JUN','SOX2','CD24','TLR4')
p2 <- plotEnrich(ego, plot_type = "geneheat", show_gene = show_gene)
p3 <- plotEnrich(ego, plot_type = "geneheat", show_gene = show_gene, fold_change = logfc)
p1 / p2 / p3 + plot_annotation(tag_levels = "A")
```

(ref:heatmapScap) Heatmap plot of enrichment analysis.

(ref:heatmapCap) **Heatmap plot of enrichment analysis.** all genes (A, default), selected genes (B), selected genes with logFC value (C).

```{r heatmap, fig.height=8, fig.width=12, fig.cap="(ref:heatmapCap)", fig.scap="(ref:heatmapScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 20))
```

## Chord Plot

Inspired by [GOplot](https://wencke.github.io/), chord plot is reproduced using `r CRANpkg("ggplot2")` and it shows similar content with heatmap plot.

- `gene_space`: The space between the gene labels and the chord.

```{r include= F}
logfc = geneList
library(patchwork)
show_gene = c('JUN','SOX2','CD24','TLR4')
p1 <- plotEnrich(ego, plot_type = "genechord", show_gene = show_gene) +
  ggplot2::theme(legend.position = "none")
p2 <- plotEnrich(ego,
  plot_type = "genechord", show_gene = show_gene, fold_change = logfc,
  remove_legend_text = T,gene_space = 0.5
)
pp <- p1 + p2
```

```{r eval= F}
library(patchwork)
show_gene = c('JUN','SOX2','CD24','TLR4')
p1 <- plotEnrich(ego, plot_type = "genechord", 
                 show_gene = show_gene) +
  ggplot2::theme(legend.position = "none")
p2 <- plotEnrich(ego, plot_type = "genechord", 
                 show_gene = show_gene, 
                 fold_change = logfc,
                 gene_space = 0.5)
p1 + p2 + plot_annotation(tag_levels = "A")
```

(ref:chordScap) Chord plot of enrichment analysis.

(ref:chordCap) **Heatmap plot of enrichment analysis.** selected genes (A), selected genes with logFC value (B).

```{r chord, fig.height=8, fig.width=15, fig.cap="(ref:chordCap)", fig.scap="(ref:chordScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 20))
```

## Wordcloud Plot

Wordcloud plot shows term emphasis based on text frequency.

(ref:wordcloudScap) Wordcloud plot of enrichment analysis.

(ref:wordcloudCap) **Wordcloud plot of enrichment analysis.** 

```{r wordcloud, fig.height=5, fig.width=8, fig.cap="(ref:wordcloudCap)", fig.scap="(ref:wordcloudScap)",results='hide'}
plotEnrich(ego, plot_type = "wordcloud")
```

## Upset Plot

Inspired by `r CRANpkg("ComplexUpset")`, the upset plot shows the association between genes and enriched terms. Unlike common Venn diagram, it could emphasize complex relationship among many gene sets.

(ref:upsetScap) Upset plot of enrichment analysis.

(ref:upsetCap) **Upset plot of enrichment analysis.** 

```{r upset, fig.height=5, fig.width=10, fig.cap="(ref:upsetCap)", fig.scap="(ref:upsetScap)"}
plotEnrich(ego, plot_type = "upset",main_text_size = 15,legend_text_size = 8)
```


## Network Plot

Inspired by `enrichplot::emapplot`, `genekitr` reproduced this plot to enhance modification.
Enriched terms are the nodes in the network, the overlapping gene sets are edges. Terms with more overlapped genes are tend to cluster together and the edge will be thicker.

- GO analysis result could use five methods to calculate the similarity between nodes: "Resnik", "Lin", "Rel", "Jiang", "Wang" and "JC"(Jaccard’s similarity index) methods.
- KEGG only supports "JC" method.

User could define `layout` argument derived from `r CRANpkg("ggraph")`, including "nicely" (default), "circle", "dh", "drl", "fr", "graphopt", "grid","lgl", "kk", "mds", "randomly", "star" etc.

> For more information about the `layout`, you could refer to: ["Introduction to ggraph: Layouts"](https://www.data-imaginist.com/2017/ggraph-introduction-layouts/)

```{r include= F}
library(patchwork)
library(igraph)
library(ggraph)
p1 <- plotEnrich(ego, plot_type = "network", scale_ratio = 0.5)
p2 <- plotEnrich(ego, plot_type = "network", layout = "circle", scale_ratio = 0.5)
p3 <- plotEnrich(ego, plot_type = "network", 
                 layout = "grid", sim_method = "Wang",
                 up_color = "#a32a31", down_color = "#f7dcca")
pp <- (p1 + p2) / p3
```

```{r eval= F}
library(patchwork)
library(igraph)
library(ggraph)
p1 <- plotEnrich(ego, plot_type = "network", scale_ratio = 0.5)
p2 <- plotEnrich(ego, plot_type = "network", 
                 layout = "circle", scale_ratio = 0.5)
p3 <- plotEnrich(ego, plot_type = "network", 
                 layout = "grid", sim_method = "Wang",
                 up_color = "#a32a31", down_color = "#f7dcca")
(p1 + p2) / p3 + plot_annotation(tag_levels = "A")
```

(ref:networkScap) Network plot of enrichment analysis.

(ref:networkCap) **Network plot of enrichment analysis.** JC method and nicely layout(A, default), circle layout (B), grid layout and Want method (C).

```{r network, fig.height=10, fig.width=15, fig.cap="(ref:networkCap)", fig.scap="(ref:networkScap)",echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 20))
```


## GO-specific: WEGO Plot

> To visulize more than one ontology of GO in one plot, user could utilize `wego plot`

Inspired by [WEGO](https://wego.genomics.cn/), `genekitr` utilized `r CRANpkg("ggplot2")` to reproduce this plot.

Here we generate two ontologies (MF and CC) result. 

```{r}
# 1st step: prepare input IDs
data(geneList, package = "genekitr")
id <- names(geneList)[abs(geneList) > 2]

# 2nd step: prepare CC and MF gene sets
go_cc <- geneset::getGO(org = "human",ont = "cc")
go_mf <- geneset::getGO(org = "human",ont = "mf")

# 3rd step: analysis
ego_cc <- genORA(id, geneset = go_cc)
ego_mf <- genORA(id, geneset = go_mf)

# 4th step: merge two data frames
# Note: each data frame should add new column "Ontology"
ego_cc <- ego_cc %>% dplyr::mutate(Ontology = "cc") %>% dplyr::rename(ID = 1)
ego_mf <- ego_mf %>% dplyr::mutate(Ontology = "mf") %>% dplyr::rename(ID = 1)
all_ego <- rbind(ego_cc,ego_mf)
```

**NOTICE:** `plotEnrich` has a parameter `n_term` for WEGO plot, which specify number of terms. If you want to plot all terms, just set `n_term` higher.

(ref:wegoplotScap) WEGO plot of enrichment analysis.

(ref:wegoplotCap) **WEGO plot of enrichment analysis.** 

```{r wegoplot, fig.height=3, fig.width=8, fig.cap="(ref:wegoplotCap)", fig.scap="(ref:wegoplotScap)"}
plotEnrich(all_ego, plot_type = "wego", n_term = 5)
```

## GO-specific: Map Plot

GO terms are built in a directed acyclic graph with a parent-child relationship.
Here the map plot utilized `r Biocpkg("GOSemSim")` to extract parent and child terms, also utilized `r CRANpkg("ggraph")` and `r CRANpkg("igraph")` to draw with default layout "sugiyama".

To avoid too much unrelated terms messed up, `genekitr` **only shows the closest parent and child with selected terms**. Besides, the top three parent term with more edges will be plotted.

```{r include= F}
library(igraph)
library(ggraph)
pp <- plotEnrich(ego, plot_type = "gomap", wrap_length = 25,
                 up_color = '#a32a31',down_color = '#3665a6')
```

```{r eval= F}
library(igraph)
library(ggraph)
plotEnrich(ego, plot_type = "gomap", wrap_length = 25,
           up_color = '#a32a31',down_color = '#3665a6')
```

(ref:mapScap) Map plot of enrichment analysis.

(ref:mapCap) **Map plot of enrichment analysis.** 

```{r map, fig.height=5, fig.width=5, fig.cap="(ref:mapCap)", fig.scap="(ref:mapScap)",echo = F, dpi=300}
pp
```

## GO-specific: Terms Heatmap Plot

Inspired by `r Biocpkg("rrvgo")`, `genekitr` scratches main codes to cluster GO terms. Also it could use five methods to calculate the similarity between terms.

```{r eval= F}
plotEnrich(ego, plot_type = "goheat", sim_method = "Rel")
```

(ref:goheatScap) GO heatmap plot of enrichment analysis.

(ref:goheatCap) **GO heatmap plot of enrichment analysis.** 

```{r goheat, fig.height=5, fig.width=10, fig.cap="(ref:goheatCap)", fig.scap="(ref:goheatScap)",echo = F, dpi=300}
plotEnrich(ego,
  plot_type = "goheat", sim_method = "Rel",
  main_text_size = 10
)
```


## GO-specific: Terms Tangram Plot

According to [rrvgo vignettes](http://bioconductor.org/packages/release/bioc/vignettes/rrvgo/inst/doc/rrvgo.html), tangram plot is space-filling visualization of hierarchical structures. The terms are grouped and colored based on their parent, and the space used by the term is proportional to the score.  

(ref:gotangramScap) GO tangram plot of enrichment analysis.

(ref:gotangramCap) **GO tangram plot of enrichment analysis.** 

```{r gotangram, fig.height=5, fig.width=8, fig.cap="(ref:gotangramCap)", fig.scap="(ref:gotangramScap)",results='hide'}
plotEnrich(ego, plot_type = "gotangram", sim_method = "Rel")
```



## Plot Theme

`Genekitr` provides a function `plot_theme()` to modify all plots theme including text size, border, legend, color, etc.

```{r include= F}
library(patchwork)
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "dot")
p2 <- plotEnrich(ego, plot_type = "dot", main_text_size = 13, legend_text_size = 13)
p3 <- plotEnrich(ego, plot_type = "dot", border_thick = 3, remove_grid = F)
p4 <- plotEnrich(ego,
  plot_type = "dot",
  remove_main_text = T,
  remove_legend_text = T,
  remove_legend = T
)
pp <- p1 + p2 + p3 + p4
```


```{r eval = F}
library(patchwork)
p1 <- plotEnrich(ego, plot_type = "dot")
p2 <- plotEnrich(ego,
  plot_type = "dot",
  main_text_size = 10,
  legend_text_size = 10
)

p3 <- plotEnrich(ego,
  plot_type = "dot",
  border_thick = 3,
  remove_grid = F
)

p4 <- plotEnrich(ego,
  plot_type = "dot",
  remove_main_text = T,
  remove_legend_text = T,
  remove_legend = T
)
p1 + p2 + p3 + p4 + plot_annotation(tag_levels = "A")
```

(ref:themeScap) Plot theme.

(ref:themeCap) **Plot theme.** default theme (A), modify text size (B),  modify grid line and border size (C) and remove all text and legend (D).


```{r theme, fig.width=15, fig.height=10, fig.cap="(ref:themeCap)", fig.scap="(ref:themeScap)", echo = F, dpi=300}
pp + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 30))
```

## Advanced Plot

### Two-group barplot for up/down regulated pathways

Up and down regulated genes could pass to `genORA` separately and visualize together as two-group barplot.

Here we take GO result as an example:

```{r}
# 1st step: prepare input IDs
# Since the geneList is logFC decreasing ordered, we could take first 100 as up-regulated genes and vice versa.
data(geneList, package = "genekitr")
up_genes <- head(names(geneList), 100)
down_genes <- tail(names(geneList), 100)

# 2nd step: prepare gene set
hg_gs <- geneset::getGO(org = "human",ont = "bp")

# 3rd step: ORA analysis separately
up_go <- genORA(up_genes, geneset = hg_gs)
down_go <- genORA(down_genes, geneset = hg_gs)

dim(up_go)
dim(down_go)

# here we only take 10 terms of each
up_go <- head(up_go, 10)
down_go <- head(down_go, 10)
```

There are two visulization types: 

- `plot_type = "one"`: both up and down-regulated pathways are plotted together
- `plot_type = "two"`: up and down-regulated pathways are plotted separately. **Recommended if the number of both groups is similar.**

(ref:twogroup1Scap) Visualize two groups together when `plot_type = "one"`

(ref:twogroup1Cap) **Visualize two groups together when `plot_type = "one"`.** 

```{r twogroup1, fig.width=7, fig.height=5, fig.cap="(ref:twogroup1Cap)", fig.scap="(ref:twogroup1Scap)", dpi=300}
plotEnrichAdv(up_go, down_go,
              plot_type = "one",
              term_metric = "FoldEnrich",
              stats_metric = "p.adjust",
              xlim_left = 25, xlim_right = 15) +
  theme(legend.position = c(0.15, 0.9))
```


```{r include= F}
library(patchwork)
color <- c("#a32a31", "#f7dcca", "#3665a6", "#d5e4ef")
left <- suppressMessages(plotEnrich(up_go,
  plot_type = "bar",
  term_metric = "FoldEnrich",
  stats_metric = "qvalue",
  up_color = color[1], down_color = color[2],
  main_text_size = 15,
  legend_text_size = 10
) +
  scale_y_discrete(limits = rev) +
  scale_x_reverse() +
  theme(
    axis.title.y = element_blank(),
    legend.position = c(0.2, 0.8)
  ))

right <- suppressMessages(plotEnrich(down_go,
  plot_type = "bar",
  term_metric = "FoldEnrich",
  stats_metric = "qvalue",
  up_color = color[3], down_color = color[4],
  main_text_size = 15,
  legend_text_size = 10
) +
  scale_y_discrete(position = "right") +
  theme(
    axis.title.y = element_blank(),
    legend.position = c(0.8, 0.2)
  ))

pp <- left + right
```


```{r eval = F}
plotEnrichAdv(up_go, down_go,
              plot_type = "two",
              term_metric = "FoldEnrich",
              stats_metric = "qvalue"
)
```

(ref:twogroup2Scap) Visualize two groups seperately when `plot_type = "two"`.

(ref:twogroup2Cap) **Visualize two groups seperately when `plot_type = "two"`.** 

```{r twogroup2, fig.width=15, fig.height=5, fig.cap="(ref:twogroup2Cap)", fig.scap="(ref:twogroup2Scap)", dpi=300, echo = F}
pp
```