1.JTK

使用metacycle包进行JTK的分析，其实metacycle包中可选的分析方法有ARSER(Yang, 2010)，JTK_CYCLE( Hughes, 2010)和Lomb-Scargle(Glynn, 2006) 三种
下面只使用JTK进行分析，输入数据长这样

library(MetaCycle)
meta2d(infile="./metacycle/clock_input.csv", 
       filestyle="csv", 
       outdir="./metacycle",
       outIntegration = "onlyIntegration",
       timepoints=rep(seq(1,21,4),6), #时间点1,5,9,13,17,21,1,5,9,...
       #timepoints=rep(seq(1,21,4),each=6), #时间点1,1,1,1,1,1,5,5,5,...
       cycMethod = "JTK")

下面是JTK结果，分别是p值，校正q值，周期，调整相位，振幅

2.cosinor分析

使用cosinor2包进行余弦分析，输入数据长这样

library(cosinor2)
xx <- read.csv("cosinor/clock_cosinor.csv",check.names = F)
x1 <- xx %>% 
  select(1:8) %>%
  arrange(gene,group)
temp=table(x1$gene)[1]
index=length(x1$gene)/temp

analysis <- function(num){
  x=x1
  max=num*temp-0.5*temp
  min=(num-1)*temp+1
  (x=x[min:max,])
  name=x$gene[1]
  group=x$group[1]
  x=x[,3:8]
  fit1 =population.cosinor.lm(data = x, time = as.integer(colnames(x)), period = 24)
  det1=cosinor.detect(fit1)
  det2=cosinor.PR(fit1)
  a = cbind(fit1$coefficients,det1,det2)
  a$gene = name
  a$group=group
  
  x=x1
  max=num*temp
  min=num*temp-0.5*temp+1
  (x=x[min:max,])
  name=x$gene[1]
  group=x$group[1]
  x=x[,3:8]
  fit2 =population.cosinor.lm(data = x, time =  as.integer(colnames(x)), period = 24,plot = F)
  det1=cosinor.detect(fit2)
  det2=cosinor.PR(fit2)
  b = cbind(fit2$coefficients,det1,det2)
  b$gene = name
  b$group=group
  
  
  c=rbind(a,b)
  contrast <- as.data.frame(cosinor.poptests(fit1, fit2))
  contrast$gene=name
  contrast$group="C_vs_N" #fit1来自group C，fit2来自group N
  list1 <- list()
  list1$a <- c
  list1$b <- contrast
  return(list1) 
}
res0 <- data.frame()
res1 <- data.frame()
for (i in seq(1,index)) {
  res = analysis(i)
  res0 = rbind(res0,res$a)
  res1 = rbind(res1,res$b)
}

res0长这样，p表示振荡的显著性，p-value表示观测数据和估计数据之间的相关性是否显著
Tips: Acrophase必小于等于0，其是相对于参考时间0°的弧度制，用(负)度表示，360°(2π)等于周期

res1长这样， 可以看两组的mesor，amp，acr三项的具体平均值和差异是否显著

关于如何画图，我使用的是cosinor包的函数，两组输入数据是一样的，整理稍有不同

xx <- read.csv("cosinor/clock_cosinor.csv",check.names = F)
x1 <- xx %>% 
  select(1:8) %>%
  gather(time,value,3:8) %>% 
  mutate(time=as.numeric(time)) %>%
  arrange(gene) %>%
  mutate(group=ifelse(group=="N",0L,1L))
temp= table(x1$gene)[1]
index =length(x1$gene)

analysis <- function(num){
  x=x1
  max=num*temp
  min=(num-1)*temp+1
  x=x[min:max,]
  fit = cosinor.lm(value ~ time(time)+group+amp.acro(group), data = x, period = 24)
 
  x <- x %>% mutate(levels=ifelse(group=="0"," group = 0"," group = 1")) %>%
    select(gene,levels,time,value)#为了使图例一致而进行的变形
  p <- ggplot.cosinor.lm(fit,x_str = "group")+
    geom_point(aes(time,value,colour=factor(levels)),data = x)+
    theme_classic(base_size = 22)+
    theme(axis.text = element_text(colour = "black"),
          plot.margin=unit(rep(0.3,4),'lines'))+
    scale_color_discrete(name="Group",labels=c("N","C"))+
    scale_x_continuous(limits = c(0,24),breaks = c(1,5,9,13,17,21))+
    labs(x="Time",y=paste0(str_to_title(x$gene[1])," ","expression"))
  p
  ggsave(filename = paste0(str_to_title(x$gene[1]),".png"),plot = p,width = 7,height = 7)
}

for (i in seq(1,index/temp)) {analysis(i)}

其实cosinor包也能做振荡检测和差异分析，但是结果读取不太友好，细细比较cosinor和cosinor2两个包，cosinor出图好看一点，cosinor2结果更易读，说到底两者的结论其实是没有什么差异的

ps：我是初学者，如有错误或遗漏，敬请批评指正