limour/limour-blog.github.io @ 3ed496751516045820e5d63460672f6b334eb106

title: 【迁移】使用 JTK_CYCLE 算法分析生物节律 urlname: shi-yong--JTK-CYCLE--suan-fa-fen-xi-sheng-wu-jie-lv date: 2022-09-19 20:17:29

tags: JTK_CYCLE

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

conda安装包

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矩阵

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

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

分析周期节律

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')

绘制热图

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)

绘制某个基因的拟合曲线

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()

2022-09-19-【迁移】使用-JTK-CYCLE-算法分析生物节律.md 2.6 KB Istoric Crud