基因表达分析
此示例说明如何使用神经网络寻找面包酵母的基因表达谱模式。
问题:面包酵母(酿酒酵母)中的基因表达分析
这一分析的目的是了解酿酒酵母(通常称为面包酵母或啤酒酵母)的基因表达。这种真菌用于烤面包和基于葡萄发酵制造葡萄酒。
当酿酒酵母引入富含葡萄糖的培养基时,葡萄糖可以转化为乙醇。最初,酵母通过一种称为“发酵”的代谢过程将葡萄糖转化为乙醇。然而,一旦葡萄糖供应耗尽,酵母就从葡萄糖的厌氧发酵转向乙醇的有氧呼吸。此过程称为双峰转换。这需要着重关注,因为它伴随着基因表达的重大变化。
本示例使用 DNA 微阵列数据来研究酿酒酵母在双峰转换过程中几乎所有基因的瞬时基因表达。
您需要 Bioinformatics Toolbox™ 才能运行此示例。
if ~nnDependency.bioInfoAvailable errordlg('This example requires Bioinformatics Toolbox.'); return; end
数据
此示例使用的数据来源:DeRisi, JL, Iyer, VR, Brown, PO."Exploring the metabolic and genetic control of gene expression on a genomic scale."Science.1997 Oct 24;278(5338):680-6.PMID:9381177
完整的数据集可以从基因表达综合网站 https://www.yeastgenome.org 下载
首先将数据加载到 MATLAB® 中。
load yeastdata.mat
基因表达水平在双峰转换期间的七个时间点测量而得的。变量 times
包含在试验中测量表达水平的时间。变量 genes
包含测量其表达水平的基因的名称。变量 yeastvalues
包含试验中七个时间步的 "VALUE" 数据或 LOG_RAT2N_MEAN,即 CH2DN_MEAN 与 CH1DN_MEAN 之比的 log2。
要了解数据的大小,您可以使用 numel(genes)
来显示数据集中有多少基因。
numel(genes)
ans = 6400
genes 是一个由基因名称组成的元胞数组。您可以使用 MATLAB 元胞数组索引来访问条目:
genes{15}
ans = 'YAL054C'
这表明变量 yeastvalues
的第 15 行包含 ORF YAL054C
的表达水平。
过滤基因
数据集相当大,许多信息对应于在试验过程中没有显示任何要关注的变化的基因。为了更容易找到关注的基因,首先要通过删除未显示任何感兴趣的点的基因表达谱来缩小数据集的大小。共有 6400 个表达谱。您可以使用多种方法将其缩减为包含最重要基因的子集。
如果您浏览基因列表,会看到几个标记为 'EMPTY' 的点。这些是数组上的空点,虽然它们可能有相关联的数据,但出于此示例的目的,您可以将这些点视为噪声。这些点可以使用 strcmp
函数找到,并使用索引命令从数据集中删除。
emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)
ans = 6314
在 yeastvalues 数据中,您还会看到表达式级别标记为 NaN 的几个位置。这表明在特定的时间步没有收集到该点的数据。处理这些缺失值的一种方法是使用特定基因在某段时间的均值或中位数来估算它们。此示例使用一种不太严格的方法,即直接丢弃未测量到一个或多个表达水平的任何基因的数据。
函数 isnan
用于识别具有缺失数据的基因,索引命令用于删除具有缺失数据的基因。
nanIndices = any(isnan(yeastvalues),2); yeastvalues(nanIndices,:) = []; genes(nanIndices) = []; numel(genes)
ans = 6276
如果您绘制所有其余谱的表达廓,会看到大多数谱是平坦的,并且与其他谱没有显著差异。此平坦数据显然是有用的,因为它表明与这些谱相关联的基因没有受到双峰转换的显著影响;然而,在此示例中,您感兴趣的是伴随双峰转换而在表达中发生巨大变化的基因。您可以使用 Bioinformatics Toolbox™ 中的过滤函数来删除各种其表达谱不能提供基因受代谢变化影响的有用信息的基因。
您可以使用 genevarfilter
函数过滤出随时间变化很小的基因。该函数返回与可变基因大小相同的逻辑数组,其中 1 对应于方差大于第 10 个百分位数的 yeastvalues 行,0 对应于低于阈值的 yeastvalues 的行。
mask = genevarfilter(yeastvalues);
% Use the mask as an index into the values to remove the filtered genes.
yeastvalues = yeastvalues(mask,:);
genes = genes(mask);
numel(genes)
ans = 5648
函数 genelowvalfilter
删除绝对表达值非常低的基因。请注意,基因过滤函数也可以自动计算过滤后的数据和名称。
[mask, yeastvalues, genes] = ... genelowvalfilter(yeastvalues,genes,'absval',log2(3)); numel(genes)
ans = 822
使用 geneentropyfilter
删除熵低的谱的基因:
[mask, yeastvalues, genes] = ... geneentropyfilter(yeastvalues,genes,'prctile',15); numel(genes)
ans = 614
主成分分析
现在您已获得可处理的基因列表,可以寻找这些谱之间的关系。
归一化数据的标准差和均值使网络可以基于每个输入值的范围以同等权重看待它们。
主成分分析 (PCA) 是一种有用的方法,可用于降低大型数据集(如微阵列分析中的数据集)的维度。该方法分离数据集的主成分,消除那些对数据集变化贡献极小的成分。
这两个设置变量可用于将 mapstd
和 processpca
应用于新数据以保持一致。
[x,std_settings] = mapstd(yeastvalues'); % Normalize data [x,pca_settings] = processpca(x,0.15); % PCA
首先使用 mapstd
对输入向量进行归一化,使它们具有零均值和单位方差。processpca
是实现 PCA 算法的函数。传递给 processpca
的第二个参量是 0.15。这意味着 processpca
会消除那些对数据集总变异贡献不到 15% 的主成分。变量 pc
现在包含 yeastvalues 数据的主成分。
使用 scatter
函数可以可视化主成分。
figure scatter(x(1,:),x(2,:)); xlabel('First Principal Component'); ylabel('Second Principal Component'); title('Principal Component Scatter Plot');
聚类分析:自组织映射
现在可以使用自组织映射 (SOM) 聚类算法对主成分进行聚类。
selforgmap
函数创建一个自组织映射网络,然后可以使用 train
函数训练该网络。
输入大小为 0,因为网络尚未配置成与我们的输入数据相匹配。这将在训练网络时进行。
net = selforgmap([5 3]); view(net)
现在网络已准备就绪,可以开始训练。
神经网络训练工具显示正在接受训练的网络和用于训练该网络的算法。它还显示训练过程中的训练状态以及停止训练的条件(以绿色突出显示)。
底部的按钮用于打开有用的绘图,这些图可以在训练期间和训练后打开。算法名称和绘图按钮旁边的链接可打开有关这些主题的文档。
net = train(net,x);
使用 plotsompos
在数据的前两个维度的散点图上显示网络。
figure plotsompos(net,x);
您可以通过查找数据集中每个点的最近节点,使用 SOM 来分配聚类。
y = net(x); cluster_indices = vec2ind(y);
使用 plotsomhits
查看向映射中的每个神经元分配了多少向量。
figure plotsomhits(net,x);
您还可以使用 Statistics and Machine Learning Toolbox™ 中提供的其他聚类算法(如分层聚类和 K 均值)进行聚类分析。
术语表
ORF - 开放阅读框 (ORF) 是基因序列的一部分,包含一个碱基序列,不被终止序列打断,可能对蛋白质编码。