1.JTK

使用metacycle包進(jìn)行JTK的分析，其實(shí)metacycle包中可選的分析方法有ARSER(Yang, 2010)，JTK_CYCLE( Hughes, 2010)和Lomb-Scargle(Glynn, 2006) 三種
下面只使用JTK進(jìn)行分析，輸入數(shù)據(jù)長這樣

library(MetaCycle)
meta2d(infile="./metacycle/clock_input.csv", 
       filestyle="csv", 
       outdir="./metacycle",
       outIntegration = "onlyIntegration",
       timepoints=rep(seq(1,21,4),6), #時間點(diǎn)1,5,9,13,17,21,1,5,9,...
       #timepoints=rep(seq(1,21,4),each=6), #時間點(diǎn)1,1,1,1,1,1,5,5,5,...
       cycMethod = "JTK")

下面是JTK結(jié)果，分別是p值，校正q值，周期，調(diào)整相位，振幅

2.cosinor分析

使用cosinor2包進(jìn)行余弦分析，輸入數(shù)據(jù)長這樣

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表示觀測數(shù)據(jù)和估計數(shù)據(jù)之間的相關(guān)性是否顯著
Tips: Acrophase必小于等于0，其是相對于參考時間0°的弧度制，用(負(fù))度表示，360°(2π)等于周期

關(guān)于如何畫圖，我使用的是cosinor包的函數(shù)，兩組輸入數(shù)據(jù)是一樣的，整理稍有不同

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)#為了使圖例一致而進(jìn)行的變形
  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)}

其實(shí)cosinor包也能做振蕩檢測和差異分析，但是結(jié)果讀取不太友好，細(xì)細(xì)比較cosinor和cosinor2兩個包，cosinor出圖好看一點(diǎn)，cosinor2結(jié)果更易讀，說到底兩者的結(jié)論其實(shí)是沒有什么差異的

ps：我是初學(xué)者，如有錯誤或遺漏，敬請批評指正