查看: 412|回复: 20

[R语言] 再分享一种分析基因表达趋势的方法

  [复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
461
注册时间
2020.6.16
在线时间
304 小时

发表于 6 天前 | 显示全部楼层 |阅读模式
趋势分析(Series Test of Cluster)是一种挖掘相同表达模式基因的重要方法,之前在《师兄,师兄,趋势分析怎么做?》一文中已经介绍过,趋势分析本质上是将表达趋势相同的基因进行归类,如下,详情可阅读原文。


下面再为大家介绍一种采用模糊聚类算法的趋势分析方法,聚类效果如下,与上文提到基于STEM的方法(一般3~5个时间点)相比,可支持更多时间节点的数据,当然聚类效果也不可避免有一点“模糊”。


以下是一些高分文章的应用范例:

Cell, 2015)

Nature, 2019)

下面就一起看看如何分析吧!

1. 读入数据

[AppleScript] 纯文本查看 复制代码
#安装 Mfuzz包;
#if (!requireNamespace("BiocManager", quietly = TRUE))
#install.packages("BiocManager")
BiocManager::install("Mfuzz")
#载入Mfuzz包;
library(Mfuzz)

#设置工作目录;
setwd("C:/Users/86136/Desktop/Mfuzz")
#读入测试数据;
#输入数据一般为差异基因或蛋白的表达矩阵;
#推荐使用归一化(Normalisation)后的数值,比如FPKM值;
df <- read.csv("testdata.csv",header = T)
head(df)



注意
Normalisation是为了使不同的样本具有可比性;而Mfuzz中的standardisation是为了使不同的基因或蛋白质具有可比性。

[AppleScript] 纯文本查看 复制代码
#这里的范例数据每组有3个重复,计算组内均值;
T0 <- apply(df[,2:4], 1, mean)
T1 <- apply(df[,5:7], 1, mean)
T2 <- apply(df[,8:10], 1, mean)
T3 <- apply(df[,11:13], 1, mean)
#生成新的表达量矩阵;
df1 <- data.frame(T0,T1,T2,T3)
row.names(df1) <- df$gene
#转成矩阵;
mat <- as.matrix(df1)
#查看前6行;
head(mat)



[AppleScript] 纯文本查看 复制代码
#创建用于Mfuzz的对象;
dt <- new("ExpressionSet",exprs = mat)
dim(dt)
#Features Samples
#9345 4
#过滤缺失值超过25%的基因;
dt.r <- filter.NA(dt, thres=0.25)
#以均值的方式填充缺失值;
dt.f <- fill.NA(dt.r,mode="mean")
#过滤低表达或表达量变化不大的基因;
#由于是差异基因,这里不做过滤;
tmp <- filter.std(dt.f,min.std=0)

#对数据进行标准化;
dt.s <- standardise(tmp)
#查看标准化后的数据;
df.s <- dt.s@assayData$exprs
head(df.s)



2.模糊聚类

这里的模糊聚类(Soft clustering)是指Mfuzz包用到的fuzzy c-means clustering算法,思路其实与常规的hard clustering (如k-means)类似。

[AppleScript] 纯文本查看 复制代码
#使用mestimate函数估计m值;
m1 <- mestimate(dt.s)
#执行模糊聚类;
set.seed(007)
cl <- mfuzz(dt.s,c=8,m=m1)


这里的参数c即centers,指期望得到的cluster个数(比较主观,多尝试),参数m指fuzzifier m值,由上述函数估计得出。

3.图表绘制

[AppleScript] 纯文本查看 复制代码
#对聚类结果进行可视化;
mfuzz.plot(dt.s,cl,mfrow=c(2,4),
new.window= FALSE,
time.labels=colnames(dt.s))



[AppleScript] 纯文本查看 复制代码
#自定义配色;
library(RColorBrewer)
mycol <- c("cyan","yellow","orangered")
mycolor <- colorRampPalette(mycol)(100)
#使用自定义配色批量绘图;
mfuzz.plot(dt.s,cl,mfrow=c(2,4),
new.window= FALSE,
time.labels=colnames(dt.s),
colo = mycolor)



[AppleScript] 纯文本查看 复制代码
#查看每个cluster的基因数量;
cl$size
#[1] 910 1276 778 994 2032 1408 918 1025
#查看前6个基因所属簇;
head(cl$cluster)
#Ldhb Hsp90aa1 Uchl1 Zar1 Prdx2 Bcl2l10
#5 5 5 5 5 2
#查看基因在不同cluster的membership值;
head(cl$membership)



这里得membership值有点类似相关性系数,利用基因的membership值可对基因进行分簇(例如Ldhb基因的最大membership值为0.39226651,对应cluster5),此外,作图时也可用于建立渐变色映射。

[AppleScript] 纯文本查看 复制代码
#整合标准化表达量和聚类结果;
gene_cluster <- data.frame(df.s,cluster=cl$cluster)
#预览结果;
head(gene_cluster)
#提取指定簇下的基因;
cl_5 <- subset(gene_cluster,cluster==5)
#预览结果;
head(cl_5)



[AppleScript] 纯文本查看 复制代码
#也可以绘制单个cluster(如cluster1)的中心线折线图;
plot(cl$centers[1,],type="l",main="cluster1",
col="red",ylab="Expression")



[AppleScript] 纯文本查看 复制代码
#导出数据;
write.csv(gene_cluster, 'gene_cluster.csv')


趋势分析完成后,接下来可对感兴趣的Cluster进一步做基因功能分析,如GO、KEGG富集分析等,直到挖掘到感兴趣的基因。好啦,今天就分享到这里啦~

▼参考资料▼
Walther D M, Kasturi P, Zheng M, et al. Widespread proteome remodeling and aggregation in aging C. elegans[J]. Cell, 2015, 161(4): 919-932.
Zhou W, Sailani M R, Contrepois K, et al. Longitudinal multi-omics of host–microbe dynamics in prediabetes[J]. Nature, 2019, 569(7758): 663-671.
http://mfuzz.sysbiolab.eu/



本文作者:基迪奥-莫北

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
31
注册时间
2016.1.8
在线时间
556 小时

发表于 6 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 6 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2021.6.4
在线时间
5 小时

发表于 6 天前 | 显示全部楼层
求下示例数据,谢谢!
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
31
注册时间
2016.1.8
在线时间
556 小时

发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
148 小时

灌水之王


发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
148 小时

灌水之王


发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
148 小时

灌水之王


发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
148 小时

灌水之王


发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2017.11.3
在线时间
234 小时

发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2021.6.22
在线时间
16 小时

发表于 5 天前 | 显示全部楼层
学习
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
81 小时

发表于 5 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2016.1.13
在线时间
270 小时

发表于 4 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 4 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 4 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2018.7.5
在线时间
20 小时

发表于 4 天前 | 显示全部楼层
漠北老师厉害了,这个图真好看
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2022.1.10
在线时间
2 小时

发表于 4 天前 | 显示全部楼层
哈哈哈哈哈哈哈哈哈
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2022.1.10
在线时间
2 小时

发表于 4 天前 | 显示全部楼层
新的一天加油!
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表