这边文章的代码主要源于网上已有的精品推文,根据自身数据和作图时出现的报错“Error in plot.new() : figure margins too large”进行了个性化修改
值得注意的是:maSigPro软件包主要用于差异基因分析,因此其内置的函数和方法都是基于差异基因的分析。它可以通过比较多组样本,识别出在多个时间点或不同组之间显著变化的基因,并对这些差异基因进行进一步分析。因此,它并不适用于对整体基因进行聚类分析。若您需要对整体基因进行聚类分析,建议使用其他聚类分析软件,例如R语言中的hclust()函数、kmeans()函数等(chatgpt给出的建议)
下面是输入文件格式示例
NC代表对照组,SC代表卒中组;01/03/07/14/28为五个不同时间点,R1/2/3代表三个生物学重复
下面是具体代码部分
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")