cosinor分析之前已经介绍过两个了（cosinor，cosinor2），其实R中还有三个包可以进行cosinor分析（card::cosinor，season::cosinor和psych::cosinor）,下面详细介绍一下

1.psych::cosinor

#构造测试数据集
set.seed(1)
sim1 <- data.frame(time=rep(as.POSIXlt(seq(1,25,4)*3600,origin = "2021-01-01 00:00:00"),2),
                   zt_time=rep(seq(1,25,4),2),
                   group=rep(c("WT","KO"),each=7),
                   value=c(seq(1,4),seq(3,1),seq(8,2,-2),seq(4,8,2))+rnorm(14,mean = 0,sd = 0.5))
sim1
#     time                zt_time group     value
# 1  2021-01-01 09:00:00       1    WT 0.6867731
# 2  2021-01-01 13:00:00       5    WT 2.0918217
# 3  2021-01-01 17:00:00       9    WT 2.5821857
# 4  2021-01-01 21:00:00      13    WT 4.7976404
# 5  2021-01-02 01:00:00      17    WT 3.1647539
# 6  2021-01-02 05:00:00      21    WT 1.5897658
# 7  2021-01-02 09:00:00      25    WT 1.2437145
# 8  2021-01-01 09:00:00       1    KO 8.3691624
# 9  2021-01-01 13:00:00       5    KO 6.2878907
# 10 2021-01-01 17:00:00       9    KO 3.8473058
# 11 2021-01-01 21:00:00      13    KO 2.7558906
# 12 2021-01-02 01:00:00      17    KO 4.1949216
# 13 2021-01-02 05:00:00      21    KO 5.6893797
# 14 2021-01-02 09:00:00      25    KO 6.8926501

用法也非常简单

cosinor(angle = "zt_time",#时间点
        x="value",#检测基因的表达值
        code = "group",#分组
        data = sim1)
######或者#######
cosinor(angle = 2,
        x=4,
        code = "group",#code参数不能为列号
        data = sim1)

1到6列分别为相位，拟合性能，振幅，标准差，均值（基线），截距

可惜的没有震荡性检测的p值和振荡差异总体p值比较，也无法比较相位，振幅，基线差异的p值

注意：经过比较，该包拟合的cosinor曲线参数和下面两个包的拟合结果相差稍大，不太建议使用

绘图

cosinor.plot(angle ="zt_time",#时间点
      x="value",#检测基因的表达值
      IDloc = "group",#分组
      ID="WT"，#组中变量名
      data=sim1)
######或者#######
cosinor.plot(angle = 2,
      x=4,
      IDloc = 3,
      ID="KO",
      data=sim1)

2.card::cosinor

沿用上面构造的模拟数据集，该函数用法如下：

m <- cosinor(value ~ zt_time, #公式
             data = sim1[1:7,], #数据
             tau = 24)#欲拟合的周期
summary(m)

红线第一处为个体cosinor模型，还有一种是总体cosinor模型，我还没弄清楚到底啥意思，运行时总是报错，暂且不表，
第二处为拟合的公式，
第三处为拟合的周期，
第四处为拟合的参数：

分别为均值（基线），振幅，而相位要换算一下，abs(phi1)/2π×24，此处为13.055≈13.1

注意：和上面psych::cosinor拟合的结果不太一样，接下来看看下一个函数

3.season::cosinor

沿用上面构造的模拟数据集，该函数用法如下：

res1 = cosinor(value~1, #拟合公式
      date='time', #时间列名，类型必须为POSIXlt类型
      type='hourly', #拟合类型为小时，昼夜节律为24小时
      data=sim1[1:7,])
summary(res1)###单个
plot(res1)

从上到下依次为观测数，振幅（和card::season结果接近），峰值相位，波谷相位，振荡显著性

注意：sim1数据集中time和zt_time要换算一下，此处的时间为time，前面两个函数使用的都是zt_time，所以time=21.1即为zt_time=13.1，这点来看，也是和card::cosinor结果较为接近，与psych::cosinor相差一点点

回归系数矩阵中第一行第一列为均值（基线，同样与card::cosinor结果较为接近），
cosw和sinw两行最后一列 两个p值只要有一个即说明振荡具有显著性

比较的话，这样做：

res = cosinor(value~group, date='time', type='hourly', data=sim1)
# summary(res)###比较
# plot(res)
summary(res$glm)

差异比较的话有一个值比较有用，图中红框位置，但是只有一个总体p值，该函数无法比较相位，振幅，基线差异的p值

总结：这三个包中的cosinor函数都稍稍有点鸡肋，两个连p值都没有，统计学上没什么太大意义，一个检测单个基因单组好使，但是对于检测单个基因组间差异就爱莫能助了，因为它只有振荡差异总体的p值，就少那么一点点功能（组间相位，振幅，基线差异的p值）就完美了，此外，两个包（psych和season）都是基于base系统作图，作图功能聊胜于无，另外一个包card倒是提供了基于ggplot2的函数ggcosinor，可是不支持时间点样本无重复画法。不过，用来评估单个基因的cosinor曲线，season::cosinor倒是个不错的选择，参数齐全且准确（和card::cosinor结果几乎一致），同时提供p值。