使用 Mfuzz 包聚类分析并自定义绘图

天下没有不散的宴席

1.引言

对于多个时间点的 RNA-seq 数据做差异分析将会有多个比较组，对于分析显然也会更加复杂一些。这种情况我们可以采用 聚类分析 ，将每个特定时期具有表达模式类似的基因聚类到一个 cluster，后续我们可以针对这些感兴趣的 cluster 进行后续分析。

今天介绍的 Mfuzz 包采用 模糊 c 均值聚类分析(fuzzy c-means algorithm) ，用于识别相似的基因表达谱，一定程度上降低了噪声对聚类结果的干扰，而且这种算法 有效的定义了基因和 cluster 之间的关系，即基因是否属于某个 cluster，对应的值为 memebership 。

2.安装

BiocManager::install("Mfuzz")

3.使用

我们需要提供归一化的表达量数据，TPM/FPKM/RPKM 均可，对于含有生物学重复的样本，可以取均值。

这里我们使用示例数据集进行演示:

# 读取数据
exps <- read.delim('c:/Users/admin/Desktop/protein_exp.txt',row.names = 1,header = T) %>%
  as.matrix()

# 查看数据
head(exps,3)
          zygote  X2.cell  X4.cell  X8.cell    morula blastocyst
Oog4   1.3132282 1.237078 1.325978 1.262073 0.6549312  0.2067114
Psmd9  1.0917337 1.315989 1.174417 1.064756 0.8685598  0.4845448
Sephs2 0.9859232 1.201026 1.123076 1.084673 0.8878931  0.7174088

行名为基因，列是受精卵不同发育的时期。

Mfuzz 聚类时要求是一个 ExpressionSet 类型 的对象，所以先对我们的数据构建这样一个对象:

# 构建对象
myset <- new("ExpressionSet",exprs = exps)

# 根据标准差去除样本间差异太小的基因
myset <- filter.std(myset,min.std=0)

# 标准化
myset <- standardise(myset)

# 聚类,人为指定聚类个数
cluster_number <- 10

# 第二个参数称之为fuzzifier值，用小写字母m表示，可以通过函数评估一个最佳取值
m <- mestimate(myset)

# 设置种子,以便可以重复上次运行的结果
set.seed(123)
mfuzz_res <- mfuzz(myset, c = cluster_number, m = m)

# 绘图
mfuzz.plot2(myset, # classExpressionSet
            cl = mfuzz_res, # 聚类结果
            mfrow = c(2, 5), # 图形排列
            time.labels = colnames(exps), # 横坐标标签
            x11 = F # 关闭弹出显示窗口
            )

自己换个颜色，mfuzz.plot2 函数还有其它可以调节的参数，大家自己取探索调试:

# 换颜色
library(RColorBrewer)
color <- colorRampPalette(rev(c("#E02401", "Yellow","#0F52BA", "#3E7C17")))(1000)

mfuzz.plot2(myset, # classExpressionSet
            cl = mfuzz_res, # 聚类结果
            mfrow = c(3, 4), # 图形排列
            time.labels = colnames(exps), # 横坐标标签
            x11 = F, # 关闭弹出显示窗口
            colo = color
)

这样调整的限度非常有限，而且没有颜色的注释，不如自己整理数据用 ggplot 绘图。

4.提取并匹配数据

查看每个 cluster 基因数量:

# 查看每个聚类群中基因数量
cluster_size <- mfuzz_res$size
names(cluster_size) <- 1:cluster_number
cluster_size
  1   2   3   4   5   6   7   8   9  10
479 209 244 829 402 231 598 342 211 222

查看基因所属 cluster :

# 查看每个基因所属的cluster
head(mfuzz_res$cluster)
Oog4   Psmd9  Sephs2  Nhlrc2 Trappc4   Ywhah
 5       5       5      10       5      10

查看每个基因与 cluster 的 membership :

# 查看基因与每个cluster的membership
head(mfuzz_res$membership)
                  1           2           3           4         5          6           7          8
Oog4    0.002222202 0.002859622 0.003910825 0.002207282 0.6473041 0.02534048 0.002399312 0.23212880
Psmd9   0.002165038 0.002822132 0.003648413 0.002132773 0.6470210 0.02620424 0.002303071 0.23416681
Sephs2  0.005289313 0.006941265 0.008328799 0.005172955 0.5082228 0.05106428 0.005593873 0.20919426
Nhlrc2  0.002830995 0.004295915 0.003497748 0.002629085 0.2424768 0.01256882 0.002704833 0.03946422
Trappc4 0.009214473 0.013404199 0.012850647 0.008740730 0.4311962 0.04180379 0.009005823 0.14050977
Ywhah   0.002951047 0.004880814 0.003381517 0.002677187 0.2043156 0.00999969 0.002661010 0.03215145
                  9         10
Oog4    0.014978531 0.06664882
Psmd9   0.011944222 0.06759231
Sephs2  0.019884453 0.18030803
Nhlrc2  0.007428461 0.68210311
Trappc4 0.042353455 0.29092096
Ywhah   0.008630061 0.72835165

查看归一化的值:

# 查看归一化的值
head(myset@assayData$exprs,3)
            zygote   X2.cell   X4.cell   X8.cell     morula blastocyst
Oog4    0.67469689 0.5106688 0.7021608 0.5645078 -0.7432820  -1.708752
Psmd9   0.31432997 1.0827515 0.5976492 0.2218889 -0.4503866  -1.766233
Sephs2 -0.07985995 1.1404532 0.6982307 0.4803662 -0.6360019  -1.603188

1.把 cluster 信息与原始表达值整合

# 与原始数据整合
raw_cluster_anno <- cbind(exps,cluster = mfuzz_res$cluster)

# 查看
head(raw_cluster_anno,3)
          zygote  X2.cell  X4.cell  X8.cell    morula blastocyst cluster
Oog4   1.3132282 1.237078 1.325978 1.262073 0.6549312  0.2067114       5
Psmd9  1.0917337 1.315989 1.174417 1.064756 0.8685598  0.4845448       5
Sephs2 0.9859232 1.201026 1.123076 1.084673 0.8878931  0.7174088       5

# 保存
write.csv(raw_cluster_anno,file = 'raw_cluster_anno.csv',row.names = T)

2.把 cluster 信息与归一化表达值整合

# 与归一化数据整合
norm_cluster_anno <- cbind(myset@assayData$exprs,cluster = mfuzz_res$cluster)

# 查看
head(norm_cluster_anno,3)
            zygote   X2.cell   X4.cell   X8.cell     morula blastocyst cluster
Oog4    0.67469689 0.5106688 0.7021608 0.5645078 -0.7432820  -1.708752       5
Psmd9   0.31432997 1.0827515 0.5976492 0.2218889 -0.4503866  -1.766233       5
Sephs2 -0.07985995 1.1404532 0.6982307 0.4803662 -0.6360019  -1.603188       5

# 保存
write.csv(norm_cluster_anno,file = 'norm_cluster_anno',row.names = T)

5.提取归一化数据重新绘图

# 提取归一化数据绘图

# membership矩阵处理
mem <- cbind(mfuzz_res$membership,cluster2 = mfuzz_res$cluster) %>%
  as.data.frame() %>%
  mutate(gene = rownames(.))

# 查看
head(mem,3)
                 1           2           3           4         5          6           7         8
Oog4   0.002222202 0.002859622 0.003910825 0.002207282 0.6473041 0.02534048 0.002399312 0.2321288
Psmd9  0.002165038 0.002822132 0.003648413 0.002132773 0.6470210 0.02620424 0.002303071 0.2341668
Sephs2 0.005289313 0.006941265 0.008328799 0.005172955 0.5082228 0.05106428 0.005593873 0.2091943
                9         10 cluster2   gene
Oog4   0.01497853 0.06664882        5   Oog4
Psmd9  0.01194422 0.06759231        5  Psmd9
Sephs2 0.01988445 0.18030803        5 Sephs2

然后整理每个基因属于具体 cluster 的 membership 的值: