---
title: "Introduction to MetaCycle"
author: "Gang Wu, Xavier Li, Matthew Carlucci, Ron Anafi, Michael Hughes, Karl Kornacker, and John Hogenesch"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to MetaCycle}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
---
The **MetaCycle** package is mainly used for detecting rhythmic signals from large scale time-series data. Depending on features of each time-series data, **MetaCycle** incorporates ARSER(ARS), JTK_CYCLE(JTK), and Lomb-Scargle(LS) properly for periodic signal detection, and it could also output integrated analysis results if required. 

This vignette introduces implementation of method selection and integration steps of **MetaCycle**, which are not explained in detail in help files. For looking at how to use two main functions--**meta2d** and **meta3d** of this package, please see 'Examples' part of each function's help file. 

The **MetaCycle** source code is available on [github](https://github.com/gangwug/MetaCycle).

## DATASETS
### Two main categories
A typical time-series dataset from a non-human organism is a two-dimensional matrix. Each row indicates one molecular's profile along time, and all molecule at any one time point are detected from the same sample. It is usually not necessary to keep track of which individual organism does a sample come from. For ease of explanation, we named this kind of dataset as 2D time-series dataset. Take the time-series transcriptome dataset from mouse liver as an example.

```{r}
library(MetaCycle)
head(cycMouseLiverRNA[,1:5])
```

For time-series datasets from human, it is usually essential to keep track of the subject identity and group for each sample. In this case, one matrix stores experimental values of detected molecule from all samples, while another matrix is necessary to store subject identity and group for each sample. This kind of dataset is named as 3D time-series dataset. For example, a time-series dataset from human blood is shown as below. 

The subject information matrix:

```{r}
set.seed(100)
row_index <- sample(1:nrow(cycHumanBloodDesign), 4)
cycHumanBloodDesign[row_index,]
```

The corresponding experimental values:

```{r}
sample_id <- cycHumanBloodDesign[row_index,1]
head(cycHumanBloodData[,c("ID_REF", sample_id)])
```

A 3D time-series dataset could be divided into multiple 2D time-series datasets, with all experimental values for an individual under the same treatment in one 2D time-series dataset. For example, we could extract all experimental values from "AF0004" under "SleepExtension" into one 2D time-series dataset.

```{r}
group_index <- which(cycHumanBloodDesign[, "group"] == "SleepExtension")
cycHumanBloodDesignSE <- cycHumanBloodDesign[group_index,]
sample_index <- which(cycHumanBloodDesignSE[, "subject"] == "AF0004")
sample_AF0004 <- cycHumanBloodDesignSE[sample_index, "sample_library"]
cycHumanBloodDataSE_AF0004 <- cycHumanBloodData[, c("ID_REF", sample_AF0004)]
head(cycHumanBloodDataSE_AF0004)
```

### Detail types of 2D time-series dataset
One kind of usual 2D time-series dataset is evenly sampled once at each time point, and the interval between neighbour time points is integer. Not all datasets are as simple as this. There are datasets with replicate samples, or with missing values, or un-evenly sampled, or sampled with a non-integer interval. Examples of these types of dataset are shown in the below table. 

Data Type                    | Point 1| Point 2| Point 3| Point 4| Point 5| Point 6
---------------------------  | -------| -------| -------| -------| -------| -------
The usual data               | CT0    | CT4    | CT8    | CT12   | CT16   | CT20  
With missing value           | CT0    | NA     | CT8    | CT12   | CT16   | CT20  
With replicates              | CT0    | CT0    | CT8    | CT8    | CT16   | CT16  
With un-even interval        | CT0    | CT2    | CT8    | CT10   | CT16   | CT20
With non-integer interval    | CT0    | CT4.5  | CT9    | CT13.5 | CT18   | CT22.5  

Of course, some datasets may be a combination of two or more of above types of data.

Data Type                             | Point 1| Point 2| Point 3| Point 4| Point 5| Point 6
------------------------------------  | -------| -------| -------| -------| -------| -------
With replicates and missing value     | CT0    | CT0    | CT8    | NA     | CT16   | CT16  
With un-even interval and replicates  | CT0    | CT2    | CT2    | CT10   | CT16   | CT20

## METHOD SELECTION
### The Pros and Cons of meta2d methods
Each method (i.e. JTK, ARS, LS) has its pros and cons. It is ideal to know these trade-offs prior to decisions on study design. For example, the features of meta2d methods are listed in the below figure. 

```{r, echo=FALSE, warning=FALSE, out.width = "640px"}
knitr::include_graphics("./images/image1.png")
```

### Method selection based on the sampling pattern
As shown in the above figure, the accuracy of rhythmic analysis is influenced a lot by the sampling pattern. Although higher sampling resolution usually gives more accurate results, the ideal experimental design depends on the finial research aim. For example, biological replicates are important for identifying differentially expressed genes between two time points (e.g. midnight and noon). The high sampling resolution (2h/2days) is necessary for accurately reporting the percentage of rhythmic genes at genome-wide level. If the goal is to get accurate phase or period length of the rhythmic gene, much higher sampling resolution is required. A brief summary of method selection based on the sampling pattern is shown in the below figure. 

```{r, echo=FALSE, warning=FALSE, out.width = "640px"}
knitr::include_graphics("./images/image2.png")
```

**We suggest to think about method selection during the experimental design step, which will make the latter analysis much easier.** 

**Experimental design >>> statistical methods.**

### The automatic method selection strategy used by meta2d 
The **meta2d** function in **MetaCycle** is designed to analyze 2D time-series datasets, and it could automatically select proper methods to analyze different types of input datasets. The implementation strategy used for **meta2d** is shown in the flow chart (drawn with "diagram" package).

```{r, echo=FALSE, warning=FALSE, out.width = "640px"}
knitr::include_graphics("./images/image3.png")
```

For analyzing 3D time-series datasets, the **meta3d** function in **MetaCycle** is suggested. It firstly divides the input dataset into multiple 2D time-series datasets based on individual information, and then use the defined method through calling **meta2d** to analyze each divided dataset. 

## INTEGRATION
In addition to selecting proper methods to analyze different kinds of datasets, **MetaCycle** could also output integrated results. In detail, **meta2d** integrates analysis results from multiple methods and **meta3d** integrates analysis results from multiple individuals.

### Pvalue
[Fisher's method](https://en.wikipedia.org/wiki/Fisher%27s_method) is implemented in both **meta2d** and **meta3d** for integrating multiple p-values. The below formula is used to combine multiple p-values into one test statistic (X^2^). 

$$X^2_{2k} \sim -2\sum_{i=1}^k ln(p_i)$$

X^2^ has a chi-squared distribution with 2k degrees of freedom (k is the number of p-values), when all the null hypotheses are true, and each p-value is independent. The combined p-value is determined by the p-value of X^2^.

### Notes on the Fisher's method
We acknowledge that Hutchison A.L., et al. pointed out the p-value integration issue using Fisher's method in **meta2d** function in two papers (*Correcting for Dependent P-values in Rhythm Detection* and *Bootstrapping and empirical Bayes methods improve rhythm detection in sparsely sampled data*. There is a similar issue mentioned in the MetaCycle paper discussing about the performance of N-version programming (NVP). *"In some cases, NVP may not give better results than the single most suitable method, but it will rarely give the worst results. For example, in analyzing time-series datasets with high temporal resolution (every 1 h over 2 days; Supplementary Fig.S3), NVP does better than ARS but not as well as LS or JTK (Under these conditions, all three methods do relatively well but have similar mode failure, while ARS has relatively more false positive observations.)."*

Hutchison A.L., et al. also mentioned in the paper, *"While the combinations of the methods integrated improperly with the Fisher method appear to outperform the individual methods at p-values below typical significance cutoffs (Figs. 4C andD, dark gray), once the p-values are accurately calculated using the Brown method, the combined methods (Figs. 4C and D, light gray) underperform the individual methods for low p-values."*. The advantage of Fisher method comparing with single method is compreshensively evaluated in the MetaCycle paper using different sampling patterns and rhythmic curves. We agree that Fisher's method is not ideal. But it shows its power in practice. So in the current version of **MetaCycle**, Fisher's method is still the default selection. We will update it immediately whenever we find an integration method that is beautiful in theory and powerful in practice.

As we know, each method has its pros and cons. There is not a perfect rhythmic detection method for all time-series datasets. As more rhythmic detection methods are published in coming years, more comprehensive evaluation works are required in this field. In our opinion, one of the most important thing is to generate a series of gold standard evaluation datasets. Multiple methods should be evaluated on the same datasets. Otherwise, it is not an apples to apples comparison when claiming one method is better than others. 

### Period and phase
The integrated period from **MetaCycle** is an arithmetic mean value of multiple periods, while phase integration based on mean of circular quantities is implemented in **meta2d** and **meta3d**. The detailed steps are as below.

* convert phase values to polar coordinates $\alpha_j$
* convert polar coordinates to cartesian coordinates ($cos\alpha_j$, $sin\alpha_j$)
* compute the arithmetic mean of these points and its corresponding polar coordinate $\bar{\alpha}$
$$\bar{\alpha} = atan2(\frac{\sum_{j=1}^n sin\alpha_j}{n}, \frac{\sum_{j=1}^n cos\alpha_j}{n})$$
* convert the resulting polar coordinate to a integrated phase value

```{r, warning=FALSE}
# given three phases
pha <- c(0.9, 0.6, 23.6)
# their corresponding periods
per <- c(23.5, 24, 24.5)
# mean period length
per_mean <- mean(per)
# covert to polar coordinate
polar <- 2*pi*pha/per
# get averaged ploar coordinate
polar_mean <- atan2(mean(sin(polar)), mean(cos(polar)))
# get averaged phase value
pha_mean <- per_mean*polar_mean/(2*pi)
pha_mean
```

### Amplitude calculation
**meta2d** recalculates the amplitude with the following model:

$$Y_i = B + TRE*(t_i - \frac{\sum_{i=1}^n t_i}{n}) + A*cos(2*\pi*\frac{t_i - PHA}{PER})$$

where $Y_i$ is the observed value at time $t_i$; B is baseline level of the time-series profile; TRE is trend level of the time-series profile; A is the amplitude of the waveform. PHA and PER are integrated period and phase mentioned above. In this model, only B, TRE and A are unknown parameters, which could be calculated with the ordinary least squares (OLS) method. The baseline and trend level are explained in the below example.

```{r, echo=FALSE, warning=FALSE, fig.width=6.65, fig.height=5}
getAMP <- function(expr, per, pha, tim=18:65)
{ 
    trendt <- tim - mean(tim[!is.na(tim) & !is.nan(tim)])
    cost <- cos(2*pi/per*(tim - pha))
    fit <- lm(expr~trendt + cost)
    fitcoef <- fit$coefficients
    basev <- fitcoef[1]
    trendv <- fitcoef[2]
    ampv <- fitcoef[3]
    fitexp <- basev + trendv*trendt + ampv*cost
    outL <- list("base"=basev, "trend"=trendv, "amp"=ampv, "fit"=fitexp)
    return(outL)
}

cirD <- cycVignettesAMP
ampL <- getAMP(expr=as.numeric(cirD[1,24:71]), per=cirD[1, "meta2d_period"], pha=cirD[1, "meta2d_phase"])

lay<-layout(cbind(1, 2), widths=c( lcm(cm(4.5)), lcm(cm(1.5)) ), heights=lcm(cm(4.5)) )
par(mai=c(0.65,0.6,0.4,0.05),mgp=c(2,0.5,0),tck=-0.01)
xrange <- c(18, 65)
yrange <- c(200, 2350)
plot(18:65, cirD[1,24:71], type="b", xlim=xrange, ylim=yrange, xlab="Circadian time(CT)", ylab="Expression value",  main=cirD[1,1], cex.main=1.2)
par(new=T)
plot(18:65, ampL[[4]], type="b", xlim=xrange, ylim=yrange, col="red", xlab="", ylab="", main="")
abline(h=ampL[[1]], lty=3, col="purple", lwd=1.5)
lines(18:65, 500+ampL[[2]]*(18:65-mean(18:65)), lty=4, col="orange", lwd=1.5)
legend("topleft", legend=c("Raw value", "OLS fitted value"), col=c("black", "red"), pch=1, bty="n")
legend("topright", legend=c("Baseline", "Trend"), col=c("purple", "orange"), lty=c(3, 4), lwd=1.5, bty="n" )

par(mai=c(0.5,0.05,0.4,0.1),mgp=c(2,0.3,0),tck=-0.01);
plot(x=NULL,y=NULL,xlim=c(0,10),ylim=c(0,10),type="n", xaxt="n",yaxt="n",bty="n",xlab="",ylab="",main="")
text(rep(1,3), c(8, 5, 2), c("Base = ", "Trend = ", "AMP = "), adj=0)
text_value <- unlist(ampL)
text(rep(6,6), c(8, 5, 2), round(text_value[1:3], 1), adj=0)
```

In addition, **meta2d** also outputs a relative amplitude value (rAMP), which could be easily taken as the ratio between amplitude and baseline (if |B| >= 1). The amplitude value is associated with the general expression level, which indicates highly expressed genes may always have larger amplitude than lowly expressed genes. The rAMP may be used to compare the amplitude values among genes with different expression levels. For example, *Ugt2b34* has a larger amplitude than *Arntl*, but its rAMP is smaller than *Arntl*.

```{r, echo=FALSE, warning=FALSE, fig.width=6.65, fig.height=5}
cirD <- cycVignettesAMP
cirM <- as.matrix(cirD[2:3, 24:71])
expmax <- apply(cirM, 1, max)
cirM <- cirM/expmax

lay<-layout(cbind(1, 2), widths=c( lcm(cm(4.5)), lcm(cm(1.5)) ), heights=lcm(cm(4.5)) )
par(mai=c(0.65,0.6,0.4,0.05),mgp=c(2,0.5,0),tck=-0.01)
xrange <- c(18, 65)
yrange <- c(0, 1)
colname <- c("red", "blue")
grey_trans <- rgb(191/255,191/255,191/255,0.65);

par(mai=c(0.65,0.6,0.2,0.05),mgp=c(2,0.5,0),tck=-0.01)
plot(NULL,NULL,xlim=xrange,ylim=yrange,xaxt="n",yaxt="n",xlab="Circadian time(CT)",ylab="Exp/Max", main="");
rect_xL <- c(18, 36, 60)
rect_yL <- c(24, 48, 65)
rect(rect_xL, rep(-0.1, 3), rect_yL, rep(1.1,3), col=grey_trans, border=NA, bty="n")

for (i in 1:2)
{
  loessD <- data.frame(expd=as.numeric(cirM[i,]),tp=18:65);
  exploess <- loess(expd~tp, loessD, span = 0.2);
  expsmooth <- predict(exploess, data.frame(tp=18:65));
  lines(18:65,expsmooth,lwd=1.2,col=colname[i]);
}

xpos <- c(seq(18,60,by=6), 65)
axis(side=1,at=xpos,labels=xpos,mgp=c(0,0.2,0),tck=-0.01,cex.axis=0.8)
ypos <- seq(0, 1, by=0.2)
axis(side=2,at=ypos,labels=ypos,mgp=c(0,0.2,0),tck=-0.01,cex.axis=0.8)

par(mai=c(0.5,0.05,0.2,0.1),mgp=c(2,0.3,0),tck=-0.01)
plot(x=NULL,y=NULL,xlim=c(0,10),ylim=c(0,10),type="n", xaxt="n",yaxt="n",bty="n",xlab="",ylab="",main="")
lines(c(0.2,2.3), c(9.5,9.5), col="blue", lwd=1.5)
text(c(5,0,0), c(9.5, 8, 6.5), c("Ugt2b34", "AMP = ", "rAMP = "), col="blue", adj=0)
text(c(5,5), c(8, 6.5), round(as.numeric(cirD[3, 22:23]), 2), col="blue", adj=0)

lines(c(0.2,2.3), c(4.5,4.5), col="red", lwd=1.5)
text(c(5,0,0), c(4.5, 3, 1.5), c("Arntl", "AMP = ", "rAMP = "), col="red", adj=0)
text(c(5,5), c(3,1.5), round(as.numeric(cirD[2, 22:23]), 2), col="red", adj=0)
```

Based on the calculated baseline, amplitude and relative amplitude values by **meta2d**, **meta3d** calculates the corresponding integrated values with arithmetic mean of multiple individuals in each group. 