1.JTK
使用metacycle包進(jìn)行JTK的分析穿香，其實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), #時間點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結(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π)等于周期
res1長這樣乏盐， 可以看兩組的mesor佳窑，amp，acr三項的具體平均值和差異是否顯著
關(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)}
其實cosinor包也能做振蕩檢測和差異分析何吝，但是結(jié)果讀取不太友好溉委，細(xì)細(xì)比較cosinor和cosinor2兩個包鹃唯，cosinor出圖好看一點，cosinor2結(jié)果更易讀瓣喊，說到底兩者的結(jié)論其實是沒有什么差異的
ps：我是初學(xué)者俯渤，如有錯誤或遺漏，敬請批評指正
