```{r setup, include=FALSE}

# knitr setup options
knitr::opts_chunk$set(
  warning=FALSE, 
  message=FALSE, 
  echo = TRUE)

# load all the packages
library(tidyverse)
library(GGally)
library(ggExtra)
library(ggthemes)
# library(hrbrthemes)
library(ggpubr)
library(pheatmap)

```


```{r data}
library(tidyverse)
gapminder2011 <- read_csv("../data/Gapminder_vars_2011.csv")
```


# `ggmarginal`

https://cran.r-project.org/web/packages/ggExtra/vignettes/ggExtra.html

```{r margins_FoodvsLifeExp, fig.width=10, fig.height=5}

library(ggExtra)

p <- ggplot(data = gapminder2011,
        aes(x = FoodSupplykcPPD, 
            y = LifeExpectancyYrs,
            color = four_regions)
        ) +
  geom_point(alpha = .4) +
  scale_color_discrete(
    name = "Regions",
    labels = c("Africa", "Americas", 
               "Asia", "Europe")
    ) +
  theme(legend.position="bottom") +
  labs(
    x = "Daily Food Supply PP (kc)",   
    y = "Life Expectancy (years)",   
    title = "Scatterplot"
    )
```



```{r margins_FoodvsLifeExp_out, fig.width=10, fig.height=5}
ggMarginal(p,
  type = "density",
  margins = "both",
  groupColour = TRUE,
  groupFill = TRUE
)
```


## `GGally::ggpairs()`
https://ggobi.github.io/ggally/index.html

```{r fig.width=10, fig.height=5}
library(GGally)
gapminder2011 %>% 
  select(FoodSupplykcPPD:WaterSourcePrct) %>% # specifying which columns to use
  ggpairs()
```


# Genomic data

```{r, echo=FALSE}
#Install ggbio 
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ggbio")
```

```{r}
browseVignettes("ggbio")
```


```{r}
library("ggbio")
data("hg19IdeogramCyto", package = "biovizBase")
Ideogram(hg19IdeogramCyto, subchr = "chr1")

```

```{r}
library("GenomicRanges")
data("darned_hg19_subset500", package = "biovizBase")
autoplot(darned_hg19_subset500, layout = "karyogram",
         aes(color = exReg, fill = exReg))
```

```{r}
data("ideoCyto", package = "biovizBase")
dn = darned_hg19_subset500
seqlengths(dn) = seqlengths(ideoCyto$hg19)[names(seqlengths(dn))]
dn = keepSeqlevels(dn, paste0("chr", c(1:22, "X")))
autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg))
```



# Gene Expression



# Pasilla Data

```{r}
pasilla_data <- read_csv("../data/gene_expr_pasilla_results.csv")
glimpse(pasilla_data)
```


# Volcano Plot

```{r volcanoplot_nice, include = TRUE}
# Create subset for labeling
pasilla_data_top = pasilla_data %>%
  filter(-log10(padj) > 10, 
         abs(log2FoldChange) > 2.5)

ggplot(data = pasilla_data,
       aes(x = log2FoldChange,
           y = log10(padj))) +
  geom_point() + 
  scale_y_reverse() +
  aes(color = padj < 0.05) +
  ggrepel::geom_text_repel(
    data = pasilla_data_top, 
    aes(label = gene), color = "black",
    box.padding = 0.5, min.segment.length = 0) +
  xlim(c(-7,7)) + 
  geom_vline(xintercept = c(-2.5, 2.5), 
             lty = "dashed", color="grey") + 
  ggthemes::theme_clean() + 
  labs(
    x = bquote(~Log[2]~ "fold change"),
    y = bquote(~Log[10]~adjusted~italic(P)),
    title = "Volcano Plot",
    subtitle = "Gene Expression of Pasilla Data"
  )
```


```{r volcanoplot, include=FALSE}
ggplot(data = pasilla_data,
       aes(x = log2FoldChange,
           y = log10(padj))) +
  geom_point() + 
  scale_y_reverse() +
  aes(color = padj < 0.05) +
  ggrepel::geom_text_repel(
    data = pasilla_data_top, 
    aes(label = gene), color = "black",
    box.padding = 0.5,
    min.segment.length = 0) +
  xlim(c(-7,7)) + 
  geom_vline(xintercept = c(-2.5, 2.5), 
             lty = "dashed", color="grey") + 
  ggthemes::theme_clean() + 
  labs(
    x = bquote(~Log[2]~ "fold change"),
    y = bquote(~Log[10]~adjusted~italic(P)),
    title = "Volcano Plot",
  subtitle="Gene Expression of Pasilla Data"
  )
```




# Heatmap with `pheatmap::pheatmap()`

It's possible to make heatmaps in ggplot2 with `geom_tile()`, but there are many other better functions using base R that cluster and annotate the data. This is using `pheatmap` package.


```{r pasilla_heat, cache=FALSE}
# select expression data
pasilla_heat <- pasilla_data %>%
  select(treated1:untreated4)
# subtract off gene-specific means
pasilla_heat <- pasilla_heat - rowMeans(pasilla_heat)
# calculate standard deviation of each centered gened
sd_gene <- apply(pasilla_heat,1,sd)
# select top 500 most variable
pasilla_heat <- 
  pasilla_heat[order(sd_gene, decreasing = TRUE)[1:500],]

# create annotation data
pasilla_col <- data.frame(
  trt = factor(c(rep("trt",3), rep("untrt",4))),
  id = 1:7,
  row.names=colnames(pasilla_heat))
```


```{r}
head(pasilla_heat, n = 3)

pasilla_col
```

# Heatmap with `pheatmap::pheatmap()`


```{r heatmap, eval=FALSE}
pheatmap::pheatmap(
  mat = pasilla_heat,
  show_rownames = FALSE,
  annotation_col = pasilla_col
)
```


# Side by side plot with [`ggpubr`](https://rpkgs.datanovia.com/ggpubr/)

```{r ggpubr}
p1 <- ggplot(data = pasilla_data,
       aes(x = log2FoldChange,
           y = -log10(padj),
           color = log10(baseMean))) +
  geom_point() + 
  geom_vline(xintercept = c(-2.5, 2.5), 
             lty = 2, color="grey") + 
  theme_few() + scale_color_viridis_c() + 
  labs(x = bquote(~Log[2]~ "fold change"),
       y = bquote(~Log[10]~adjusted~italic(P)),
       title = "Volcano Plot")
p2 <- ggplot(data = pasilla_data,
       aes(x = baseMean,
           y = log2FoldChange,
           color = log10(baseMean))) +
  geom_point() + 
  scale_x_log10() + 
  geom_hline(yintercept = 0, color = "red") + 
  theme_few() + scale_color_viridis_c() + 
  labs(y = bquote(~Log[2]~ "fold change"),
       x = bquote(~Log[10]~ "mean expression"),
       title = "MA Plot")
```


```{r ggpubr_out,out.height="80%", out.width="80%", fig.width=6, fig.height=6}
ggpubr::ggarrange(p1, p2, labels = "AUTO",
  common.legend = TRUE, legend = "bottom")
```

