tidy-genomics-talk.Rmd

---
title: "Tidy Analysis of Genomic Data"
author: |
  | Michael Love
  | Dept of Genetics &
  | Dept of Biostatistics
  | UNC-Chapel Hill
date: "UVA ~ October 2023"
output: beamer_presentation
urlcolor: blue
---

```{r setup, echo=FALSE}
suppressPackageStartupMessages(library(tidyverse))
knitr::opts_chunk$set(cache = TRUE)
```

# Data organization depends on purpose

![](non-tidy.png)

# "Tidy data" is organized for programming

One row per observation, one column per variable

```{r include=FALSE}
dat <- read_delim("data.tsv")
dat$value <- runif(nrow(dat))
dat$drug <- factor(dat$drug)
```

```{r echo=FALSE}
head(dat)
```

# The pipe

```
command | command | command > output.txt
```

\vspace{2em}

> "Pipes rank alongside the hierarchical file system and regular expressions as one of the most powerful yet elegant features of Unix-like operating systems."

<http://www.linfo.org/pipe.html>

\vspace{2em}

In R we use `%>%` or `|>` instead of `|` to chain operations.

# Verb-based operations

In the R package *dplyr*:

\small
* `mutate()` adds new variables that are functions of existing variables.
* `select()` picks variables based on their names.
* `filter()` picks cases based on their values.
* `slice()` picks cases based on their position.
* `summarize()` reduces multiple values down to a single summary.
* `arrange()` changes the ordering of the rows.
* `group_by()` perform any operation by group.

<https://dplyr.tidyverse.org/>
\normalsize

# Summarize after grouping

A useful paradigm is to *group* data and then *summarize*:

```{r eval=FALSE}
dat %>%
  filter(!outlier) %>%
  group_by(drug, genotype) %>%
  summarize(mu_hat = mean(value))
```

# Summarized output

```{r echo=FALSE, message=FALSE}
dat %>%
  filter(!outlier) %>%
  group_by(drug, genotype) %>%
  summarize(mu_est = mean(value))
```

# Piping directly into plots facilitates data exploration

```{r fig.dim=c(5,2)}
dat %>%
  mutate(newvalue = value^2) %>%
  ggplot(aes(genotype, newvalue)) + 
  geom_boxplot() + 
  facet_wrap(~drug)
```

# Summary I

* I teach both base R and "tidy"
* Both are wrappers, choose based on 1) efficiency 2) flow
* I use the former for writing software, latter for scripting
* Students know dplyr/ggplot2 already
* Next:
  - tidy for genomic ranges
  - tidy for matrix data (scRNA-seq)

# Genomic range data is already tidy

![](narrowpeak.png)

# Great packages in Bioconductor to work with ranges

* [LOLA](https://code.databio.org/LOLA/) - facilitates testing overlaps, fast, useful databases
* [COCOA](https://code.databio.org/COCOA/) - explore sample variation along genome
* [GenomicDistributions](http://code.databio.org/GenomicDistributions/) - annotate, visualize distribution with respect to other features (genes)
* [regioneR]( https://bioconductor.org/packages/regioneR/) - permutation testing
* [ChIPpeakAnno](https://bioconductor.org/packages/ChIPpeakAnno/) - facilitates downstream analysis

Going to talk now about data exploration

# Exploring data with tidy syntax

\large
Helps avoid intermediate variables, and tucks away control code

\vspace{1em}

```{r eval=FALSE}
dat3 <- dat2[dat2$signal > 5]

# vs.

dat %>%
  filter(signal > 5)
```

\normalsize

```{r echo=FALSE, fig.align="center", out.width="25%"}
knitr::include_graphics("plyranges.png")
```

This is *plyranges* from Stuart Lee, Michael Lawrence and Di Cook

# Bringing range data into R

ENCODE mouse embryonic fibroblast, H3K4me1:

\vspace{1em}

```{r echo=FALSE}
suppressPackageStartupMessages(library(plyranges))
```

```{r}
library(plyranges)
pks <- read_narrowpeaks("ENCFF231UNV.bed.gz")
```

or equivalently:

```{r eval=FALSE}
pks <- read.csv("file.csv") %>% 
  rename(seqnames = chr) %>%
  as_granges()
```

```{r echo=FALSE}
#library(GenomeInfoDb)
#si <- Seqinfo(genome="mm10")
#si <- keepStandardChromosomes(si)
#save(si, file="si.rda")
load("si.rda")
seqlevels(pks) <- seqlevels(si)
seqinfo(pks) <- si
```

# Another common paradigm, separating single column

```{r eval=FALSE}
pks <- read.delim("file.tsv") %>%
  tidyr::separate_wider_delim(
    location, 
    delim=":|-", # e.g. chr1:123-456
    into=c("seqnames","start","end")
  ) %>%
  as_granges()
```

# Ranges are rows, metadata are columns

\footnotesize
```{r}
pks %>% 
  slice(1:3) %>% # first 3 ranges
  select(signalValue) # just one metadata column
```
\normalsize

# Example use of *plyranges*

\Large

* Suppose query ranges, `tiles` (e.g. ~1 Mb)
* Find all overlaps between `pks` and `tiles`
* Perform computation on the overlaps
* Many other choices in Bioc for enrichment (e.g. LOLA)

\normalsize

# Example use of *plyranges*

```{r echo=FALSE}
tile0 <- data.frame(seqnames="chr1", 
                    start=51e6 + 1, 
                    width=3e6) %>%
  as_granges()
tiles <- tile0 %>%
  tile_ranges(1e6) %>%
  select(-partition) %>%
  mutate(tile_id = 1:3)
seqinfo(tiles) <- si
```

Created with `tile_ranges` (see also `tileGenome`):

\vspace{1em}

\footnotesize
```{r}
tiles
```
\normalsize

# Consider genomic overlaps as a `join`

```{r echo=FALSE, fig.align="center", out.width="50%"}
# https://www.flickr.com/photos/hellothomas/5073821890
knitr::include_graphics("woodjoin.png")
```

* We are joining two sources of information by match
* How would you then pick top scoring peak (`pks`) per `tile`?
* What verbs would be involved?

# Consider overlaps as a `join`

\footnotesize
```{r}
pks %>%
  select(score) %>% # just `score` column
  join_overlap_inner(tiles) %>% # overlap -> add cols from tiles
  group_by(tile_id) %>% # group matches by which tile
  slice(which.max(score)) # take the top scoring peak
```
\normalsize

# Counting overlaps

* Use "`.`" to specify self within a command
* Add number of overlaps to each entry in `tiles`:
* Can specify `maxgap` and/or `minoverlap`

\vspace{1em}

\footnotesize
```{r}
tiles %>% 
  mutate(n_overlaps = count_overlaps(., pks))
```
\normalsize

# More complex cases

* For peaks near genes, compute correlation of cell-type-specific accessibility and expression (Wancen Mu) → similar to COCOA
* For regulatory variants falling in open chromatin peaks, visualize their distribution stratified by SNP and peak categories (Jon Rosen)
* For looped and un-looped enhancer-promoter pairs, compare average ATAC and RNA time series, while controlling for genomic distance and contact frequency (Eric Davis)

# Nest $\rightarrow$ map $\rightarrow$ unnest

```{r eval=FALSE}
library(purrr)
library(broom)
pks %>%
  join_overlap_inner(tiles) %>%
  as_tibble() %>%
  select(tile_id, signalValue, qValue) %>%
  nest(data = -tile_id) %>%
  mutate(fit = map(data, 
                   ~lm(signalValue ~ qValue, data=.)
                   ),
         stats = map(fit, glance)) %>%
  unnest(stats)
```

# Nest $\rightarrow$ map $\rightarrow$ unnest

```{r echo=FALSE}
library(purrr)
library(broom)
pks %>%
  join_overlap_inner(tiles) %>%
  as_tibble() %>%
  select(tile_id, signalValue, qValue) %>%
  nest(data = -tile_id) %>%
  mutate(fit = map(data, ~lm(signalValue ~ qValue, data=.)),
         stats = map(fit, glance)) %>%
  unnest(stats) %>%
  select(tile_id, data, fit, r.squared)
```

# More *plyranges*-based tutorials online

* *plyranges* vignettes (on Bioc and GitHub)
* Enrichment of peaks and genes: "Fluent Genomics" workflow
* Null regions: *nullranges* vignettes (on Bioc and GitHub)
* Other examples, incl. bootstrap: "Tidy Ranges Tutorial"
* `#tidiness_in_bioc` and `#nullranges` Slack channels

# Summary: tidy analysis for genomic range data

```{r echo=FALSE, fig.show="hold", fig.align="center", out.width="25%"}
knitr::include_graphics(c("dplyr.png","GenomicRanges.png"))
```

```{r echo=FALSE, fig.show="hold", fig.align="center", out.width="25%"}
knitr::include_graphics(c("plyranges.png","nullranges.png"))
```

\small
*nullranges* development sponsored by CZI EOSS ![](czi.png){width=50px}
\normalsize

# Tidy analysis of matrix data

```{r echo=FALSE, fig.align="center", out.width="50%"}
knitr::include_graphics("tt_roadmap.png")
```

tidy-* from Stefano Mangiola (WEHI) *et al.*

# Example use of tidySingleCellExperiment

```{r message=FALSE, echo=FALSE}
library(tidySingleCellExperiment)
sce <- tidySingleCellExperiment::pbmc_small
library(scran)
var_genes <- sce %>%
    modelGeneVar() %>%
    getTopHVGs(prop=0.1)
library(scater) # for next chunk
library(ggplot2) # for next chunk
```

```{r fig.dim=c(4,3), fig.align="center", out.width="50%"}
sce %>%
  scater::runPCA(ncomp=2, subset_row=var_genes) %>%
  ggplot(aes(PC1, PC2, color=groups)) + 
  geom_point()
```

# Example use of tidySingleCellExperiment

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(ggforce)
colLabels(sce) <- sce %>%
    buildSNNGraph(use.dimred="PCA") %>%
    igraph::cluster_walktrap() %$%
    membership %>%
    as.factor()
```

```{r fig.dim=c(4,3), fig.align="center", out.width="50%", message=FALSE}
sce %>%
  join_features(c("CCL5","CST3")) %>%
  ggplot(aes(label, .abundance_logcounts)) + 
  geom_violin() +
  geom_sina() +
  facet_wrap(~.feature)
```

# More complex cases

* Join extra cell-level data
* Perform nested analyses per cell population
* Create a custom expression signature from subset of genes
* Find genes near ChIP-seq peaks, convert to pseudobulk, plot

See [our Bioc2023 workshop](https://tidyomics.github.io/tidyomicsWorkshopBioc2023/articles/tidyGenomicsTranscriptomics.html)
and [tidyseurat](https://stemangiola.github.io/tidyseurat/) / [tidySCE](https://stemangiola.github.io/tidySingleCellExperiment/)

# Altogether, "tidyomics"

<https://github.com/tidyomics>

```{r echo=FALSE, fig.show="hold", fig.align="center", out.width="30%"}
knitr::include_graphics(c("tidyomics1.png", "tidyomics2.png"))
```

# Reading

\small
* Hutchison, WJ, Keyes, TJ, *et al.* The tidyomics ecosystem: Enhancing omic data analyses *bioRxiv* (2023) [10.1101/2023.09.10.557072](https://doi.org/10.1101/2023.09.10.557072)
* Lee, S, Cook, D, Lawrence, M. plyranges: a grammar of genomic data transformation. *Genome Biology* (2019) [10.1186/s13059-018-1597-8](https://doi.org/10.1186/s13059-018-1597-8)
* Lee S, Lawrence M, Love MI. Fluent genomics with plyranges and tximeta. *F1000Research* (2020) [10.12688/f1000research.22259.1](https://doi.org/10.12688/f1000research.22259.1)

Tidy analysis for matrix data:

* Mangiola, S, Molania, R, Dong, R et al. tidybulk: an R tidy framework for modular transcriptomic data analysis. *Genome Biology* (2021) [10.1186/s13059-020-02233-7](https://doi.org/10.1186/s13059-020-02233-7)
* tidySE, tidySCE, tidyseurat
  [stemangiola.github.io/tidytranscriptomics](https://stemangiola.github.io/tidytranscriptomics)

# Extra slides

# plyranges pointers

* TSS: `anchor_5p() %>% mutate(width=1)`
* Overlaps can specify `*_directed` or `*_within`
* Flatten/break up ranges: `reduce_ranges`, `disjoin_ranges`
* Concatenating ranges: `bind_ranges` with `.id` argument
* Overlaps are handled often with "joins": `join_overlap_*`, 
  `join_nearest`, `join_nearest_downstream`, etc.
* Also `add_neareast_distance`
* Load *plyranges* last to avoid name masking with *AnnotationDbi*
  and *dplyr*