【迁移】使用 JTK_CYCLE 算法分析生物节律

Last updated on January 20, 2024 pm

JTK是一种非参数检测程序,能从芯片数据中检测循环转录本。除了计算每个转录本最佳的相位(LAG)、振幅(AMP)和周期(PER)外,JTK还计算了置换检验P值(ADJ.P)和Benjamini-Hochberg q值 (BH.Q)。与常规的周期检测算法相比,JTK具有更好的检验效能、更高的计算效率和更强的鲁棒性。R语言的metacycle包实现了ARSERJTK_CYCLELomb-Scargle三种分析方法。

conda安装包

1
2
3
4
5
6
7
8
conda create -n metacycle -c conda-forge r-base=4.1.3
conda activate metacycle
conda install -c conda-forge r-metacycle -y
conda install -c conda-forge r-cosinor -y
conda install -c conda-forge r-tidyverse -y
conda install -c conda-forge r-irkernel -y
Rscript -e "IRkernel::installspec(name='metacycle', displayname='r-metacycle')"
install.packages('cosinor2')

准备数据

通过比对,得到的counts矩阵

1
2
3
4
5
6
require(MetaCycle)
require(tidyverse)
tmp <- read.csv('zctcount.csv', row.names = 1)
head(cycMouseLiverRNA[,1:5])
tmp <- tmp[c('ID','ZT16.con.1', 'ZT16.con.2', 'ZT16.con.3', 'ZT16.con.4', 'ZT28.con.1', 'ZT28.con.2', 'ZT28.con.3', 'ZT28.con.4')]
write.csv(tmp, file="tmp.csv", row.names=FALSE)

f_counts2TMM

1
2
3
4
tmp <- read.csv('tmp.csv', row.names = 1)
tmp <- f_counts2TMM(tmp)
tmp <- log2(tmp + 1)
write.csv(tmp, file="tmp.csv")

分析周期节律

1
2
3
meta2d(infile="tmp.csv",filestyle="csv",outdir="example", cycMethod="JTK", timepoints=c(16,16,16,16,28,28,28,28),outRawData=TRUE)
r1 <- read.csv('example/JTKresult_tmp.csv')
r2 <- read.csv('example/meta2d_tmp.csv')

绘制热图

1
2
3
4
5
6
require(pheatmap)
r1 <- read.csv('example/JTKresult_tmp.csv')
r1 <- subset(r1, ADJ.P < 0.1)
dat <- r1[-(1:6)]
p3 <- pheatmap(dat, scale = "row", cluster_cols = F,
border_color = NA, show_rownames = F,treeheight_row = 0)

绘制某个基因的拟合曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
get_sin_lm <- function(PER, LAG, AMP, mean=0){
function(xvar){
-(AMP/2) * cos(2 * pi* ((xvar + LAG)/PER)) + mean
}
}
dat <- r1[-(1:6)]
dat <- as.data.frame(t(dat))
colnames(dat) <- r1$CycID
index = 5
a <- get_sin_lm(r1$PER[[index]], r1$LAG[[index]], r1$AMP[[index]], mean = mean(unlist(r1[index,-(1:6)])))
mean(unlist(r1[index,-(1:6)]))
dat1 <- dat[index]
dat1$x <- c(16,16,16,16,28,28,28,28)
names(dat1) <- c('y', 'x')
ggplot(dat1, aes(x=x, y = y)) +
stat_function(fun = a) + geom_point()

【迁移】使用 JTK_CYCLE 算法分析生物节律
https://hexo.limour.top/shi-yong--JTK-CYCLE--suan-fa-fen-xi-sheng-wu-jie-lv
Author
Limour
Posted on
September 19, 2022
Updated on
January 20, 2024
Licensed under