呂永春 趙曉芳 李 華曾祥睿* 徐 旻*

(*中國科學(xué)院計(jì)算技術(shù)研究所 北京 100190) (**中國科學(xué)院大學(xué) 北京 100190) (***中國科學(xué)院智能信息處理重點(diǎn)實(shí)驗(yàn)室 北京 100190) (****美國卡內(nèi)基梅隆大學(xué)計(jì)算機(jī)科學(xué)學(xué)院 匹茲堡 15213)

0 引 言

近年來，電子低溫顯微鏡技術(shù)(electron cryomicroscopy，cryo-EM)在大分子復(fù)合體的3維結(jié)構(gòu)恢復(fù)中發(fā)揮了越來越重要的作用。該技術(shù)促進(jìn)了大分子和復(fù)合體的3維結(jié)構(gòu)進(jìn)一步被研究，揭示其功能可能在生物細(xì)胞機(jī)理、制藥、疾病治療等方面產(chǎn)生重大的突破。電子冷凍斷層掃描(electron cryotomography, cryo-ET)技術(shù)是電子低溫顯微鏡技術(shù)(cryo-EM)的一種應(yīng)用。電子冷凍斷層掃描技術(shù)可以實(shí)現(xiàn)生物大分子和細(xì)胞亞納米分辨率(1～4 nm)的3維成像(3維密度圖像)。Cryo-ET技術(shù)類似人體的CT掃描，由于受機(jī)械臂的限制，對(duì)樣品進(jìn)行有角度采樣時(shí)，通常只能旋轉(zhuǎn)到±70 °，導(dǎo)致部分角度無法采樣，使得重構(gòu)的3維圖像在傅里葉空間存在部分角度范圍信息的缺失。另外，為了保證電子束不損傷樣品，成像過程中會(huì)使用較低劑量的電子，這樣就會(huì)造成噪聲很大，使得信號(hào)和噪聲的比率(signal-to-noise ratio，SNR)很小，產(chǎn)生的原始2維圖像很模糊，也影響3維重構(gòu)的分辨率。為了提高cryo-ET技術(shù)重構(gòu)的分辨率，局部斷層平均(subtomogram averaging, SA)技術(shù)被使用以提高3維重構(gòu)的信噪比，從而得到高分辨的大分子復(fù)合體3維結(jié)構(gòu)。

局部斷層平均技術(shù)需要對(duì)3維顆粒進(jìn)行旋轉(zhuǎn)和平移，最小化3維顆粒與參考顆粒間不相似值，實(shí)現(xiàn)與參考顆粒的對(duì)齊。在局部斷層3維對(duì)齊過程中，每個(gè)局部斷層3維顆粒對(duì)齊都包含6個(gè)參數(shù)(3個(gè)旋轉(zhuǎn)參數(shù)，3個(gè)平移參數(shù))，計(jì)算量很大。為減少局部斷層3維對(duì)齊的計(jì)算量，Kovacs等人[1，2]在傅里葉空間利用球諧函數(shù)實(shí)現(xiàn)快速旋轉(zhuǎn)匹配，把3維數(shù)據(jù)轉(zhuǎn)換到2維空間進(jìn)行旋轉(zhuǎn)匹配，得到旋轉(zhuǎn)參數(shù)，利用傅里葉空間性質(zhì)，得到平移參數(shù)。Xu等人[3]改進(jìn)了局部斷層快速旋轉(zhuǎn)對(duì)齊方法，使用傅里葉空間的3維數(shù)據(jù)進(jìn)行旋轉(zhuǎn)參數(shù)計(jì)算，然后計(jì)算平移參數(shù)，得到6個(gè)參數(shù)更加準(zhǔn)確，但這樣局部斷層對(duì)齊是粗粒度的，為了得到全局擇優(yōu)的旋轉(zhuǎn)和平移參數(shù)，需要對(duì)局部斷層進(jìn)行細(xì)化對(duì)齊。Xu等人[4]提出利用萊文貝格-馬夸特方法(Levenberg-Marquardt algorithm)實(shí)現(xiàn)上述6個(gè)參數(shù)的更新，實(shí)現(xiàn)局部斷層的精對(duì)齊。雖然Xu的方法是較早提出的局部斷層細(xì)化對(duì)齊方法，但該方法每次計(jì)算會(huì)處理整個(gè)3維體，導(dǎo)致計(jì)算量較大。為減少計(jì)算量，本文提出利用小批量梯度下降算法(mini-batch gradient descent, MBGD)實(shí)現(xiàn)局部斷層3維顆粒細(xì)化對(duì)齊，得到局部最優(yōu)解；同時(shí)引入Spark分布式框架進(jìn)行局部斷層細(xì)化對(duì)齊，利用Spark進(jìn)行候選集的分布式并行計(jì)算，得到全局擇優(yōu)旋轉(zhuǎn)和平移參數(shù)，實(shí)現(xiàn)局部斷層的對(duì)齊。通過對(duì)仿真數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)齊，本文提出的基于MBGD和Spark的細(xì)化對(duì)齊方法比現(xiàn)有的基線對(duì)齊方法，具有明顯的優(yōu)勢(shì)，不僅計(jì)算速度快，同時(shí)對(duì)齊精度也有所提高。

本文的主要貢獻(xiàn)如下：(1) 提出一種3維版本的小批量梯度下降算法，并應(yīng)用到局部斷層3維對(duì)齊過程中。(2) 利用小批量梯度下降算法進(jìn)行局部斷層3維顆粒對(duì)齊的優(yōu)化，減少局部斷層3維對(duì)齊的計(jì)算量。(3) 首次實(shí)現(xiàn)Spark分布式局部斷層細(xì)化對(duì)齊，利用Spark進(jìn)行候選參數(shù)的分布式計(jì)算，實(shí)現(xiàn)較快得到最佳對(duì)齊參數(shù)。

本文第1節(jié)對(duì)相關(guān)研究進(jìn)行介紹。第2節(jié)對(duì)局部斷層對(duì)齊的相關(guān)工作進(jìn)行介紹，提出基于小批量梯度下降的局部斷層對(duì)齊算法，并對(duì)算法進(jìn)行詳細(xì)介紹，利用小批量梯度下降細(xì)化對(duì)齊算法實(shí)現(xiàn)局部斷層對(duì)齊，利用Spark框架實(shí)現(xiàn)對(duì)候選集的分布式并行計(jì)算，得到全局擇優(yōu)參數(shù)。第3節(jié)介紹實(shí)驗(yàn)數(shù)據(jù)，實(shí)驗(yàn)環(huán)境，并從對(duì)齊精度和速度兩方面對(duì)2種對(duì)齊算法進(jìn)行比較。最后總結(jié)全文。

1 相關(guān)研究

局部斷層對(duì)齊是局部斷層平均的前提和關(guān)鍵步驟，當(dāng)所有3維顆粒都與參考顆粒進(jìn)行對(duì)齊，才能對(duì)所有3維顆粒進(jìn)行平均，從而提高局部斷層平均結(jié)構(gòu)的信噪比，得到高分辨率3維結(jié)構(gòu)。

1.1 電子冷凍斷層掃描技術(shù)和局部斷層平均技術(shù)

電子冷凍斷層掃描(cryo-ET)技術(shù)是從一個(gè)物體的投影圖像重構(gòu)獲得物體內(nèi)部結(jié)構(gòu)的技術(shù)，通過獲取同一物體的多個(gè)連續(xù)角度下的2維投影圖來反向重構(gòu)它的3維結(jié)構(gòu)。與醫(yī)院中使用的CT掃描類似，簡(jiǎn)單地說，電子斷層掃描技術(shù)就是將一個(gè)樣品(物體)沿著一個(gè)與電子束垂直的軸旋轉(zhuǎn)，每旋轉(zhuǎn)一個(gè)角度，采集這個(gè)物體在相對(duì)應(yīng)方向上的2維投影圖，通過對(duì)這些2維投影圖的處理(如圖像配準(zhǔn)、濾波等)，然后將處理過的不同角度2維投影圖利用反向重構(gòu)技術(shù)進(jìn)行重構(gòu),如加權(quán)背投影(weighted back projection，WBP)等算法，獲得樣品整體3維結(jié)構(gòu)的技術(shù)(見圖 1)。

圖1 電子冷凍斷層掃描技術(shù)成像示意圖[5]

電子冷凍斷層掃描技術(shù)研究對(duì)象包括不具有均一性的蛋白、病毒(如包膜病毒)、細(xì)胞器、細(xì)胞等，因此電子冷凍電鏡斷層掃描技術(shù)也是目前唯一能夠研究原位生物信息的強(qiáng)有力工具。

電子冷凍斷層掃描技術(shù)的優(yōu)點(diǎn)在于旋轉(zhuǎn)角度參數(shù)已知，不需要重新求解，但存在缺失錐(missing wedge)、單張信噪比低、不同投影角度間存在非均勻信噪比等問題，其中缺失錐是電子冷凍斷層掃描技術(shù)面臨的最大問題。為了解決上述問題，需要使用局部斷層平均技術(shù)。

局部斷層平均技術(shù)是將單顆粒3維重構(gòu)與電子斷層成像技術(shù)相結(jié)合的技術(shù)。它的基本思路是：先對(duì)同構(gòu)象的分子顆粒做電子斷層重構(gòu)，再將這些同構(gòu)象的分子顆粒從重構(gòu)體中挑選出來，進(jìn)行類似單顆粒分析技術(shù)的對(duì)齊、分類、平均等步驟，這樣消除了缺失錐，提高了3維顆粒的信噪比，從而實(shí)現(xiàn)高分辨的3維結(jié)構(gòu)(見圖 2)。

1.2 局部斷層對(duì)齊

局部斷層平均技術(shù)工作流程的核心是迭代對(duì)齊和平均過程。在迭代對(duì)齊過程，N個(gè)局部斷層3維顆粒(3個(gè)旋轉(zhuǎn)參數(shù)、3個(gè)平移參數(shù))共有6N-1個(gè)維度，非凸優(yōu)化問題需要解決。局部斷層平均的前提就是局部斷層3維圖像顆粒都已經(jīng)對(duì)齊，故局部斷層3維顆粒對(duì)齊是局部斷層平均中最核心的任務(wù)。由于局部斷層3維顆粒數(shù)據(jù)存在缺失錐，故在對(duì)齊過程中應(yīng)該考慮缺失錐問題，最普遍的校正缺失錐的方法是采用約束性相關(guān)系數(shù)[6](constrained correlation coefficient， CCC)。在實(shí)空間，缺失信息造成3維結(jié)構(gòu)的變形，影響局部斷層數(shù)據(jù)的分類和平均，所以局部斷層對(duì)齊過程中通過約束性互相關(guān)解決缺失錐的影響。

為了加快局部斷層對(duì)齊的速度，Bartesaghi等人[2]將沿射線經(jīng)過傅里葉空間原點(diǎn)的所有傅里葉系數(shù)大小投影到單位球面上相應(yīng)的點(diǎn)上，然后對(duì)相應(yīng)的2幅2維球面圖像利用球諧函數(shù)進(jìn)行2維匹配，得到旋轉(zhuǎn)參數(shù)，并且利用傅里葉空間性質(zhì)，得到平移參數(shù)。Xu等人[3]改進(jìn)了該方法，實(shí)現(xiàn)了3維的快速旋轉(zhuǎn)匹配，并利用局部斷層整體數(shù)據(jù)進(jìn)行對(duì)齊。這些局部斷層對(duì)齊方法都是粗粒度的。Xu等人[4]利用了萊文貝格-馬夸特方法，迭代計(jì)算3個(gè)旋轉(zhuǎn)和3個(gè)平移參數(shù)，利用多線程進(jìn)行并行計(jì)算，實(shí)現(xiàn)局部斷層細(xì)化對(duì)齊。但這些細(xì)化對(duì)齊方法在每次迭代過程中需要計(jì)算整個(gè)局部斷層3維體，非常耗時(shí)。

1.3 Spark分布式計(jì)算框架

Apache Spark[7,8]是一個(gè)通用的開源分布式集群計(jì)算框架。Spark提供了一個(gè)在集群上實(shí)現(xiàn)數(shù)據(jù)并行和編程的接口。 Spark是基于彈性分布式數(shù)據(jù)集(resilient distributed dataset， RDD)。RDD是分布在一組計(jì)算機(jī)上并以容錯(cuò)方式維護(hù)的只讀多集數(shù)據(jù)項(xiàng)。在Spark框架里，客戶端任務(wù)轉(zhuǎn)換成RDD，然后RDD經(jīng)過一系列轉(zhuǎn)換算子(transformation)操作，通過利用行動(dòng)算子(action)觸發(fā)任務(wù)的執(zhí)行。因?yàn)镽DD具有容錯(cuò)機(jī)制的特性，Spark框架在內(nèi)存處理數(shù)據(jù)，適合快速數(shù)據(jù)處理和迭代處理。

在局部斷層對(duì)齊過程中，本文使用Spark分布式并行計(jì)算各初始旋轉(zhuǎn)和平移參數(shù)對(duì)各節(jié)點(diǎn)分別計(jì)算，得到對(duì)應(yīng)的局部最優(yōu)值，然后對(duì)各節(jié)點(diǎn)的計(jì)算結(jié)果進(jìn)行比較，得到全局擇優(yōu)的旋轉(zhuǎn)參數(shù)和平移參數(shù)，從而快速實(shí)現(xiàn)局部斷層細(xì)化對(duì)齊。

2 基于小批量梯度下降和Spark分布式框架的局部斷層細(xì)化對(duì)齊

2.1 基于實(shí)空間的局部斷層對(duì)齊優(yōu)化

局部斷層3維圖像定義為一個(gè)可積分函數(shù)，V(x)：R3→R，局部斷層3維圖像的平移操作ΛT，且T∈R3，定義為[6]：

ΛTV(x):=V(x-T)

(1)

局部斷層3維圖像旋轉(zhuǎn)操作ΛR可以表示為[6]：

ΛRV(x):=V[R-1(x)]

(2)

其中旋轉(zhuǎn)R是一個(gè)3×3旋轉(zhuǎn)矩陣。

局部斷層3維V(x)的旋轉(zhuǎn)和平移操作可以表示為[6]：

ΛT(ΛRV(x))=V(R-1(x)-T)

(3)

局部斷層3維變換參數(shù)β包括成對(duì)的旋轉(zhuǎn)和平移操作，表示為β=(R,T)=(φ,θ,ψ,τ1,τ2,τ3)T，其中旋轉(zhuǎn)參數(shù)R=(φ,θ,ψ)T可視為使用ZYZ慣例的歐拉角[9]，平移參數(shù)T=(τ1,τ2,τ3)T。

局部斷層3維圖像V1的缺失錐由傅里葉空間的3維模板M表示，其中在頻率測(cè)量可用的區(qū)域中值為1，在有限傾斜范圍導(dǎo)致沒有數(shù)據(jù)覆蓋的區(qū)域中值為0。2個(gè)局部斷層3維圖像V1和V2在傅里葉空間對(duì)齊后的重疊區(qū)域Ω:=M(V1)ΛRM(V2)。當(dāng)2個(gè)局部斷層3維圖像對(duì)齊時(shí)，約束性互相關(guān)只考慮傅里葉空間中旋轉(zhuǎn)的最佳重疊區(qū)域，并根據(jù)傅里葉空間的平移不變性，消除平移的影響，即只考慮在傅里葉空間旋轉(zhuǎn)操作。為了降低局部斷層周圍噪聲的影響，本文在實(shí)空間定義一個(gè)二值化模板函數(shù)M。

在實(shí)空間，局部斷層3維圖像V1的正則化函數(shù)表示為[6]：

(4)

(5)

局部斷層3維圖像V2的約束性函數(shù)表示為[6]：

(6)

其中，

(7)

(8)

在傅里葉空間的平移不變性，只需考慮旋轉(zhuǎn)參數(shù)，通過快速旋轉(zhuǎn)匹配算法[2，3]，得到一組初始旋轉(zhuǎn)候選集合{R1，R2，…，RN}，然后通過快速平移匹配算法[10]，獲取一組對(duì)應(yīng)的平移候選集合{T1，T2，…，TN}。這樣形成N對(duì)變換候選集合，{(R1,T1), (R2,T2)，…，(RN,TN)}，其中N是候選集合數(shù)，但這些初始成對(duì)候選集合不能實(shí)現(xiàn)局部斷層的精對(duì)齊，需要進(jìn)行局部斷層細(xì)化對(duì)齊。

給定一組參數(shù){R,T}，本文在實(shí)空間通過小批量梯度下降算法[11]對(duì)局部斷層3維圖像進(jìn)行局部細(xì)化對(duì)齊。利用該算法可以實(shí)現(xiàn)局部斷層3維圖像與參考圖像進(jìn)行局部精對(duì)齊。在局部斷層體積V中，得到一組新的旋轉(zhuǎn)參數(shù)值Rk和平移參數(shù)值Tk，使得標(biāo)準(zhǔn)化歐式距離變得越來越?。?/p>

dRk, Tk≥dRk+1, Tk+1

(9)

在局部斷層3維圖像對(duì)齊時(shí)，通常一個(gè)局部斷層3維圖像固定作為參考體積(無缺失錐)，另一個(gè)局部斷層3維圖像(有缺失錐)進(jìn)行旋轉(zhuǎn)和平移操作，與固定的參考體積進(jìn)行對(duì)齊。然而不能直接對(duì)局部斷層體積使用MBGD算法。針對(duì)局部斷層體積的特點(diǎn)，本文設(shè)計(jì)沿局部斷層3維圖像的x軸進(jìn)行小批量采樣，這樣就可以在局部斷層3維圖像上使用小批量梯度下降算法進(jìn)行變換參數(shù)的計(jì)算，實(shí)現(xiàn)局部斷層體積的細(xì)化對(duì)齊。

通過MBGD算法不斷更新β=(R,T)，定義3維局部斷層體積對(duì)齊的損失函數(shù)L為：

(10)

其中n是3維局部斷層體積沿x軸的長(zhǎng)度，Hβ(xi)∶=(V1(xi)*-ΛβV2(xi)*)2。

MBGD算法的迭代表達(dá)如下：

(11)

其中在每次迭代中，局部斷層3維體中沿x軸的截面坐標(biāo)i在{1,…,n-B}被隨機(jī)地選擇，小批量的長(zhǎng)度為B，αk是步長(zhǎng)。

通過式(10)和式(11)，描述算法的流程(算法1)。Xu的對(duì)齊方法每次迭代需要全部數(shù)據(jù)參與計(jì)算，計(jì)算量大，但能實(shí)現(xiàn)較快收斂。而隨機(jī)梯度下降算法[12](stochastic gradient descent, SGD)每次迭代只需一個(gè)樣本，雖然計(jì)算速度快，但不是每次迭代都向著整體優(yōu)化方向，收斂速度緩慢。基于MBGD細(xì)化對(duì)齊算法每次僅需要小批量的數(shù)據(jù)參與梯度計(jì)算，所以計(jì)算速度快，且能較快實(shí)現(xiàn)收斂。所以本文提出的基于MBGD細(xì)化對(duì)齊算法結(jié)合了Xu的對(duì)齊算法和SGD算法的優(yōu)點(diǎn)，既能實(shí)現(xiàn)較快的計(jì)算時(shí)間，又能保證收斂的速度，適合在局部斷層體積精對(duì)齊中應(yīng)用。

算法1基于MBGD局部斷層細(xì)化對(duì)齊方法最小化約束性不相似值12nni=1Hβ(xi) d=0且i={1, 2,…, n} eps=0.0001 oldd=mind=1000 fork=0; k≤maxIter; k++do 隨機(jī)選擇截面i, i={1, 2, …, n-B} d=B+iiHβ(xi)=B+ii(V1(xi)?-ΛβV2(xi)?)2 β=β-αBi+BiHβ(xi)′ ifd

2.2 基于Spark架構(gòu)分布式并行局部斷層細(xì)化對(duì)齊過程

為了實(shí)現(xiàn)全局擇優(yōu)的旋轉(zhuǎn)和平移參數(shù)，本文首次提出基于Spark架構(gòu)分布式并行局部斷層細(xì)化對(duì)齊過程，首先需要對(duì)不同旋轉(zhuǎn)候選參數(shù)和平移候選參數(shù)中執(zhí)行多次局部斷層細(xì)化對(duì)齊，實(shí)現(xiàn)每個(gè)候選參數(shù)(旋轉(zhuǎn)和平移)得到局部最優(yōu)，然后比較不同旋轉(zhuǎn)參數(shù)和平移參數(shù)的局部最優(yōu)值，得到全局擇優(yōu)的旋轉(zhuǎn)和平移參數(shù)。為了實(shí)現(xiàn)不同候選集同步計(jì)算，通過1.3節(jié)的分析，本文選取Spark架構(gòu)實(shí)現(xiàn)并行局部斷層細(xì)化對(duì)齊，得到全局擇優(yōu)的旋轉(zhuǎn)和平移參數(shù)。雖然Spark在內(nèi)存中處理數(shù)據(jù)，具有容錯(cuò)機(jī)制和冗余機(jī)制，適合高頻數(shù)據(jù)交換和大數(shù)據(jù)迭代處理，但Spark仍然是一個(gè)粗粒度的分布式框架，為了實(shí)現(xiàn)基于Spark架構(gòu)分布式局部斷層細(xì)化對(duì)齊，需要設(shè)計(jì)新穎的數(shù)據(jù)交換和轉(zhuǎn)移代碼。

本文設(shè)計(jì)基于MBGD和Spark框架分布式局部斷層細(xì)化對(duì)齊，具體流程如圖3，通過快速旋轉(zhuǎn)匹配和快速平移匹配，得到N對(duì)旋轉(zhuǎn)和平移集合數(shù)組，把數(shù)組轉(zhuǎn)換成RDD，并對(duì)RDD進(jìn)行分區(qū)，然后通過調(diào)用mapPartitions操作，執(zhí)行各節(jié)點(diǎn)任務(wù)的分配，調(diào)用collect操作，每個(gè)節(jié)點(diǎn)并行執(zhí)行基于MBGD局部斷層細(xì)化對(duì)齊，返回各節(jié)點(diǎn)局部最優(yōu)變換參數(shù)，通過比較，得到全局擇優(yōu)參數(shù)。

圖3 Spark框架的分布式局部斷層細(xì)化對(duì)齊算法流程

3 實(shí)驗(yàn)與分析

3.1 仿真數(shù)據(jù)

從Protein Data Bank (PDB)數(shù)據(jù)庫下載分子伴侶蛋白質(zhì)(GroEL)原子模型，設(shè)定了分辨率和空間體素值，并對(duì)GroEL數(shù)據(jù)進(jìn)行低通濾波。

本文使用Situs PDB2VOL[13]程序獲得GroEL電子密度圖，對(duì)GroEL電子密度圖進(jìn)行隨機(jī)旋轉(zhuǎn)和平移操作，利用散焦值模擬對(duì)比度傳遞函數(shù)(contrast transfer function， CTF)。在指定的傾斜范圍和角度增量條件下，對(duì)GroEL電子密度圖進(jìn)行投影，以模擬斷層數(shù)據(jù)進(jìn)行有角度采樣過程。對(duì)投影圖像添加高斯噪聲和調(diào)制傳遞函數(shù)噪聲(modulation transfer function noise， MTF)來模擬電子光學(xué)效應(yīng)。投影后的2維圖像利用加權(quán)反投影算法(weighted back projection， WBP)進(jìn)行3維重構(gòu)，這樣就產(chǎn)生模擬的局部斷層3維數(shù)據(jù)集。

GroEL蛋白質(zhì)原子模型(PDB ID:1KP8)用于產(chǎn)生網(wǎng)格點(diǎn)數(shù)為64×64×64的局部斷層3維體積，其中體素尺寸為0.6 nm×0.6 nm×0.6 nm， 散焦為-6 μm。

在傾斜范圍為±60 °、角度增量為1 °的3種不同SNR(0.01，0.03，0.003)條件下，分別模擬20個(gè)局部斷層數(shù)據(jù)，這些模擬數(shù)據(jù)進(jìn)行隨機(jī)旋轉(zhuǎn)和平移操作。在傾斜范圍為±40 °、角度增量為1 °的3種不同SNR(0.01，0.03，0.003)條件下，也進(jìn)行相同的操作。GroEL局部斷層仿真數(shù)據(jù)經(jīng)常出現(xiàn)在局部斷層分析方面的相關(guān)文獻(xiàn)中[2-4]。

沿x-z平面得到的不同傾斜范圍和SNR條件下的中心切片，如圖4所示。在圖4中，具有較小傾斜范圍(如±40 °)和較低SNR(0.003)的局部斷層數(shù)據(jù)明顯具有較大的變形，肉眼幾乎無法分辨其中所包含的結(jié)構(gòu)。

3.2 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)的局部斷層數(shù)據(jù)為分子伴侶蛋白質(zhì)復(fù)合體(GroEL和GroEL/ES)數(shù)據(jù)集[6]。

為了收集這些GroEL14GroES7復(fù)合體，F(xiàn)?rster等人[6]采用如下流程：1 μM GroEL14和5 μM GroES7需要放在5 mM MgCl2，5 mM KCL， 5 mM ADP，1 mM DTT和12.5 mM Hepes (pH 7.5)的緩沖液中進(jìn)行孵育，并在30 ℃下培育15 min。使用網(wǎng)格將3.5 μl蛋白質(zhì)溶液和0.5 μl的10 nm BSA膠體金懸浮液進(jìn)行混合。樣品用插入式冷凍法進(jìn)行玻璃化。

圖4 仿真局部斷層數(shù)據(jù)切片(x-z平面)

樣品在傾斜角度±65 °，2 °或2.5 °角度增量條件下，使用Tecnai G2 Polara顯微鏡(配備2k×2k FEI CCD相機(jī))進(jìn)行低溫采樣和投影，并利用UCSF tomography軟件在單軸傾斜下進(jìn)行2維投影數(shù)據(jù)的獲取。圖像在2k×2k 像素的CCD相機(jī)和7～4 μm的散焦水平上進(jìn)行記錄。物體像素尺寸為 0.6 nm。

3.3 實(shí)驗(yàn)環(huán)境

本文實(shí)現(xiàn)的基于MBGD和Spark分布式方法的局部斷層細(xì)化對(duì)齊算法在4臺(tái)服務(wù)器上運(yùn)行，每臺(tái)服務(wù)器配置包括2個(gè)1.7 GHz Intel Xeon Bronze 3104 CPUs，含有12個(gè)物理核，1個(gè)千兆網(wǎng)卡，30 G內(nèi)存。

Spark集群部署采用stand alone模式，1個(gè)為管理節(jié)點(diǎn)，4個(gè)為工作節(jié)點(diǎn)。代碼開發(fā)是在Spark 2.0和Python 2.7環(huán)境下進(jìn)行。

3.4 實(shí)驗(yàn)數(shù)據(jù)分類

數(shù)千個(gè)實(shí)驗(yàn)的分子伴侶蛋白質(zhì)復(fù)合體(GroEL和GroEL/ES)數(shù)據(jù)集也包含假定粒子，因此需要手動(dòng)挑選，并且與局部斷層平均顆粒進(jìn)行對(duì)齊，排除低的互相關(guān)系數(shù)(如CCC≤0.42)，剩下的顆粒被挑選進(jìn)行局部斷層對(duì)齊和分類。實(shí)驗(yàn)的～800 kDa GroEL14和GroEL14/GroES7局部斷層復(fù)合體數(shù)據(jù)集作為局部斷層對(duì)齊和分類研究的準(zhǔn)標(biāo)準(zhǔn)[6,14,15]。在數(shù)據(jù)集中786個(gè)局部斷層顆粒以任意方向和非監(jiān)督方式對(duì)所有局部斷層顆粒進(jìn)行平均而對(duì)齊。

本文使用MCO-A[13]分類算法對(duì)GroEL和GroEL/ES復(fù)合體分類，其中含有10個(gè)初始類和7倍對(duì)稱。本文用MCO-A方法得到最終3個(gè)不同的類，其結(jié)果與先前在文獻(xiàn)[6，14-16]中發(fā)布的結(jié)果一致。由MCO-A分類方法產(chǎn)生的每個(gè)分類平均的中心切片顯示在圖5中。

圖5 使用MCO-A方法分類GroEL和GroEL/ES復(fù)合體

3.5 不同局部斷層對(duì)齊算法精度比較

為了比較基于MBGD局部斷層細(xì)化對(duì)齊方法與Xu的對(duì)齊方法，本文使用無噪聲的參考體積，在不同傾斜范圍和不同信噪比條件下，20個(gè)仿真局部斷層數(shù)據(jù)與參考體積進(jìn)行對(duì)齊，然后通過假設(shè)檢驗(yàn)計(jì)算兩種方法對(duì)應(yīng)的互相關(guān)系數(shù)值，計(jì)算對(duì)應(yīng)的t值和P值，比較2種算法的對(duì)齊精度。

無噪聲參考體積通過GroEL結(jié)構(gòu)(PDB ID:1KP8)產(chǎn)生的。參考體積低通濾波至6 nm分辨率，并用作對(duì)齊過程的初始參考。

在Xu的對(duì)齊方法使用約束性互相關(guān)方法評(píng)價(jià)其對(duì)齊精度，為了科學(xué)地評(píng)價(jià)2種對(duì)齊算法精度，在基于MBGD細(xì)化對(duì)齊方法中也使用約束性互相關(guān)方法評(píng)價(jià)其對(duì)齊算法的精度。

使用基于成對(duì)數(shù)據(jù)的t檢驗(yàn)來比較2種不同對(duì)齊方法的精度。對(duì)2種對(duì)齊方法得到的一批成對(duì)的互相關(guān)值，利用基于MBGD對(duì)齊方法得到一系列互相關(guān)值減去Xu的對(duì)齊方法得到的一系列互相關(guān)值，這樣就形成了一系列差值，然后通過計(jì)算差值的P值，來比較兩者對(duì)齊方法是否存在顯著性。

如表1和圖6所示，在傾斜角范圍為±60 °下，本文提出的基于MBGD細(xì)化對(duì)齊方法與Xu的對(duì)齊方法進(jìn)行比較，各自使用模擬的信噪比為0.003的局部斷層數(shù)據(jù)進(jìn)行對(duì)齊，得到系列CCC差值的平均值為正，對(duì)應(yīng)的P值為1.01E-12(P<0.01)，說明2種方法具有非常明顯的差異，基于MBGD細(xì)化對(duì)齊方法的性能優(yōu)于Xu的對(duì)齊方法。但在傾斜角范圍為±60 °、信噪比為0.01和0.03條件下，基于MBGD細(xì)化對(duì)齊方法和Xu的對(duì)齊方法無顯著的差異(P>0.05)。

表1 在傾斜范圍±60 °下2種對(duì)齊算法間P值比較

圖6 不同傾斜范圍和信噪比下2種對(duì)齊方法互相關(guān)差值的平均值

如表2和圖6所示，在傾斜角范圍為±40 °、信噪比為0.003條件下，仿真的局部斷層數(shù)據(jù)缺失信息更多，且信噪比也更低，基于MBGD細(xì)化對(duì)齊方法與Xu的對(duì)齊方法進(jìn)行比較，得到系列CCC差值的均值為正，對(duì)應(yīng)的P值為2.15E-05(P<0.01)，表明2種對(duì)齊方法有非常顯著差異，基于MBGD細(xì)化對(duì)齊方法的對(duì)齊性能也優(yōu)于Xu的對(duì)齊方法。而在傾斜角范圍為±40 °、信噪比為0.01和0.03條件下，2種對(duì)齊方法沒有顯著差異(P>0.05)。

表2 在傾斜范圍±40 °下2種對(duì)齊算法間P值比較

實(shí)際上，利用電子冷凍斷層掃描技術(shù)重構(gòu)的斷層圖像，因?yàn)殡娮觿┝康?，存在缺失錐，造成3維圖像信噪比低，信噪比為0.003時(shí)更接近于實(shí)驗(yàn)條件下的局部斷層3維圖像。通過對(duì)2種對(duì)齊算法精度的比較，在傾斜范圍為±60 °和±40 °下，信噪比為0.003，通過對(duì)表1、表2和圖6的分析，基于MBGD細(xì)化對(duì)齊方法在對(duì)齊精度上優(yōu)于Xu的對(duì)齊方法，這也表明基于MBGD細(xì)化對(duì)齊方法更加適合對(duì)實(shí)驗(yàn)環(huán)境下局部斷層數(shù)據(jù)圖像對(duì)齊。

3.6 不同局部斷層對(duì)齊算法運(yùn)算時(shí)間比較

為了客觀公正地比較2種對(duì)齊算法的運(yùn)行時(shí)間，2種對(duì)齊算法都用Python語言進(jìn)行實(shí)現(xiàn)。

本文使用隨機(jī)方向上的局部斷層平均體積作為兩者對(duì)齊算法的初始參考，每次局部斷層對(duì)齊算法收斂時(shí)，就會(huì)得到一個(gè)新的參考和對(duì)應(yīng)的分辨率，通過在指定的最大迭代次數(shù)下，得到最佳的分辨率值對(duì)應(yīng)的運(yùn)算時(shí)間和迭代次數(shù)。這種初始參考稱為無參考策略[16]是不需要外部參考，因?yàn)橥獠繀⒖紩?huì)導(dǎo)致參考體積的偏差。

首先比較每個(gè)對(duì)齊算法使用一次的運(yùn)行時(shí)間。利用在傾斜角范圍為±60 °、信噪比為0.003條件下仿真的20個(gè)局部斷層數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。如圖7所示，本文提出的基于MBGD細(xì)化對(duì)齊算法運(yùn)行時(shí)間是63 s，而Xu的方法耗時(shí)是150 s。

圖7 2種對(duì)齊算法運(yùn)算時(shí)間

通過2種對(duì)齊算法運(yùn)行時(shí)間的比較，基于MBGD局部斷層細(xì)化對(duì)齊算法的運(yùn)行時(shí)間大體上是Xu對(duì)齊方法的一半，甚至更短。

然后利用GroEL局部斷層實(shí)驗(yàn)數(shù)據(jù)比較2種對(duì)齊方法在得到最佳分辨率時(shí)的迭代次數(shù)。把實(shí)驗(yàn)數(shù)據(jù)分成相同的兩部分，每部分獨(dú)立進(jìn)行對(duì)齊和迭代，然后對(duì)每次迭代后對(duì)齊的數(shù)據(jù)進(jìn)行平均，使用金標(biāo)準(zhǔn)的FSC(FSC=0.143)方法計(jì)算分辨率，平均全部數(shù)據(jù)為最新的參考，不斷循環(huán)，直到循環(huán)結(jié)束。利用金標(biāo)準(zhǔn)FSC為0.143的規(guī)則，基于MBGD細(xì)化對(duì)齊算法和Xu的對(duì)齊方法分別對(duì)GroEL局部斷層實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)齊和平均，記錄得到最佳分辨率對(duì)應(yīng)的迭代次數(shù)。

首先使用基于MBGD細(xì)化對(duì)齊算法對(duì)GroEL局部斷層實(shí)驗(yàn)數(shù)據(jù)進(jìn)行迭代對(duì)齊和平均，經(jīng)過5次迭代，得到最佳分辨率為30 ?的GroEL局部斷層平均(見圖8)。

圖8 2種對(duì)齊方法對(duì)GroEL局部斷層實(shí)驗(yàn)數(shù)據(jù)平均

然后，使用Xu的對(duì)齊方法對(duì)GroEL局部斷層實(shí)驗(yàn)數(shù)據(jù)進(jìn)行迭代對(duì)齊和平均，經(jīng)過9次迭代，最終GroEL局部斷層數(shù)據(jù)平均的最佳分辨率為32.5 ?(見圖8)。

通過對(duì)得到最佳分辨率時(shí)迭代次數(shù)的比較，基于MBGD局部斷層細(xì)化對(duì)齊算法得到最佳分辨率的迭代次數(shù)也近似是Xu方法得到最佳分辨率迭代次數(shù)的一半。

基于上述2種運(yùn)算時(shí)間的分析，在細(xì)化對(duì)齊得到的最佳分辨率基本不變的條件下，基于MBGD局部斷層細(xì)化對(duì)齊算法的運(yùn)算時(shí)間明顯優(yōu)于Xu的對(duì)齊方法，近似為Xu方法的一半，并且基于MBGD局部斷層細(xì)化對(duì)齊算法的迭代次數(shù)基本上為Xu的對(duì)齊方法迭代次數(shù)的1/2。

4 結(jié) 論

針對(duì)局部斷層對(duì)齊過程計(jì)算量大，且需要細(xì)化對(duì)齊的問題，本文提出基于小批量梯度下降的細(xì)化對(duì)齊算法，在實(shí)空間實(shí)現(xiàn)優(yōu)化約束性不相似值。通過對(duì)仿真局部斷層數(shù)據(jù)進(jìn)行測(cè)試，在傾斜范圍為±60 °和±40 °、信噪比為0.003條件下，通過成對(duì)數(shù)據(jù)的t檢驗(yàn)，P值都小于0.01，證明本文提出的基于MBGD細(xì)化對(duì)齊算法在對(duì)齊精度上明顯勝過Xu的對(duì)齊方法。相對(duì)于Xu的對(duì)齊方法，基于MBGD細(xì)化對(duì)齊方法更加適合對(duì)實(shí)驗(yàn)環(huán)境下(如信噪比為0.003)局部斷層數(shù)據(jù)進(jìn)行對(duì)齊。

通過對(duì)仿真局部斷層數(shù)據(jù)和實(shí)驗(yàn)GroEL數(shù)據(jù)進(jìn)行對(duì)齊，在對(duì)齊算法運(yùn)行一次的時(shí)間比較上，本文提出的基于MBGD對(duì)齊方法的運(yùn)算時(shí)間大體為Xu對(duì)齊方法運(yùn)算時(shí)間的一半；在得到最佳分辨率的迭代次數(shù)比較上，本文提出的基于MBGD局部斷層細(xì)化對(duì)齊算法迭代次數(shù)基本上為Xu的對(duì)齊方法迭代次數(shù)的1/2，同時(shí)基于MBGD對(duì)齊方法得到的最佳分辨率(30 ?)也略優(yōu)于Xu對(duì)齊方法得到的最佳分辨率(32.5 ?)。本文提出的基于MBGD局部斷層細(xì)化對(duì)齊方法有利于局部斷層對(duì)齊的結(jié)果優(yōu)化和更改。

此外，本文在局部斷層對(duì)齊過程中使用Spark框架實(shí)現(xiàn)局部斷層細(xì)化對(duì)齊，該框架可以實(shí)現(xiàn)局部斷層細(xì)化對(duì)齊的分布式計(jì)算，從而得到全局擇優(yōu)參數(shù)，相對(duì)于其他分布式框架，Spark框架實(shí)現(xiàn)更簡(jiǎn)單。