最近在进行lumicycle实验，接触到了细胞上面的衰减的昼夜节律。
关于这一点，和动物造模有所不同，在细胞上面，由于细胞过于密集+培养时长延长+培养基介质长期不更换的因素，细胞的昼夜节律会越来越衰减，关于这种情况下我们应该怎么计算呢？
首先，查阅资料，了解Lumicycle这一仪器的使用，这个仪器一般我们接触不到，在这里就不在详细描述了，只要知道它可以用来记录细胞的昼夜节律就可以了，当然动物的组织切片也是可以的，就是操作会更复杂一点。
关于我们所使用的细胞，我们使用的Bmal1::dluc U2OS细胞和Per2::dluc U2OS细胞。这两个细胞历史也比较悠久了，下面简单介绍一下它们是怎么构建的。
图中各个元件的作用
SIN-LTR 自灭活长末端复制片段，帮助慢病毒将其序列整合到宿主基因组中
RRE REV反应元件，一种约 350 个核苷酸、高度结构化的顺式作用 RNA 元件，对于病毒复制至关重要
PPT HIV-1 中央多嘌呤束，作为正链合成的引物产生一个flap元件，是所有逆转录病毒RNA基因组形成双链cDNA的关键
P SV40 Promoter SV40启动子序列，引入在Per2启动子上游，避免活性位点整合效应
T SV40 terminator SV40终止子序列，引入在Per2启动子上游，与SV40启动子序列合在一起构建“增强子陷阱诱饵”，阻断整合位点附近增强子活性
mPer2 小鼠Per2启动子序列
dluc 含有PEST序列的荧光素酶基因，在PEST的帮助下荧光素酶快速降解以保证不在细胞质内积累，实时反应上游启动子活性
IRES 内部核糖体加入位点，遇到终止密码子翻译不停滞，继续向下翻译
EGFP 增强型绿色荧光蛋白序列，用于流式分选
WPRE 土拨鼠转录后调控元件，提高基因表达
使用如上图所示的慢病毒载体感染U2OS细胞，之后使用流式分选筛查GFP荧光活性较强和振幅较大的细胞，这类细胞被称为Per2::dluc U2OS。Bmal1::dluc U2OS构建方案类似。
我们在这里使用了Bmal1::dluc U2OS细胞，其基础振幅相对于Per2::dluc U2OS要更高一点。
使用lumicycle生物节律性光度测定仪记录细胞的昼夜节律，连续记录7天以上，读入Lumicycle analysis软件，简单查看一下振荡趋势。
很容易看出，细胞的昼夜节律振幅是不断下降的，其中，第一个周期由于更换培养基介质带来的高瞬时发光所以弃之不用。
这样的图像是不易计算昼夜节律的，我们需要对其进行基准化，也就是减去一个Mesor（均值）。
程序自带的基准化算法有三种，分别是Polynomial fit，Running Average，Butterworth，之后我们可以得到基准化后的图像
程序自带的拟合算法有7种，分别是Sin fit，Sin fit (damped)，Sin fit (undamp)，LM Fit (Damped Sine)，LM Fit (Sin)，Periodogram，FFT。
对于这些基准化和拟合算法，我们查阅文献看看别人怎么使用的。
由上可知，我们选择一阶多项式进行基准化，选择指数衰减形式的正（余）弦函数进行曲线拟合。
既然是使用指数衰减形式的正（余）弦函数进行拟合，那么我们完全可以使用circacompare进行计算，见Introduction to circacompare，里面提到了它也可以用来计算指数形式的振幅衰减。
R中操作如下：
#1.收集要读入的文件和实验分组情况
library(readxl)
library(tidyverse)
library(stringr)
library(magrittr)
library(circacompare)
library(ggplot2)
file1 <- list.files(".",pattern = "subtracted.csv")
x <- str_remove(file1,"\\d_subtracted.csv")
unique(x)
group <- factor(x,labels = c("B","CYS","DEX","VNN"))#分组
要读入的数据（经过一阶多项式基准化的）长这样：
第一列日期
第二列时间
第三列相对于零点的天数
第四列基准化后的值
第五列程序衰减余弦拟合的值
当然，也完全可以读入原始数据，即未经过基准化的，然后自己一阶函数线性拟合一下，原始的数据长这样：
第一列日期
第二列时间
第三列相对于零点的天数
第四列原始值
第五列程序一阶多项式拟合的基准值，第四列减去第五列即为上面第四列中基准化后的值
#拟合得基准（自己预测的）值，并与程序输出的基准值比较
x1 <- read.csv("B51_Raw11111.csv",skip = 3,header = F)[,3:5]
colnames(x1) <- c("day","raw","predicted_program")
res <- with(x1,lm(raw~poly(day,1)))
a <-data.frame(day=x1$day)
predicted_myself <- predict(res,a)
a$predicted_myself <- predicted_myself
judge <- a$predicted_myself-x1$predicted_program
summary(judge<0.01)#拟合的基准值与程序输出基准值差异都在0.01内，误差非常小
# Mode TRUE
#logical 1186
#之后可以使用程序输出的基准值，也可以使用我们自己拟合的预测基准值，之后拿原始值减去基准值即可得基准化后的值，可以用于后期衰减cosinor拟合
x1$subtracted_program <- x1$raw-x1$predicted_program
#或者
x1$subtracted_myself <- x1$raw-a$predicted_myself
#2.读入每个样品的振荡数据并合并起来
dat1 <- data.frame()
for (i in seq(24)) {
example1 <- read.csv(file = file1[i],skip = 3,header = F)[,3:5]#跳过前三行，选取第三到五列
colnames(example1) <- c("day","subtracted","fitted")
example1$id <- ifelse(i%%6==0,6,i%%6)
example1$time <- example1$day*24
example1$group <- group[i]
example1$measure <- example1$subtracted
example1 <- example1[,c("id","group","time","measure","day","subtracted","fitted")]
example1 <- example1[example1$day>=1 & example1$day<=6, ]
dat1 <- rbind(dat1,example1)
}
# library(openxlsx)
# write.xlsx(dat1,file = "all.xlsx",overwrite = T)
dat1长这样：
#3.计算每个样本的振荡情况并写出
data2 <- data.frame(parameter=c("rhythmic_p","mesor","amplitude","alpha_decay","phase_radians","peak_time_hours","period"))
for (x in levels(group)) {
for(y in seq(6)){
res <- circa_single(x = dat1 %>% filter(group==x & id==y),
col_time = "time",
col_outcome = "measure",
period = NA,
control=list(period_param=T,decay_params=c("alpha")))
res1 <- res$summary
data2<- cbind(data2,res1$value)
colnames(data2)[ncol(data2)] <- paste0(x,y,"_value")
}
}
library(openxlsx)
write.xlsx(data2,file = "all_stat.xlsx",overwrite = T)
data2长这样：
