前不久小编接触了这样一个分析，给定基因表达矩阵后，通过自组织映射（Self-organizing map，SOM）技术识别其中的高表达基因集，以便和功能建立关联。

下文将该方法简称为SOM聚类分析，因为它就是一种基于神经网络的聚类算法。那么，SOM聚类在表达谱数据中是如何应用的呢？本篇我们就来看一下。

相关文献分析思路解读

为了帮助大家理解这种方法，首先来看文献“NLRP3 inflammasome activation drives tau pathology”中的部分内容。

作者构建小鼠模型，模拟额颞叶痴呆（FTD）病理学效应。Tau22小鼠转基因了人类tauFTD变体，并在一段时间内发展为tau病理学。获取野生型（WT）和Tau22小鼠的脑组织并提取RNA，包括2、8、11月龄的小鼠，进行RNA芯片分析基因表达，获得脑组织在tau病理学过程中显著被激活的基因。

试验流程图

作者通过SOM聚类，鉴定了参与tau蛋白的致病性基因集。分别根据野生型或Tau22小鼠在3、8、11月时间点的相似表达水平将基因分组，并根据表达水平升高将其定义为每种条件的特征基因。共识别了6个主要的表达模块，同一模块内的基因集具有相似的表达模式，它们在该时间点均处于相对高表达的状态，暗示它们在这些时间点中发挥主要活性。

基因表达谱的SOM聚类分析

为了明确这些高度活跃的基因发挥的功能，随后作者通过基因集富集分析（GSEA），比较Tau22小鼠相较于WT小鼠中哪些途径被激活。在3个月大的Tau22小鼠中，特征基因显示与免疫反应的联系，表明在疾病开始时特征基因就参与了免疫过程；而在疾病发展后期，小鼠中上调的基因参与了诸如“应激反应”等功能中，且高达73%的特征基因与干扰素相关。

R语言的SOM聚类分析

不难理解，上述文献中，作者通过SOM聚类识别高度表达的基因集，用作功能分析。

本篇模仿该文献中的思路，展示如何在R中执行基因表达谱的SOM聚类。

使用kohonen包执行SOM聚类，首先安装加载该包。

#安装及加载kohonen包
#install.packages('kohonen')
library(kohonen)

#数据集yeast是800种酵母基因的微阵列细胞周期数据
data(yeast)

dat <- yeast$elu
head(dat)  #展示其中一种类型的芯片表达谱

示例基因表达矩阵

接下来基于示例的基因表达谱，对各基因执行SOM聚类。

##1、标准化表达值数据
dat_scale <- as.matrix(scale(dat))
names(dat_scale) <- names(dat)
head(dat_scale)

##2、训练SOM模型
set.seed(123)

#定义网络的大小和形状，以10*10为例
som_grid <- somgrid(xdim = 10, ydim = 10, topo = 'rectangular')  

#多层SOM
som_model <- supersom(dat_scale, grid = som_grid, keep.data = TRUE)

##3、作图，颜色代表了模块内基因平局表达值
coolBlueHotRed <- function(n, alpha = 0.7) rainbow(n, end=4/6, alpha=alpha)[n:1]

color_by = apply(som_model$data[[1]], 1, mean)
unit_colors <- aggregate(color_by, by = list(som_model$unit.classif), FUN = mean, simplify = TRUE)
unit_dat <- data.frame(value = rep(0, 100))
unit_dat[unit_colors$Group.1,'value'] <- unit_colors$x

plot(som_model, type = 'property', property = unit_dat[[1]], palette.name = coolBlueHotRed, 
    shape = 'straight', keepMargins = TRUE, border = NA) 

基于表达值的聚类图谱

如上过程基于基因表达值进行了聚类，获得了聚类模块，并按模块内基因的平均表达值赋值了模块颜色。随后，即可从图中判断选择高表达的模块，将其中的基因挑选出来，作为发挥生物学过程的“活跃”基因集。

那么，如何获得各模块中，都包含哪些基因呢？参考以下操作。

##获取每个SOM中心点相关的基因
som_model_class <- data.frame(name=rownames(som_model$data[[1]]), code_class=som_model$unit.classif)
head(som_model_class)

write.table(som_model_class, 'gene_SOM.txt', row.names = FALSE, sep = '\t', quote = FALSE)

基因和所属聚类模块的对应关系

这样，就将基因名称和其所属模块对应起来了。

最后，识别高表达的模块，并从中进一步筛选基因集就可以了。这些基因集既然存在高表达，那么必然会和功能密不可分。了解它们的功能，可以初步执行GO或KEGG富集分析进行探索，如开篇展示的文献中思路那样，不再多说了。