刘小晴,左清曈,刘 青,刘昌灵,杨 刚,龚新奇

1(中国人民大学 信息学院,北京 100872)2(中国人民大学 数学科学研究院,北京 100872)

1 引 言

冷冻电镜成像技术是获取蛋白质等生物分子结构的重要途径之一,其流程分为生物制样、低剂量成像、图像处理三个阶段.图像处理过程中的粒子挑选是重构出生物分子三维结构的基础.由于物理性质所限,冷冻电镜原始图像信噪比极低,单幅图像尺寸较大、粒子颗粒数量多,人工挑选几乎不可行,因此需要效果好、自动化程度高的粒子挑选方法.鉴于深度学习方法在图像处理领域展现出的优势,使之处理冷冻电镜图像成为可能[1].

与传统算法相比,深度学习方法的优势如下:

1)效果好:在图像处理领域,深度学习方法具有较高的准确率,可以更加充分、准确地提取图像特征,完成物体检测等工作,对于大规模数据集也有着更好的适应性.

2)效率高:传统方法中,对于每张输入的图片,都是基于模型进行计算,算法复杂度高,运行时间长;而采用深度学习方法,在训练完模型并经过验证之后,对于新数据,只需调用模型进行预测,处理过程简便快捷.

为了便于冷冻电镜图像提取蛋白质粒子,以及在深度学习方法和应用方面的进一步探索,本文尝试提出冷冻电镜图像处理的新思路,用于图像去噪和颗粒挑选.

本文第1节介绍冷冻电镜相关知识以及相关工作.第2节介绍卷积神经网络模型以及自动编码机模型,结合两者对冷冻电镜图像进行处理,提出了EM-CAE算法(Electron Microscopy-Convolutional AutoEncoder).第3节介绍小波变换相关内容,结合小波变换和卷积自动编码机模型,进一步改进算法,提出DWT-CAE(Discrete Wavelet Transform-Convolutional AutoEncoder)方法,从而更加有效地完成冷冻电镜的去噪工作.第4节为实验设计,由于缺少纯净粒子图像,本文设计了一种蛋白质粒子投影算法,产生了一组protein-projection标准数据集.用小波分析及降噪自动编码机,在protein-projection数据集与mrc原始数据集上进行实验,并对实验结果进行分析.最后对全文进行总结.

2 相关工作

本节主要介绍冷冻电镜图像处理的相关背景,综述了目前已有的处理技术及软件.

2.1 冷冻电镜相关背景

冷冻电镜技术主要包括三项:冷冻电镜二维晶体技术、冷冻电镜单颗粒技术及冷冻电镜断层成像技术.其中冷冻电镜单颗粒技术是结构生物学的重要研究方法之一,对于解析生物大分子以及病毒结构有着重要意义[21].

2.2 冷冻电镜图像处理过程

目前主流的冷冻电镜图像处理软件的流程大致分为以下几步:CTF估计与矫正,颗粒挑选,2D图像分析,三维重构和优化,优化密度图像,确定分辨率并发布.

其中颗粒挑选[22]是冷冻电镜图像处理的基础工作,指的是从原始图像中确定生物分子所在位置和轮廓.目前,粒子自动化识别方法主要分为两大类,一类是模板匹配算法,另一类是特征识别算法.模板匹配算法的理论基础是互函数原理.算法中的模板,可以是已知结构的蛋白质的三维形状投影到二维平面的结果,也可以采用人工方法或其他方式,从冷冻电镜的原始图像中挑选出的含有粒子的局部图像.挑选恰当的模板可以通过比较目标与模板之间的差异而判断粒子.基于模板的匹配算法主要包括Roseman的算法[2]、Sigworth 的算法[3]、Ludtke的算法[4]、Bern的算法[5]、Penczek 的算法[6]等.这一类方法的主要困难之一在于模板获取,很难得到能够充分代表粒子特征的一组模板.同时,模板匹配算法对噪声较为敏感,不同区域的噪声不同,会严重影响算法的性能,如Ludtke的算法[4]中,采用对模板取平均的方式,但这种方式与原始图像相比,信息损失较多,并不能充分展示分子图像的特征.特征识别算法的主要思想是对分子图像的局部或全局特征进行识别,并做出判断.对于特征的选取,可以选择分子图像的统计特征,如Hall的算法[7]和Mallick的算法[8],也可以采用互相关函数特征,如Sorzano的算法[9]等.基于特征识别的算法的不足之处在于,特征对噪声较为敏感,算法复杂度较高且识别率不高.利用挑选好的粒子图像,可以使用较为高效的方式预测和构建3D分子结构[10].图1为含粒子的冷冻电镜原始2D图像.其中(a)为mrc格式,(b)为tiff格式,(c)为局部放大图.

图1 冷冻电镜原始图像Fig.1 Original image of CryoEM

RELION[11](REgularised LIkelihood OptimisatioN)是目前应用较为广泛的一款冷冻电镜图像处理软件,在实验中证实具有较好的效果.

RELION的理论基础是贝叶斯概率模型.在挑选粒子过程中,首先对图像进行傅立叶变换,用最大似然估计的方法计算每块区域中粒子可能出现的概率,并根据贝叶斯概率模型做出判断.RELION处理图像过程中,可以根据图像本身特征生成参数,最终重构出高质量的三维结构.由于使用最大似然方法在整幅图像上进行计算,复杂度较高,因此RELION的性能并不让人满意,处理一批图像需要的时间较长.为了提高软件效率,Sjors Scheres团队开发了可在GPU上运行的RELION2.0,提高了软件处理速度.

EMAN/EMAN2[12]是另一款基于贝叶斯方法编写的冷冻电镜图像处理软件,由Python语言开发.与RELION相比,EMAN/EMAN2的优势在于可以很方便的调用其中各个功能模块,包括图像读取、颗粒挑选、2D分类等;同时,EMAN/EMAN2也具有图形界面,便于使用者进行操作.

Ali Punjani团队开发出了CryoSpark[13],用随机蒙特卡罗方法代替了RELION中的最大似然估计,进一步提高了处理速度.

DeepPicker[14]以CNN方法为理论基础,实现对冷冻电镜图像中粒子的高效自动挑选,并取得了较好的效果.与RELION软件及专家人工判断结果相比,DeepPicker在颗粒挑选方面达到了较高的准确率,但在DeepPicker 软件的工作流程中,去噪环节只使用了简单的高斯去噪方法.本文则尝试在冷冻电镜图像处理的去噪过程中结合多种深度学习方法和手段,以期得到效果更好的粒子.

3 深度学习方法与EM-CAE模型

本节介绍深度学习方法卷积神经网络模型,以及自动编码机模型,并提出将两者结合的EM-CAE模型来处理冷冻电镜图像.通过实验检验了该方法的优缺点,并在下一节加入小波变化的处理方法,改进了模型.

卷积神经网络以及自动编码机

卷积神经网络(Convolutional Neutral Network,CNN)是人工神经网络的一种,是当前应用最广泛的深度学习模型之一,常用于解决图像处理等问题.卷积神经网络的结构具有权值共享的特点,与生物的神经网络十分相似.该结构减少了网络模型的复杂度和参数的数量.当网络的输入是图像时,这一优点表现得更加明显.

图2中给出了典型的CNN模型的示意图[15].

图2 CNN结构示意图Fig.2 Structure of CNN

图3中给出了一般的自动编码机的结构示意图.

当输入数据为图像时,可以将卷积神经网络与自动编码机结合形成卷积自动编码机(convolution autoencoder,简称CAE).卷积自动编码机中,以卷积层代替全连接层进行特征提取.在医学图像的处理上,比如在乳腺癌的图像上进行核检测的方面已经开始应用自动编码机的方法[17].

图3 自动编码机结构图Fig.3 Structure of auto-encoder

本文采用卷积自动编码机尝试去除冷冻电镜图像中的噪声,提出EM-CAE(Electron Microscopy-Convolutional AutoEncoder)模型,具体算法流程如算法1所示.通过在不同数据集上实验发现,该方法兼备同时去噪和粒子挑选的功能,但在原始图像上效果较差,其对于原始图像中较为复杂的噪声特征学习效果不好,因此本文在下一节尝试加入小波变换的处理方法来提高对于原始数据中复杂噪声的训练效果.

算法1.EM-CAE

输入:纯净的模拟数据集x,模拟粒子与噪声叠加的数据y,含有粒子的原始图片z

输出:编码结果x′,y′,z′

1)根据输入数据(x,x′)及(y,y′),按比例划分为训练集、验证集、测试集,三个集合中均含有纯净数据集对应的含噪声数据

2)训练卷积自动编码机模型,当loss函数在训练集和验证集上均稳定收敛,保存模型

3)观察模型对于protein-projection数据集的去噪效果,若效果良好,转到4);否则返回2)并调整模型

4)使用模型对含粒子的mrc原始图像数据进行预测,输出结果并观察.

4 小波变换以及DWT-CAE模型

本节介绍小波变换的相关内容,并基于上一节的EM-CAE模型,提出DWT-CAE模型用于冷冻电镜图像去噪.

小波变换与DWT-CAE模型

小波变换(wavelet transform,WT)[18]是信号和图像处理领域的经典方法之一.

小波变换的主要特点是通过变换过程,使信号或数据的某一方面特征得以充分展示,这样可以对原始数据进行局部化分析,通过伸缩、平移变换,对原始数据逐渐进行多尺度变换和细化,最后实现高频处对时间进行细分,低频处对频率进行细分,自动适应时频信号分析的要求[23],且便于对细节信息进行进一步的处理.

对于信号f(x),其离散小波(DWT,discrete wavelet transform)定义为:

(1)

其逆变换为

(2)

其中C是一个与信号量无关的常数.

降噪是小波变换的功能之一,其基本原理是小波的尺度可变性能够有效对信号进行集中,可以设置合理的阈值,对小波系数进行筛选,并将筛选后的结果进行重构,从而对图像进行降噪.

利用小波变换去噪的关键在于对小波的筛选.由于自动编码机在特征提取上效果显著,本文考虑将小波变换与自动编码机相结合,构造DWT-CAE(Discrete Wavelet Transform-Convolutional AutoEncoder)算法进行图像去噪.图4为算法流程图.图中虚线左侧为训练过程,虚线右侧为测试过程.图中各符号的物理意义与算法2一致.

算法的详细步骤如算法2所示:

算法2.DWT-CAE

输入:纯净的模拟数据集x,模拟粒子与噪声叠加的数据y,含有粒子的原始图片z

输出:重构的图像数据集z′

1) for image in x:

2) 对image进行一层二阶小波变换,选取haar小波基

3) 得到图像的小波系数cAx,(cHx,cVx,cDx)

4) for image in y:

1http://www.rcsb.org/

5) 对image进行一层二阶小波变换,选取haar小波基

6) 得到图像的小波系数cAY,(cHY,cVY,cDY)

7) 令input_x=(cHx,cVx,cDx),input_y=(cHY,cVY,cDY)

8) 以input_x,input_y为输入,训练卷积降噪自动编码机

9) 调节参数,直至模型收敛

10) 输出训练模型

11) for image in z:

12) 对image进行一层二阶小波变换,选取haar小波基

13) 得到图像的小波系数cAz,(cHz,cVz,cDz)

14) 利用10)中的模型对(cHz,cVz,cDz)进行预测

17) 根据16)中的小波系数coeffs对图像进行重构

18) z′= pywt.wavedec2(x,mode=′haar′)

19) 观察实验结果并进行相应分析

图4 DWT-CAE模型图Fig.4 Model of DWT-CAE

5 实验过程与结果

本节进行了实验设计和结果分析.由于目前无法得到纯净的白噪声及粒子图像作为最理想的训练样本,考虑到实验目的是采用降噪自动编码机去除冷冻电镜图像中的噪声,模型主要学习噪声的特征,因此利用已知结构的真实蛋白质模拟二维粒子图像,与噪音数据叠加而生成protein-projection数据集.并介绍了实验结果以及不同方法的对比分析.

5.1 实验平台

本文实验环境为Ubuntu14.04操作系统,Intel(R)Xeon(R)2.40GHz CPU,64GBCPU内存,4GB内存GPU.算法使用Python2.7实现,深度学习部分采用keras框架完成.

5.2 实验数据集

5.2.1 mrc原始数据集

γ-分泌酶(γ-secretase)是由四个亚单位组成的膜内蛋白水解酶,主要参与β-淀粉样蛋白前体(APP)和Notch等重要跨膜蛋白的切割和水解过程.γ-分泌酶是导致阿尔兹海默氏病的重要因素之一.

本文使用的γ-secretase数据集来自[19].数据集中有900张冷冻电镜图像.数据集由两部分组成:冷冻电镜照片mrc文件和与之对应的电镜图像标注文件star.其中mrc文件尺寸为3710(高)*3838(宽),像素值为0～65535,star文件中包含粒子中心点的坐标以及粒子的角度和2D分类信息,实验中主要用到粒子中心点的位置信息.

5.2.2 protein-projection数据集

本文选取来自于蛋白质数据库1的蛋白质结构文件(pdb文件)中的原子坐标信息,将蛋白质分子投影到二维平面,作为训练数据集中的纯净粒子.实验共选取了109个pdb文件,读取蛋白质分子坐标,并进行多方向投影,得到9810张图像,并以此作为训练数据的一部分.生成模拟粒子图像的步骤如下:

算法3.名称:生成模拟粒子形状图像数据集

输入:蛋白质分子结构文件(pdb)

输出:模拟粒子形状图像数据集

1) get pdbfilelist

2) for pdbfile in pdbfilelist:

3 get coordinate={x,y,z}

4) for (xi,Yi,zi) in coordinate:/*计算每个中心点与其他中心点之间的最小距离*/

7)Di=min(disi,j)

8)v=random(20,70) /*随机模拟原子大小*/

9)ri=Di×v

10) a,b,c=random(0,180) /*随机选择投影方向*/

11) Projection=(a,b,c) /*确定投影方向*/

12) Get img_projection

13) Save(img_projrction)

14) get protein-projection dataset

用{I}表示通过算法3的算法生成的粒子形状数据集,用{N}表示从mrc原始数据集提取的噪声集合.对{I}进行随机信号衰减,衰减幅度为e-1或e-2,之后与{N}叠加,得到protein-projection数据.图5给出了数据的效果图,其中图5(a)为模拟粒子,图5(b)为纯噪声,图5(c)为protein-projection数据.

图5 protein-projection数据集效果图Fig.5 Result of protein-projection dataset

5.3 数据预处理

5.3.1 读取mrc文件及star文件

本文实验中用到的冷冻电镜数据包含mrc文件和star文件两部分.

mrc文件中包含照片文件的像素信息,每个像素点的灰度值表示电子信号的强弱.首先,用EMAN2软件包中的EMData模块读入数据,并存储为等价的png文件.将像素值进行L2标准化,将像素矩阵的值压缩至0～255.

star文件是mrc文件经过relion软件完成粒子挑选后输出的文件.文件中前两列是图像中粒子中心点的横纵坐标.实验中,用python语言对文件进行读取,并保存粒子坐标,为后续处理做准备.

5.3.2 直方图均衡化

直方图均衡化是图像处理常用的方法之一.其基本思想是对图像中像素较为密集的灰度级进行展宽,而对图像中像素较为稀疏的灰度进行压缩,这样可以对像原取值的范围进行动态扩展,能够有效提高图像的对比度和灰度色调的变化,从而使图像更清晰.

5.3.3 冷冻电镜图像切割

为了提供训练和测试样本,首先需要从原始图像中获得含有粒子及只含噪声的局部区域,并保存为256*256尺寸的图像.对于mrc原始图像,根据star文件中给出的粒子中心点(x,y),将其周围大小为256*256的区域保存,得到粒子数据集.

5.4 实验结果

本文实现了EM-CAE模型和DWT-CAE模型两种方法,使用EM-CAE方法的效果如图6所示.

图6 在protein-projection数据集(a)和mrc原始数据集(b)上采用EM-CAE的结果Fig.6 Result of EM-CAE

应用DWT-CAE模型对mrc原始图像去噪的部分实验效果如图7所示.

图7 DWT-CAE实验结果图Fig.7 Result of DWT-CAE

图7第一行为mrc原始图像,第二行为去噪后图像,红色轮廓线内为粒子.

我们在研究冷冻电镜图像去噪方法的过程中,先后采用了EM-CAE和小波变换方法,并将小波变换与自动编码机结合,提出了DWT-CAE降噪方法.表1对基于这些方法的实验进行了对比总结:

峰值信噪比(Peak Signal to Noise Ratio,PSNR)和结构相似性(structural similarity index,SSIM)是两种常用的衡量图像相似性方法.由于冷冻电镜图像去噪的目的是为了获得更加清晰而且尽量保留原图像特征、与原图像相似度较高的结果,因此可以采用PSNR和SSIM对不同方法的结果进行检验.

表1 去噪实验方法及比较
Table 1 Experiments of denoising methods and comparison

方法名称EM-CAE小波变换DWT-CAE理论基础卷积去噪自动编码机小波变换小波变换,自动编码机protein-pro-jection数据集效果去噪结果与真实数据的均方误差在10-6～10-7数量级,认为可以完全去除噪声效果较差,无法去除噪声可以去除部分噪声.低频部分去噪结果与真实数据的均方误差在10-6～10-7数量级,认为可以完全去除噪声mrc原始数据集效果无法对含有粒子的图像进行去噪可以去除部分噪声,但粒子边缘仍不清晰去噪效果较好,可以得到较为清晰的粒子边缘

PSNR的定义如下:

(3)

其中n为每个像素点的最大比特数.对于标准图像n=8.PSNR值越大,表示图像间的差异越小.

SSIM的定义如下:

SSIM(X,Y)=l(X,Y)×c(X,Y)×s(X,Y)

(4)

其中l(X,Y),c(X,Y),s(X,Y)分别表示图像的亮度、对比度、结构:

(5)

(6)

(7)

这里μx与μY分别表示图像X与Y的均值,σx与σY表示图像X与Y的方差,而σXY表示图像X与Y的协方差.C1～C3为常数,可以保证分母不为0.

SSIM值可以综合反映变化前后的图像在视觉上的差异,用均值作为亮度的估计,标准差作为对比度的估计,协方差作为结构相似程度的度量.SSIM并不是准确反映出图像像素矩阵的误差,而是综合评价图像变化前后像素结构上的相似度.

表2 实验结果对比
Table 2 Results of different methods and comparison

MethodsPSNRSSIMwavelet transform3.74440.2605DWT-CAE3.79840.3504Gaussian filter3.76600.3946

表2以mrc原始数据集中的粒子局部图像为例给出了采用不同方法进行处理的结果与原图之间的PSNR、SSIM值.

从表2可以看出,用PSNR指标衡量,DWT-CAE方法优于其他几种算法;而用SSIM衡量,高斯滤波方法效果最好,DWT-CAE次之.由于高斯滤波的方法对每个像素邻域内其他像素进行平滑时,不同位置的像素被赋予不同的权值并平均,会有图像信息的丢失[20],而本文的方法通过小波变换后的去噪会根据小波系数重构图像,尽量减少了图像信息的损失.所以可以认为本文的方法效果较为突出.

6 总 结

本文针对冷冻电镜图像特点,采用CNN和自动编码机模型,对冷冻电镜图像进行粒子挑选和去噪.

提出了将小波变换与自动编码机相结合的DWT-CAE模型,用于对冷冻电镜图像去噪,通过在pdb蛋白质投影数据集和mrc原始数据集上的实验,取得良好的效果.DWT-CAE 是我们研究的针对CryoEM图像处理软件DeepConstructing中的重要模块,它为后续的粒子挑选与分类、三维重构提供高质量的数据,我们开发的CNN粒子挑选和分类模块也验证了其有效性.本文设计了蛋白质分子投影算法,利用已解析出三维结构的蛋白质,生成标准数据集,形成对真实数据的完好模拟,可作为同类研究的测试集.