上周的公众号处理了不同时间序列的数据集，但因为是内置的数据集，很多分析流程都已经被pipeline函数包装了，那如果是自己的时间序列数据集该怎么分析呢？
曾老师就让我学习一下这个包，今天咱就浅学一下吧~
Package ‘Mfuzz’
以数据集GSE120418为例，是转录组的数据集哦
- 主要内容：Transcriptome-wide analysis of gene expression using detached first-pair rosette leaves before culture (time 0), at 10min to 12h after detachment Col-0, coi1-2 and sdg8-2 seedlings.
- 课题设计：To reveal the early molecular events upon wounding during DNRR from Arabidopsis leaf explants cultured on B5 medium without exogenous hormones, we carried out an RNA-seq experiment using leaf explants before culture (i.e. time 0), at 10 and 30 min after detachment, and at 1 to 12 h after detachment.
安装包
代码语言：javascript复制# BiocManager::install("Mfuzz")
library(Mfuzz)
library(limma)
library(clusterProfiler)
# BiocManager::install("org.At.tair.db") ###一定要注意物种，这个数据集是拟南芥
library(org.At.tair.db)
library(ggplot2)
library(ggstatsplot)
library(tidyverse)
下载数据
代码语言：javascript复制d='GSE120418_RAW/'
fs = list.files(d,pattern = '_Col_') ##这里我就想看看对照组的情况，就把对照组的GSM样本提取出来
fs
library(data.table)
tmp=fread(file.path(d,fs[1])) ##先读一个文件看看
colnames(tmp)
dim(tmp)
lapply(fs, function(f){
print(dim(fread(file.path(d,f ))))
})
gene.count = do.call(cbind,
lapply(fs, function(f){
ceiling(fread(file.path(d,f ),data.table = F)[,5])#Ceiling函数返回最接近输入值但大于输入值的值
}))
head(gene.count)
rownames(gene.count) = tmp$gene_id
library(stringr)
colnames(gene.count) = str_split(fs,'_',simplify = T)[,4] ##提出时间序列
gene.count[1:4,1:4]
dat=gene.count ##便于后续差异分析
temp = str_split(fs,'_0_',simplify = T)[,2]
temp1 <- paste0("Col_",str_split(temp,"_r",simplify = T)[,1])
temp1
colnames(dat) <- temp1
dat[1:4,1:4]
dat <- log2(edgeR::cpm(dat) 1)
dat[1:4,1:4]
dim(dat)
library(limma)
avereps_df <- t(limma::avereps( t(dat) , ID = colnames(dat)))##对相同时间序列的表达值取平均
avereps_df[1:4,1:4]
colnames(avereps_df)
save(avereps_df,file = 'avereps_df.Rdata')
数据过滤
代码语言：javascript复制load(file = 'avereps_df.Rdata')
avereps_df[1:4,1:4]
colnames(avereps_df)
avereps_df = avereps_df[,c( "time0", "10min","30min" ,
"1h" , "2h", "4h" , "8h" , "12h" )]
## 2.1 Filtering----
# 去除表达量太低或者在不同时间点间变化太小的基因等步骤
# Mfuzz聚类时要求是一个ExpressionSet类型的对象，所以需要先用表达量构建这样一个对象。
eset <- new("ExpressionSet",exprs = avereps_df)
# 根据标准差去除样本间差异太小的基因
eset <- filter.std(eset,min.std=0)
# 10818 genes excluded. ，不同的数据集去除的基因数量不一样
eset
## 2.2 Standardisation----
# 聚类时需要用一个数值来表征不同基因间的距离，Mfuzz中采用的是欧式距离，
# 由于普通欧式距离的定义没有考虑不同维度间量纲的不同，所以需要先进行标准化
eset <- standardise(eset)
## 2.3 Setting of parameters for FCM clustering----
# Mfuzz中的聚类算法需要提供两个参数，
# 第一个参数为希望最终得到的聚类的个数，这个参数由我们直接指定
# 第二个参数称之为fuzzifier值，用小写字母m表示，可以通过函数评估一个最佳取值
c <- 9
m <- mestimate(eset) # 评估出最佳的m值
cl <- mfuzz(eset, c = c, m = m) # 聚类
## 2.4 glimpse results----
# 在cl这个对象中就保存了聚类的完整结果，对于这个对象的常见操作如下
cl$size # 查看每个cluster中的基因个数
cl$cluster[cl$cluster == 1] # 提取某个cluster下的基因
## cluster cores
# membership values can also indicate the similarity of vectors to each other.
eset
cl.thres <- acore(eset,cl,min.acore=0.5) ## extracts genes forming the alpha cores of soft clusters
head(cl.thres[[1]])
table(cl$cluster)
unlist(lapply(cl.thres, nrow))
查看集群之间的耦合
代码语言：javascript复制# coupling between clusters
overlap.cl <- overlap(cl)
pdf('mfuzz_overlap_plot.pdf',height = 4,width = 5)
p.overlaps <- overlap.plot(cl, over = overlap.cl, thres = 0.05)
p.overlaps
dev.off()
可视化
代码语言：javascript复制## 2.5 visualise----
library(RColorBrewer)
color.2 <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000)
pdf('mfuzz_clusters_plot.pdf',height = 7,width = 12)
mfuzz.plot(eset,cl,mfrow=c(3,3),
new.window= FALSE,
time.labels= colnames(eset) ,
colo = color.2)
dev.off()
pdf('mfuzz_clusters_plot01.pdf',height = 7,width = 12)
mfuzz.plot2(eset, cl, mfrow = c(3, 3),
centre = T, x11 = F, centre.lwd = 0.2)
批量输出
代码语言：javascript复制# 3. 批量输出聚类所含的基因 ----------------------------------------------------------
getwd()
dir.create(path = "mfuzzGenes", recursive = T)
for (i in 1:9) {
potname = names(cl$cluster[unname(cl$cluster) == i])
write.csv(cl[[4]][potname, i], paste0("mfuzzGenes", "/mfuzz_", i, ".csv"))
}
然而
这个包只能进行聚类，是找不了有处理对照组的差异基因的
所以得引进它⬇
Package ‘maSigPro’
安装
代码语言：javascript复制# BiocManager::install('maSigPro')
library(maSigPro)
读取原始文件
代码语言：javascript复制# 1.读取原始数据 ------------------------------------------------------------------
d='GSE120418_RAW/'
fs = list.files(d) ##注意这时候我是把整个数据集都读进来了
fs
library(data.table)
tmp=fread(file.path(d,fs[1])) ##先读一个文件看看
colnames(tmp)
dim(tmp)
lapply(fs, function(f){
print(dim(fread(file.path(d,f ))))
})
处理count矩阵
代码语言：javascript复制# 2.处理count矩阵 ---------------------------------------------------------------
gene.count = do.call(cbind,
lapply(fs, function(f){
ceiling(fread(file.path(d,f ),data.table = F)[,5])
#Ceiling函数返回最接近输入值但大于输入值的值，和floor是相对的
}))
rownames(gene.count) = tmp$gene_id
library(stringr)
colnames(gene.count) = str_split(fs,'_',simplify = T)[,4] ##提出时间序列
gene.count[1:4,1:4]
dat=gene.count ##便于后续差异分析
temp <- str_split(fs,"_r",simplify = T)[,1]
temp1 <- gsub('^GSM[0-9]*_','',temp)
temp1
colnames(dat) <- temp1
dat[1:4,1:4]
dat <- log2(edgeR::cpm(dat) 1)
dat[1:4,1:4]
dim(dat)
###按照官方文档改一下列名
colnames(dat)=ifelse(seq(1,ncol(dat),1) %%2 !=0,
gsub("_S[0-9]*","_1",colnames(dat)),
gsub("_S[0-9]*","_2",colnames(dat))
)
colnames(dat)
保存矩阵信息
代码语言：javascript复制save(gene.count,dat,file = "readR_count_matrix.Rdata")
gene.count[1:4,1:4]
dat[1:4,1:4]
colnames(dat)
建立回归模型
代码语言：javascript复制dat[1:4,1:4]
library(stringr)
##3.1 构建design矩阵
edesign <- read.xlsx("edesign.xlsx") ###把10分钟和30分钟写成了1/6和1/2,我是手动用excel创建的
rownames(edesign) <- colnames(dat)
head(edesign)
colnames(edesign) <- str_to_title(colnames(edesign))
head(edesign)
#### CREATE EXPERIMENTAL DESIGN（官方文档是直接用R创建的，怎么顺手怎么来）
# Time <- rep(c(rep(c(1:3), each = 3)), 4)
# Replicates <- rep(c(1:12), each = 3)
# Control <- c(rep(1, 9), rep(0, 27))
# Treat1 <- c(rep(0, 9), rep(1, 9), rep(0, 18))
# Treat2 <- c(rep(0, 18), rep(1, 9), rep(0,9))
# Treat3 <- c(rep(0, 27), rep(1, 9))
# edesign <- cbind(Time, Replicates, Control, Treat1, Treat2, Treat3)
# rownames(edesign) <- paste("Array", c(1:36), sep = "")
##3.2 生成回归矩阵(makeDesignMatrix)
table(edesign$Time)
design <- make.design.matrix(edesign, degree = 6) ###这里degree的值我在文档里没有找到确切的依据，因为我的分组是7，定了6
design$groups.vector
关于degree参数，官方说明是这样的
代码语言：javascript复制This example has three time points, so we can consider up to a quadratic regression model (degree = 2). **Larger number of time points would potentially allow a higher polynomial degree.**大家见仁见智
##3.3 寻找重要基因(p.vector)
fit<-p.vector(dat,#标准化的表达矩阵
design,#实验设计的矩阵make.design.matrix生成
Q=0.01,#显著性水平
MT.adjust="BH",
min.obs=22
#最低表达样本数不应小于(degree 1)xGroups 1
)
save(fit,edesign,design,file = "fit.Rdata")
差异基因分析
代码语言：javascript复制load("fit.Rdata")
fit$i ##显著性基因的数量
fit$SELEC ##差异表达矩阵
fit$i # returns the number of significant genes
fit$alfa # gives p-value at the Q false discovery control level
# null??
fit$SELEC # is a matrix with the significant genes and their expression values
开始我是这样写的：
代码语言：javascript复制##3.4 寻找显著性差异(T.fit)
tc.tstep <- T.fit(data = fit , alfa = 0.01)
tc.tstep <- T.fit(data =fit , alfa = 0.05)
然后就报错了
代码语言：javascript复制Error in terms.formula(formula, data = data) :
公式里有'.'，而没有'data'这一参数
上网查也没找到答案，然后
代码语言：javascript复制##3.4 寻找显著性差异(T.fit)
tstep <- T.fit(fit, # p.vector结果
step.method = "two.ways.forward",
alfa = 0.01) # 在逐步回归中用于变量选择的显著性水平
就OK啦
那就是step.method参数的选择，因为默认是“backward”,但为啥我的数据没法用这个方法呢？如果有统计学方向的同学，欢迎评论区答疑呀~~~~
官方说明在这里：
代码语言：javascript复制The step.method can be ”backward” or ”forward” indicating whether the step procedure starts from the model with all or none variables. Use method ”two.ways.backward” and ”two.ways.forward” to allow variables to both get in and out. At each regression step the p-value of each variable is computed and variables get in/out the model when this p-value is lower or higher than the given cut-off value alfa.
# 3.5 获取显著性基因列表
sigs <- get.siggenes(tstep, # T.fit结果
rsq = 0.80, # 逐步回归中的R-squared截至值
vars ="groups")
# sigs1 <- get.siggenes(tstep, # T.fit结果
# rsq = 0.6, # 逐步回归中的R-squared截至值
# vars ="each") ##会给出时间点和实验条件的所有组合对应差异基因列表
save(sigs,file = 'maSigPro_deg_results.Rdata')
可视化
代码语言：javascript复制# 3.6 韦恩图
dev.off()
suma2Venn(sigs$summary[, c(2:3)])
suma2Venn(sigs$summary[, c(1:3)])
# colnames(sigs1$summary)
# suma2Venn(sigs1$summary[, c(4:6)])
# suma2Venn(sigs1$summary[, c(8:10)])
# 3.7 聚类
sigs$sig.genes$CoilvsControl$g
see.genes(sigs$sig.genes$CoilvsControl, # 差异基因表达矩阵
show.fit = T, # 是否显示回归拟合线(虚线)
dis =design$dis, # 回归设计矩阵
cluster.method="hclust" , # 聚类方法
cluster.data = 1,
k = 9) # 聚类数目
#换一种聚类方式
# install.packages("mclust")
# 加载
library(mclust)
see.genes(sigs$sig.genes$CoilvsControl, # 差异基因表达矩阵
show.fit = T, # 是否显示回归拟合线(虚线)
dis =design$dis, # 回归设计矩阵
cluster.method="Mclust" , # 聚类方法
cluster.data = 1,
k.mclust=TRUE) # 聚类数目
我作出来的图太丑了，不便展示，
