数据挖掘 机器学习 人工智能 算法 RNA 34. SCI 文章基因蛋白表达时间趋势分析及划分聚类群 (Mfuzz)

转录组生信分析教程

桓峰基因公众号推出转录组分析教程，有需要生信的老师可以联系我们！转录分析教程整理如下：

RNA 1. 基因表达那些事--基于 GEO

RNA 2. SCI文章中基于GEO的差异表达基因之 limma

RNA 3. SCI 文章中基于T CGA 差异表达基因之 DESeq2

RNA 4. SCI 文章中基于TCGA 差异表达之 edgeR

RNA 5. SCI 文章中差异基因表达之 MA 图

RNA 6. 差异基因表达之-- 火山图 (volcano)

RNA 7. SCI 文章中的基因表达——主成分分析 (PCA)

RNA 8. SCI文章中差异基因表达--热图 (heatmap)

RNA 9. SCI 文章中基因表达之 GO 注释

RNA 10. SCI 文章中基因表达富集之--KEGG

RNA 11. SCI 文章中基因表达富集之 GSEA

RNA 12. SCI 文章中肿瘤免疫浸润计算方法之 CIBERSORT

RNA 13. SCI 文章中差异表达基因之 WGCNA

RNA 14. SCI 文章中差异表达基因之 蛋白互作网络 (PPI)

RNA 15. SCI 文章中的融合基因之 FusionGDB2

RNA 16. SCI 文章中的融合基因之可视化

RNA 17. SCI 文章中的筛选 Hub 基因 (Hub genes)

RNA 18. SCI 文章中基因集变异分析 GSVA

RNA 19. SCI 文章中无监督聚类法 （ConsensusClusterPlus）

RNA 20. SCI 文章中单样本免疫浸润分析 （ssGSEA）

RNA 21. SCI 文章中单基因富集分析

RNA 22. SCI 文章中基于表达估计恶性肿瘤组织的基质细胞和免疫细胞(ESTIMATE)

RNA 23. SCI文章中表达基因模型的风险因子关联图(ggrisk)

RNA 24. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER)

RNA 25. SCI文章中估计组织浸润免疫细胞和基质细胞群的群体丰度(MCP-counter)

RNA 26. SCI文章中基于转录组数据的基因调控网络推断 (GENIE3)

RNA 27 SCI文章中转录因子结合motif富集到调控网络 (RcisTarget)

RNA 28 SCI 文章中基于RNA-seq数据反褶积揭示肿瘤免疫结构的分子和药理学 (quanTIseq)

RNA 29. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER2.0)

RNA 30. SCI文章中基于TCGA和GTEx数据挖掘神器（GEPIA2）

RNA 31. SCI文章临床蛋白质组肿瘤在线数据挖掘神器（CPTAC）

RNA 32. SCI文章临床多组学肿瘤在线数据挖掘神器（UALCAN）

RNA 33. SCI文章肿瘤在线数据挖掘神器（cBioportal）

这期介绍基因表达，或蛋白表达时间趋势以及聚类，这项分析同样非常常见，本期同时复现 CELL 期刊文章中的图，利用一个简单的表达矩阵即可完成绘制，所以好多老师觉得文章空洞，没内涵，可以多角度的分析基因表达，或蛋白表达的时间趋势，还是蛮有意思，下面看细节吧！

简介

在微阵列数据分析中，经常使用聚类技术。这些方法大多是基于数据的硬聚类，其中一个基因(或样本)被精确地分配到一个聚类。然而，硬集群遭受一些缺点，如对噪声和信息丢失的敏感性。相比之下，软聚类方法可以将基因分配给几个群。可以克服传统硬聚类技术的缺点，并提供进一步的优势。因此，Mfuzz的R包实现了微阵列数据分析的聚类工具。

软件包安装

if (!require("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("Mfuzz")

BiocManager::install("marray")

数据读取

library(Mfuzz)

data(yeast)

head(yeast@assayData$exprs[, 1:5])

## cdc28_0 cdc28_10 cdc28_20 cdc28_30 cdc28_40

## YDR132C 0.19 0.30 -0.29 0.29 -0.31

## YMR012W -0.15 -0.15 -0.04 -0.28 -0.39

## YLR214W 0.38 0.30 -0.68 -0.52 -0.43

## YLR116W 0.17 0.06 -0.21 0.19 0.33

## YDR203W 0.85 -0.10 -0.56 -0.31 -0.43

## YEL059C-A 0.45 0.20 0.06 0.10 -0.21

实例操作

软件包自带例子

第一步：我们排除丢失超过25％的测量值的基因。注意缺失值应在基因表达矩阵中用NA表示。

# Data pre-processing

yeastF <- filter.NA(yeast)

## 49 genes excluded.

yeastF <- fill.NA(yeastF)

标准化

yeastF <- standardise(yeastF)

软聚类和可视化

# Soft clustering and visualisation

cl <- mfuzz(yeastF, c = 20, m = 1.25)

mfuzz.plot(yeastF, cl = cl, mfrow = c(4, 5))

绘制颜色条

mfuzzColorBar(col = "fancy", main = "Membership", cex.main = 1)

黄色或绿色的线条对应于成员价值较低的基因；红色和紫色线对应具有高隶属价值的基因。

mfuzz.plot2(yeastF, cl = cl, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, mfrow = c(4,

5))

查看每个cluster中的基因个数，提取某个cluster下的基因并查看基因和cluster之间的membership

cl$size

## [1] 135 128 208 157 151 147 135 171 140 148 121 151 166 160 172 179 156 130 68

## [20] 128

cluster_gene <- as.data.frame(cl$cluster)

membership <- cl$membership

聚类稳定性

FCM参数m的变化还允许研究聚类的稳定性。

这里，我们将稳定聚类定义为仅随参数m的变化而显示出结构上微小变化的聚类。稳定的聚类通常比其他聚类紧凑。与之形成鲜明对比的是，如果m增加，弱聚类内部结构失去或消失了 。

cl1 <- mfuzz(yeastF, c = 16, m = 1.35) #m增加。

mfuzz.plot(yeastF, cl = cl1, mfrow = c(4, 4), time.labels = seq(0, 160, 10))

全局聚类结构

软聚类的一个有趣的特征是聚类之间的重叠或耦合。偶联表明两个簇共有多少个基因。耦合度低的集群总体上表现出明显的差异模式。如果耦合较大，则聚类模式更相似。因此，耦合为成对的聚类定义相似性度量。

O <- overlap(cl)

Ptmp <- overlap.plot(cl, over = O, thres = 0.05)

由于定义了聚类之间的关系，因此可以分析通过软聚类获得的全局聚类结构。与分层聚类相似，全局聚类结构可以按聚类数确定的不同分辨率进行检查C。

cl2 <- mfuzz(yeastF, c = 10, m = 1.25)

mfuzz.plot(yeastF, cl = cl2, mfrow = c(3, 4))

对于较小的C，仅获得数据中存在的主要聚类。

O1 <- overlap(cl2)

overlap.plot(cl2,over=O1,P=Ptmp,thres=0.05)

## PC1 PC2 PC3 PC4 PC5

## cdc28_0 -0.07145877 0.15989427 -0.56689231 -0.100418973 -0.09290026

## cdc28_10 -0.21737594 0.24894738 -0.19707752 -0.351395955 -0.13959866

## cdc28_20 -0.21574121 0.32785374 0.24854966 -0.004234176 -0.19422055

## cdc28_30 -0.32132410 0.11799740 0.19446405 0.204994971 -0.20601231

## cdc28_40 -0.24691558 -0.20801026 0.36161114 0.139134032 -0.05127488

## cdc28_50 -0.23778553 -0.34762911 -0.01038037 0.193852214 -0.03796400

## cdc28_60 -0.20049605 -0.37982795 0.15314325 -0.083137917 0.13182695

## cdc28_70 -0.15958646 -0.41183114 -0.04027852 -0.120196236 0.21016136

## cdc28_80 0.06431964 -0.37707227 -0.37256435 -0.080180386 0.09474066

## cdc28_90 0.23384862 0.19212590 0.30044761 -0.089656247 0.41540001

## cdc28_100 0.19869316 0.14378420 -0.03818160 0.352203488 0.57637272

## cdc28_110 0.29491424 0.10937444 -0.02967616 0.400575137 -0.16504514

## cdc28_120 0.18681339 -0.16733554 -0.21213307 0.494426873 -0.31485432

## cdc28_130 0.29303823 -0.15463362 0.23357183 0.008812839 -0.36773211

## cdc28_140 0.31318385 -0.20057985 0.15888131 -0.172334380 -0.06216872

## cdc28_150 0.29509320 -0.06744186 0.16339264 -0.367004607 -0.19679768

## cdc28_160 0.35137129 -0.06806756 -0.03669656 -0.195443200 -0.09038667

如果c增加，则会出现具有不同模式的子集群，因为它们共享整个表达模式，所以从主集群派生的子集群通常会紧密耦合。

cl3 <- mfuzz(yeastF, c = 25, m = 1.25)

mfuzz.plot(yeastF, cl = cl3, mfrow = c(5, 5))

最后，软聚类产生空的聚类，以进一步增加c。

O2 <- overlap(cl3)

overlap.plot(cl3, over = O2, P = Ptmp, thres = 0.05)

## PC1 PC2 PC3 PC4 PC5

## cdc28_0 -0.07145877 0.15989427 -0.56689231 -0.100418973 -0.09290026

## cdc28_10 -0.21737594 0.24894738 -0.19707752 -0.351395955 -0.13959866

## cdc28_20 -0.21574121 0.32785374 0.24854966 -0.004234176 -0.19422055

## cdc28_30 -0.32132410 0.11799740 0.19446405 0.204994971 -0.20601231

## cdc28_40 -0.24691558 -0.20801026 0.36161114 0.139134032 -0.05127488

## cdc28_50 -0.23778553 -0.34762911 -0.01038037 0.193852214 -0.03796400

## cdc28_60 -0.20049605 -0.37982795 0.15314325 -0.083137917 0.13182695

## cdc28_70 -0.15958646 -0.41183114 -0.04027852 -0.120196236 0.21016136

## cdc28_80 0.06431964 -0.37707227 -0.37256435 -0.080180386 0.09474066

## cdc28_90 0.23384862 0.19212590 0.30044761 -0.089656247 0.41540001

## cdc28_100 0.19869316 0.14378420 -0.03818160 0.352203488 0.57637272

## cdc28_110 0.29491424 0.10937444 -0.02967616 0.400575137 -0.16504514

## cdc28_120 0.18681339 -0.16733554 -0.21213307 0.494426873 -0.31485432

## cdc28_130 0.29303823 -0.15463362 0.23357183 0.008812839 -0.36773211

## cdc28_140 0.31318385 -0.20057985 0.15888131 -0.172334380 -0.06216872

## cdc28_150 0.29509320 -0.06744186 0.16339264 -0.367004607 -0.19679768

## cdc28_160 0.35137129 -0.06806756 -0.03669656 -0.195443200 -0.09038667

文章复习

Gao等（2017）基于蛋白质谱的方法，研究了小鼠胚胎着床前发育过程中的蛋白质组。

本文共涉及了6个发育阶段，受精卵（Zygote）、二分胚（2-cell）、四分胚（4-cell）、八分胚（8-cell）、桑葚胚（Morula）和囊胚（Blastocyst）。为了将蛋白质功能与胚胎发育相结合，作者首先表征了蛋白质丰度与胚胎发育阶段的时间关系，根据所有蛋白质在每个阶段的丰度信息，通过Mfuzz包对这些蛋白质执行了时间序列的聚类。最终获得了10组聚类群（即10组蛋白群），代表了胚胎蛋白质的10种动力学模式，并在随后对这10组蛋白群的丰度变化与其功能展开了更细致的讨论（如基因集富集分析，蛋白网络分析等）。

数据来源

Gao等（2017）的蛋白质表达矩阵表格可以在原文献的补充材料Table S1中自行下载（Excel表格中，Sheet名称为“union_all_protein_exp_cluster”），或者直接下载：

https://www.cell.com/cms/10.1016/j.celrep.2017.11.111/attachment/bc1eedf3-89d1-473f-9b66-b44a17994029/mmc2.xlsx

蛋白表达数据读取并构建对象：

protein <- read.delim("uniq_all_protein_exp.txt", row.names = 1, check.names = FALSE,

sep = "\t")

protein <- as.matrix(protein)

head(protein)

## zygote 2-cell 4-cell 8-cell morula blastocyst

## Oog4 1.3132282 1.2370781 1.325978 1.262073 0.6549312 0.2067114

## Psmd9 1.0917337 1.3159888 1.174417 1.064756 0.8685598 0.4845448

## Sephs2 0.9859232 1.2010257 1.123076 1.084673 0.8878931 0.7174088

## Nhlrc2 0.9856354 1.0387869 1.061926 1.076825 0.9716945 0.8651322

## Trappc4 1.0775310 0.9757542 1.065544 1.080973 0.9732145 0.8269832

## Ywhah 1.0485306 1.0212216 1.117839 1.199569 1.0384096 0.5744298

mfuzz_class <- new("ExpressionSet", exprs = protein)

预处理缺失值或者异常值：

mfuzz_class <- filter.NA(mfuzz_class, thres = 0.25)

## 0 genes excluded.

mfuzz_class <- fill.NA(mfuzz_class, mode = "mean")

mfuzz_class <- filter.std(mfuzz_class, min.std = 0)

## 0 genes excluded.

标准化：

mfuzz_class <- standardise(mfuzz_class)

Mfuzz 基于 fuzzy c-means 的算法进行聚类，需手动定义目标聚类群的个数，例如这里我们为了重现原作者的结果，设定为 10，即期望获得 10 组聚类群。

需要设定随机数种子，以避免再次运行时获得不同的结果

set.seed(123)

cluster_num <- 10

mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))

绘制时间趋势图：

mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(protein))

查看每个聚类群中各自包含的蛋白数量

cluster_size <- mfuzz_cluster$size

names(cluster_size) <- 1:cluster_num

cluster_size

## 1 2 3 4 5 6 7 8 9 10

## 479 209 244 829 402 231 598 342 211 222

查看耦合性即成对的聚类定义相似性度量：

O <- overlap(mfuzz_cluster)

Ptmp <- overlap.plot(mfuzz_cluster, over = O, thres = 0.05)

References:

Kumar L, E Futschik M. Mfuzz: a software package for soft clustering of microarray data. Bioinformation. 2007;2(1):5-7. Published 2007 May 20. doi:10.6026/97320630002005Gao Y, Liu X, Tang B, et al. Protein Expression Landscape of Mouse Embryos during Pre-implantation Development. Cell Rep. 2017;21(13):3957-3969. doi:10.1016/j.celrep.2017.11.111

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

桓峰基因和投必得合作，文章润色优惠85折，需要文章润色的老师可以直接到网站输入领取桓峰基因专属优惠券码：KYOHOGENE，然后上传，付款时选择桓峰基因优惠券即可享受85折优惠哦！https://www.topeditsci.com/

有想进生信交流群的老师可以扫最后一个二维码加微信，备注“单位+姓名+目的”，有些想发广告的就免打扰吧，还得费力气把你踢出去！

好文推荐

 您阅读本篇文章共花了：