知识分享| 转录组个性化分析(1)——时序分析

       近年来，随着测序成本的不断降低，转录组测序已经成为解析基因组的功能成分、揭示细胞和组织的分子成分以及理解发育和疾病必不可少的技术手段，已广泛应用于基础研究、临床诊断等领域。发表论文作为检验研究成果的一条重要途径，大家都在思考高分文章到底该从何入手，一方面需要全面、细致的实验设计，但另一方面，区别于常规的分析内容，清晰明了又好看的结果图片也都会为文章添砖加瓦！那转录组测序除了常见的差异基因通路富集及结构分析之外，还有哪些适用的分析内容呢？

       本文给大家带来“转录组时序分析”实操，希望给各位老师数据分析带来一些帮助。更多转录组个性化分析内容，后续将陆续推送，大家敬请期待！

1.时间序列分析的意义

       1）将表达相似的基因聚集为1类，表达相似的基因可能具有相似的生物学功能或者共同的细胞起源；

       2）基因调控的顺序，基因间存在诱导或者阻遏；

       3）同一种处理在不同的时间下，生物学动态变化。

2.Mfuzz软件包

       Mfuzz是专门的做转录变化的时间趋势分析的方法，核心算法基于模糊c均值聚类（Fuzzy C-Means Clustering，FCM），是基于R编程的软件包。

       聚类方法：如果一个聚类方法假定一个样本只能属于一个类，或类的交集为空集，那么该方法称为硬聚类（hard clustering）方法；如果一个样本可以属于多个类，或类的交集不为空集，那么该方法称为软聚类（soft clustering）方法。
3.Mfuzz包安装及加载

       使用 Bioconductor 安装 Mfuzz 包

install.packages('BiocManager')

BiocManager::install('Mfuzz')    ##当出现 package(s) not installed when version(s) same as current; use `force = TRUE` to re-install:

BiocManager::install('Mfuzz',force = TRUE) ##具体命令

library(Mfuzz)    #加载 Mfuzz 包
4.数据准备

       all.fpkm.xls ##基因表达矩阵

        #该示例中，第一列为基因或蛋白名称，第一行为时间样本序号（按时间顺序提前排列好，若存在生物学重复需提前取均值）                                                       

5.读取矩阵表格

gene <- read.delim('all.fpkm.xls', row.names = 1, check.names = FALSE)

gene <- as.matrix(mRNA)
6.构建对象

mfuzz_class <- new('ExpressionSet',exprs = gene)
7.缺失值或者异常值过滤

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

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

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

#通常来说，比较受欢迎的就是样本的标准差作为阈值
8.标准化数据

mfuzz_class <- standardise(mfuzz_class)

#Mfuzz 基于 fuzzy c-means 的算法进行聚类

#需手动定义目标聚类群的个数，如望获得 10 组聚类群

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

set.seed(123)

cluster_num <- 10

mfuzz_cluster<- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))
9.作图

#time.labels 参数设置时间轴，需要和原基因表达数据集中的列对应

#颜色、线宽、坐标轴、字体等细节也可以添加其他参数调整，此处略，详见函数帮助

mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(3, 5), time.labels = colnames(gene))
10.查看相关结果

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

cluster_size <- mfuzz_cluster$size

names(cluster_size) <- 1:cluster_num

cluster_size

查看每个蛋白所属的聚类群

head(mfuzz_cluster$cluster)

#Mfuzz 通过计算 membership 值判断基因所属的聚类群，以最大的 membership 值为准

查看各基因的 membership 值

head(mfuzz_cluster$membership)

提取所有基因所属的聚类群，并和它们的原始表达值整合在一起

gene_cluster <- mfuzz_cluster$cluster

gene_cluster <- cbind(gene[names(gene_cluster), ], gene_cluster)

head(gene_cluster)

write.table(gene_cluster, 'gene_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)
11.结果解读

图片来自官网示例（部分图）

如何评判聚类的好坏？

因为该图由基因的变化趋势聚类得到的，所以如果聚类越好，意味着每个cluster

中线段的走向一致。

如何判断聚类的个数？

Dmin图