论文部分内容阅读
摘要：随着基因测序技术的不断发展和完善，基因芯片技术逐渐成熟，产生了大量可供分析的基因表达数据，使得研究和比较大量基因成为了可能。聚类分析能够检测出表达谱相似的基因群，将基因表达谱相似程度高的基因划分到同一类中，从而识别出未知功能的基因，是目前研究基因表达数据使用的主要技术之一。本文在对基因表达数据分析方法的研究基础上，采用了一种改进的、基于矩阵变换的基因表达数据层次聚类（NHC）方法，用于改善聚类效果。
关键词：层次聚类；基因芯片；基因表达数据
中图分类号：Q786 文献标识码：A 文章编号：1007-9599 （2012） 24-0046-03
1 引言
生物信息学是20世纪80年代末兴起的新学科，他是涉及数学、生物学、计算机科学和工程学的交叉学科，是一门理论和实践相结合的学科。
基因芯片技术为生物信息学中的一项非常重要的前沿技术，该技术产生于20世纪90年代。随着基因测序技术的发展，果蝇、小鼠、线虫等生物基因测序的完成，以及人类基因组计划的顺利进行，基因序列数据正在以前所未有的速度迅速增长。生物信息学已经进入了一个全新的时代，确定每条基因在整个生物体的发育和活动中所起的调控功能，揭示成千上万基因间的相互作用关系已经为人类基因组计划研究的重点，人类基因组计划从此进入了后基因时代。此时，需要一种能够监测基因行为变化的技术，使研究者能够同时监测不同发育阶段、不同组织、不同健康状态下的基因表达情况。传统的实验技术此时已经无法满足这种需求了，于是基因芯片技術应运而生，为分析基因功能提供了一个十分理想的技术平台。它能同步研究大规模基因在各种不同生理状态下的表达状况，使得同时研究和比较大量基因的功能和特性成为可能。
近年来，数据挖掘作为知识发现的重要手段已经在越来越多的领域得到广泛地应用。在生物信息学领域，也已有多种数据挖掘技术应用于基因表达数据分析。由于聚类分析能够根据基因之间的相似程度将基因表达谱相似程度高的基因划分到同一子类，识别出表达谱相似的基因群，从而使研究者可以通过抑制基因或者结合特征数据来对功能未知的基因进行研究，在基因表达数据分析中得到的广泛应用。已经成为目前分析研究基因表达数据采用的主要技术之一。
2 基因表达数据聚类分析
2.1 基因表达数据的表示
基因芯片作为一种先进的、大规模、高通量的检测技术，能够同步对大规模的基因进行快速检测，即使是对少量试验样品的几次重复的小规模试验，也能得到成千上万条基因的表达数据。例如，只包含5个实验样品和2次重复的小规模实验就可以得到近10万个数据点。所谓基因表达数据，是指通过基因芯片技术得到的样本中各种基因的转录产物mRNA在细胞中的丰度。这些基因表达数据中蕴含着丰富的基因活动信息，为研究基因之间的相关性、基因表达的变化、不同条件对基因活动的影响提供了数据基础。同时，基因表达数据还可以反映细胞当前的活动状态，药物对肿瘤细胞的效果等。这些信息的获得可以帮助我们深入地认识诸如基因功能、发育、衰老等生物过程的本质，因此，基因芯片技术在医学临床诊断、揭示疾病发生机制等方面都得到了应用。
基因表达数据反映的是基因转录产物mRNA在细胞中的丰度，目前基因表达数据主要通过寡核苷酸微阵列（又称基因芯片）和cDNA微阵列两种检测技术获得。基因芯片技术主要分为4个步骤：芯片方阵的构建、样品的制备、生物分子反应和信号的检测。（1）芯片制备，先将玻璃片或硅片进行表面处理，然后将待检测样本的寡核苷酸片段或DNA片段按顺序排列在片芯上形成微阵列。（2）样品制备，由于生物样品的复杂性，一边样品不能直接与芯片反应。可将样品进行生物处理，并用荧光素标记cDNA分子，以提高检测的灵敏度。（3）生物分子反应，该步骤是芯片检测的关键一步，需要选择合适的反应条件使生物分子间反应处于最佳状况中，减少生物分子之间的错配比率。（4）芯片信号检测，最后将芯片置入芯片扫描仪中，通过扫描和计算机分析技术获得样品中基因序列的表达信息。在对基因表达数据进行聚类分析时，所得的数据是在m种不同的样本测得的n个基因的表达水平，这些基因表达数据以m×n维矩阵表示，矩阵元素aij表示第i个基因在第j个样本中测得的基因表达水平。
2.2 基因表达数据的预处理
由于基因表达数据中含有噪声和异常值，因此，在对基因表达数据进行聚类分析之前需要对原始的基因表达数据进行预处理。预处理需要进行填补丢失数据、对不完整数据进行清理等工作，同时，需要根据分析的目从原始数据中筛选出符合条件的数据，并选择合适的数据转换方法，将数据转换成满足分析要求的表示形式。
2.3 聚类分析
目前，聚类分析是基因表达数据分析中最常用的方法之一，常用的基因表达数据聚类分析算法有层次聚类，K-均值聚类，自组织映射。聚类分析通过相似度来衡量基因之间的相似度，并将相似度大的基因划分到同一个子类，使得同一子类中的基因之间的差异尽可能小，同时使得不同子类的基因之间的差异尽可能大。衡量数据间相似程度的方法有很多，比如皮尔逊相关系数，欧氏距离信息熵等，基因表达数据聚类分析中常采用皮尔逊相关系数作为相似度测度。
在基因表达数据分析中，通过对基因表示数据进行聚类分析，将具有相同功能的基因划分到同一个子类中，这些处于同一子类中的基因通常都是表达模式相似的基因，即为共表达基因。因此，我们可以通过聚类分析由已知基因的功能来推断同一子类中未知基因的功能，或者结合临床特性数据来推断某一子类中的基因的功能。
3 基于矩阵变换的层次聚类算法
3.1 相似度
聚类分析中，需要一种方法来表示各个体之间的相似度，相似度的定义是聚类分析的基础，本文采用皮尔逊相关系数作为相似度测度，若两个向量完全相同，则相似度为1，若两个向量没有相似之处，则为0。 3.2 凝聚度
对基于Pearson相关系数度量的簇的凝聚度可用簇内相似度均值来度量
（1）
其中，Ck为第k个簇，uij为Ck簇内第i个对象和第j个对象的余弦相似度，M= n（n- 1）/2，n为Ck簇内对象数。
3.3 熵
本文采用熵来度量聚类效果，熵值表示同一类对象在聚类簇集中的分散程度，处于同一类中的对象熵值越高越分散，相应的聚类质量就越差，其公式为
（2）
其中，pij=mij/mi，mi是簇i中对象的个数，mij是簇i中类j的个数，k是类的个数。簇集合的总熵用每个簇的熵的加权和表示，即：
（3）
其中，m为数据集中对象的总数。
3.4 算法描述
层次聚类又称系统聚类， 具有很好的通用性， 使用非常广泛。层次聚类方法对给定的数据集进行层次的分解，直到某种条件满足为止。具体又可分为聚合式层次聚类，分裂式层次聚类的两种方案。两个过程互为逆过程，它们通过一系列合并或者分裂操作来得到聚类结果。
基于矩阵变换的层次聚類方法通过重排相似度矩阵X的行和列，得到矩阵的划分点，然后通过划分点将重排后的相似度矩阵划分为分块矩阵，比如矩阵Y.
（4）
为了得到划分点d，定义如下函数：
（5）
其中，d为划分点所在的行数（或列数），当
（6）
达达到最大时，d的值即为划分点所在的行数（或列数），然后对子矩阵X11 和X22再进行分解，直到各个子矩阵的凝聚度大于给定值。
本文采用R语言作为算法描述语言，R是一个有着强大统计分析和作图功能的免费、开源软件系统，在统计分析和数据挖掘中有着广泛的应用，在生物信息学中也得到相应的应用。基因矩阵变换的层次聚类的R语言描述如下：
算法：divide
输入：相似度矩阵X[n][n]
输出：分块矩阵X11，X22
divide = function（X）{
for（i in 1：（n-1）{
#找出A中第i行从第i+1列到第n列元素的最大值所在的列k
k=which.max（X[i，（i+1）：n]）
swap_row（X，i+1，k） #交换第i+1行和第k行
swap_col（X，i+1，k） #交换第i+1列和第k列
}
Fd = c（）
for（d in 1：（n-1））{
#计算： Fd=M（ d）（A11）* M（ d）（A22）-M（ d）（A12）* M（ d）（A21）
tmpfd=Md（X，1，1，d）*Md（X，2，2，d）-Md（X，1，2，d）*Md（X，2，1，d）
Fd=c（Fd，tmpfd）
}
d=which.max（Fd） #找出Fd 最大值时d的值
m=dim（X）[1] #X中对象的数目
list（X11=X[1：d，1：d]，X22=X[（d+1）：m，（d+1）：m]） #输出分块矩阵X11，X22
}
算法：cluster
输入：X，Tck（给定值）
输出：聚类结果
cluster=function（X，Tck）{
if（Tc（X）>Tck）{
res=divide（X）
tc11=Tc（res$A11）
tc22=Tc（res$A22）
if（tc11 if（tc22 }else{
output_cluster_result（X）
}
}
4 算法评估
为了测试基于矩阵变换的层次聚类在基因表达数据聚类分析中的效果 ，我们采用了斯坦福大学分子生物系网站上的酵母基因片段（Yeast Cell Cycle）数据集。该数
（下转第58页）
关键词：层次聚类；基因芯片；基因表达数据
中图分类号：Q786 文献标识码：A 文章编号：1007-9599 （2012） 24-0046-03
1 引言
生物信息学是20世纪80年代末兴起的新学科，他是涉及数学、生物学、计算机科学和工程学的交叉学科，是一门理论和实践相结合的学科。
基因芯片技术为生物信息学中的一项非常重要的前沿技术，该技术产生于20世纪90年代。随着基因测序技术的发展，果蝇、小鼠、线虫等生物基因测序的完成，以及人类基因组计划的顺利进行，基因序列数据正在以前所未有的速度迅速增长。生物信息学已经进入了一个全新的时代，确定每条基因在整个生物体的发育和活动中所起的调控功能，揭示成千上万基因间的相互作用关系已经为人类基因组计划研究的重点，人类基因组计划从此进入了后基因时代。此时，需要一种能够监测基因行为变化的技术，使研究者能够同时监测不同发育阶段、不同组织、不同健康状态下的基因表达情况。传统的实验技术此时已经无法满足这种需求了，于是基因芯片技術应运而生，为分析基因功能提供了一个十分理想的技术平台。它能同步研究大规模基因在各种不同生理状态下的表达状况，使得同时研究和比较大量基因的功能和特性成为可能。
近年来，数据挖掘作为知识发现的重要手段已经在越来越多的领域得到广泛地应用。在生物信息学领域，也已有多种数据挖掘技术应用于基因表达数据分析。由于聚类分析能够根据基因之间的相似程度将基因表达谱相似程度高的基因划分到同一子类，识别出表达谱相似的基因群，从而使研究者可以通过抑制基因或者结合特征数据来对功能未知的基因进行研究，在基因表达数据分析中得到的广泛应用。已经成为目前分析研究基因表达数据采用的主要技术之一。
2 基因表达数据聚类分析
2.1 基因表达数据的表示
基因芯片作为一种先进的、大规模、高通量的检测技术，能够同步对大规模的基因进行快速检测，即使是对少量试验样品的几次重复的小规模试验，也能得到成千上万条基因的表达数据。例如，只包含5个实验样品和2次重复的小规模实验就可以得到近10万个数据点。所谓基因表达数据，是指通过基因芯片技术得到的样本中各种基因的转录产物mRNA在细胞中的丰度。这些基因表达数据中蕴含着丰富的基因活动信息，为研究基因之间的相关性、基因表达的变化、不同条件对基因活动的影响提供了数据基础。同时，基因表达数据还可以反映细胞当前的活动状态，药物对肿瘤细胞的效果等。这些信息的获得可以帮助我们深入地认识诸如基因功能、发育、衰老等生物过程的本质，因此，基因芯片技术在医学临床诊断、揭示疾病发生机制等方面都得到了应用。
基因表达数据反映的是基因转录产物mRNA在细胞中的丰度，目前基因表达数据主要通过寡核苷酸微阵列（又称基因芯片）和cDNA微阵列两种检测技术获得。基因芯片技术主要分为4个步骤：芯片方阵的构建、样品的制备、生物分子反应和信号的检测。（1）芯片制备，先将玻璃片或硅片进行表面处理，然后将待检测样本的寡核苷酸片段或DNA片段按顺序排列在片芯上形成微阵列。（2）样品制备，由于生物样品的复杂性，一边样品不能直接与芯片反应。可将样品进行生物处理，并用荧光素标记cDNA分子，以提高检测的灵敏度。（3）生物分子反应，该步骤是芯片检测的关键一步，需要选择合适的反应条件使生物分子间反应处于最佳状况中，减少生物分子之间的错配比率。（4）芯片信号检测，最后将芯片置入芯片扫描仪中，通过扫描和计算机分析技术获得样品中基因序列的表达信息。在对基因表达数据进行聚类分析时，所得的数据是在m种不同的样本测得的n个基因的表达水平，这些基因表达数据以m×n维矩阵表示，矩阵元素aij表示第i个基因在第j个样本中测得的基因表达水平。
2.2 基因表达数据的预处理
由于基因表达数据中含有噪声和异常值，因此，在对基因表达数据进行聚类分析之前需要对原始的基因表达数据进行预处理。预处理需要进行填补丢失数据、对不完整数据进行清理等工作，同时，需要根据分析的目从原始数据中筛选出符合条件的数据，并选择合适的数据转换方法，将数据转换成满足分析要求的表示形式。
2.3 聚类分析
目前，聚类分析是基因表达数据分析中最常用的方法之一，常用的基因表达数据聚类分析算法有层次聚类，K-均值聚类，自组织映射。聚类分析通过相似度来衡量基因之间的相似度，并将相似度大的基因划分到同一个子类，使得同一子类中的基因之间的差异尽可能小，同时使得不同子类的基因之间的差异尽可能大。衡量数据间相似程度的方法有很多，比如皮尔逊相关系数，欧氏距离信息熵等，基因表达数据聚类分析中常采用皮尔逊相关系数作为相似度测度。
在基因表达数据分析中，通过对基因表示数据进行聚类分析，将具有相同功能的基因划分到同一个子类中，这些处于同一子类中的基因通常都是表达模式相似的基因，即为共表达基因。因此，我们可以通过聚类分析由已知基因的功能来推断同一子类中未知基因的功能，或者结合临床特性数据来推断某一子类中的基因的功能。
3 基于矩阵变换的层次聚类算法
3.1 相似度
聚类分析中，需要一种方法来表示各个体之间的相似度，相似度的定义是聚类分析的基础，本文采用皮尔逊相关系数作为相似度测度，若两个向量完全相同，则相似度为1，若两个向量没有相似之处，则为0。 3.2 凝聚度
对基于Pearson相关系数度量的簇的凝聚度可用簇内相似度均值来度量
（1）
其中，Ck为第k个簇，uij为Ck簇内第i个对象和第j个对象的余弦相似度，M= n（n- 1）/2，n为Ck簇内对象数。
3.3 熵
本文采用熵来度量聚类效果，熵值表示同一类对象在聚类簇集中的分散程度，处于同一类中的对象熵值越高越分散，相应的聚类质量就越差，其公式为
（2）
其中，pij=mij/mi，mi是簇i中对象的个数，mij是簇i中类j的个数，k是类的个数。簇集合的总熵用每个簇的熵的加权和表示，即：
（3）
其中，m为数据集中对象的总数。
3.4 算法描述
层次聚类又称系统聚类， 具有很好的通用性， 使用非常广泛。层次聚类方法对给定的数据集进行层次的分解，直到某种条件满足为止。具体又可分为聚合式层次聚类，分裂式层次聚类的两种方案。两个过程互为逆过程，它们通过一系列合并或者分裂操作来得到聚类结果。
基于矩阵变换的层次聚類方法通过重排相似度矩阵X的行和列，得到矩阵的划分点，然后通过划分点将重排后的相似度矩阵划分为分块矩阵，比如矩阵Y.
（4）
为了得到划分点d，定义如下函数：
（5）
其中，d为划分点所在的行数（或列数），当
（6）
达达到最大时，d的值即为划分点所在的行数（或列数），然后对子矩阵X11 和X22再进行分解，直到各个子矩阵的凝聚度大于给定值。
本文采用R语言作为算法描述语言，R是一个有着强大统计分析和作图功能的免费、开源软件系统，在统计分析和数据挖掘中有着广泛的应用，在生物信息学中也得到相应的应用。基因矩阵变换的层次聚类的R语言描述如下：
算法：divide
输入：相似度矩阵X[n][n]
输出：分块矩阵X11，X22
divide = function（X）{
for（i in 1：（n-1）{
#找出A中第i行从第i+1列到第n列元素的最大值所在的列k
k=which.max（X[i，（i+1）：n]）
swap_row（X，i+1，k） #交换第i+1行和第k行
swap_col（X，i+1，k） #交换第i+1列和第k列
}
Fd = c（）
for（d in 1：（n-1））{
#计算： Fd=M（ d）（A11）* M（ d）（A22）-M（ d）（A12）* M（ d）（A21）
tmpfd=Md（X，1，1，d）*Md（X，2，2，d）-Md（X，1，2，d）*Md（X，2，1，d）
Fd=c（Fd，tmpfd）
}
d=which.max（Fd） #找出Fd 最大值时d的值
m=dim（X）[1] #X中对象的数目
list（X11=X[1：d，1：d]，X22=X[（d+1）：m，（d+1）：m]） #输出分块矩阵X11，X22
}
算法：cluster
输入：X，Tck（给定值）
输出：聚类结果
cluster=function（X，Tck）{
if（Tc（X）>Tck）{
res=divide（X）
tc11=Tc（res$A11）
tc22=Tc（res$A22）
if（tc11
output_cluster_result（X）
}
}
4 算法评估
为了测试基于矩阵变换的层次聚类在基因表达数据聚类分析中的效果 ，我们采用了斯坦福大学分子生物系网站上的酵母基因片段（Yeast Cell Cycle）数据集。该数
（下转第58页）