利用masigpro包进行时间序列分析

 这边文章的代码主要源于网上已有的精品推文，根据自身数据和作图时出现的报错“Error in plot.new() : figure margins too large”进行了个性化修改

值得注意的是：maSigPro软件包主要用于差异基因分析，因此其内置的函数和方法都是基于差异基因的分析。它可以通过比较多组样本，识别出在多个时间点或不同组之间显著变化的基因，并对这些差异基因进行进一步分析。因此，它并不适用于对整体基因进行聚类分析。若您需要对整体基因进行聚类分析，建议使用其他聚类分析软件，例如R语言中的hclust()函数、kmeans()函数等(chatgpt给出的建议)

下面是输入文件格式示例

gene.fpkm.xlsx（实际基因应该有成百上千个，太少了结果不可信）

NC代表对照组，SC代表卒中组；01/03/07/14/28为五个不同时间点，R1/2/3代表三个生物学重复

my_group_design.xlsx

 下面是具体代码部分

rm(list = ls())
setwd("D:/AAA/【R】/maSigPro")
getwd()
dir()

library("maSigPro")
library("openxlsx")

# 载入并整理自己的数据
my_data <- read.xlsx("gene.fpkm.xlsx")
my_group_design <- read.xlsx("my_group_design.xlsx")
rownames(my_data)=my_data[,1]
my_data=my_data[,-1]
rownames(my_group_design)=my_group_design[,1]
my_group_design=my_group_design[,-1]

# 列出时间点和对应样本数
table(as.data.frame(my_group_design)[,"Time"])
head(my_group_design)

# 根据时间点进行回归模型中degree的设置
design <- make.design.matrix(my_group_design, degree = 4) # degree代表自由度，取值为时间点减去1。
design$groups.vector

# 一般作为差异基因的初步筛选
fit <- p.vector(my_data,
                design)
                Q = 0.05,
                MT.adjust = "BH",
                min.obs = 20)

# 差异基因数目
fit$i

# 差异基因表达矩阵
fit$SELEC

# 差异检验的p值
fit$p.vector
fit$p.adjusted

# 写出差异结果
write.csv(fit$SELEC,"my_degs_matrix.csv")

# 对fit中的结果进行组间比较，计算在组间表达显著差异的基因
# 有两种拟合方法：backward，forward，这里用backward
tstep <- T.fit(fit, step.method = "backward", alfa = 0.05)

# 提取组间差异表达分析结果
sigs <- get.siggenes(tstep, rsq = 0.6, vars = "groups")
# 进行可视化 Stroke vs control
png("plot02.png", width = 800, height = 800)
see.genes(sigs$sig.genes$StrokevsControl,
          show.fit = F,
          dis =design$dis,
          cluster.method="hclust" ,
          cluster.data = 1, k = 6)
dev.off()

#提取cluster列表
names(sigs$sig.genes)
png("plot02.png", width = 800, height = 800)
cl = see.genes(sigs$sig.genes$StrokevsControl,newX11 = F,alfa = 0.05, k = 6)## 对差异结果进行聚类分析
dev.off()
cluster = as.data.frame(cl$cut)  ## 提取聚类信息
head(cluster)
write.csv(cluster,"DElncRNA_cluster.csv")