蛋白质冷冻电镜图像的分类算法论文大学论文.doc

论文

蛋白质冷冻电镜图像的分类算法研究

Classification of Protein Cryo-EM Images

蛋白质冷冻电镜图像的分类算法研究
摘要
本篇论文的主要工作是对模拟的蛋白质二维投影图像进行分类。这个课题源于冷冻电镜蛋白质的三维重构，对蛋白质投影图像的分类平均是蛋白质三维重构的重要前提。我们首先根据蛋白质冷冻电镜图像的形成原理，模拟产生不同投影角度的蛋白质二维图像。然后提取图像的平移旋转不变量，为之后的分类做准备。由于模拟产生的图像中蛋白质的投影图像都是居中的，所以可近似忽略其平移变量。旋转不变量是通过对图像做极坐标变换和傅里叶变换而提取到的。接着运用k-medoids算法,一种典型的基于划分的聚类方法,对处理后的图像进行分类。最后采用了由Fowlkes和Mallows提出的一种评估聚类结果的方法，用来验证本课题的分类结果。在结果分析中，我们给输入图像分别加入不同的高斯白噪声，然后对比不同信噪比图片的聚类结果，得出该算法能够处理的信噪比的范围。此外，我们还对比了在输入图像总数不同或类别数不同时的分类结果。以此来全面评估k-medoids算法对蛋白质投影图像的分类效果。

关键词：蛋白质，图像，分类，旋转不变性，k-medoids算法
CLASSIFICATION OF PROTEIN CRYO-EM IMAGES
Abstract
In this paper our main work is to classify simulated protein two-dimensional projection images. This problem is related with three-dimensional reconstruction of protein structure from cryo-EM data, the process of classification and averaging of the projected image is an important prerequisite for three-dimensional reconstruction of protein. First, we produce simulated projection images of protein from various directions, in accordance with the actual collection of cryo-EM images . Then extract the translation and rotation invariant from these images, to be used as features for the following classification. The projected images are centered when we produce simulation images, so it can be approximated by ignoring its translational variables. We use the polar coordinate transformation and Fourier transform to obtain rotational invariances of images. Next use classical clustering methods, namely k-medoids algorithm, to classify the images. Finally, in order to assess the quality of the classification results, we use the index proposed by the Fowlkes and Mallows. In the experiments, we added different Gaussian white noise to the input image, and compared the results of clustering with different signal to noise ratio,. In order to fully assess the k-medoids algorithm for projected image classification, we also assessed the quality of the classification by comparing the results of different number of input images and different number of clasters.

Keywords: protein, images, classification, rotational invariance, k-medoids algorithm
目录蛋白质冷冻电镜图像的分类算法研究 I
摘要 I
Abstract II
目录 III
第1章绪论 1
1.1 课题背景和意义 1
1.2 国内外研究现状 1
1.3 论文组织结构 2
第2章 聚类 3
2.1 聚类的定义 3
2.1.1 聚类定义 3
2.1.2 相似性度量 3
2.2 聚类的步骤 4
2.3 聚类算法的分类 4
2.3.1 基于层次的方法 4
2.3.2 基于分割的方法 5
2.3.3 基于密度的方法 5
2.3.4 基于网格的方法 5
2.3.5 基于模型的方法 5
2.4 小结 5
第3章 k-medoids 算法 6
3.1 k-medoids 算法的基本原理 6
3.2 k-medoids 算法的流程 6
第4章 对蛋白质冷冻电镜下图像的分类处理 8
4.1 模拟产生蛋白质冷冻电镜下的图像 8
4.2 对冷冻电镜下蛋白质图片的分类 8
4.2.1 提取特征 8
4.2.2 相似性度量 10
4.2.3 k-medoids 算法聚类过程 10
4.3 对聚类结果的分析评估 10
第五章 总结与展望 13

IV

毕业论文 蛋白质冷冻电镜图像的分类算法研究
第1章绪论
1.1 课题背景和意义
蛋白质作为生命体的主要组成成分，生命活动的主要执行者，是生物有机体生命现象的直接体现者，而仅凭目前已知的蛋白质根本无法阐明各种复杂的生命活动过程。而以蛋白质为主体的生物大分子的功能主要取决于它们的三维结构。因此，研究蛋白质结构对于了解各种蛋白质的性质和功能就非常重要。
研究蛋白质结构的方法主要有三种：X射线晶体学(X-ray)、核磁共振技术(NMR)和冷冻电子显微镜(Cryo-EM)。目前，由于前两种技术具有不易测定大分子复合物结构等局限性，冷冻电镜三维重构已经成为研究生物大分子结构和功能的强有力手段。冷冻电镜三维重构首先利用快速冷冻技术对蛋白质样品进行冷冻固定，然后利用冷冻电镜对样品进行电子成像，接着将底片数字化，对数字化的图像进行二维图像分析——选点、分类、校正和平均，最后完成该蛋白质的三维重构[1]。
由于蛋白质结构的整个三维重构过程极其复杂且会耗费大量的时间，本课题就截取过程中的关键步骤，即二维图像的分类过程。我们期望改进图像的分类算法，以取得良好的分类效果，致使下一步三维重构时，能够提高重构模型的正确性，并降低蛋白质三维重构的难度，减少重构消耗的时间。
1.2 国内外研究现状
目前，对于单粒子重构过程中的二维图像分类平均，主要有两种方法，IMAGIC和SPIDER。IMAGIC使用多元变量统计分析法(MSA)和多参考对准法(MRA)进行二维图像分类、平均。MSA先对较大的图像集进行压缩和去噪，然后用基于分层的分类方法对图像集进行高效的分类。聚类后的图像为MRA进行类平均操作提供参考。因为原本是同一角度的投影图像，有时只是经过了小幅度的平移和旋转，所以为了实现不同角度的投影图像的分类，我们需要预先提取一些不变特征，例如，用自相关函数(ACF)或双自相关函数(DACF)进行提取。SPIDER先用无参考对准法(RFA)，接着用旋转不变特征向量作为K-means聚类的输入，对二维图像分类。RFA试着在全局范围内去对准图像。这种优化方法主要是在所有图像与其平均值的偏差平方总和最小(即最小方差)时，找到它们的旋转和变换的对准参数[11]。另外，如今很多三维重构的软件中也包含二维图像的分类平均过程，不同的软件采用了各种不同的算法。
二维图像分类平均过程大部分采用RFA，但是当图像来自大量不同角度的投影时，该方法就会出错。而且如果要实现全局对准，这样就需要所有图片两两比较一次，对准就需要n(n-1)/2次操作。这种操作耗时且没必要，因为大部分时间都用在完全不同的图片的对准上了[11]。
1.3 论文组织结构
本课题主要模拟蛋白质冷冻电镜图像三维重构中的二维图像的分类过程。采用了聚类中的k-medoids算法，对图像分类，并分析了该算法的分类效果，统计评估了分类结果。本文各章安排如下：
第一章 为绪论，简要介绍了蛋白质冷冻电镜下三维重构的背景，并说明了选题的意义。接着还详述了现如今二维图像分类平均过程在国内外的研究现状。
第二章 介绍了聚类的相关概念及其分类。
第三章 介绍了k-medoids 算法的基本原理。
第四章 介绍对蛋白质冷冻电镜下图像的整个处理过程，阐述了整个过程的步骤以及相关的原理。
第五章 为总结与展望。
第2章 聚类
2.1 聚类的定义
2.1.1 聚类定义
聚类是根据所提取的样本特征的相似性，将输入的数据集划分为几种不同的子数据集，分类的结果是使得同一类中的数据相似度较高，而不同类间的数据差别较大。特征提取就是将样本数据转化为若干个特征的过程。从数学上讲，特征提取相当于把一个物理模式变为一个随机向量。提取的特征必须能够唯一标识该样本，且一般具有伸缩、旋转、平移不变性。数据间的相似性一般用最短距离或相似系数来度量。
聚类与分类不同。聚类属于无监督式学习，即开始对数据集分类时，并没有数据相关的类信息，只是从所有数据中随机选取样本作为类参考，然后根据相似性度量规则将数据划分入各个类中。相反，分类属于监督式学习，分类模型中存在数据样本，且这些数据的类标号已知，分类就是从训练样本集中提取出分类的规则，用于标识其他类标号未知的数据。
2.1.2 相似性度量
正如2.1.1中所述，聚类的依据是数据间的相似性，即从样本中所提取的特征的相似程度。通常描述样本相似性的方法主要有两种：距离和相似系数。
(1) 距离：假设用n个特征变量来描述样本，那么该样本X即为n维特征向量，这样我们就可以将样本看成n维欧式空间中的一个点，然后选用某种距离公式，求与其他各点间的距离。距离越小，则表示两样本间相似性越高，反之，则越低。一般采用欧几里得距离作为相似性度量。
(2) 相关系数：这是衡量两个随机变量间线性相关程度的指标。X、Y表示两个样本，其相似性计算公式如下：

若两个变量的相关系数越接近1，则表明它们越相似；反之，越接近0，则差异越大。
2.2 聚类的步骤
(1) 特征提取。这一步的输入是原始样本数据，然后根据具体情况决定提取哪些特征来刻画样本的属性和结构。特征提取的输出是一个向量或矩阵。有时，提取的特征变量过多，不利于后续的分析处理，这时可以进行降维处理，降低空间维数，方便运算。另外，提取的特征变量之间，必须是相互独立的，所以可用主成成分分析法消除特征变量间的相关性。
(2) 相似性度量。这一步主要是对提取出的特征向量进行相似性计算，它将上一步输出作为输入，用选定的方法计算，如欧几里得距离或相似性系数法。计算后得到一个距离矩阵或相似系数矩阵。
(3) 聚类。这一步选用适当的聚类算法，并利用上一步的输出，对样本聚类，最后得到分类结果，一般输出是一个标有类标号的一维向量。
(4) 结果评估。即评估并验证该聚类结果的准确度。
2.3 聚类算法的分类
聚类算法主要分为基于层次的方法、基于分割的方法、基于密度的方法、基于网格的方法和基于模型的方法。
2.3.1 基于层次的方法
基于层次的聚类方法是将数据按层次划分为不同的类，同时形成一个相应的树状图。该方法又可以分为两类，即聚合和分裂。聚合的方法又称为自下而上的方法，一开始将每一个样本作为单独的类，然后自下而上逐渐合并相近的样本或类，直到所有样本合并为一个，或满足某个终止条件。分裂的方法恰好与前者相反，它又称为自上而下的方法，即先将所有的数据集看作一个类，之后开始迭代，把一个大类分为两个差别很大的小类，直到每个样本单独为一类，或者满足某一终止条件。常见的基于层次的聚类方法有BIRCH，CURE和CHAMELEON等。
基于层次的聚类，在类的合并和分裂过程中，需要估计类间的相似性。通常以距离来衡量。常用的方法有：最短距离法，最长距离法，中间距离法，类平均法等。
2.3.2 基于分割的方法
基于分割的聚类方法，处理思路如下，给n个样本，用户输入要划分的类数目k，且k≤n。开始将n个样本随机划分为k个类，然后重新选取各个类的中心以调整划分，重复前两步，最后在整个分类过程收敛时结束操作。目前，这种方法的典型代表是k-means算法和k-medoids算法。k-means与k-medoids主要的区别是，k-means的聚类中心是类中所有对象的均值或中心点，而k-medoids的聚类中心只是类中的某个对象，该对象与类中其他各个对象距离之和最小。另外，常见的k-medoids算法有PAM，CLARA和CLARANS。
2.3.3 基于密度的方法
很多分割方法是基于对象间的距离展开聚类操作的，但这样只能发现球状型的簇，而很难发现任意其他形状的簇，因此提出基于密度的聚类方法。其主要的思想是：从数据对象的分布密度出发，当临近区域的密度（数据或对象的数目）超出某个阈值时，就可以继续聚类。这样不仅可以发现任意形状的类，还可以有效去除噪声。常见的基于密度的算法有DBSCAN，DBCLASD，OPTICS和DENCLUE等。
2.3.4 基于网格的方法
基于网络的方法是把空间量化成为有限个单元，形成一个网络结构。所有的聚类操作都在这个网络结构上进行。这种方法的优点是处理速度快，其时间复杂度只与网格单元的数目有关，而与输入数据对象的数目无关。典型算法代表有STING，WaveCluster和CLIQUE等。
2.3.5 基于模型的方法
基于模型的聚类方法为每一个簇假定一个模型，然后寻找能够很好的满足这个模型的数据集。该算法主要分两类：统计学方法和网络神经方法。其中，统计学方法有COBWEB算法，网络神经方法有SOM算法。
2.4 小结
本章开头详细讲述了聚类的定义，并对比了聚类与分类两者间的细微差别，然后介绍了聚类中常见的相似性度量的方法。接着，概括了聚类的一般步骤。最后，简单介绍了各种聚类方法，并列举其对应的典型算法。这些内容都为之后详细介绍k-medoids算法做了充足的基础理论准备。
第3章 k-medoids 算法
K-medoids算法是一种基于划分的聚类方法。通常该算法先将输入数据集映射到d维欧几里得空间中，以欧式距离为相似性度量指标，两个数据对象间的距离越小，就越相似。该算法的目标是找到k个聚类中心点，并划分数据集，使得所有数据到其最近中心点的距离的平方之和最小。这种算法是一种局部最优收敛的启发式算法。相对于k-means 算法，当处理数据中存在“噪声”或孤立点数据时，k-medoids 算法更为健壮，因为k-medoids 算法的“中心点”不像k-means中的“均值”那样易受极端数据的影响。但k-medoids算法找聚类中心点时消耗较大。
3.1 k-medoids 算法的基本原理
假设使用‖p,q‖表示点p和点q之间的欧氏距离，那么d维欧几里得空间中的集合Q就可表示为： ，其中，p为集合Q的中心点。Q的中心点C(Q)必须满足集合Q中各点到中心点距离的平方之和最小，可表示为：
。
K-medoids 算法想要找到k个类的中心点集合M={}，并将输入的数据集划分为k类，就必须满足所有点到其最近的中心点的距离的平方之和最小，即值最小。
, (i=1,2,…,k)
3.2 k-medoids 算法的流程
K-medoids 算法的处理流程如下：
输入：聚类的类别数k和包含n个对象的数据集。
输出：是含有类标识的k个聚类。
(1) 从n个对象中任意选取k个对象作为初始的聚类中心点。
(2) 分别计算每个对象到k个中心点的距离，然后将其划分入距离最近的类中。
(3) 所有对象分类完成后，重新调整每个类的中心点，使每个类的中心点与类中其他点的距离的平方之和最小。
(4) 若J变化，则跳至(2)；否则，表示J达到最小值，跳转至(5)。
(5) 输出聚类结果。
K-medoids算法的流程图如下：
把其他对象分入与其距离最近的类中
任选k个聚类中心
开始
输入聚类个数k和n个对象

否
结束
输出聚类结果
J是否最小
重新调整每类的中心

是

图3.1 k-medoids算法流程图
第4章 对蛋白质冷冻电镜下图像的分类处理
整个处理过程首先模拟产生蛋白质冷冻电镜下的投影图片，然后用k-medoids 算法对其进行分类，最后评估在不同的噪声背景下该算法的分类效果，即评估图片信噪比在某一范围内时，应用k-medoids算法进行分类的效果。
4.1 模拟产生蛋白质冷冻电镜下的图像
事实上，蛋白质冷冻电镜下的图像是从显微镜下拍摄的冷冻蛋白质的投影，要想完整的重构三维立体蛋白质结构，就必须有数量充足且从不同角度投影的蛋白质二维图像。由于实验器材的限制，依据实际投影原理，编写matlab程序，最终模拟随机产生10类不同投影角度的图像，且每类图像又包含20张在平面上具有不同旋转角的图片。接着，为了评估该分类算法对现实中含有噪声的图像的处理能力，利用matlab软件中的awgn()函数，为投影图像加不同程度的高斯白噪声。如图4.1所示：

(a) Clean (b) SNR=1 (c) SNR=7/10 (d) SNR=4/10 (e) SNR=1/10
图4.1 加不同高斯白噪声的投影图像
4.2 对冷冻电镜下蛋白质图片的分类
4.2.1 提取特征
由于同一角度投影的图像有时也会因为在平面上具有不同的旋转角度，或者图像偏离图片中心位置的不同，而显得不同。所以对蛋白质图像分类时，首先要提取图片的旋转不变量和平移不变量作为特征。但是由于模拟的图像都是居中的，所以可以忽略平移量。故只需考虑图像的旋转不变量。输入的图像是一个二维矩阵，可将其看作直角坐标系下的函数，现有两张同样的图片，只是在平面上有一定的旋转角，两张图像如图4.2中a、b所示，分别用函数分别表示为f (x,y)和h(u,v)，然后对其进行极坐标变换，这样图像又可表示为f(r,θ) 和 h(w,φ)，如图4.2中c、d所示，即得到以r，θ为行列的二维矩阵。当r=w时，同一半径上的圆只是经过了某个角度的旋转，即二维矩阵的同一行只是经过了一个平移。接着，只对极坐标变换后的图像的行进行傅里叶变换，即只沿半径方向作变换，由于原图像与旋转后的图像只是在半径方向上经过了平移，依据傅里叶变换的空间平移
(a) (b)

(c) (d)
图4.2 极坐标变换
性质，这两张图像的幅值相等。由以上所述原理，即可得出图片的旋转不变特性。

4.2.2 相似性度量
该过程选用距离作为对象间的相似性度量。根据4.2.1 所提取的特征向量，利用相似性系数计算公式计算所有图像两两间的相似度，得到一个相似系数矩阵。若两对象相似，则所得的值越接近1，反之，则越接近0。之后，用1减去该矩阵则得到一个新矩阵，若矩阵中的值越接近0，则表示相关的两个对象越相似，反之越接近1，则相差越大。这个新矩阵即可看作所有图像的距离矩阵，作为该过程分类的相似性度量。
4.2.3 k-medoids 算法聚类过程
将4.2.2所得的距离矩阵和分类的类别数10作为输入，利用k-medoids算法对其进行聚类，最后输出聚类结果。部分聚类结果如下图4.4所示：
图4.4中共有5类图片，每类3个，即图(a)、(b)、(c)、(d)、(e) 各为一类。

(a) (b) (c) (d) (e)
图4.4 部分分类结果
4.3 对聚类结果的分析评估
运用Fowlkes和Mallows提出的评估聚类结果的方法，编写程序对4.2节分类结果进行统计评估，我们用Fowlkes-Mallows Index(FM)表示评估结果。若FM=1，则表示分类结果完全正确。也就是说，FM越接近1，表示分类效果越好。图4.5和图4.6即为评估统计的结果，横坐标为不同的信噪比，纵坐标为FM。N表示输入图像的总数，K为分类的类别数。图中最后一个点横坐标为11，表示图像在没有噪声的情况下的分类结果。
单独观察图4.5(a)或图4.5(b)可知，信噪比大于1/2时，FM都大于0.7，逼近0.8，表示分类效果都很好。也就是说，用该算法处理实际蛋白质冷冻电镜图像时，图像降噪处理过程只需将图像的信噪比提升至1/2以上即可。另外，从蛋白质三维重构的操作实践中得出，一般在二维图像的分类平均过程中，处理的图像数目大约为200至500张，对比图4.5中的图(a)和图(b)，可得出，输入图像分别为200、500时，分类结果并无太大差异，即该算法处理图像分类的能力满足蛋白质二维分类平均过程的要求。

(a) K=10 (b) K=15
图4.6 对比不同类别数的分类结果

(a) N=100 (b) N=500
图4.5 对不同信噪比图像的分类结果评估

对比图4.6(a)和4.6(b),当分类的类别数分别为10和15时，由图像可得，信噪比大于1/2时，分类效果比较相似，都大于0.7，逼近0.8。虽然K=15的分类效果普遍略低于K=10的分类结果，但是，它们之间的差别不大。由实践经验得知，蛋白质三维重构时，分类数目一般为10到15，因此，该算法可满足要求，并且性能稳定。
第五章 总结与展望
本论文模拟了蛋白质三维重构中，对二维图像的分类过程。运用k-medoids算法对其进行分类，并评估了分类结果，最终证实k-medoids算法可以很好的处理该分类过程，并且适用于蛋白质二维投影图像的分类处理。之前，我已经学习了EMAN2软件的使用，它是对蛋白质投影进行单粒子三维重构的软件。并了解了EMAN2中对蛋白质二维投影图像的分类平均过程。接下来，我们可以在此基础上模拟整个分类平均过程，最终熟悉并掌握该过程的所有细节。最后对该算法进行改进，进一步提高分类结果的准确度，并最大限度节省处理时间。之后将该算法用于EMAN2中二维图像分类平均的处理过程，争取能够发布改进后的EMAN2新版本。
参考文献
[1]林亚静，刘志杰，龚为民.蛋白质结构研究[J].生命科学，2007，19(03).
[2]郭雨珍，冯恩民.蛋白质结构研究现状与展望[J].生物信息学，2005，5(04).
[3]张建辉.K-means聚类算法研究及应用[D].武汉：武汉理工大学，2007.
[4]冯超.K-means聚类算法的研究[D].大连：大连理工大学，2008.
[5]张建萍.基于聚类分析的K_means算法研究及应用[J].计算机应用研究， 2007,24(5):166.
[6]蔡元萃.聚类算法研究综述[J].科技情报开发与经济，2007,17(1):145.
[7]贺玲等.数据挖掘中的聚类算法综述[J].计算机应用研究，2007:10-11.
[8]高小梅，冯云，冯兴杰.增量式K_Medoids聚类算法[J].计算机工程，2005,31.
[9]沈庭芝，方子文.数字图像处理及模式识别[M].北京：北京理工大学出版社， 1998.11-13.
[10]夏良正.数字图像处理[M].南京：东南大学出版社，1999.48-55.
[11]Zhao Z,Singer A. Rotationally Invariant Image Representation for Viewing Direction Classification in Cryo-EM[J].arXiv preprint arXiv:1309.7643,2013.
[12]G. Frahling, C. Sohler.A fast k-means implementation using Coresets[J]. Proceedings of the Twenty-second Annual Symposium on Computational Geometry, 2006, 18(6):1-3.

14