Introduction to MetaCycle

In MetaCycle: Evaluate Periodicity in Large Scale Data

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.

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.

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:

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

The corresponding experimental values:

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.

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.

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.

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).

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 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

# 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}) + Acos(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.

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.

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.

MetaCycle documentation built on May 2, 2019, 9:14 a.m.