今天在研究一篇關(guān)于硬葉蘭（Cymbidium mannii）的基因組學(xué)的文章的時(shí)候偶然發(fā)現(xiàn)了一個(gè)很好的分組轉(zhuǎn)錄組表達(dá)水平的R包Mfuzz，所以來(lái)記錄一下基于示例數(shù)據(jù)的代碼運(yùn)行過(guò)程。
圖1 硬葉蘭在24小時(shí)周期內(nèi)基因表達(dá)模式的Mfuzz簇。每個(gè)簇中的生物鐘(黑色)和核心CAM基因紅色)在右側(cè)被標(biāo)記
該包最早被開發(fā)應(yīng)用在不同的小鼠細(xì)胞內(nèi)蛋白組的變化，隨后逐漸被應(yīng)用于越來(lái)越多的多組數(shù)據(jù)的聚類表達(dá)分析。在介紹該包的使用之前先簡(jiǎn)單的介紹一下該包的優(yōu)點(diǎn)。以硬葉蘭基因組文章的舉例，分別將7個(gè)時(shí)段（7個(gè)組）的轉(zhuǎn)錄組的差異表達(dá)水平進(jìn)行整理，以往的分析方法通常是對(duì)差異表達(dá)變化相似的基因進(jìn)行書面的整理，并不能給讀者直觀的體現(xiàn)，很難瞬間明白變化的規(guī)律。而Mfuzz可以通過(guò)將7個(gè)時(shí)段的所有基因的表達(dá)量變化進(jìn)行聚類和整理，從而篩選出若干組（由人為指定）表達(dá)量變化相似的差異基因集，并且將每組基因的表達(dá)量基因的變化繪制在一張圖上。
接下來(lái)是Mfuzz的運(yùn)行全程的代碼記錄。
#切換到工作路徑
setwd('D:/生物信息學(xué)數(shù)據(jù)庫(kù)集合/Mfuzz/')
#首次使用需要下載Mfuzz包，如果沒安裝過(guò)BioManager也要安裝BioManager
install.packages('BiocManager')
BiocManager::install("Mfuzz")
library("Mfuzz")
#讀入示例數(shù)據(jù)，示例數(shù)據(jù)在另外一個(gè)博主的網(wǎng)盤中，直接下載即可，以txt結(jié)尾的就是示例數(shù)據(jù)
protein <- read.delim('mmc2.union_all_protein_exp.txt', row.names = 1, check.names = FALSE)
protein <- as.matrix(protein)
mfuzz_class <- new('ExpressionSet',exprs = protein)
#預(yù)處理缺失值或者異常值，這一步主要是對(duì)你的數(shù)據(jù)進(jìn)行初步的處理，將一些可能存在誤差或者錯(cuò)誤的數(shù)據(jù)剔除
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)
#標(biāo)準(zhǔn)化數(shù)據(jù)
mfuzz_class <- standardise(mfuzz_class)
#需手動(dòng)定義目標(biāo)聚類群的個(gè)數(shù)，例如這里我們?yōu)榱酥噩F(xiàn)原作者的結(jié)果，設(shè)定為 10，即期望獲得 10 組聚類群
#需要設(shè)定隨機(jī)數(shù)種子，以避免再次運(yùn)行時(shí)獲得不同的結(jié)果
set.seed(123)
cluster_num <- 10
mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))
#作圖
#time.labels 參數(shù)設(shè)置時(shí)間軸，需要和原基因表達(dá)數(shù)據(jù)集中的列對(duì)應(yīng)
#顏色、線寬、坐標(biāo)軸、字體等細(xì)節(jié)也可以添加其他參數(shù)調(diào)整，此處略，詳見函數(shù)幫助
mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(protein))
#到了上一步就已經(jīng)分析結(jié)束，并且得到了結(jié)果，但是一般我們還要知道那些基因具有相似的表達(dá)量變化所以還要把結(jié)果以文本的形式提取出來(lái)
#查看每個(gè)聚類群中各自包含的蛋白數(shù)量
cluster_size <- mfuzz_cluster$size
names(cluster_size) <- 1:cluster_num
cluster_size
#查看每個(gè)蛋白所屬的聚類群
head(mfuzz_cluster$cluster)
#Mfuzz 通過(guò)計(jì)算 membership 值判斷蛋白質(zhì)所屬的聚類群，以最大的 membership 值為準(zhǔn)
#查看各蛋白的 membership 值
head(mfuzz_cluster$membership)
#最后，提取所有蛋白所屬的聚類群，并和它們的原始表達(dá)值整合在一起
protein_cluster <- mfuzz_cluster$cluster
protein_cluster <- cbind(protein[names(protein_cluster), ], protein_cluster)
head(protein_cluster)
write.table(protein_cluster, 'protein_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)
參考文獻(xiàn)：
Gao Yawei, Liu Xiaoyu, Tang Bin et al. Protein Expression Landscape of Mouse Embryos during Pre-implantation Development. Cell Rep, 2017, 21: 3957-3969.
Fan W, He ZS, Zhe M, Feng JQ, Zhang L, Huang Y, Liu F, Huang JL, Ya JD, Zhang SB, Yang JB, Zhu A, Li DZ. High-quality Cymbidium mannii genome and multifaceted regulation of crassulacean acid metabolism in epiphytes. Plant Commun. 2023 Sep 11;4(5):100564.
示例數(shù)據(jù)下載，https://pan.baidu.com/s/1wXYop3wql8SZqcM2CzpJIg，密碼4mtu.
