---

歡迎關(guān)注天下博客：http://blog.genesino.com/2017/11/gene-cluster/
共表達(dá)基因的尋找是轉(zhuǎn)錄組分析的一個(gè)部分，樣品多可以使用WGCNA，樣品少可直接通過(guò)聚類(lèi)分析如K-means、K-medoids (比K-means更穩(wěn)定)或Hcluster或設(shè)定pearson correlation閾值來(lái)選擇共表達(dá)基因。

下面將實(shí)戰(zhàn)演示K-means、K-medoids聚類(lèi)操作和常見(jiàn)問(wèn)題：如何聚類(lèi)分析，如何確定合適的cluster數(shù)目，如何繪制共表達(dá)密度圖、線圖、熱圖、網(wǎng)絡(luò)圖等。

獲得模擬數(shù)據(jù)集

MixSim是用來(lái)評(píng)估聚類(lèi)算法效率生成模擬數(shù)據(jù)集的一個(gè)R包。

library(MixSim) # 獲得5個(gè)中心點(diǎn)，8維屬性的數(shù)據(jù)模型 data = MixSim(MaxOmega=0, K=5, p=8, ecc=0.5, int=c(10, 100)) # 根據(jù)模型獲得1000次觀察的數(shù)據(jù)集 A <- simdataset(n=1000, Pi=data$Pi, Mu=data$Mu, S=data$S, n.out=0) data <- A$X # 數(shù)據(jù)標(biāo)準(zhǔn)化 data <- t(apply(data, 1, scale)) # 定義行列名字 rownames(data) <- paste("Gene", 1:1000, sep="_") colnames(data) <- letters[1:8] head(data) ## a b c d e f ## Gene_1 -0.04735251 -0.7147304 0.3836436 1.322786 0.9718988 -0.5468517 ## Gene_2 0.09276733 -0.8066507 0.5476909 1.351780 0.8679073 -0.6019107 ## Gene_3 -0.08751894 -0.6461075 0.4371506 1.522767 0.8031865 -0.6904081 ## Gene_4 0.11065988 -0.7327674 0.4550544 1.379773 0.9304277 -0.5532253 ## Gene_5 -0.02722127 -0.7833089 0.6700604 1.448916 0.7128284 -0.6266295 ## Gene_6 0.15119823 -0.7468292 0.4859932 1.351159 0.9179421 -0.5625206 ## g h ## Gene_1 0.4127370 -1.782130 ## Gene_2 0.2852284 -1.736813 ## Gene_3 0.3420581 -1.681128 ## Gene_4 0.1808146 -1.770737 ## Gene_5 0.2936467 -1.688292 ## Gene_6 0.1821925 -1.779136

K-means

K-means稱(chēng)為K-均值聚類(lèi)；k-means聚類(lèi)的基本思想是根據(jù)預(yù)先設(shè)定的分類(lèi)數(shù)目，在樣本空間隨機(jī)選擇相應(yīng)數(shù)目的點(diǎn)做為起始聚類(lèi)中心點(diǎn)；然后將空間中到每個(gè)起始中心點(diǎn)距離最近的點(diǎn)作為一個(gè)集合，完成第一次聚類(lèi)；獲得第一次聚類(lèi)集合所有點(diǎn)的平均值做為新的中心點(diǎn)，進(jìn)行第二次聚類(lèi)；直到得到的聚類(lèi)中心點(diǎn)不再變化或達(dá)到嘗試的上限，則完成了聚類(lèi)過(guò)程。

聚類(lèi)模擬如下圖：

聚類(lèi)過(guò)程需要考慮下面3點(diǎn)：

1.需要確定聚出的類(lèi)的數(shù)目。可通過(guò)遍歷多個(gè)不同的聚類(lèi)數(shù)計(jì)算其類(lèi)內(nèi)平方和的變化，并繪制線圖，一般選擇類(lèi)內(nèi)平方和降低開(kāi)始趨于平緩的聚類(lèi)數(shù)作為較優(yōu)聚類(lèi)數(shù), 又稱(chēng)elbow算法。下圖中拐點(diǎn)很明顯，5。

tested_cluster <- 12 wss <- (nrow(data)-1) * sum(apply(data, 2, var)) for (i in 2:tested_cluster) {wss[i] <- kmeans(data, centers=i,iter.max=100, nstart=25)$tot.withinss } plot(1:tested_cluster, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")

2.K-means聚類(lèi)起始點(diǎn)為隨機(jī)選取，容易獲得局部最優(yōu)，需重復(fù)計(jì)算多次，選擇最優(yōu)結(jié)果。

library(cluster) library(fpc) # iter.max: 最大迭代次數(shù) # nstart: 選擇的隨機(jī)集的數(shù)目 # centers: 上一步推測(cè)出的最優(yōu)類(lèi)數(shù)目 center = 5 fit <- kmeans(data, centers=center, iter.max=100, nstart=25) withinss <- fit$tot.withinss print(paste("Get withinss for the first run", withinss)) ## [1] "Get withinss for the first run 44.381088289378" try_count = 10 for (i in 1:try_count) {tmpfit <- kmeans(data, centers=center, iter.max=100, nstart=25)tmpwithinss <- tmpfit$tot.withinssprint(paste(("The additional "), i, 'run, withinss', tmpwithinss))if (tmpwithinss < withinss){withins <- tmpwithinssfit <- tmpfit} } ## [1] "The additional 1 run, withinss 44.381088289378" ## [1] "The additional 2 run, withinss 44.381088289378" ## [1] "The additional 3 run, withinss 44.381088289378" ## [1] "The additional 4 run, withinss 44.381088289378" ## [1] "The additional 5 run, withinss 44.381088289378" ## [1] "The additional 6 run, withinss 44.381088289378" ## [1] "The additional 7 run, withinss 44.381088289378" ## [1] "The additional 8 run, withinss 44.381088289378" ## [1] "The additional 9 run, withinss 44.381088289378" ## [1] "The additional 10 run, withinss 44.381088289378" fit_cluster = fit$cluster

簡(jiǎn)單繪制下聚類(lèi)效果

clusplot(data, fit_cluster, shade=T, labels=5, lines=0, color=T, lty=4, main='K-means clusters')

3.預(yù)處理：聚類(lèi)變量值有數(shù)量級(jí)上的差異時(shí)，一般通過(guò)標(biāo)準(zhǔn)化處理消除變量的數(shù)量級(jí)差異。聚類(lèi)變量之間不應(yīng)該有較強(qiáng)的線性相關(guān)關(guān)系。(最開(kāi)始模擬數(shù)據(jù)集獲取時(shí)已考慮)

K-medoids聚類(lèi)

K-means算法執(zhí)行過(guò)程，首先需要隨機(jī)選擇起始聚類(lèi)中心點(diǎn)，后續(xù)則是根據(jù)聚類(lèi)結(jié)點(diǎn)算出平均值作為下次迭代的聚類(lèi)中心點(diǎn)，迭代過(guò)程中計(jì)算出的中心點(diǎn)可能在觀察數(shù)據(jù)中，也可能不在。如果選擇的中心點(diǎn)是離群點(diǎn) (outlier)的話，后續(xù)的計(jì)算就都被帶偏了。而K-medoids在迭代過(guò)程中選擇的中心點(diǎn)是類(lèi)內(nèi)觀察到的數(shù)據(jù)中到其它點(diǎn)的距離最小的點(diǎn)，一定在觀察點(diǎn)內(nèi)。兩者的差別類(lèi)似于平均值和中值的差別，中值更為穩(wěn)健。

圍繞中心點(diǎn)劃分（Partitioning Around Medoids，PAM）的方法是比較常用的(cluster::pam)，與K-means相比，它有幾個(gè)特征: 1.
接受差異矩陣作為參數(shù)；2. 最小化差異度而不是歐氏距離平方和，結(jié)果更穩(wěn)健；3. 引入silhouette plot評(píng)估分類(lèi)結(jié)果，并可據(jù)此選擇最優(yōu)的分類(lèi)數(shù)目; 4. fpc::pamk函數(shù)則可以自動(dòng)選擇最優(yōu)分類(lèi)數(shù)目，并根據(jù)數(shù)據(jù)集大小選擇使用pam還是clara (方法類(lèi)似于pam，但可以處理更大的數(shù)據(jù)集) 。

fit_pam <- pamk(data, krange=2:10, critout=T)

不同的分類(lèi)書(shū)計(jì)算出的silhouette值如下，越趨近于1說(shuō)明分出的類(lèi)越好。

## 2 clusters 0.5288058 ## 3 clusters 0.6915659 ## 4 clusters 0.8415226 ## 5 clusters 0.8661989 ## 6 clusters 0.7415207 ## 7 clusters 0.5862313 ## 8 clusters 0.4196284 ## 9 clusters 0.2518583 ## 10 clusters 0.116984

獲取分類(lèi)的數(shù)目

fit_pam$nc ## [1] 5 layout(matrix(c(1, 2), 1, 2)) plot(fit_pam$pamobject) layout(matrix(1)) #改回每頁(yè)一張圖

獲取分類(lèi)信息

fit_cluster <- fit_pam$pamobject$clustering

數(shù)據(jù)提取和可視化

以pam的輸出結(jié)果為例 (上面兩種方法的輸出結(jié)果都已處理為了同一格式，后面的代碼通用)。

1.獲取每類(lèi)數(shù)值的平均值，利用線圖一步畫(huà)圖法獲得各個(gè)類(lèi)的趨勢(shì)特征。

cluster_mean <- aggregate(data, by=list(fit_cluster), FUN=mean) write.table(t(cluster_mean), file="ehbio.pam.cluster.mean.xls", sep='\t',col.names=F, row.names=T, quote=F)

2.獲取按照分類(lèi)排序的數(shù)據(jù)，繪制熱圖 (點(diǎn)擊查看)。

dataWithClu <- cbind(ID=rownames(data), data, fit_cluster) dataWithClu <- as.data.frame(dataWithClu) dataWithClu <- dataWithClu[order(dataWithClu$fit_cluster),] write.table(dataWithClu, file="ehbio.pam.cluster.xls", sep="\t", row.names=F, col.names=T, quote=F)

3.獲取每類(lèi)數(shù)據(jù)，繪制多線圖和密度圖

cluster3 <- data[fit_cluster==3,] head(cluster3) ## a b c d e f ## Gene_413 -1.2718728 0.6957162 0.9963399 -0.1895966 0.1786798 1.225407 ## Gene_414 -1.1705230 0.6765085 0.8689340 -0.2155533 0.4176178 1.251818 ## Gene_415 -0.9545339 0.5635188 0.8437158 -0.1360588 0.2084771 1.316728 ## Gene_416 -1.0888687 0.8269888 0.7590209 -0.3090701 0.4478664 1.275057 ## Gene_417 -1.1230295 0.8282559 0.9112640 -0.2524612 0.3966905 1.104951 ## Gene_418 -1.1291253 0.9574904 0.8405449 -0.1200131 0.1964983 1.155913 ## g h ## Gene_413 -0.11021165 -1.524462 ## Gene_414 -0.22266636 -1.606135 ## Gene_415 -0.03562464 -1.806222 ## Gene_416 -0.32047268 -1.590522 ## Gene_417 -0.21171923 -1.653951 ## Gene_418 -0.27792177 -1.623386 cluster3 <- t(cluster3)

多線圖，繪制見(jiàn)線圖繪制。

data_rownames <- rownames(cluster3) data_mean <- data.frame(id=data_rownames, data_mean=rowMeans(cluster3)) data_mean ## id data_mean ## a a -1.1708878 ## b b 0.7811888 ## c c 0.8443212 ## d d -0.2008149 ## e e 0.2382650 ## f f 1.2601860 ## g g -0.1612029 ## h h -1.5910553 library(reshape2) library(ggplot2) data_melt <- melt(cluster3) colnames(data_melt) <- c("id", "Gene", "Expr") ggplot(data_melt, aes(id, Expr)) + stat_density_2d(aes(alpha=..level.., group=1)) + stat_smooth(data=data_mean, mapping=aes(x=id, y=data_mean, colour="red", group=1), se=F) + theme(legend.position='none') ## `geom_smooth()` using method = 'loess'

等高線的顏色越深表示對(duì)應(yīng)的Y軸的點(diǎn)越密，對(duì)平均值的貢獻(xiàn)越大；顏色淺的點(diǎn)表示分布均勻。不代表點(diǎn)的多少。等高線的變化趨勢(shì)與平均值曲線一致。

參考

GIF來(lái)自 https://commons.wikimedia.org/wiki/File:Kmeans_animation.gif

PAM更多描述見(jiàn)http://shiyanjun.cn/archives/1398.html

總結(jié)

以上是生活随笔為你收集整理的基因共表达聚类分析及可视化的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。