基因表达谱分析

探索数据集

本例使用了DeRisi等人于1997年发表的酵母基因表达微阵列研究数据。作者使用DNA微阵列研究了几乎所有基因的时间基因表达酿酒酵母在从发酵到呼吸的代谢转变过程中。在语音转换期间的7个时间点测量表达水平。完整的数据集可以从Gene Expression Omnibus网站下载，https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28．

的MAT-fileyeastdata.mat包含表达式值(log2)的比率CH2DN_MEAN和CH1DN_MEAN)，从实验中的七个时间步骤，基因的名称，以及测量表达水平的时间序列。

负载yeastdata.mat

来了解您可以使用的数据的大小元素个数(基因)显示数据集中包含了多少个基因。

元素个数(基因)

Ans = 6400

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

基因{15}

ans = 'YAL054C'

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

Plot (times, yeastvalues(15，:)) xlabel(的时间(小时));ylabel (“Log2相对表达水平”);

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

Plot (times, 2.^yeastvalues(15，:)) xlabel(的时间(小时));ylabel (“相对表达水平”);

与这种ORF相关的基因ACS1在双氧转换过程中似乎被强烈上调。你可以通过在同一图上绘制多条线来比较这个基因的表达与其他基因的表达。

持有在Plot (times, 2.^yeastvalues(16:26，:)') xlabel(的时间(小时));ylabel (“相对表达水平”);标题(“剖面表达水平”);

筛选基因

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

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

emptySpots = strcmp(“空”,基因);yeastvalues(emptySpots，:) = [];基因(emptySpots) = [];元素个数(基因)

Ans = 6314

在数据集中还有几个地方，表达式级别被标记为南．这表明在特定的时间步没有为这个点收集数据。处理这些缺失值的一种方法是使用特定基因随时间的平均值或中位数数据来推算它们。这个例子使用了一种不太严格的方法，简单地扔掉了没有测量到一个或多个表达水平的任何基因的数据。这个函数isnan用于识别缺失数据的基因，索引命令用于删除缺失数据的基因。

nanIndices = any(isnan(yeastvalues)，2);yeastvalues(nanIndices，:) = [];基因(nanIndices) = [];元素个数(基因)

得分= 6276

如果您要绘制所有剩余概要文件的表达式概要文件，您将看到大多数概要文件是平坦的，并且与其他概要文件没有显著的不同。这个平面数据显然是有用的，因为它表明与这些谱相关的基因不受双氧转移的显著影响;然而，在本例中，您感兴趣的是伴随双音转换的表达大变化的基因。您可以使用Bioinformatics Toolbox™中的过滤功能来删除具有各种类型的配置文件的基因，这些配置文件不能提供有关受代谢变化影响的基因的有用信息。

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

Mask = genevarfilter(酵母值);Yeastvalues = Yeastvalues (mask，:);基因=基因(mask);元素个数(基因)

Ans = 5648

这个函数genelowvalfilter去除绝对表达值非常低的基因。请注意，这些过滤函数也可以自动计算过滤后的数据和名称，因此没有必要使用掩码对原始数据进行索引。

[mask,yeastvalues,genes] = genelowvalfilter(酵母值，基因)“absval”log2 (3));元素个数(基因)

答案= 822

最后，您可以使用该函数geneentropyfilter去除那些熵值较低的基因，例如在数据的第15百分位的熵值水平。

[mask,yeastvalues,genes] = geneentropyfilter(酵母值，基因)“prctile”15);元素个数(基因)

Ans = 614

聚类分析

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

cordist = pdist(酵母值，“相关系数”);clusterTree = linkage(corrDist，“平均”);

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

cluster = cluster(clusterTree，“maxclust”16);

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

数字为C = 1:16 subplot(4,4, C);Plot (times,yeastvalues((clusters == c)，:)');轴紧结束sgtitle (“配置文件的分层聚类”);

统计和机器学习工具箱也有一个K-means聚类函数。同样，找到了16个聚类，但由于算法不同，这些聚类不一定与通过分层聚类找到的聚类相同。

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

rng (“默认”);[cidx, ctrs] = kmeans(yeastvalues,16)“距离”，“相关系数”，“代表”5,“disp”，“最后一次”);数字为C = 1:16 subplot(4,4, C);Plot (times,yeastvalues((cidx == c)，:)');轴紧结束sgtitle (“k -均值聚类”);

重复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 -均值聚类”);

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

cgObj = clustergram(yeastvalues(:，2:end)，“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命令查看方差的累积和。

sum(pcvars./sum(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, hsv(8)包含(“第一主成分”);ylabel (“第二主成分”);标题(带彩色聚类的主成分散点图);

自组织映射

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

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

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

P = zscores(:，1:2)';Net = selforgmap([4 4]);

使用默认参数训练网络。

net = train(net,P);

使用plotsom在数据的散点图上显示网络。请注意，SOM算法使用随机起始点，因此每次运行的结果都会有所不同。

图绘制(P(1:)、P (2:)“.g”，“markersize”, 20)在net.layers plotsom (net.iw {1}, {1} .distances)从

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

距离= dist(P'，net. iw{1}');[d,cndx] = min(距离，[]，2);% CNDX包含集群索引图gscatter (P(1:)、P (2:), cndx, hsv(元素个数(独特(cndx))));传说从;持有在net.layers plotsom (net.iw {1}, {1} .distances);持有从

关闭所有图形。

关闭(“所有”);删除(cgObj);删除(h);

参考文献

[10]王晓明，王晓明，王晓明，“基因表达的遗传调控研究”，中国生物医学工程学报，2009,31(4):387 - 398。