---
title: "Quantifying differential rhythmicity between conditions"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quantifying differential rhythmicity between conditions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = '#>',
  message = FALSE,
  fig.align = 'center',
  fig.retina = 2)
```

## Introduction

Here we show how to use `limorhyde2` to quantify rhythmicity and differential rhythmicity in data from multiple conditions. The data are based on liver samples from wild-type and Rev-erb$\alpha/\beta$ double-knockout mice ([Cho et al. 2012](https://doi.org/10.1038/nature11048) and [GSE34018](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34018)).

## Load packages

```{r load_packages}
library('data.table')
library('ggplot2')
library('limorhyde2')
library('qs')

# doParallel::registerDoParallel() # register a parallel backend to minimize runtime
theme_set(theme_bw())
```

## Load the data

The expression data are in a matrix with one row per gene and one column per sample. The metadata are in a table with one row per sample. To save time and space, the expression data include only a subset of genes.

```{r load_data}
y = GSE34018$y
y[1:5, 1:5]

metadata = GSE34018$metadata
metadata
```

## Fit linear models and compute posterior fits

Because the samples were acquired at relatively low temporal resolution (every 4 h), we use three knots instead of the default four, which reduces the flexibility of the spline curves. We specify `condColname` so `getModelFit()` knows to fit a differential rhythmicity model.

```{r fit}
fit = getModelFit(y, metadata, nKnots = 3L, condColname = 'cond')
fit = getPosteriorFit(fit)
```

## Get rhythm statistics

Next, we use the posterior fits to compute rhythm statistics for each gene in each condition.

```{r rhythm_stats}
rhyStats = getRhythmStats(fit)
print(rhyStats, nrows = 10L)
```

## Get differential rhythm statistics

We can now calculate the rhythmic differences for each gene between any two conditions, here between wild-type and knockout.

```{r diff_rhythm_stats}
diffRhyStats = getDiffRhythmStats(fit, rhyStats)
print(diffRhyStats, nrows = 10L)
```

We can examine the distributions of the statistics in various ways, such as ranking genes by difference in peak-to-trough amplitude (no p-values necessary) or plotting difference in peak-to-trough amplitude vs. difference in mean expression.

```{r plot_stats, fig.width = 4, fig.height = 4}
print(diffRhyStats[order(diff_peak_trough_amp)], nrows = 10L)

ggplot(diffRhyStats) +
  geom_point(aes(x = diff_mesor, y = diff_peak_trough_amp), alpha = 0.2) +
  labs(x = bquote(Delta * 'mesor (norm.)'), y = bquote(Delta * 'amplitude (norm.)'))
```

## Get observed and fitted time-courses

We can compute the expected measurements for one or more genes at one or more time-points in each condition, which correspond to the fitted curves. Here we plot the posterior fits and observed expression for three genes (converting from gene id to gene symbol).

```{r expected_meas}
genes = data.table(
  id = c('13170', '12686', '26897'),
  symbol = c('Dbp', 'Elovl3', 'Acot1'))

measFit = getExpectedMeas(fit, times = seq(0, 24, 0.5), features = genes$id)
measFit[genes, symbol := i.symbol, on = .(feature = id)]
print(measFit, nrows = 10L)
```

Next we combine the observed expression data and metadata. The curves show how `limorhyde2` is able to fit non-sinusoidal rhythms.

```{r plot_timecourse, fig.width = 7, fig.height = 2.75}
measObs = mergeMeasMeta(y, metadata, features = genes$id)
measObs[genes, symbol := i.symbol, on = .(feature = id)]
print(measObs, nrows = 10L)

ggplot() +
  facet_wrap(vars(symbol), scales = 'free_y', nrow = 1) +
  geom_line(aes(x = time, y = value, color = cond), data = measFit) +
  geom_point(aes(x = time %% 24, y = meas, color = cond, shape = cond),
             size = 1.5, data = measObs) +
  labs(x = 'Zeitgeber time (h)', y = 'Expression (norm.)',
       color = 'Condition', shape = 'Condition') +
  scale_x_continuous(breaks = seq(0, 24, 4)) +
  scale_color_brewer(palette = 'Dark2') +
  scale_shape_manual(values = c(21, 23)) +
  theme(legend.position = 'bottom')
```