---
title: "hidecan"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{hidecan}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(tibble)
library(dplyr)
library(purrr)
library(stringr)
library(viridis)
```
```{r setup}
library(hidecan)
```
## Input data
The `hidecan` package takes as input tibbles (data-frames) of GWAS and DE results or candidate genes. The input data-frames should contain some mandatory columns, depending on the type of data.
A list of example input datasets can be obtained via the `get_example_data()` function:
```{r get-example-data}
x <- get_example_data()
str(x, max.level = 1)
```
### GWAS results
GWAS results should be provided as a tibble or data-frame, with one row per genetic marker. The data-frame should contain at least the following columns:
- `chromosome`: character column, giving the ID of the chromosome on which each marker is located;
- `position`: numeric column, the physical position along the chromosome (in base pairs - bp) of the marker;
- either `score` or `padj`: numeric column, providing either the score (i.e. `-log10(p-value)`) or the adjusted p-value of the marker. If a `score` column is provided, the `padj` column will be ignored. If only a `padj` column is provided, a `score` column will be constructed as `-log10(padj)`.
Any other column present in the data-frame will be ignored. An example of valid input is shown below:
```{r example-data-gwas}
head(x[["GWAS"]])
```
### Differential expression results
DE results should be provided as a tibble or data-frame, with one row per gene. The data-frame should contain at least the following columns:
- `chromosome`: character column, giving the ID of the chromosome on which each gene is located;
- `start` and `end`: numeric columns, giving the starting and end position of the gene along the chromosome (in bp). These two columns will be used to calculate the position of the gene as the half-way point between the start and end of the gene.
- either `score` or `padj`: numeric column, providing either the score (i.e. `-log10(p-value)`) or the adjusted p-value of the gene. If a `score` column is provided, the `padj` column will be ignored. If only a `padj` column is provided, a `score` column will be constructed as `-log10(padj)`.
- either `foldChange` or `log2FoldChange`: numeric column, giving either the fold-change or log2(fold-change) of the gene. If a `log2FoldChange` column is provided, the `foldChange` column will be ignored. If only a `foldChange` column is provided, a `log2FoldChange` will be constructed as `log2(foldChange)`.
Any other column present in the data-frame will be ignored. An example of valid input is shown below:
```{r example-data-de}
head(x[["DE"]])
```
(Note that in the example dataset, some genes have missing values in the `padj` column; this corresponds to genes that have been filtered out via [independent filtering in the DESeq2 package](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pvaluesNA)).
### Candidate genes
A list of candidate genes (e.g. genes previously found associated with a trait of interest based on literature search) can be provided as a tibble or data-frame, with one row per gene. This data-frame can also contain variants of interest (see below). The data-frame should contain at least the following columns:
- `chromosome`: character column, giving the ID of the chromosome on which each gene is located;
- `start` and `end`: numeric columns, giving the starting and end position of the gene along the chromosome (in bp). These two columns will be used to calculate the position of the gene as the half-way point between the start and end of the gene. For genomic variants or markers, simply set both the `start` and `end` columns to the physical position of the marker.
- `name`: character column, giving the name of the candidate gene that will be displayed in the HIDECAN plot. Set to `NA` for any subset of genes to remove their label in the plot (can help if many genes are very close, to avoid cluttering the plot).
Any other column present in the data-frame will be ignored. An example of valid input is shown below:
```{r example-data-can}
head(x[["CAN"]])
```
## Creating a HIDECAN plot
The `hidecan_plot()` function creates a HIDECAN plot. It takes as input the data-frames presented above, as well as the score and log2(fold-change) thresholds used to select the significant markers and genes.
In this example, we only show markers with a score above 4, which corresponds to a p-value below $1\times10^{-4}$, and genes with a score above 1.3, which corresponds to a p-value of 0.05. We don't place any threshold on the log2(fold-change) of the genes:
```{r hidecan-plot, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]], ## data-frame of GWAS results
de_list = x[["DE"]], ## data-frame of DE results
can_list = x[["CAN"]], ## data-frame of candidate genes
score_thr_gwas = -log10(0.0001), ## sign. threshold for GWAS
score_thr_de = -log10(0.05), ## sign. threshold for DE
log2fc_thr = 0 ## log2FC threshold for DE
)
```
Note that it is possible to provide only a subset of the possible input data, e.g. only GWAS results and a list of candidate genes:
```{r hidecan-plot-gwas-can-only, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
can_list = x[["CAN"]],
score_thr_gwas = 4
)
```
## Removing empty chromosomes
By default, the HIDECAN plot shows all chromosomes present in the input data. However, it is possible that some of the chromosomes appear empty, as they do not contain any significant gene or marker, nor any candidate gene. In this case, it is possible to exclude such "empty" chromosomes from the HIDECAN plot, through the `remove_empty_chrom` argument.
We will demonstrate that by increasing the score threshold applied to the GWAS results, in order to get fewer significant markers. In this case, chromosomes 0, 6, 9 and 10 do not contain any significant marker of gene of interest:
```{r hidecan-plot-with-empty-chrom, fig.width = 10, fig.height = 10}
## Chromosomes 0, 6, 9 and 10 are empty
hidecan_plot(
gwas_list = x[["GWAS"]],
can_list = x[["CAN"]],
score_thr_gwas = 5
)
```
By setting the `remove_empty_chrom` argument to `TRUE`, these chromosomes will be removed from the plot:
```{r hidecan-plot-without-empty-chrom, fig.width = 10, fig.height = 8}
hidecan_plot(
gwas_list = x[["GWAS"]],
can_list = x[["CAN"]],
score_thr_gwas = 5,
remove_empty_chrom = TRUE
)
```
## Selecting chromosomes and genomic positions
It is possible to specify which chromosomes should be represented in the HIDECAN plot, via the `chroms` argument. For example, with the following command we restrict the plot to chromosomes 7 and 8:
```{r hidecan-select-chroms, fig.width = 10, fig.height = 3}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
chroms = c("ST4.03ch07", "ST4.03ch08")
)
```
We can also "zoom in" on some or all chromosomes, through the `chrom_limits` argument. To zoom in on all chromosomes at once, we pass to the `chrom_limits` argument an integer vector of length 2, which gives the lower and upper limits in bp to use. For example here, we focus on the 10-20Mb region of each chromosome:
```{r hidecan-chrom-limits-all-chrom, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
chrom_limits = c(10e6, 20e6)
)
```
Alternatively, we can apply different limits to some of the chromosomes, by passing a named list to the argument. The names of the list should match the chromosomes name, and each element should be an integer vector of length 2 giving the lower and upper limits in bp to use for the corresponding chromosome. For example, we will focus on the 10-20Mb region for chromosome 1, and the 30-40Mb region for chromosome 5, and leave all other chromosomes as is:
```{r hidecan-chrom-limits-some-chrom, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
chrom_limits = list("ST4.03ch01" = c(10e6, 20e6),
"ST4.03ch05" = c(30e6, 40e6))
)
```
The two options `chroms` and `chrom_limits` can be used together:
```{r hidecan-select-chroms-and-chrom-lims, fig.width = 10, fig.height = 3}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
chroms = c("ST4.03ch07", "ST4.03ch08"),
chrom_limits = list("ST4.03ch07" = c(50e6, 55e6),
"ST4.03ch08" = c(45e6, 50e6))
)
```
## Colour genes by log2(fold-change)
By default, in a HIDECAN plot, the points representing both significant markers and DE genes are coloured according to their GWAS/DE score. However, it is possible to colour the DE genes by their log2(fold-change) value instead, by setting the `colour_genes_by_score` argument to `FALSE`:
```{r hidecan-genes-colour-log2fc, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
colour_genes_by_score = FALSE
)
```
Genes with a negative log2(fold-change) will be represented with a shade of blue, and genes with a positive log2(fold-change) will be represented with a shade of red.
## More than one GWAS, DE or candidate gene list
The `hidecan_plot()` function can take as an input lists of data-frames for GWAS results, DE results or candidate genes. This way, it is possible to visualise more than one GWAS or DE analyses at once, for example if investigating several traits at once or comparing more than two treatment groups.
For this example, we'll focus on chromosomes 7 and 8 (only for clarity of the plot):
```{r making-small-dataset}
library(dplyr)
library(purrr)
library(stringr)
## Retaining only markers and genes on chromosomes 7 and 8
x_small <- x |>
map(~ filter(.x, str_detect(chromosome, "(07|08)")))
```
We'll create a second data-frame of GWAS results by shuffling the marker scores in the example dataset:
```{r creating-second-gwas-tibble}
## Creating a second GWAS result tibble by shuffling
## the marker scores from the original data
gwas_1 <- x_small[["GWAS"]]
gwas_2 <- gwas_1 |>
mutate(score = sample(score))
```
We can pass both GWAS results data-frames to the `hidecan_plot()` function as a list:
```{r hidecan-multiple-gwas-input, fig.width = 10, fig.height = 3}
hidecan_plot(
gwas_list = list(gwas_1, gwas_2),
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0
)
```
By default, the two GWAS tracks will be given unique y-axis labels, as can be seen above. It is possible to customise this by naming the elements in the input list:
```{r multiple-gwas-input-named, fig.width = 10, fig.height = 3}
hidecan_plot(
gwas_list = list("Trait 1" = gwas_1,
"Trait 2" = gwas_2),
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0
)
```
## Defining chromosomes length
By default, the `hidecan_plot()` function calculates the length of the different chromosomes based on the input data, by looking at the maximum position of genes and markers on each chromosome. However, it is also possible to pass on a tibble of chromosome length (in bp) through the `chrom_length` argument.
```{r make-potato-chrom-length}
library(tibble)
## Chromosomes length as recorded in Ensembl Plants
potato_chrom_length <- c(
ST4.03ch00 = 45813526,
ST4.03ch01 = 88663952,
ST4.03ch02 = 48614681,
ST4.03ch03 = 62190286,
ST4.03ch04 = 72208621,
ST4.03ch05 = 52070158,
ST4.03ch06 = 59532096,
ST4.03ch07 = 56760843,
ST4.03ch08 = 56938457,
ST4.03ch09 = 61540751,
ST4.03ch10 = 59756223,
ST4.03ch11 = 45475667,
ST4.03ch12 = 61165649
) |>
## turn a named vector into a tibble
enframe(name = "chromosome",
value = "length")
head(potato_chrom_length)
```
```{r hidecan-chrom-length, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
chrom_length = potato_chrom_length
)
```
Note that in this case we can't really see the difference with the computed chromosome length values.
## Controlling the plot properties
The `hidecan_plot()` function offers several arguments to control different aspects of the HIDECAN plot. For example, it is possible to specify the number of rows or columns the plot should have, through the `n_rows` and `n_cols` arguments. Note that only one of these arguments will be considered (`n_rows` takes precedence):
```{r hidecan-nrows, fig.width = 10, fig.height = 8}
## Specifying the number of rows
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
n_rows = 3
)
```
```{r hidecan-ncols, fig.width = 10, fig.height = 10}
## Specifying the number of columns
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.005),
log2fc_thr = 0,
n_cols = 3
)
```
In addition, it is possible to:
- add a title and subtitle to the plot (`title` and `subtitle` arguments);
- control the position of the legend (`legend_position` argument);
- control the size of the points (`point_size` argument);
- control the size of the labels and padding around the text for the candidate genes labels (`label_size` and `label_padding` arguments).
## Creating a Manhattan plot
It is possible to visualise GWAS results with a Manhattan plot through the `manhattan_plot()` function. The function takes as input a data-frame (or list of data-frames) of GWAS results, which must follow the same rules as for `plot_hidecan()` function (see [this section](hidecan.html#gwas-results) of the vignette). A horizontal line representing the significance threshold can be added to the plot, by passing the threshold value to the `score_thr` argument.
```{r manhattan-plot, fig.width = 11, fig.height = 6}
manhattan_plot(x[["GWAS"]], score_thr = 4)
```
As for the `hidecan_plot()` function, it is possible to pass a list of GWAS results to the function, in order to draw several Manhattan plots in one figure. The names of the list will be used to identify the different plots:
```{r manhattan-plot-facets, fig.width = 11, fig.height = 12}
manhattan_plot(list("Trait A" = x[["GWAS"]],
"Trait B" = x[["GWAS"]]),
score_thr = 4,
ncol = 1)
```
## Controlling the aesthetics (colours, point shapes, labels) of the tracks
The default look of the different data types (GWAS peaks, DE genes and candidate genes) in a HIDECAN plots is controlled by the arguments returned by `hidecan_aes()`. More specifically, the function controls:
- the label of each data type on the x-axis (e.g. 'GWAS peaks' for GWAS results),
- the colour of the vertical lines used to showcase the position of markers or genes,
- the shape of the points,
- whether the name of the markers/genes are displayed (by default, will only show for candidate genes),
- the colour palette used for the score or log2 fold-change of the points.
The default values used for these parameters have been chosen to ensure that the resulting plot is easily understandable. However, it is possible to change them. To do so, we'll start by saving the defaults returned by `hidecan_aes()`:
```{r hidecan-aes}
default_aes <- hidecan_aes()
str(default_aes, max.level = 1)
default_aes[[1]]
```
As we can see above, the function returns a list of the default aesthetics to use for GWAS results (the `GWAS_data_thr` element), for DE results (`DE_data_thr` element), for candidate genes (`CAN_data_thr` element), and for custom tracks (`CUSTOM_data_thr` element - see the [Adding custom tracks](web_only/custom_tracks.html) article for more information). Each element is itself a list, containing the value for each attribute mentioned above.
Let's modify the label given to the GWAS peaks, and the colour palette used for DE score:
```{r modify-hidecan-aes}
library(viridis)
default_aes$GWAS_data_thr$y_label <- "GWAS signif. markers"
default_aes$DE_data_thr$fill_scale <- scale_fill_viridis(
"DE gene score",
option = "inferno",
guide = ggplot2::guide_colourbar(
title.position = "top",
title.hjust = 0.5,
order = 3 # ensures this legend is shown after the one for GWAS scores
)
)
```
We can use these modified values by passing this list to `hidecan_plot()` through its `custom_aes` argument:
```{r use-custom-aes, fig.width = 10, fig.height = 10}
hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
custom_aes = default_aes
)
```
## Viewport error
If you are working on RStudio, you may encounter the following error:
```{r show-window-size-error, echo = FALSE, error = TRUE, fig.width = 1, fig.height = 1}
hidecan_plot(gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
label_size = 2)
```
This is caused by the plotting window being too small. Try increasing the size of the plotting window in the RStudio console. Alternatively, you can save the plot into an R object, then use `ggplot2::ggsave()` to save it into a file:
```{r use-ggsave, eval = FALSE}
p <- hidecan_plot(
gwas_list = x[["GWAS"]],
de_list = x[["DE"]],
can_list = x[["CAN"]],
score_thr_gwas = -log10(0.0001),
score_thr_de = -log10(0.05),
log2fc_thr = 0,
label_size = 2
)
ggplot2::ggsave("hidecan_plot.pdf", p, width = 10, height = 10)
```