探索数据集

本例使用的数据来自DeRisi等人1997[1]发表的酵母基因表达的微阵列研究。作者使用DNA微阵列研究了几乎所有基因的瞬时基因表达酿酒酵母酿酒酵母在从发酵到呼吸的代谢过程中。在diauxic移位期间的7个时间点测量表达水平。完整的数据集可以从基因表达综合网站下载，http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28．

的MAT-fileyeastdata.mat包含表达式值(log2的ratioCH2DN_MEAN和ch1dn_mean.）从实验中的七个时间步骤，基因的名称和测量表达水平的时间阵列。

负载yeastdata.mat

要了解您可以使用的数据的大小numel（基因）来显示数据集中包含了多少基因。

numel（基因）

ans = 6400.

你可以通过索引变量来访问与实验相关的基因名称基因，一个代表基因名称的细胞阵列。例如，第15个元素基因YAL054C。这表示变量的第15行yeastvalues包含YAL054C的表达水平。

基因{15}

ans = ' YAL054C '

可以使用一个简单的图来显示这个ORF的表达式概要文件。

情节（次，yeastvalues（15，:)）xlabel（的时间(小时)）;ylabel（'log2相对表达水平'）;

您还可以绘制实际表达式比率，而不是log2变换的值。

情节(*,2。^ yeastvalues(15日:))包含(的时间(小时)）;ylabel（“相对表达水平”）;

与该ORF，ACS1相关的基因似乎在辅助偏移期间强烈上调。您可以通过在同一图中绘制多条线来将该基因的表达与其他基因的表达进行比较。

抓住在情节(*,2。^ yeastvalues(: 16:26)”)包含(的时间(小时)）;ylabel（“相对表达水平”）;标题(概要文件表达水平的）;

过滤的基因

通常情况下，基因表达数据集包含与基因相对应的信息，这些基因在实验过程中没有显示出任何有趣的变化。为了更容易地找到感兴趣的基因，可以将数据集的大小减少到只包含最重要基因的某个子集。

如果你仔细看一下基因列表，你会看到一些标记为“空”．这些是数组上的空点，虽然它们可能有与之关联的数据，但对于本例，您可以将这些点视为噪声。这些点可以用比较字符串函数，并使用索引命令从数据集中删除。

emptyspots = strcmp（“空”，基因）;yeastvalues（空虚，:) = [];基因（Everyspots）= [];numel（基因）

ans = 6314

还有在数据集中看到的几个位置，其中表达级别被标记为南．这表明在特定时间步骤中没有收集数据。处理这些缺失值的一种方法是使用时间随时间使用特定基因的均值或数据的平均值或中位数来施加它们。此示例使用较少的严格方法，即仅丢弃未测量一个或多个表达水平的任何基因的数据。功能isnan.用于识别缺失数据的基因，索引命令用于删除缺失数据的基因。

nanIndices =任何(isnan (yeastvalues), 2);yeastvalues (nanIndices:) = [];基因(nanIndices) = [];numel（基因）

ANS = 6276.

如果您要绘制所有其余配置文件的表达式配置文件，您将看到大多数配置文件是扁平的，与其他配置文件没有显著差异。这个平坦的数据显然是有用的，因为它表明与这些谱图相关的基因不受双auxic转移的显著影响;然而，在这个例子中，您感兴趣的是伴随双auxic转移的表达变化较大的基因。您可以使用生物信息学工具箱™中的过滤功能来删除具有不同类型配置文件的基因，这些配置文件不能提供有关受代谢变化影响的基因的有用信息。

你可以使用GenevarFilter.滤除随时间变化小的基因的功能。该函数返回与变量大小相同的逻辑数组(即掩码)基因与相对应的行yeastvalues方差大于第10个百分位数，零对应于低于阈值的。您可以使用掩码对值建立索引，并删除过滤的基因。

面具= Genevarfilter（yeastvalues）;yeastvalues = yeastvalues（面具，:);基因=基因（面膜）;numel（基因）

ans = 5648.

功能genelowvalfilter去除绝对表达值很低的基因。注意，这些筛选函数还可以自动计算筛选后的数据和名称，因此没有必要使用掩码索引原始数据。

[掩盖，酵母，基因] = Genelowvalfilter（酵母，基因，基因，“absval”，log2（3））;numel（基因）

ans = 822

最后，您可以使用该函数geneentropyfilter移除那些拥有低熵值的基因，例如数据中15%的熵值水平。

(面具,yeastvalues,基因)= geneentropyfilter (yeastvalues,基因,'proctile'15);numel（基因）

ans = 614

聚类分析

现在您已经有了一个可管理的基因列表，您可以使用Statistics和Machine Learning Toolbox™中的一些不同的聚类技术来查找配置文件之间的关系。对于分层聚类，函数Pdist.计算配置文件和之间的成对距离连锁创建分层群集树。

corrDist = pdist (yeastvalues,'corr'）;clustertree =链接（corrdist，“平均”）;

这簇函数根据截断距离或最大簇数计算簇。在这种情况下maxclust选项用于识别16个不同的集群。

clusters = cluster（clustertree，“maxclust”16);

这些簇中的基因的谱可以使用简单的环路绘制在一起次要情节命令。

数字为C = 1:16 subplot(4,4, C);情节(次,yeastvalues ((: = = c)集群)');轴紧结束sgtitle（“配置文件的分层聚类”）;

统计和机器学习工具箱也有一个k -均值聚类函数。同样，我们发现了16个簇，但由于算法不同，这些簇不一定与分层聚类发现的簇相同。

初始化随机数生成器的状态，以确保这些命令生成的数字与本示例的HTML版本中的数字匹配。

rng (“默认”）;[cidx, ctrs] = kmeans(酵母值，16，“距离”那'corr'那“代表”5,“disp”那“最后一次”）;数字为C = 1:16 subplot(4,4, C);图（次数，yeastvalues（（cidx == c），:)'）;轴紧结束sgtitle（“概况的K-Means聚类”）;

复制1,21迭代，距离总和= 23.4699。复制2,22迭代，距离总和= 23.5615。复制3,10次迭代，距离总和= 24.823。复制4,28迭代，距离总和= 23.4501。复制5,19迭代，距离总和= 23.5109。最佳距离总和= 23.4501

不用画所有的剖面图，你可以只画质心。

数字为C = 1:16 subplot(4,4, C);情节(次,点击率数据(c:) ');轴紧轴从结束sgtitle（“概况的K-Means聚类”）;

你可以使用clustergram函数创建表达式级别的热图和从分层聚类的输出生成的树状图。

cgObj = clustergram (yeastvalues(:, 2:结束),“RowLabels”，基因，'columnlabels'，次（2：结束））;

主成分分析

主组件分析（PCA）是一种有用的技术，可用于降低大数据集的维度，例如微阵列。PCA也可用于在嘈杂数据中找到信号。功能mapcaplot.计算数据集的主要成分，并创建结果的散点图。您可以交互式地从其中一个图中选择数据点，这些点在另一个图中自动高亮显示。这让您可以同时可视化多个维度。

H = MapCaplot（酵母值，基因）;

注意，前两个主成分的分数的散点图显示有两个不同的区域。这并不意外，因为过滤过程删除了许多低变异或低信息的基因。这些基因会出现在散点图的中间。

如果要查看主组件的值，则主成分分析统计信息和机器学习工具箱中的功能用于计算数据集的主组件。

[PC，Zscores，PCVARS] = PCA（yeastValues）;

第一个输出，个人电脑，是由主分量组成的矩阵yeastvalues数据。矩阵的第一列是第一个主分量，第二列是第二个主分量，以此类推。第二个输出,ZScores.，包括主成分分数，即主要成分空间中的酵母值的代表。第三个输出，PCVARS.，包含主成分方差，它给出了数据方差中有多少由每个主成分所占的衡量。

很明显，第一个主要成分占模型中的大部分方差。您可以计算每个组件的方差的确切百分比，如下所示。

pcvars./sum(pcvars）* 100

Ans = 79.8316 9.5858 4.0781 2.6486 2.1723 0.9747 0.7089

这意味着几乎90%的方差是由前两个主成分引起的。你可以使用cumsum命令查看差异的累积总和。

cumsum (pcvars. /笔(pcvars) * 100)

Ans = 79.8316 89.4174 93.4955 96.1441 98.3164 99.2911 100.0000

如果您希望在绘制主组件的绘图方面有更多的控制，可以使用散射函数。

图散射(zscores (: 1), zscores (:, 2));包含('第一个主要成分'）;ylabel（'第二个主要成分'）;标题('主成分散点图'）;

创建散点图的替代方法是函数gscatter来自统计学和机器学习工具箱。gscatter创建一个分组散点图，其中每一组的点都有不同的颜色或标记。您可以使用clusterdata，或任何其他聚类功能，以对点进行分组。

图pcclusters = clusterdata(zscores(:，1:2)，“maxclust”，8，'连锁'那“影音”）;gscatter (zscores (: 1) zscores (:, 2), pcclusters)包含('第一个主要成分'）;ylabel（'第二个主要成分'）;标题(“主要成分散射块与彩色簇”）;

自组织映射

如果您有深度学习工具箱™，您可以使用自组织映射(SOM)来聚类数据。

%检查是否安装了深度学习工具箱如果~((存在“selforgmap”），'文件') disp ('此示例的自组织地图部分需要深度学习工具箱。）返回结束

这selforgmap函数创建一个新的SOM网络对象。本例将使用前两个主组件生成SOM。

p = zscores（：，1：2）';net = selforgmap（[4 4]）;

使用默认参数训练网络。

净=火车（网，P）;

使用Plotsom.在数据的散点图上显示网络。请注意，SOM算法使用随机起始点，因此结果将从运行运行时变化。

图绘制(P(1:)、P (2:)“.g”那“markersize”，20）持有在Plotsom（net.iw {1,1}，net.layers {1} .distans）持有从

您可以通过找到数据集中每个点最近的节点来使用SOM分配集群。

距离= dist (P’,net.IW{1}”);[d, cndx] = min(距离,[],2);% CNDX包含集群指数图gscatter (P(1:)、P (2:), cndx);传说从；抓住在Plotsom（net.iw {1,1}，net.layers {1} .distans）;抓住从

关闭所有数据和应用程序。

关闭（“所有”）;删除（cgobj）;删除（h）;nntraintool（'关闭'）;

参考文献

[1] DeRisi J.L.， Iyer, V.R. and Brown, P.O.，“在基因组规模上探索基因表达的代谢和遗传控制”，科学，278(5338):680- 6,1997。