上周的暑期生信黑馬培訓有老師提出要做SOM分析，最後卡在code plot只能出segment plot卻出不來line plot。查了下，沒看到解決方案。今天看了下源碼，設置了一個參數，得到趨勢圖。也順便學習了SOM分析的整個過程，整理下來，以備以後用到。
SOM分析基本理論
SOM (Self-Organizing Feature Map,自組織特徵圖)是基於神經網絡方式的數據矩陣和可視化方式。與其它類型的中心點聚類算法如K-means等相似，SOM也是找到一組中心點 (又稱為codebook vector)，然後根據最相似原則把數據集的每個對象映射到對應的中心點。在神經網絡術語中，每個神經元對應於一個中心點。
與K-means類似，數據集中的每個對象每次處理一個，判斷最近的中心點，然後更新中心點。與K-means不同的是，SOM中中心點之間存在拓撲形狀順序，在更新一個中心點的同時，鄰近的中心點也會隨著更新，直到達到設定的閾值或中心點不再有顯著變化。最終獲得一系列的中心點 (codes)隱式地定義多個簇，與這個中心點最近的對象歸為同一個簇。
SOM強調簇中心點之間的鄰近關係，相鄰的簇之間相關性更強，更有利於解釋結果，常用於可視化網絡數據或基因表達數據。
Even though SOM is similar to K-means, there is a fundamental difference. Centroids used in SOM have a predetermined topographic ordering relationship. During the training process, SOM uses each data point to update the closest centroid and centroids that are nearby in the topographic ordering. In this way, SOM produces an ordered set of centroids for any given data set. In other words, the centroids that are close to each other in the SOM grid are more closely related to each other than to the centroids that are farther away. Because of this constraint, the centroids of a two-dimensional SOM can be viewed as lying on a two-dimensional surface that tries to fit the n-dimensional data as well as possible. The SOM centroids can also be thought of as the result of a nonlinear regression with respect to the data points. At a high level, clustering using the SOM technique consists of the steps described in Algorithm below:
1: Initialize the centroids.
2: repeat
3: Select the next object.
4: Determine the closest centroid to the object.
5: Update this centroid and the centroids that are close, i.e., in a specified neighborhood.
6: until The centroids don't change much or a threshold is exceeded.
7: Assign each object to its closest centroid and return the centroids and clusters.SOM分析實戰
下面是R中用kohonen包進行基因表達數據的SOM分析。
加載或安裝包
### LOAD LIBRARIES - install with:
#install.packages(c("kohonen")
library(kohonen)讀入數據並進行標準化
data <- read.table("ehbio_trans.Count_matrix.xls", row.names=1, header=T, sep="\t")
# now train the SOM using the Kohonen method
# 標準化數據
data_train_matrix <- as.matrix(t(scale(t(data))))
names(data_train_matrix) <- names(data)
head(data_train_matrix)untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011 trt_N61311
ENSG00000223972 1.6201852 -0.5400617 -0.5400617 -0.5400617 -0.5400617
ENSG00000227232 -1.0711639 1.0274429 0.6776751 0.8525590 -1.2460478
ENSG00000278267 -1.6476479 1.3480756 0.1497862 0.7489309 -0.4493585
ENSG00000237613 2.4748737 -0.3535534 -0.3535534 -0.3535534 -0.3535534
ENSG00000238009 -0.3535534 -0.3535534 -0.3535534 -0.3535534 2.4748737
ENSG00000268903 -0.7020086 0.9025825 -0.7020086 -0.7020086 -0.7020086
trt_N052611 trt_N080611 trt_N061011
ENSG00000223972 1.6201852 -0.5400617 -0.5400617
ENSG00000227232 -1.2460478 0.5027912 0.5027912
ENSG00000278267 0.7489309 0.1497862 -1.0485032
ENSG00000237613 -0.3535534 -0.3535534 -0.3535534
ENSG00000238009 -0.3535534 -0.3535534 -0.3535534
ENSG00000268903 0.9025825 -0.7020086 1.7048781訓練SOM模型
# 定義網絡的大小和形狀
som_grid <- somgrid(xdim = 10, ydim=10, topo="hexagonal")
# Train the SOM model!
som_model <- supersom(data_train_matrix, grid=som_grid, keep.data = TRUE)可視化SOM結果
# Plot of the training progress - how the node distances have stabilised over time.
# 展示訓練過程，距離隨著疊代減少的趨勢，判斷疊代是否足夠；最後趨於平穩比較好
plot(som_model, type = "changes")計量每個SOM中心點包含的基因的數目
## custom palette as per kohonen package (not compulsory)
coolBlueHotRed <- function(n, alpha = 0.7) {
rainbow(n, end=4/6, alpha=alpha)[n:1]
}
# shows the number of objects mapped to the individual units.
# Empty units are depicted in gray.
plot(som_model, type = "counts", main="Node Counts", palette.name=coolBlueHotRed)計量SOM中心點的內斂性和質量
# map quality
# shows the mean distance of objects mapped to a unit to
# the codebook vector of that unit.
# The smaller the distances, the better the objects are
# represented by the codebook vectors.
plot(som_model, type = "quality", main="Node Quality/Distance", palette.name=coolBlueHotRed)鄰居距離-查看潛在邊界點
# 顏色越深表示與周邊點差別越大，越是分界點
# neighbour distances
# shows the sum of the distances to all immediate neighbours.
# This kind of visualization is also known as a U-matrix plot.
# Units near a class boundary can be expected to have higher average distances to their neighbours.
# Only available for the "som" and "supersom" maps, for the moment.
plot(som_model, type="dist.neighbours", main = "SOM neighbour distances", palette.name=grey.colors)查看SOM中心點的變化趨勢
#code spread
plot(som_model, type = "codes", codeRendering="lines")獲取每個SOM中心點相關的基因
table(som_model$unit.classif)# 只顯示一部分
1 2 3 4 5 6
197 172 434 187 582 249
95 96 97 98 99 100
168 919 226 419 193 241# code是從左至右，從下至上進行編號的
som_model_code_class = data.frame(name=rownames(data_train_matrix), code_class=som_model$unit.classif)
head(som_model_code_class) name code_class
1 ENSG00000223972 81
2 ENSG00000227232 37
3 ENSG00000278267 93
4 ENSG00000237613 51
5 ENSG00000238009 11
6 ENSG00000268903 4SOM結果進一步聚類
# 選擇合適的聚類數目
# show the WCSS metric for kmeans for different clustering sizes.
# Can be used as a "rough" indicator of the ideal number of clusters
mydata <- as.matrix(as.data.frame(som_model$codes))
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
for (i in 2:15) wss[i] <- sum(kmeans(mydata, centers=i)$withinss)
par(mar=c(5.1,4.1,4.1,2.1))
plot(1:15, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares", main="Within cluster sum of squares (WCSS)")# Form clusters on grid
## use hierarchical clustering to cluster the codebook vectors
som_cluster <- cutree(hclust(dist(mydata)), 6)
# Colour palette definition
cluster_palette <- function(x, alpha = 0.6) {
n = length(unique(x)) * 2
rainbow(n, start=2/6, end=6/6, alpha=alpha)[seq(n,0,-2)]
}
cluster_palette_init = cluster_palette(som_cluster)
bgcol = cluster_palette_init[som_cluster]
#show the same plot with the codes instead of just colours
plot(som_model, type="codes", bgcol = bgcol, main = "Clusters", codeRendering="lines")
add.cluster.boundaries(som_model, som_cluster)
有一些類的模式不太明顯，以後再看怎麼優化。
SOM獲取基因所在的新類
som_model_code_class_cluster = som_model_code_class
som_model_code_class_cluster$cluster = som_cluster[som_model_code_class$code_class]
head(som_model_code_class_cluster) name code_class cluster
1 ENSG00000223972 81 2
2 ENSG00000227232 37 8
3 ENSG00000278267 93 8
4 ENSG00000237613 51 7
5 ENSG00000238009 11 4
6 ENSG00000268903 4 3映射某個屬性到SOM圖
# 此處選擇一個樣本作為示例，可以關聯很多信息，
# 比如基因通路，只要在矩陣後增加新的屬性就可以。
color_by_var = names(data_train_matrix)[1]
color_by = data_train_matrix[,color_by_var]
unit_colors <- aggregate(color_by, by=list(som_model$unit.classif), FUN=mean, simplify=TRUE)
plot(som_model, type = "property", property=unit_colors[,2], main=color_by_var, palette.name=coolBlueHotRed)