In thomazbastiaanssen/kronos: Microbiome Oriented Circadian Rhythm Analysis Toolkit

Preparation

Here, we will compare the performance of Kronos to that of three common R-based tools to assess rhythmic data. All code available at https://github.com/thomazbastiaanssen/kronos/blob/main/benchmark/

First, we will load the required packages.

#Wrangling
library(tidyverse)

#Plotting
library(ggforce)
library(ggvenn)
library(UpSetR)

#Load circadian rhythm packages
library(kronos) #devtools::install_github("thomazbastiaanssen/kronos") 
library(limorhyde)
library(cosinor)
library(cosinor2)

Now we will load and prepare our high dimensional microbiome data for assessment.

source("00_load_and_clean_WGS.R")

Next, we apply the four packages to the data. Code to run JTK_CYCLE, Limorhyde and Cosinor2 was adapted from their respective vignettes.

source("01a_run_kronos.R")
source("01b_run_jtk_cycle.R")
source("01c_run_limorhyde.R")
source("01d_run_cosinor2.R")

Here, we collate all outcome parameters into a single object to compare them.

res_tot = do.call(cbind, list(res_kronos, res_jtk, res_limo, res_cosinor2))

colnames(res_tot) = paste(colnames(res_tot), c(rep("kronos", ncol(res_kronos)),
                                               rep("jtk",    ncol(res_jtk)),
                                               rep("limo",    ncol(res_limo)),
                                               rep("cosinor2",    ncol(res_cosinor2))),
                          sep = "_")

Comparing p-value distributions between programs

Theoretically, we expect a uniform distribution between 0-1 (for the cases where H~0~ is true) as well as an inflation of low p-values (for the cases where H~1~ is true).

par(mfcol = c(2, 2))

hist(res_kronos$p.val,   breaks = 20)
hist(res_jtk$p.val,      breaks = 20)
hist(res_cosinor2$p.val, breaks = 20)
hist(res_limo$p.val,     breaks = 20)
mtext("P-value histograms", side = 3, line = -1, outer = TRUE)

Both Limorhyde and Kronos show a desirable distribution of p-values, whereas JTK_CYCLE has an inflation on the higher end and Cosinor2 has a skewed, non-uniform distribution in the H~0~ region.

Comparing acrophases between programs

res_tot %>%
  rownames_to_column("ID") %>%
  dplyr::select(!gene_id_limo) %>% 
  group_by(ID) %>% 
  filter(min(p.val_kronos,p.val_limo,p.val_jtk,p.val_cosinor2) < 0.05) %>% 
  ungroup() %>% 

  ggplot() +
  aes(x = .panel_x, y = .panel_y) +
  geom_point(aes(x = .panel_x, y = .panel_y)) +
  ggforce::facet_matrix(vars(acro_kronos, acro_jtk, acro_cosinor2)) +
  theme_bw() +
  ggtitle("Concordance in acrophase detection between Kronos, JTK_CYCLE and Cosinor2")

All three programs that produce acrophases are in agreement as to where the acrophases of features foudn to be rhythmic at p < 0.05 at least once are. Limorhyde does not have a clear method to return acrophases.

Plot agreement in hits between programs

In the venn diagram we see how many features are considered hits (p < 0.05) by the respective programs Alternatively, the UpSetR plot below shows the same.

venn_plot <- list(kronos    = row.names(res_tot[res_tot$p.val_kronos   < 0.05,]),
     jtk_cycle = row.names(res_tot[res_tot$p.val_jtk      < 0.05,]),
     limorhyde = row.names(res_tot[res_tot$p.val_limo     < 0.05,]),
     cosinor2  = row.names(res_tot[res_tot$p.val_cosinor2 < 0.05,])) %>% 
  ggvenn::ggvenn(text_size = 2.5) + 
  ggtitle("Venn diagram showing agreement between\nfeatures rhythmic at p < 0.05")

venn_plot$layers[[3]]$data$x <- c(-1, -0.8, 0.8, 1) #adjust margins
venn_plot

list(kronos    = row.names(res_tot[res_tot$p.val_kronos   < 0.05,]),
     jtk_cycle = row.names(res_tot[res_tot$p.val_jtk      < 0.05,]),
     limorhyde = row.names(res_tot[res_tot$p.val_limo     < 0.05,]),
     cosinor2  = row.names(res_tot[res_tot$p.val_cosinor2 < 0.05,])) %>% 
  UpSetR::fromList() %>% 
  UpSetR::upset(., order.by = "freq", 
                mainbar.y.label = "Number of p-values < 0.05\n per intersection", 
                sets.x.label = "Number of p-values < 0.05\n per program")

We can see that Kronos never reports a significntly rhythmic feature that at least one other program doesn't also see. Furthermore, Kronos finds the second highest total number of rhythmic features at p < 0.05. Only Cosinor2 finds more hits that Kronos, but 303 of those are not corroborated by any other tool and thus may be spurious.

From this, we can conclude that Kronos performs as well as (or better than) the best comparable program in each category. Furthermore, Kronos allows for convenient testing of arbitrarily complex models and comes with pre-designed ggplot-compatible plotting functions.