超低信噪比冷凍電鏡圖像的深度學習去噪算法—DWT-CAE

劉小晴,左清曈,劉 青,劉昌靈,楊 剛,龔新奇

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粒子挑選和分類模塊也驗證了其有效性.本文設計了蛋白質分子投影算法,利用已解析出三維結構的蛋白質,生成標準數據集,形成對真實數據的完好模擬,可作為同類研究的測試集.