转录组中的基因表达模式聚类分析

实验设计对于转录组数据的分析是非常重要的，对于常规的case/control实验设计，通过两组间的差异检验就可以得到不同条件下的差异基因；对于多组的实验设计，可以每两组之间进行差异分析，也可以通过annova的检验，得到差异基因。

在多组的实验设计中，有一种类型非常特殊，就是时间序列。这里的时间序列不仅仅指的是单纯的不同时间点取样，也包括生成发育的不同阶段，疾病治疗的不同阶段等。生命活动是一个动态变化的过程，通过时间序列的实验设计可以有效的挖掘基因的表达规律。

对于时间序列的实验而言，通常会有多个时间点的设计，当然我们也可以两两之间进行差异分析，或者所有时间点进行方差分析，但是这样得到的差异基因并不能有效的代表整个时间序列中的变化，而且两两分析会得到很多的差异基因列表，也增加了后续生物学规律的挖掘的难度。

为了有效利用时间序列的信息，科学家们提出了一种新的分析策略，就是基因表达模式聚类。所谓的表达模式，其实就是基因表达量随着时间点变化的一个规律，比如对于一个疾病的临床治疗而言，选择患病，药物治疗，治愈后3个时间段分析，我们希望看到药物治疗对疾病是有效的，与此相关的表达模式包括以下两种

这种先上调后下调的变化规律就是一个特定的表达模式，符合某种特定模式的基因可能是参与相同的代谢通路，也可能是受到了相同分子的调控。对于特定表达模式下的基因进行GO, KEGG等功能富集分析，有助于挖掘潜在规律。

不同于传统的差异分析，基因表达模式聚类分析中更关键的是筛选感兴趣的表达模式，即表达量的变化规律，然后对给模式下的基因进行后续的功能富集分析。

能够进行基因表达模式聚类分析的软件有很多，STEM自带图形界面，操作简单，是使用的较为广泛的软件之一，官网如下

http://www.cs.cmu.edu/~jernst/stem/

该软件适用于处理时间节点较少的数据，比如8个或者更少，因为时间点越多，计算量越大，计算时间会特别长，占用的硬件资源也会更多。同时该软件集成了Gene Ontology数据库，可以方便的对特定表达模式下的基因进行GO富集分析。

该软件的安装非常简单，下载安装包，双击jar文件即可启动。界面如下

界面分为4个部分，第一部分为Expression Data Info, 用于加载基因的表达量数据，第二部分为Gene Info, 用于加载基因对应的注释信息，第三部分为Options, 用于设置聚类的算法和参数，第四部分为Execute, 用于执行程序。

1. Expression Daa Info

该软件支持读取芯片和NGS两种数据，通过Data File加载表达量文件，表达量文件的格式如下

第一列为SPOT, 代表探针ID, 每个探针ID必须是唯一的，对于NGS数据，可以不要这一列信息；第二列为Gene Symbol, 如果没有gene symbol信息，对应的单元格为空或者用0填充，其他列为每个时间点对应的表达量，时间点按照顺序排列，表达量允许为空。

如果每个时间点有生物学重复的话，可以通过Repeat Data加载生物学重复的表达量文件，格式和上述格式完全一样，要注意的一点就是行列必须和Data File文件中一致，这意味着每个时间点的生物学重复个数必须一致。

在芯片中会出现多个探针对应同一个基因的情况，该软件会将同一个基因的多个探针的表达量取中位数，作为该基因的表达量；对于生物学重复，也是取中位数作为最终的表达量。

对于生物学重复个数不同的情况，可以自己先合并生物学重复，即选取中位数作为表达量，然后在进行分析。

上传成功后，通过View Data File按钮进行查看，示意如下

在进行聚类前，该软件首先会对原始的表达量进行转化，有以下3种策略

1. Log normalize data

2. Normalize data

3. No normalization/add 0

具体的转换规则如下

将第一个时间点作为control, 将其表达量定义为0， 其他的时间点依次与第一个时间点原本的表达量进行比较。通常我们选择log转换的方法。

从这里可以看出，STEM并不会对输入的原始数据进行样本间的归一化，所以我们的输入数据必须是可以直接在样本间比较的表达量，对于raw count的定量方式而言，需要输入归一化之后的表达量。

2. Gene  Info

STEM支持自定义基因注释，格式如下

第一列为基因或者探针ID，第二列为对应的注释，多个注释用分号:，逗号，或者 | 分隔。

3. Options

STEM支持以下两种聚类方法：

1. STEM

2. K-means

默认选择STEM聚类算法，该算法需要调整以下两个参数

每一个特定的表达模式称之为model profile, STEM通过第二个参数控制两个时间点间基因变化趋势的类别，这个参数可以看做差异的倍数，默认值为2，示意图如下

对于两个连续的时间点而言，STEM在判断变化趋势时不是简单的上调和下调两种，而是根据差异的倍数进行了细分，在上图中，根据差异的倍数可以划分出5个趋势，第一个为上调倍数2倍以上，第二个为上调倍数在1倍到2倍之间，第三个为差异倍数在1倍以内，第四个为下调倍数在1倍到2倍之间，第五个为下调倍数在2倍以上。

用C代表上面的倍数值，用n代表时间点的个数，则所有可能的profile的个数为

当时间点较多时，所有的profiile太多了，根本看不过来，所以需要对其过滤，最终挑选出多少个profile 则有第一个参数控制。STEM根据profile之间的距离，从所有的profile中挑选出距离最大的N个profile, 任意两个profile间的距离都很大，意味着它们是完全不同的profile。

在profile中，有一部分是由于生物学规律的影响而出现的特定表达模式，有一部分是随机出现的，为了剔除这部分随机性的profike, 通过特定的统计模型计算每个profile的p值，p值小于0.001的认为是真实的, 有研究价值的profile。

运行完成后，会产生如下的结果

每个profile对应一个矩形，在矩形左上角的数字为profile的编号，从0开始，其中的折线为表达量随时间变化的趋势，有颜色标记的是p值小于0.001的profile。
每个profile都是可以点击的，点击之后可以看到该profile下所有基因的变化趋势图

Profile Gene Table可以看到该profiel下具体的基因列表和转换后的表达量，示意如下

这个列表可以导出，对于感兴趣的profile, 可以导出该profile下的基因，去做pathway等富集分析。

以上就是该软件的基本用法，更加详细的用法请参考官方文档。

·end·