转录组数据的时间序列分析，你学会了吗

生信菜鸟团

发布于 2022-10-31 17:44:21

2.4K0

发布于 2022-10-31 17:44:21

文章被收录于专栏：生信菜鸟团生信菜鸟团

上周的公众号处理了不同时间序列的数据集，但因为是内置的数据集，很多分析流程都已经被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.

安装包

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

下载数据

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

数据过滤

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

查看集群之间的耦合

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

可视化

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

批量输出

# 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’

安装

# BiocManager::install('maSigPro')
library(maSigPro)

读取原始文件

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

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

保存矩阵信息

save(gene.count,dat,file = "readR_count_matrix.Rdata")

gene.count[1:4,1:4]
dat[1:4,1:4]
colnames(dat)

建立回归模型

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参数，官方说明是这样的

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

差异基因分析

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

开始我是这样写的：

##3.4 寻找显著性差异(T.fit) 
tc.tstep <- T.fit(data = fit , alfa = 0.01)
tc.tstep <- T.fit(data =fit , alfa = 0.05)

然后就报错了

Error in terms.formula(formula, data = data) : 
  公式里有'.'，而没有'data'这一参数

上网查也没找到答案，然后

##3.4 寻找显著性差异(T.fit) 
tstep <- T.fit(fit, # p.vector结果
               step.method = "two.ways.forward", 
               alfa = 0.01) # 在逐步回归中用于变量选择的显著性水平

就OK啦

那就是step.method参数的选择，因为默认是“backward”,但为啥我的数据没法用这个方法呢？如果有统计学方向的同学，欢迎评论区答疑呀~~~~

官方说明在这里：

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

可视化

# 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) # 聚类数目

我作出来的图太丑了，不便展示，😳

这个包还有很多函数值得好好探索，感兴趣的小伙伴不妨拿一个下午的时间看看，小编也是抛砖引玉，希望大家可以利用这个包写出更精彩的文章，嘿嘿

叠个buff

加载包包

rm(list = ls())
options(stringsAsFactors = F)
# 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)

加载数据

load('maSigPro_deg_results.Rdata')
G1 <- sigs$sig.genes$CoilvsControl$sig.pvalues[,1:2]
G2 <- sigs$sig.genes$Control$sig.pvalues[,1:2]
G3 <- sigs$sig.genes$Sdg8vsControl$sig.pvalues[,1:2]

load(file = "readR_count_matrix.Rdata")
dat[1:4,1:4]

查看差异基因列表

# 1. all deg gene----
g1 <- rownames(subset(G1, G1$`p-value`< 0.05))
g2 <- rownames(subset(G2, G2$`p-value` < 0.05))
g3 <- rownames(subset(G3, G3$`p-value` < 0.05))
deg.gene <- unique(c(g1,g2,g3))

# dat <- log2(edgeR::cpm(gene.count)+1)
dim(dat)
dat.deg <- dat[deg.gene,]
dim(dat.deg)

数据过滤

# 2. Mfuzz----
# limma::avereps
# 时间节点表达谱
library(stringr)
colnames(dat.deg) <- gsub('.{2}$', '',colnames(dat.deg))  
##https://www.delftstack.com/zh/howto/r/remove-last-character-in-r/
table(colnames(dat.deg))
group <- colnames(dat.deg)
DEGs_exp_averp <- t(limma::avereps(t(dat.deg), ID=group))
head(DEGs_exp_averp)
# Samle_Gene <- sample(rownames(DEGs_exp_averp),size = 10)  ##随机抽取10个基因
# DEGs_exp_averp <- DEGs_exp_averp[Samle_Gene,]
# head(DEGs_exp_averp)

## 2.1 Filtering----
# 去除表达量太低或者在不同时间点间变化太小的基因等步骤
# Mfuzz聚类时要求是一个ExpressionSet类型的对象，所以需要先用表达量构建这样一个对象。
eset <- new("ExpressionSet",exprs = DEGs_exp_averp)
# 根据标准差去除样本间差异太小的基因
eset <- filter.std(eset,min.std=0)

## 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.
cl.thres <- acore(eset,cl,min.acore=0.5)  ## a posteriori
head(cl.thres[[1]])
table(cl$cluster)
unlist(lapply(cl.thres, nrow))

dir.create('output')
lapply(seq_along(cl.thres), function(i){
  # i = 1
  pro <- paste0('cluster',i)
  print(pro)
  cluster.gene <- cl.thres[[i]]$NAME
  write.table(cluster.gene, quote = F,row.names = F,col.names = F,
              file = file.path('output/',
                               paste0(pro,'_deg_mfuzz.9-threshold_cluster.txt')))
})

可视化

# coupling between clusters
overlap.cl <- overlap(cl)
pdf('output/deg_mfuzz_overlap_plot.pdf',height = 4,width = 5)
p.overlaps <- overlap.plot(cl, overlap = overlap.cl, thres = 0.05)

dev.off()

library(RColorBrewer)
color.2 <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000)
pdf('output/deg_mfuzz_clusters_plot.pdf',height = 7,width = 12)
mfuzz.plot(eset,cl,mfrow=c(3,3),
           new.window= FALSE,
           time.labels=colnames(DEGs_exp_averp),
           colo = color.2)
dev.off()

# 3. genes unregulated in response to time----
# cluster-2 genes were upregulated within R0 then downregulated
# cluster-5 genes were upregulated at around 12 h, then downregulated
# cluster-8 genes were upregulated at around R0 to 1w, then downregulated
names(cl.thres) <- paste0('cluster',c(1:9))
# choose.cluster <- c(2,5,8)

lapply(seq_along(cl.thres), function(i){
   # i = 3
  pro <- names(cl.thres)[i]
  print(pro)
  choose.genes <- cl.thres[[i]]$NAME
  choose_matrix <- DEGs_exp_averp[choose.genes,]
  choose_matrix <- t(scale(t(choose_matrix)))
  choose_matrix[choose_matrix>2] =2
  choose_matrix[choose_matrix< -2] = -2
  colD=data.frame(Groups=colnames(choose_matrix))
  rownames(colD)=colnames(choose_matrix)
  palette = RColorBrewer::brewer.pal(12,"Set3")
  names(palette) <- colnames(choose_matrix)

  pdf(file = file.path('output/',paste0(pro, '_threshold_gene.pdf')),
      width = (ncol(choose_matrix)*0.3+2.2) *2 ,
      height = 550/100*3.5*nrow(choose_matrix)/100
      )
  pheatmap::pheatmap(choose_matrix,
                     cluster_cols = F,
                     annotation_col = colD,
                     annotation_colors = list(Groups = palette),
                     fontsize = 12,
                     border_color = NA,
                     main = paste0(pro, '_threshold_gene'))
  dev.off()
})

你有几个cluster，就会出几个图哦

基因注释

# 4. annotation----
## Convert gene ID into entrez genes
#可参考CSDN文章：「Bioconductor」再次提醒，研究植物的不要轻易相信你用的注释包
tmp <- toTable( org.At.tairENTREZID)  ###记住你的物种啊！！！
head(tmp)
colnames(tmp) <- c("NAME","entrez")
gcSample <- lapply(seq_along(cl.thres), function(i){
   # i = 1
  pro <- names(cl.thres)[i]
  print(pro)
  gene.cluster <- cl.thres[[i]]
  head(gene.cluster)
  head(tmp)
  s2e <- merge(gene.cluster,tmp,by='NAME') %>%
    dplyr::select(NAME) %>% unlist(use.names = FALSE)
  head(s2e)
  return(s2e)
})

names(gcSample) <- names(cl.thres)
unlist(lapply(gcSample, length))

通路富集

# kk <- enrichKEGG(gcSample$cluster1, organism="ath",keyType = "kegg") ,
#                  pvalueCutoff=0.05, 
#                  pAdjustMethod="BH",qvalueCutoff=0.1)
# dotplot(kk)

cm_KEGG <- compareCluster(gcSample,
                     fun="enrichKEGG",
                     organism="ath",
                     pvalueCutoff=0.05)

p1 <- dotplot(cm_KEGG) +
  scale_color_continuous(low="#fe6f5e",high = "#0095b6") + ## colors for p.adjust
  theme_bw()+
  labs(title = 'KEGG Pathway enrichment compared') + xlab(label = '') +
  theme(axis.text  = element_text(size = 15,colour = 'black',face = 'bold'),
        axis.title = element_text(size = 15),
        plot.title = element_text(size=15,face = 'bold',hjust = 0.5))
p1+ aes(shape = GeneRatio > 0.3)
ggsave('output/kegg_compareCluster_threshold_gene.pdf',height = 7,width = 12)

保存数据

save(cl.thres, gcSample,file = 'output/step3_mfuzz_results.Rdata')

听说还有一个包也可以做转录组数据的时间序列分析，留个悬念吧，以后有空追更一下

Package'ImpulseDE2'

参考文档在此：https://rdrr.io/bioc/ImpulseDE2/

小结

这两期介绍这些包包的主要目的，颇有些“工欲善其事，必先利其器”的意味。

Anyway，希望大家都能在生信菜鸟团这里找到问题的答案，当然如果没有，也欢迎评论区或后台留言，说不定下一期就是为你量身定做的呢

希望我们共同进步，早日上树！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2022-10-18，如有侵权请联系 cloudcommunity@tencent.com 删除

model

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度