51、基于主成分分析和聚类分析的基因表达分析(matlab)

发布于:2024-07-01 ⋅ 阅读:(15) ⋅ 点赞:(0)

1、主成分分析和聚类分析简介

主成分分析(Principal Component Analysis, PCA)和聚类分析(Cluster Analysis)是两种常用的数据分析方法,用于降维和数据分类。

1)主成分分析(PCA)

主成分分析是一种常用的多元统计数据分析方法,旨在通过找到数据中最重要的变量(主成分),将数据从高维空间降维到低维空间,同时保留尽可能多的信息。其基本原理如下:

  • 首先,通过协方差矩阵或相关系数矩阵计算数据间的相关性;
  • 然后,通过特征值分解或奇异值分解等方法,找到数据中最重要的主成分;
  • 最后,使用主成分来表示原始数据,实现降维。

PCA常用于特征提取、数据可视化和降维处理,帮助揭示数据中的模式和结构,发现数据之间的关系。

2)聚类分析(Cluster Analysis)

聚类分析是一种无监督学习技术,旨在将数据对象组织成类或簇,使得同一簇内的数据对象相互之间相似,而不同簇之间的数据对象差异较大。其基本原理如下:

  • 首先,通过定义一个相似性度量标准(如欧氏距离、余弦相似度等),计算数据对象之间的相似性;
  • 然后,将数据对象划分为若干个簇,使得同一簇内的数据对象之间相似度高,不同簇之间相似度低;
  • 最后,评估聚类结果的质量和有效性,调整聚类算法的参数来优化聚类效果。

聚类分析常用于数据分类、模式识别和群体分析等领域,帮助发现数据对象之间的隐藏结构和规律。

3)总结

主成分分析主要用于降维和特征提取,聚类分析用于数据分类和群体分析。这两种方法在数据分析、机器学习和模式识别等领域具有广泛的应用,有助于理解和挖掘数据背后的规律和关联。

2、基于主成分分析和聚类分析的基因表达分析说明 

解决问题

使用神经网络寻找面包酵母的基因表达谱模式

3、实验数据

数据来源

来源基因表达综合网站 https://www.yeastgenome.org

加载数据

代码

load yeastdata.mat

 4、使用 numel(genes) 显示数据集中有的基因

1)说明

基因表达水平在双峰转换期间的七个时间点测量而得的。变量 times 包含在试验中测量表达水平的时间。变量 genes 包含测量其表达水平的基因的名称。变量 yeastvalues 包含试验中七个时间步的 "VALUE" 数据或 LOG_RAT2N_MEAN,即 CH2DN_MEAN 与 CH1DN_MEAN 之比的 log2。

代码

numel(genes)

ans =

        6400

 2)genes 是一个由基因名称组成的元胞数组。

说明

变量 yeastvalues 的第 15 行包含 ORF YAL054C 的表达水平

代码

genes{15}

ans =

    'YAL054C'

 5、过滤基因

1)删除'EMPTY'点

说明:为了更容易找到关注的基因,首先要通过删除未显示任何感兴趣的点的基因表达谱来缩小数据集的大小。共有 6400 个表达谱。您可以使用多种方法将其缩减为包含最重要基因的子集。
几个标记为 'EMPTY' 的点,这些是数组上的空点,将这些点视为噪声,使用 strcmp 函数找到,并使用索引命令从数据集中删除。

代码

emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)

ans =

        6314

 2)删除NaN点

说明:标记为 NaN 的几个位置,这表明在特定的时间步没有收集到该点的数据。处理这些缺失值的一种方法是使用特定基因在某段时间的均值或中位数来估算它们。此示例使用一种不太严格的方法,即直接丢弃未测量到一个或多个表达水平的任何基因的数据。
函数 isnan 用于识别具有缺失数据的基因,索引命令用于删除具有缺失数据的基因。

nanIndices = any(isnan(yeastvalues),2);
yeastvalues(nanIndices,:) = [];
genes(nanIndices) = [];
numel(genes)

ans =

        6276

3) 过滤基因

说明:使用 genevarfilter 函数过滤出随时间变化很小的基因。该函数返回与可变基因大小相同的逻辑数组,其中 1 对应于方差大于第 10 个百分位数的 yeastvalues 行,0 对应于低于阈值的 yeastvalues 的行。


代码

mask = genevarfilter(yeastvalues);
yeastvalues = yeastvalues(mask,:);
genes = genes(mask);
numel(genes)

ans =

        5648

 4)函数 genelowvalfilter 删除绝对表达值非常低的基因。

 代码

[mask, yeastvalues, genes] = ...
    genelowvalfilter(yeastvalues,genes,'absval',log2(3));
numel(genes)

ans =

   822

5)使用 geneentropyfilter 删除熵低的谱的基因

代码
 

[mask, yeastvalues, genes] = ...
    geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)

ans =

   614

6、主成分分析

1)说明

归一化数据的标准差和均值使网络可以基于每个输入值的范围以同等权重看待它们。
主成分分析 (PCA) 是一种有用的方法,可用于降低大型数据集(如微阵列分析中的数据集)的维度。
该方法分离数据集的主成分,消除那些对数据集变化贡献极小的成分。
两个设置变量可用于将 mapstd 和 processpca 应用于新数据以保持一致。

使用 mapstd 对输入向量进行归一化,使它们具有零均值和单位方差。

代码

[x,std_settings] = mapstd(yeastvalues');  % Normalize data

2)processpca 是实现 PCA 算法的函数

说明:传递给 processpca 的第二个参量是 0.15。这意味着 processpca 会消除那些对数据集总变异贡献不到 15% 的主成分。变量 pc 现在包含 yeastvalues 数据的主成分。

代码

[x,pca_settings] = processpca(x,0.15);    % PCA

3)使用 scatter 函数可以可视化主成分

代码

figure(1)
scatter(x(1,:),x(2,:));
xlabel('First Principal Component');
ylabel('Second Principal Component');
title('主成分');

视图效果

75d6435dc47e459386665fae166b9e9b.png

 7、聚类分析:自组织映射

1)说明

使用自组织映射 (SOM) 聚类算法对主成分进行聚类

selforgmap 函数创建一个自组织映射网络,然后使用 train 函数训练该网络。

代码

net = selforgmap([5 3]);
view(net)

视图效果

a22f6f99271e46dab94c1aa3f0e01039.png

 8、训练网络

代码

net = train(net,x);

1)使用 plotsompos 在数据的前两个维度的散点图上显示网络。

 代码

figure(2)
plotsompos(net,x);

视图效果

308d1c645c7d4733b4b17a64773934c7.png

2) 通过查找数据集中每个点的最近节点,使用 SOM 来分配聚类。

代码

y = net(x);
cluster_indices = vec2ind(y);

视图效果

fe548af6f5fe4eaa9c7c54c29029fd35.png

3)使用 plotsomhits 查看向映射中的每个神经元分配了多少向量 

代码

figure(3)
plotsomhits(net,x);

视图效果

 65eba89cfd9b4830b10865e0241f749f.png

 9、总结

基于主成分分析(PCA)和聚类分析的基因表达数据分析在生物信息学和生物医学研究中具有重要意义。下面是基于MATLAB进行基因表达分析的主要步骤和总结:

  1. 数据准备:首先,需要获取基因表达数据集,包含多个基因在不同条件下的表达值。数据应进行预处理,如去除噪声、归一化处理等,以保证数据质量和可靠性。

  2. 主成分分析(PCA): a. 使用MATLAB的PCA函数对预处理后的基因表达数据进行主成分分析。 b. 通过计算数据的协方差矩阵或相关系数矩阵,找到数据中最重要的主成分,实现数据的降维和特征提取。 c. 利用MATLAB的函数进行数据可视化,绘制主成分贡献度图和主成分载荷图,帮助解释主成分的含义和解释数据方差结构。

  3. 聚类分析: a. 使用MATLAB的聚类分析函数,如kmeans或者hierarchical clustering等,对基因表达数据进行聚类。 b. 通过定义聚类的簇数、相似性度量标准和聚类算法等参数,将基因样本划分为若干个类别。 c. 利用MATLAB的聚类结果可视化工具,如热图或散点图,展示不同基因样本之间的相似性和差异性,帮助对基因表达数据进行分类和分析。

  4. 结果解释:根据主成分分析和聚类分析的结果,可以发现基因表达数据中的模式、结构和规律,识别不同基因之间的相关性和差异性,有助于推断基因调控网络和功能通路。进一步结合生物实验和文献信息,可以揭示基因在生物学过程中的重要作用和生物学意义。

通过MATLAB对基因表达数据进行主成分分析和聚类分析,可以深入理解基因表达谱的特征和变化,为基因功能研究、疾病诊断和治疗等领域提供重要支持和指导。

10、源代码

代码

%% 基于主成分分析和聚类分析的基因表达分析
%解决问题:使用神经网络寻找面包酵母的基因表达谱模式
%% 数据
%来源基因表达综合网站 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);