circacompare 是一个专为分析生物节律数据而设计的 R 包。它的主要功能是比较不同条件下的节律参数，例如振幅、周期和相位。circacompare 使用非线性混合效应模型来拟合节律数据，这使得它在处理具有重复测量和复杂实验设计的数据时表现出色。与 circacompare 相比，MetaCycle 是另一个流行的 R 包，用于生物节律分析。MetaCycle 提供了多种算法（如 ARSER、JTK_CYCLE 和 Lomb-Scargle）来检测时间序列数据中的周期性信号。它的优势在于能够处理大规模数据集，并且适用于各种不同的实验条件。

conda安装包

conda create -n zct conda-forge::r-tidyverse conda-forge::r-irkernel
conda run -n zct Rscript -e "IRkernel::installspec(name='zct', displayname='zct')"
conda run -n zct Rscript -e "install.packages('circacompare')"

导入包和数据

library(tidyverse)
library(circacompare)
library(ggplot2)

dt <- readr::read_csv('./circacompare.CSV') %>% 
mutate(group = factor(group)) %>% 
mutate(organ = factor(organ)) %>% 
mutate(project = factor(project))

两组比较

s_symbol = 'Bmal1' 
s_organ = 'Heart'
s_project = 'compare1'

# 根据参数选择数据
dt_s <- dt %>% 
subset(symbol == s_symbol & organ == s_organ & project == s_project)  %>% 
mutate(group = factor(group))
# 进行比较
result <- circacompare(x = dt_s, col_time = "time", col_group = "group", col_outcome = "measure", alpha_threshold = 1)

# 查看绘图
save_plot <- result$plot + 
theme_minimal() +  
ggtitle(paste(c(s_symbol, s_organ), collapse = '_')) +
theme(plot.title = element_text(hjust = 0.5))

# 查看统计汇总
result$summary
tmp <- circacompare:::extract_model_coefs(result$fit)
tmp
tmp['phi', c('estimate', 'std_error')] <- (tmp['phi', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp['phi1', c('estimate', 'std_error')] <- (tmp['phi1', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp

save_plot

# 保存图为 pdf
{pdf(file = paste0('pdf/', paste(c(s_symbol, s_organ, s_project), collapse = '_'), '.pdf'), width = 6, height = 6)
 print(save_plot)
dev.off()}

# 导出拟合的数据
save_plot
x <- seq(0, 24, by = 0.01)
newdf <- data.frame(time = x, g1 = save_plot$plot_env$eq_1(x), g2 = save_plot$plot_env$eq_2(x))
write.csv(newdf, paste0('csv/', paste(c(s_symbol, s_organ, s_project), collapse = '_'), '.csv'), row.names = F)

单组绘图

s_symbol = 'Rev-erbα' 
s_organ = 'Kidney'
s_project = 'KO-AL'

# 根据参数选择数据
dt_s <- dt %>% 
subset(symbol == s_symbol & organ == s_organ & project == s_project)  %>% 
mutate(group = factor(group))

# 进行统计分析
options(show.error.messages = F, warn = -1)
result <- try({
  circa_single(
    x = dt_s, col_time = "time", col_outcome = "measure", period = 24, alpha_threshold = 1,
    timeout_n = 100000,
    control = list(
      main_params = c("k", "alpha", "phi")
    )
  )
}, silent = TRUE)
options(show.error.messages = T, warn = 1)
# “k”表示中值，“alpha”表示振幅，“phi”表示相位。引入的额外参数是“tau”表示周期。

# 查看绘图
save_plot <- result$plot + 
theme_minimal() +  
ggtitle(paste(c(s_symbol, s_organ), collapse = '_')) +
theme(plot.title = element_text(hjust = 0.5))

# 查看统计汇总
result$summary
tmp <- circacompare:::extract_model_coefs(result$fit)
tmp
tmp['phi', c('estimate', 'std_error')] <- (tmp['phi', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp['phi1', c('estimate', 'std_error')] <- (tmp['phi1', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp

save_plot

# 保存图为 pdf
{pdf(file = paste0('pdf/', paste(c(s_symbol, s_organ, s_project), collapse = '_'), '.pdf'), width = 6, height = 6)
 print(save_plot)
dev.off()}

周期和衰减参数

s_symbol = 'Bmal1' 
s_organ = 'Heart'
s_project = 'compare1'

# 根据参数选择数据
dt_s <- dt %>% 
subset(symbol == s_symbol & organ == s_organ & project == s_project)  %>% 
mutate(group = factor(group))

# 进行统计分析
options(show.error.messages = F, warn = -1)
result <- try({
  circacompare(
    x = dt_s, col_time = "time", col_group = "group", col_outcome = "measure", period = 24, alpha_threshold = 1,
    timeout_n = 100000,
    control = list(
      main_params = c("k", "alpha", "phi"),
      decay_params = c("alpha"),
      grouped_params = c("alpha", "alpha_decay")
    )
  )
}, silent = TRUE)
options(show.error.messages = T, warn = 1)
# “k”表示中值，“alpha”表示振幅，“phi”表示相位。引入的额外参数是“tau”表示周期。

# 查看绘图
save_plot <- result$plot + 
theme_minimal() +  
ggtitle(paste(c(s_symbol, s_organ), collapse = '_')) +
theme(plot.title = element_text(hjust = 0.5))

# 查看统计汇总
result$summary
tmp <- circacompare:::extract_model_coefs(result$fit)
tmp
tmp['phi', c('estimate', 'std_error')] <- (tmp['phi', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp['phi1', c('estimate', 'std_error')] <- (tmp['phi1', c('estimate', 'std_error')] / (2 * pi)) * save_plot$plot_env$V['tau']
tmp

save_plot

# 保存图为 pdf
{pdf(file = paste0('pdf/', paste(c(s_symbol, s_organ, s_project), collapse = '_'), '_decay.pdf'), width = 6, height = 6)
 print(save_plot)
dev.off()}