%% 基于主成分分析和聚类分析的基因表达分析
%解决问题：使用神经网络寻找面包酵母的基因表达谱模式
%% 数据
%来源基因表达综合网站 https://www.yeastgenome.org
%加载数据
load yeastdata.mat
%% 使用 numel(genes) 显示数据集中有的基因
%基因表达水平在双峰转换期间的七个时间点测量而得的。变量 times 包含在试验中测量表达水平的时间。变量 genes 包含测量其表达水平的基因的名称。变量 yeastvalues 包含试验中七个时间步的 "VALUE" 数据或 LOG_RAT2N_MEAN，即 CH2DN_MEAN 与 CH1DN_MEAN 之比的 log2。
numel(genes)
%genes 是一个由基因名称组成的元胞数组。
%变量 yeastvalues 的第 15 行包含 ORF YAL054C 的表达水平
genes{15}
%% 过滤基因
%删除'EMPTY'点
%为了更容易找到关注的基因，首先要通过删除未显示任何感兴趣的点的基因表达谱来缩小数据集的大小。共有 6400 个表达谱。您可以使用多种方法将其缩减为包含最重要基因的子集。
%几个标记为 'EMPTY' 的点，这些是数组上的空点，将这些点视为噪声，使用 strcmp 函数找到，并使用索引命令从数据集中删除。
emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)
%删除NaN点
%标记为 NaN 的几个位置，这表明在特定的时间步没有收集到该点的数据。处理这些缺失值的一种方法是使用特定基因在某段时间的均值或中位数来估算它们。此示例使用一种不太严格的方法，即直接丢弃未测量到一个或多个表达水平的任何基因的数据。
%函数 isnan 用于识别具有缺失数据的基因，索引命令用于删除具有缺失数据的基因。
nanIndices = any(isnan(yeastvalues),2);
yeastvalues(nanIndices,:) = [];
genes(nanIndices) = [];
numel(genes)
%过滤基因
%使用 genevarfilter 函数过滤出随时间变化很小的基因。该函数返回与可变基因大小相同的逻辑数组，其中 1 对应于方差大于第 10 个百分位数的 yeastvalues 行，0 对应于低于阈值的 yeastvalues 的行。
mask = genevarfilter(yeastvalues);
yeastvalues = yeastvalues(mask,:);
genes = genes(mask);
numel(genes)
%函数 genelowvalfilter 删除绝对表达值非常低的基因。
[mask, yeastvalues, genes] = ...
genelowvalfilter(yeastvalues,genes,'absval',log2(3));
numel(genes)
%使用 geneentropyfilter 删除熵低的谱的基因
[mask, yeastvalues, genes] = ...
geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)
%% 主成分分析
%归一化数据的标准差和均值使网络可以基于每个输入值的范围以同等权重看待它们。
%主成分分析 (PCA) 是一种有用的方法，可用于降低大型数据集（如微阵列分析中的数据集）的维度。
%该方法分离数据集的主成分，消除那些对数据集变化贡献极小的成分。
%两个设置变量可用于将 mapstd 和 processpca 应用于新数据以保持一致。
%使用 mapstd 对输入向量进行归一化，使它们具有零均值和单位方差。
[x,std_settings] = mapstd(yeastvalues'); % Normalize data
%processpca 是实现 PCA 算法的函数。
%传递给 processpca 的第二个参量是 0.15。这意味着 processpca 会消除那些对数据集总变异贡献不到 15% 的主成分。变量 pc 现在包含 yeastvalues 数据的主成分。
[x,pca_settings] = processpca(x,0.15); % PCA
%使用 scatter 函数可以可视化主成分。
figure(1)
scatter(x(1,:),x(2,:));
xlabel('First Principal Component');
ylabel('Second Principal Component');
title('主成分');
%% 聚类分析：自组织映射
%使用自组织映射 (SOM) 聚类算法对主成分进行聚类
%selforgmap 函数创建一个自组织映射网络，然后使用 train 函数训练该网络。
net = selforgmap([5 3]);
view(net)
%% 训练网络
net = train(net,x);
%使用 plotsompos 在数据的前两个维度的散点图上显示网络。
figure(2)
plotsompos(net,x);
%通过查找数据集中每个点的最近节点，使用 SOM 来分配聚类。
y = net(x);
cluster_indices = vec2ind(y);
%使用 plotsomhits 查看向映射中的每个神经元分配了多少向量。
figure(3)
plotsomhits(net,x);
