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：我是初学者，如有错误或遗漏，敬请批评指正
