呂襲明1)2) 李輝1)? 尤菁1)2) 李偉1) 王鵬業(yè)1)2) 李明1)2) 奚緒光竇碩星

1)(中國(guó)科學(xué)院物理研究所,軟物質(zhì)物理重點(diǎn)實(shí)驗(yàn)室,北京凝聚態(tài)物理國(guó)家實(shí)驗(yàn)室,北京 100190)

2)(中國(guó)科學(xué)院大學(xué)物理科學(xué)學(xué)院,北京 100049)

3)(西北農(nóng)林科技大學(xué)生命科學(xué)院,楊凌712100)

單分子熒光共振能量轉(zhuǎn)移數(shù)據(jù)處理的優(yōu)化算法?

呂襲明1)2) 李輝1)? 尤菁1)2) 李偉1) 王鵬業(yè)1)2) 李明1)2) 奚緒光3)竇碩星1)2)?

1)(中國(guó)科學(xué)院物理研究所,軟物質(zhì)物理重點(diǎn)實(shí)驗(yàn)室,北京凝聚態(tài)物理國(guó)家實(shí)驗(yàn)室,北京 100190)

2)(中國(guó)科學(xué)院大學(xué)物理科學(xué)學(xué)院,北京 100049)

3)(西北農(nóng)林科技大學(xué)生命科學(xué)院,楊凌712100)

(2016年12月8日收到;2017年3月14日收到修改稿)

單分子熒光共振能量轉(zhuǎn)移(smFRET)技術(shù)是當(dāng)今單分子生物物理研究領(lǐng)域的重要實(shí)驗(yàn)手段,該技術(shù)通過測(cè)量供體、受體熒光光強(qiáng)以及二者間的共振能量轉(zhuǎn)移效率,揭示標(biāo)記位點(diǎn)間的距離,用于研究DNA、蛋白質(zhì)等生物大分子的構(gòu)象變化.然而,當(dāng)前傳統(tǒng)數(shù)據(jù)處理方法大量依賴人工干預(yù),噪音大,嚴(yán)重影響了實(shí)驗(yàn)效率和數(shù)據(jù)的可靠性.本文提出了一種針對(duì)smFRET數(shù)據(jù)的自動(dòng)分析算法.該算法主要包括三個(gè)部分:基于計(jì)算供體與受體熒光光強(qiáng)的相關(guān)系數(shù)來確定受體與供體對(duì)應(yīng)熒光點(diǎn)的自動(dòng)匹配算法、甄別錯(cuò)誤點(diǎn)的篩選算法以及基于隱馬爾可夫模型的全局?jǐn)M合算法.經(jīng)改進(jìn)后的算法大大簡(jiǎn)化了傳統(tǒng)算法中需要人工干預(yù)的步驟,而且自動(dòng)篩除了實(shí)驗(yàn)數(shù)據(jù)中主要的幾類噪音.將改進(jìn)的算法應(yīng)用于人類端粒重復(fù)序列G-四聯(lián)體(G4)DNA折疊動(dòng)力學(xué)的數(shù)據(jù)分析,結(jié)果顯示優(yōu)化算法比傳統(tǒng)算法能夠更快地得到更高信噪比的數(shù)據(jù),而且該數(shù)據(jù)結(jié)果清晰地表明G4的折疊體現(xiàn)出多態(tài)性并受到鉀離子濃度的影響.

單分子熒光共振能量轉(zhuǎn)移,數(shù)據(jù)處理算法,G-四聯(lián)體DNA,折疊動(dòng)力學(xué)

1 引言

隨著生物物理學(xué)和單分子生物技術(shù)的發(fā)展,更多的生命科學(xué)問題可以在單分子尺度上加以研究.單分子熒光共振能量轉(zhuǎn)移(single-molecule fl uorescence resonance energy transfer,smFRET)技術(shù)作為當(dāng)前單分子生物物理學(xué)領(lǐng)域的熱門技術(shù),可以實(shí)時(shí)地觀測(cè)研究許多生物大分子的構(gòu)像變化、生物大分子之間的相互作用等過程,從而直觀地說明各種生物物理學(xué)中的結(jié)構(gòu)與動(dòng)力學(xué)問題[1?5].然而,由于實(shí)驗(yàn)中使用的標(biāo)記熒光分子的不穩(wěn)定性和單分子生物實(shí)驗(yàn)流程的復(fù)雜性,smFRET實(shí)驗(yàn)發(fā)展之初往往難以獲得高信噪比的實(shí)驗(yàn)數(shù)據(jù),實(shí)驗(yàn)的重復(fù)性也較差.為了解決上述問題,前人進(jìn)行了種種嘗試,解決了熒光分子容易淬滅和閃爍的問題[6?8]、規(guī)范了實(shí)驗(yàn)操作的步驟[9]、提出了數(shù)據(jù)處理的理論基礎(chǔ)并開發(fā)了相應(yīng)的數(shù)據(jù)處理程序[10].然而,這些已有的smFRET數(shù)據(jù)處理算法依賴人工建立坐標(biāo)匹配關(guān)系,容易引入錯(cuò)誤配對(duì)及不活躍的點(diǎn),精度差且人為主觀因素大,而且擬合EFRET-t曲線時(shí)也必須對(duì)每條數(shù)據(jù)曲線分別進(jìn)行操作,自動(dòng)化程度低.

針對(duì)以上問題,我們對(duì)傳統(tǒng)處理算法提出了優(yōu)化.首先,利用smFRET數(shù)據(jù)中供體(donor)與受體(acceptor)熒光變化曲線的負(fù)相關(guān)特性實(shí)現(xiàn)了供體與受體熒光點(diǎn)坐標(biāo)對(duì)應(yīng)關(guān)系的自動(dòng)確立.其次,利用該負(fù)相關(guān)性對(duì)數(shù)據(jù)進(jìn)行自動(dòng)篩選,可有效篩除雜質(zhì)發(fā)光或熒光提前淬滅等情況.最后,通過改進(jìn)隱馬爾可夫模型(hidden Markov model,HMM)在數(shù)據(jù)擬合過程中的運(yùn)用方式,實(shí)現(xiàn)一次性全局?jǐn)M合所有EFRET-t曲線,省去人工分別擬合再整合的麻煩.經(jīng)過改進(jìn)的數(shù)據(jù)處理算法大大提高了實(shí)驗(yàn)數(shù)據(jù)的可靠性,自動(dòng)化的處理方式幫助處理大批量的實(shí)驗(yàn)結(jié)果,有效提高了實(shí)驗(yàn)效率以及實(shí)驗(yàn)的可重復(fù)性.

2 算法原理

2.1 smFRET原理及熒光分子光強(qiáng)的負(fù)相關(guān)性

供體熒光分子在受到外源激光激發(fā)后發(fā)出熒光,該熒光的波譜與受體熒光分子的激發(fā)波譜交疊,若此時(shí)受體與供體的距離足夠接近,則供體熒光光能會(huì)部分傳輸至受體,使其受激發(fā)出熒光.假設(shè)供體熒光光強(qiáng)為ID,受體光強(qiáng)為IA,則能量傳遞的效率.考慮到供體熒光泄漏、受體熒光分子直接發(fā)光等因素的影響,該公式最終被修正為,γ為修正因子[9,11,12].根據(jù)已有研究,該效率與兩熒光分子的距離R具有如下關(guān)系:,其中常數(shù)R0為F?rster半徑[9,13].由此可知,通過觀測(cè)標(biāo)記在同一生物大分子上的一對(duì)供體受體熒光分子實(shí)時(shí)的熒光光強(qiáng)變化,計(jì)算EFRET便可獲知兩個(gè)熒光分子標(biāo)記位點(diǎn)間距離R的變化信息,從而反映出生物大分子的結(jié)構(gòu)變化.

需要注意的是,在實(shí)驗(yàn)中外源激發(fā)光強(qiáng)度相對(duì)恒定,根據(jù)能量守恒,ID與IA應(yīng)滿足負(fù)相關(guān)關(guān)系,ID-t曲線與IA-t曲線的相關(guān)系數(shù)

(xi,yi分別為ID與IA的序列)應(yīng)當(dāng)為負(fù),且該負(fù)相關(guān)關(guān)系符合得越好,r應(yīng)當(dāng)越接近?1.

2.2 HMM擬合原理

馬爾可夫過程是指在已知當(dāng)前狀態(tài)的情況下,未來的演變不依賴于其過去演變的隨機(jī)過程.處于平衡狀態(tài)下生物大分子的不同狀態(tài)間的跳轉(zhuǎn)概率恒定,下一狀態(tài)僅與當(dāng)前狀態(tài)有關(guān),故其態(tài)跳轉(zhuǎn)過程屬于馬爾可夫過程[10].HMM是用來描述含有隱含參量的馬爾可夫過程,對(duì)于smFRET實(shí)驗(yàn)來說,EFRET值是顯性變化的參量,而大分子實(shí)際所處狀態(tài)則是隱含變化的參量,需要通過擬合分析得到.如圖1中的original曲線為原始EFRET-t曲線,由于實(shí)驗(yàn)中不可避免的噪音,EFRET值在一定范圍內(nèi)波動(dòng),某時(shí)刻的EFRET值到底對(duì)應(yīng)大分子的哪一個(gè)結(jié)構(gòu)具有不確定性.HMM擬合算法通過計(jì)算所有可能路徑的概率的最大值,將實(shí)驗(yàn)中在0到1之間連續(xù)變化的EFRET-t曲線擬合為在若干分離EFRET值之間跳轉(zhuǎn)的EFRET-t方波(如圖1中的squared曲線),從而確定平衡狀態(tài)下生物大分子處于各種結(jié)構(gòu)時(shí)EFRET的準(zhǔn)確高度以及EFRET變化的起止時(shí)間.

圖1 EFRET-t曲線(original)與其方波擬合曲線(squared)Fig.1.EFRET-t curve and the square-function fi tting curve.

下面以窮舉法HMM擬合兩態(tài)EFRET-t曲線為例來說明HMM擬合的原理.首先對(duì)EFRET-t曲線中EFRET做柱狀分布圖并使用多峰高斯擬合,此步需要人工輸入態(tài)數(shù)目并給出態(tài)位置的初值,設(shè)擬合后兩態(tài)S1與S2(這兩態(tài)被稱為隱狀態(tài))的高斯峰中心位置(稱為態(tài)位置)分別為E1與E2,相應(yīng)標(biāo)準(zhǔn)差為σ1與σ2,EFRET-t曲線上A點(diǎn)EFRET為EA,此值被稱為觀測(cè)值,則隱狀態(tài)為S1時(shí)觀測(cè)值為EA的概率,同理可得q(AS2),這些概率稱為測(cè)量概率.另外定義pij為相鄰兩幀之間從Si態(tài)跳轉(zhuǎn)到Sj態(tài)的概率(i,j可相等),這樣可得一個(gè)態(tài)間跳轉(zhuǎn)概率矩陣,稱為轉(zhuǎn)移矩陣.然后窮舉EFRET-t曲線所有可能的隱狀態(tài)路徑并通過觀測(cè)概率與轉(zhuǎn)移矩陣計(jì)算每種路徑的總概率大小,其中總概率最大值的情況即為最終擬合路徑.由于轉(zhuǎn)移矩陣未知,所以需要通過擬合結(jié)果來逐步迭代獲得.

以上使用窮舉法來簡(jiǎn)單解釋HMM擬合過程.在實(shí)際運(yùn)用中,通常使用一種HMM簡(jiǎn)化算法,即Baum-Welch算法,來擬合得到轉(zhuǎn)移矩陣、高斯峰中心、標(biāo)準(zhǔn)差以及該算法引入的初始概率.使用Baum-Welch算法可以對(duì)一系列具有相同隱狀態(tài)的觀察序列進(jìn)行統(tǒng)一的參數(shù)學(xué)習(xí),然后再使用Viterbi算法即可對(duì)每一條觀察序列的路徑進(jìn)行擬合[14].

3 實(shí)驗(yàn)條件

smFRET光路如圖2所示.波長(zhǎng)為532 nm的激光激發(fā)反應(yīng)池中被標(biāo)記物上的供體熒光分子(Cy3),供體熒光分子轉(zhuǎn)移部分能量到受體熒光分子(Cy5),兩者發(fā)出的熒光經(jīng)分光光路分為中心波長(zhǎng)在568和675 nm的平行兩束,后被EMCCD(DU897,Andor)收集.供體與受體熒光點(diǎn)的信號(hào)分別成像于CCD視野的左右半?yún)^(qū),通常需要實(shí)驗(yàn)后根據(jù)其坐標(biāo)對(duì)應(yīng)關(guān)系將來自同一被標(biāo)記物上的供體與受體熒光點(diǎn)匹配.

圖2 smFRET系統(tǒng)示意圖Fig.2.Schematic of the smFRET system.

為探究人類端粒上的DNA重復(fù)序列G-四聯(lián)體(G4)結(jié)構(gòu)在K+溶液中的折疊動(dòng)力學(xué).在G4的首尾段分別標(biāo)記供體(Cy3)和受體(Cy5)熒光分子.然后在10 mM/100 mM KCl(M=1 mol/L),10 mM Tris-HCl,pH 7.5緩沖液(下稱K+緩沖液)環(huán)境中測(cè)量G4上熒光分子的smFRET信號(hào),并使用優(yōu)化算法處理實(shí)驗(yàn)結(jié)果.已知G4在K+溶液存在條件下有多個(gè)折疊狀態(tài)[15?18](如圖3),統(tǒng)計(jì)EFRET柱狀分布圖與轉(zhuǎn)移矩陣后可得到更加詳細(xì)的結(jié)構(gòu)與動(dòng)力學(xué)信息.實(shí)驗(yàn)操作為:使用K+緩沖液沖洗固定有biotin-PEG與mPEG的反應(yīng)池,加入10μg/mL的鏈霉親和素(streptavidin)溶液,等待其與生物素(biotin)連接5 min,再次用K+緩沖液沖洗反應(yīng)池后加入含有G4的K+緩沖液,其中G4濃度為10 pM.待G4與反應(yīng)池表面連接30 s后用含有除氧抗淬滅系統(tǒng)[19,20](2.3 mg/mL葡萄糖(glucose),0.1 mg/mL葡萄糖過氧化酶(glucose oxidase,Sigma),0.02 mg/mL過氧化氫酶(catalase,Sigma),1 mM Trolox(Sigma))的K+緩沖液沖洗反應(yīng)池后放置在顯微鏡載物臺(tái)上開始熒光觀察錄制.

圖3 人類端粒DNA重復(fù)序列G4在K+溶液中的多種折疊結(jié)構(gòu)(a)3+1混合1型結(jié)構(gòu);(b)3+1混合2型結(jié)構(gòu);(c)椅式結(jié)構(gòu)Fig.3.Di ff erent folding conformations of human telomeric G4 in K+solution:(a)3+1 hybrid form 1;(b)3+1 hybrid form 2;(c)chair conformation.

4 smFRET數(shù)據(jù)處理方法

當(dāng)前,smFERT數(shù)據(jù)處理方法主要采用人工確定坐標(biāo)對(duì)應(yīng)關(guān)系,缺少對(duì)錯(cuò)誤點(diǎn)的篩選過程,使用HMM算法擬合單條EFRET-t曲線后再匯總擬合全局.其中Ha作為smFRET領(lǐng)域的先驅(qū),提出了使用HMM算法擬合EFRET-t曲線的方法,開拓了一種較為普遍使用的分析方法[10].這些傳統(tǒng)方法在數(shù)據(jù)提取、篩選與擬合步驟上都過于依賴人工處理,具有精度低、工作量大、人為主觀因素影響大等缺陷,需要改進(jìn).

下面將在數(shù)據(jù)提取、數(shù)據(jù)篩選以及數(shù)據(jù)擬合三個(gè)方面分別介紹優(yōu)化的自動(dòng)分析算法,并使用優(yōu)化算法處理G4的smFRET實(shí)驗(yàn)數(shù)據(jù),分析實(shí)驗(yàn)結(jié)果.

4.1數(shù)據(jù)提取

經(jīng)過顯微鏡成像的分光光路后,smFRET的供體與受體熒光點(diǎn)會(huì)分置于所錄制圖像的左右半?yún)^(qū).為了將同一生物大分子上供體與受體的熒光點(diǎn)對(duì)應(yīng)起來,需要建立兩者的坐標(biāo)對(duì)應(yīng)矩陣.傳統(tǒng)方法是通過在錄像疊加圖像上根據(jù)左右半?yún)^(qū)形貌的相似性,人工選取三對(duì)對(duì)應(yīng)的供體、受體熒光點(diǎn)來計(jì)算坐標(biāo)變換矩陣,該方法精度差且效率低,光點(diǎn)密集時(shí)尤其如此.

根據(jù)原理中介紹的ID與IA的負(fù)相關(guān)性,我們提出以下算法以便自動(dòng)準(zhǔn)確確定參考點(diǎn):以512×512像素實(shí)驗(yàn)錄像為例,首先疊加原始實(shí)驗(yàn)錄像獲得平均圖像并對(duì)其降噪,找出平均圖像上的所有光點(diǎn),記錄其坐標(biāo),并從原始實(shí)驗(yàn)錄像中提取相應(yīng)位置光強(qiáng)隨時(shí)間變化的曲線.任取左側(cè)某一供體分子,坐標(biāo)為(xd,yd),求其ID-t曲線與在右側(cè)大致對(duì)應(yīng)范圍(xd+256±5,yd±5)內(nèi)各個(gè)受體分子IA-t曲線的相關(guān)系數(shù)r.若發(fā)現(xiàn)其中出現(xiàn)r小于一定閾值r0時(shí)(如800幀錄像時(shí)r0可取?0.5),我們認(rèn)為該對(duì)熒光信號(hào)高度負(fù)相關(guān),為同一生物大分子上正常發(fā)生smFRET的一對(duì)供體與受體分子發(fā)出,記錄該供體與受體坐標(biāo)作為一對(duì)參考點(diǎn).重復(fù)找到三對(duì)不同參考點(diǎn)后(三對(duì)點(diǎn)最好在圖像上分散)即可得到坐標(biāo)變換矩陣,多次重復(fù)以上過程可確保準(zhǔn)確.上述過程可由程序自動(dòng)完成,而且即使不同實(shí)驗(yàn)中重新調(diào)整分光光路,坐標(biāo)對(duì)應(yīng)關(guān)系發(fā)生改變,也可自動(dòng)生成新的坐標(biāo)對(duì)應(yīng)矩陣而不用重新人工校準(zhǔn).利用該坐標(biāo)變換矩陣即可計(jì)算匹配所有光點(diǎn)的對(duì)應(yīng)點(diǎn),供下一步數(shù)據(jù)篩選.

使用上述優(yōu)化算法前后的結(jié)果比較如圖4.傳統(tǒng)方法計(jì)算坐標(biāo)對(duì)應(yīng)矩陣時(shí)三對(duì)參考點(diǎn)由實(shí)驗(yàn)者人工選點(diǎn)確立,但由于smFRET圖像光點(diǎn)密集,單從圖像左右半?yún)^(qū)形貌相似性上判斷參考點(diǎn),如果出現(xiàn)誤差將會(huì)得到錯(cuò)誤的坐標(biāo)對(duì)應(yīng)矩陣,從而在后續(xù)的匹配中引入錯(cuò)配的情況.如圖4(a)中,供體(Cy3)與受體(Cy5)的光強(qiáng)變化并不呈現(xiàn)FRET應(yīng)有的負(fù)相關(guān)關(guān)系.而使用我們改進(jìn)的方法,依照相關(guān)系數(shù)r自動(dòng)尋找參考點(diǎn)并計(jì)算對(duì)應(yīng)關(guān)系后,可得到同一供體分子(42,172)的正確受體分子位置(297,176),如圖4(b).可以發(fā)現(xiàn),正確與錯(cuò)誤受體坐標(biāo)僅在橫坐標(biāo)相差3個(gè)像素,使用傳統(tǒng)算法在人工點(diǎn)選時(shí)若出現(xiàn)誤差很容易導(dǎo)致坐標(biāo)對(duì)應(yīng)矩陣計(jì)算錯(cuò)誤,最終導(dǎo)致錯(cuò)誤匹配.采用優(yōu)化后的算法則可根據(jù)發(fā)生高質(zhì)量smFRET的數(shù)據(jù)點(diǎn)自動(dòng)計(jì)算坐標(biāo)對(duì)應(yīng)矩陣,即使在實(shí)驗(yàn)光路頻繁調(diào)整、坐標(biāo)對(duì)應(yīng)關(guān)系變動(dòng)的情況下依然能自動(dòng)保證匹配準(zhǔn)確.

圖4 坐標(biāo)對(duì)應(yīng)算法結(jié)果對(duì)比(a)傳統(tǒng)方法可能產(chǎn)生的錯(cuò)配,供體坐標(biāo)為(42,172),受體坐標(biāo)為(300,176);(b)優(yōu)化方法修正后正確配對(duì),供體坐標(biāo)為(42,172),受體坐標(biāo)為(297,176)Fig.4.Comparison of fi lm mapping algorithms:(a)Mismatch caused by traditional algorithm,donor coordinate(42,172),acceptor coordinate(300,176);(b)correction made by the advanced algorithm,donor coordinate(42,172),acceptor coordinate(297,176).

4.2 數(shù)據(jù)篩選

采用正確坐標(biāo)對(duì)應(yīng)關(guān)系提取的數(shù)據(jù)并不都是高質(zhì)量的smFRET數(shù)據(jù),其中可能包含雜質(zhì)發(fā)光、供體熒光分子過早淬滅等情況,相應(yīng)的EFRET-t曲線不能真實(shí)反映生物大分子的狀態(tài)變化,需要進(jìn)行篩選.傳統(tǒng)方法對(duì)于獲得的熒光信號(hào)不能進(jìn)行高效準(zhǔn)確的篩選,導(dǎo)致數(shù)據(jù)質(zhì)量不高或者人為因素影響大.經(jīng)我們優(yōu)化后的算法將對(duì)數(shù)據(jù)進(jìn)行自動(dòng)篩選.同樣依據(jù)ID與IA的負(fù)相關(guān)性,我們計(jì)算ID-t曲線與IA-t曲線的相關(guān)系數(shù)r,通過設(shè)置閾值來去掉r過大的情況.圖5為同一實(shí)驗(yàn)錄像中的三種典型熒光光強(qiáng)數(shù)據(jù),錄像幀數(shù)均為800幀.圖5(a)中r=0.09,供體與受體光強(qiáng)無相關(guān)性且光強(qiáng)在長(zhǎng)時(shí)間內(nèi)幾乎不變,可能為溶液中的雜質(zhì)受激發(fā)光或被標(biāo)記物失去生物活性.圖5(b)中r=?0.3,供體與受體光強(qiáng)在約80 s前具有顯著負(fù)相關(guān)性,說明該處確實(shí)有正確標(biāo)記的生物大分子.但在80 s后,供體分子光強(qiáng)曲線下降(綠色)但受體分子光強(qiáng)(紅色)并未上升,說明并未發(fā)生能量傳遞,而是供體分子發(fā)生了熒光淬滅,80 s之后的EFRET-t曲線不能反映生物大分子的狀態(tài)變化,該數(shù)據(jù)也需要篩除.圖5(c)中r=?0.9,ID與IA的負(fù)相關(guān)性全程顯著,為優(yōu)質(zhì)數(shù)據(jù),需要保留.經(jīng)過測(cè)試,對(duì)于800幀的數(shù)據(jù),將閾值r0設(shè)為?0.5,自動(dòng)篩選保留r

圖5 不同相關(guān)系數(shù)r的熒光光強(qiáng)數(shù)據(jù)(a)r=0.09;(b)r=?0.3;(c)r=?0.9Fig.5.Fluorescence intensity traces with di ff erent correlation coefficient r:(a)r=0.09;(b)r=?0.3;(c)r=?0.9.

圖6 EFRET柱狀分布圖對(duì)比(a)不使用相關(guān)系數(shù)r進(jìn)行數(shù)據(jù)篩選;(b)使用相關(guān)系數(shù)r進(jìn)行自動(dòng)數(shù)據(jù)篩選Fig.6.Comparison of histograms of EFRET:(a)Output without data sifting by the correlation coefficient r;(b)Output using data sifting by the correlation coefficient r.

為證明使用相關(guān)系數(shù)r自動(dòng)篩選熒光光強(qiáng)數(shù)據(jù)的算法不僅高效、客觀,并可有效提高實(shí)驗(yàn)結(jié)果的信噪比,充分利用實(shí)驗(yàn)數(shù)據(jù),我們給出對(duì)同一次實(shí)驗(yàn)錄像不使用與使用該優(yōu)化算法得到的EFRET柱狀分布圖的結(jié)果對(duì)比(EFRET-t曲線通過ID-t與IA-t曲線可得).該實(shí)驗(yàn)中,首尾熒光標(biāo)記的人類端粒G4在100 mM KCl,10 mM Tris-HCl,pH 7.5溶液的環(huán)境中測(cè)量smFRET信號(hào).圖6(a)為使用傳統(tǒng)方法而不使用相關(guān)系數(shù)r進(jìn)行數(shù)據(jù)篩選的結(jié)果.傳統(tǒng)方法為防止出現(xiàn)如圖5(b)中供體過早淬滅的情況,對(duì)所有數(shù)據(jù)均只選取前若干幀的平均EFRET值做柱狀圖,浪費(fèi)了圖5(c)情況下后半段大量有效的EFRET值信息,使得結(jié)果數(shù)據(jù)量少.同時(shí),如果傳統(tǒng)算法不對(duì)所有的光強(qiáng)曲線進(jìn)行人工篩選,則無法排除圖5(a)的情況.綜合以上因素,傳統(tǒng)方法所得到的圖6(a)中柱狀圖噪音大,數(shù)據(jù)量少.圖6(b)為使用相關(guān)系數(shù)r進(jìn)行數(shù)據(jù)篩選的結(jié)果.由于彌補(bǔ)了傳統(tǒng)方法的不足,排除了圖5(a)和圖5(b)所示的兩種主要噪音來源,同時(shí)充分利用了高質(zhì)量數(shù)據(jù),最終得到數(shù)據(jù)量大、高斯峰清晰的EFRET值分布圖,提高了實(shí)驗(yàn)效率與準(zhǔn)確性.

我們經(jīng)優(yōu)化算法處理之后的實(shí)驗(yàn)結(jié)果與他人同樣的實(shí)驗(yàn)結(jié)論一致[21],證明該改進(jìn)算法所得結(jié)果準(zhǔn)確且高效.關(guān)于前人結(jié)果中EFRET=0的受體熒光淬滅與閃爍的數(shù)據(jù)如何去除,我們將在后面的數(shù)據(jù)擬合過程中討論.

4.3 數(shù)據(jù)擬合

得到ID-t與IA-t曲線數(shù)據(jù)后可計(jì)算得到EFRET-t曲線,為了進(jìn)一步分析平衡態(tài)下生物大分子的各種結(jié)構(gòu)以及結(jié)構(gòu)間相互轉(zhuǎn)換的規(guī)律,需要將連續(xù)變化的EFRET-t曲線擬合為在幾個(gè)分離EFRET值之間跳轉(zhuǎn)的方波曲線.前人已經(jīng)提出并使用HMM模型來擬合EFRET-t曲線,然而傳統(tǒng)方法為先對(duì)每條EFRET-t曲線單獨(dú)進(jìn)行擬合,再將這些擬合結(jié)果進(jìn)行匯總再擬合.由于每一條曲線可能會(huì)呈現(xiàn)不同的態(tài)數(shù)目和態(tài)位置,所以單獨(dú)擬合時(shí)需要分別輸入態(tài)數(shù)目和態(tài)位置.這樣擬合完成同一平衡條件下所有EFRET-t曲線后需要再次整合態(tài)位置等參數(shù).整個(gè)過程需要大量的人工參數(shù)輸入,費(fèi)時(shí)費(fèi)力.

使用我們優(yōu)化的方法可以避免這個(gè)問題:首先用Baum-Welch算法對(duì)同一平衡狀態(tài)下的所有EFRET-t曲線進(jìn)行全局的參數(shù)學(xué)習(xí)擬合.由于同一平衡條件下轉(zhuǎn)移矩陣恒定,且所有分子對(duì)應(yīng)相同的隱狀態(tài),因此對(duì)所有EFRET-t曲線進(jìn)行統(tǒng)一的參數(shù)學(xué)習(xí)擬合可以直接得到分子所共有的不同態(tài)的高斯峰中心位置以及不同狀態(tài)間的跳轉(zhuǎn)概率.然后將得到的全局態(tài)位置和跳轉(zhuǎn)概率配合Viterbi算法對(duì)每一條EFRET-t曲線進(jìn)行路徑的擬合.這樣的擬合過程只需要對(duì)態(tài)數(shù)目和態(tài)位置的初值進(jìn)行一次設(shè)定,就可完成對(duì)所有同一平衡狀態(tài)下EFRET-t曲線的全局?jǐn)M合,不僅高效而且準(zhǔn)確.

我們使用優(yōu)化后的擬合程序擬合了G4在100 mM K+溶液條件下的229條EFRET-t曲線,輸入一次態(tài)數(shù)目與位置初值數(shù)據(jù)(4態(tài),EFRET初值分別為0,0.3,0.6,0.8)得到4態(tài)的最終位置分別為:0,0.33,0.59,0.77.

需要注意的是,經(jīng)過數(shù)據(jù)篩選步驟得到的EFRET柱狀圖中,EFRET=0的峰實(shí)際對(duì)應(yīng)供體熒光淬滅或閃爍的狀態(tài).在數(shù)據(jù)擬合完成后,我們可根據(jù)擬合結(jié)果進(jìn)一步篩除該態(tài)的EFRET數(shù)據(jù),得到最終的柱狀分布圖如圖7(a).再對(duì)剩下的3態(tài)跳轉(zhuǎn)概率重新歸一化,就可得到實(shí)際的態(tài)跳轉(zhuǎn)概率從而進(jìn)一步計(jì)算平衡相關(guān)參數(shù).篩除供體熒光淬滅或閃爍狀態(tài)后,G4在100 mM K+溶液條件下的擬合結(jié)果共有3態(tài)(如圖7(a)),其中未折疊態(tài)(EFRET=0.33)平均壽命為5.4 s,兩種混合型折疊構(gòu)型(如圖3(a)和(b),EFRET=0.59)平均態(tài)壽命為10 s,反平行折疊構(gòu)型(如圖3(c),EFRET=0.77)平均態(tài)壽命為51 s.同樣方法處理G4在10 mM K+溶液條件下的119條EFRET-t曲線數(shù)據(jù)(如圖7(b)),也得到3態(tài),其中未折疊態(tài)(EFRET=0.25)平均壽命為15 s,兩種混合型折疊構(gòu)型(如圖3(a)和(b),EFRET=0.57)平均態(tài)壽命為10 s,反平行折疊構(gòu)型(如圖3(c),EFRET=0.77)平均態(tài)壽命為31 s.我們的結(jié)果表明,K+對(duì)于G4的折疊能起到穩(wěn)定的作用.以上結(jié)論與通過傳統(tǒng)方法擬合得到的結(jié)果一致,但是我們的算法獲得了更加清晰的狀態(tài)特征峰[21,22].

圖7 去除受體熒光淬滅與閃爍后的EFRET柱狀分布圖(a)100 mM K+溶液條件;(b)10 mM K+溶液條件Fig.7.Histograms of EFRET.EFRET-t curves that exhibited acceptor bleaching and blinking are excluded:(a)100 mM K+;(b)10 mM K+.

我們優(yōu)化后的算法使用Baum-Welch算法對(duì)所有EFRET-t曲線進(jìn)行全局的參數(shù)學(xué)習(xí),再用全局參數(shù)對(duì)每條EFRET-t曲線進(jìn)行路徑擬合,既減少了人工操作的步驟又保證了擬合結(jié)果準(zhǔn)確.而且,通過擬合結(jié)果可將受體熒光分子淬滅與閃爍的狀態(tài)進(jìn)行有效篩除,讓G4的未折疊態(tài)與各折疊態(tài)更加明顯,最終結(jié)果更加精確.

5 結(jié)論

針對(duì)傳統(tǒng)smFRET數(shù)據(jù)處理算法在數(shù)據(jù)提取、篩選和擬合方面的不足,我們提出了一種優(yōu)化算法.運(yùn)用smFRET的原理,首先引入相關(guān)系數(shù)r來判定圖像上兩熒光點(diǎn)是否為同一生物大分子上的供體、受體的熒光標(biāo)記,自動(dòng)尋找實(shí)驗(yàn)數(shù)據(jù)兩通道的坐標(biāo)對(duì)應(yīng)關(guān)系,使得該步驟不再依賴人工,提高了精度與效率.然后,在數(shù)據(jù)篩選過程中,通過相關(guān)系數(shù)r來自動(dòng)判定并篩除雜質(zhì)發(fā)光、被標(biāo)記物無活性以及供體熒光分子淬滅等錯(cuò)誤數(shù)據(jù)干擾,提高了數(shù)據(jù)質(zhì)量的同時(shí)避免了人工篩選的低效和人的主觀傾向性.最后,通過使用Baum-Welch算法先對(duì)同一平衡條件下的所有數(shù)據(jù)進(jìn)行全局?jǐn)M合,然后使用全局?jǐn)M合參數(shù)對(duì)單條數(shù)據(jù)單獨(dú)擬合路徑,從而獲得態(tài)位置、態(tài)平均壽命,減少了HMM擬合EFRET-t曲線中的人工操作環(huán)節(jié),提高了數(shù)據(jù)處理效率.由擬合路徑可對(duì)受體熒光的淬滅與閃爍數(shù)據(jù)進(jìn)行篩除,進(jìn)一步提高了數(shù)據(jù)的信噪比,為后續(xù)分析實(shí)驗(yàn)結(jié)果提供了便利.

我們將該算法應(yīng)用于G4折疊動(dòng)力學(xué)數(shù)據(jù)的處理,得到了清晰的EFRET柱狀分布圖,展現(xiàn)出G4的未折疊態(tài)與兩個(gè)折疊態(tài),并獲得了不同狀態(tài)下的平均壽命,進(jìn)一步揭示了G4的折疊的多態(tài)性以及各態(tài)在不同K+溶液條件下的穩(wěn)定性變化.日漸成熟的smFRET實(shí)驗(yàn)技術(shù)已經(jīng)向著高精度、高通量的方向發(fā)展,相信更多自動(dòng)化的數(shù)據(jù)處理算法將會(huì)在其中起到重要的作用.

[1]Zhou R,Kozlov A G,Roy R,Zhang J,Korolev S,Lohman T M,Ha T 2011 Cell 146 222

[2]Honda M,Park J,Pugh R A,Ha T,Spies M 2009 Mol.Cell 35 694

[3]Liu C,Mckinney M C,Chen Y H,Earnest T M,Shi X,Lin L J,Ishino Y,Dahmen K,Cann I K,Ha T 2011 Biophy.J.100 1344

[4]Wu J Y,Stone M D,Zhuang X 2010 Nucl.Acids Res.38 e16

[5]Hengesbach M,Kim N K,Feigon J,Stone M D 2012 Angew.Chem.51 5876

[6]Ha T,Tinnefeld P 2012 Annu.Rev.Phys.Chem.63 595

[7]He Z C,Li F,Li M Y,Wei L 2015 Acta Phys.Sin.64 046802(in Chinese)[何志聰,李芳,李牧野,魏來2015物理學(xué)報(bào)64 046802]

[8]Li M Y,Li F,Wei L,He Z C,Zhang J P,Han J B,Lu P X 2015 Acta Phys.Sin.64 108201(in Chinese)[李牧野,李芳,魏來,何志聰,張俊佩,韓俊波,陸培祥2015物理學(xué)報(bào)64 108201]

[9]Roy R,Hohng S,Ha T 2008 Nat.Methods 5 507

[10]Mckinney S A,Joo C,Ha T 2006 Biophys.J.91 1941

[11]Lee N K,Kapanidis A N,Wang Y,Michalet X,Mukhopadhyay J,Ebright R H,Weiss S 2005 Biophys.J.88 2939

[12]Sabanayagam C R,Eid J S,Meller A 2005 J.Chem.Phys.122 061103

[13]Deniz A A,Dahan M,Grunwell J R,Ha T J,Faulhaber A E,Chemla D S,Weiss S,Schultz P G 1999 PNAS 96 3670

[14]Rabiner L R 1989 Proc.IEEE 77 257

[15]Ambrus A,Chen D,Dai J,Bialis T,Jones R A,Yang D 2006 Nucl.Acids Res.34 2723

[16]Gray R D,Trent J O,Chaires J B 2014 J.Mol.Biol.426 1629

[17]Tippana R,Xiao W,Myong S 2014 Nucl.Acids Res.42 8106

[18]Li Y,Liu C,Feng X,Xu Y,Liu B F 2014 Anal.Chem.86 4333

[19]Cordes T,Vogelsang J,Tinnefeld P 2009 J.Am.Chem.Soc.131 5018

[20]Hubner C G,Renn A,Renge I,Wild U P 2001 J.Chem.Phys.115 9619

[21]Lee J Y,Okumus B,Kim D S,Ha T 2005 PNAS 102 18938

[22]Noer S L,Preus S,Gudnason D,Aznauryan M,Mergny J L,Birkedal V 2016 Nucl.Acids Res.44 464

PACS:87.80.Nj,87.15.ad,87.15.Cc,87.15.B–DOI:10.7498/aps.66.118701

An optimization algorithm for single-molecule fl uorescence resonance(smFRET)data processing?

Lü Xi-Ming1)2)Li Hui1)?You Jing1)2)Li Wei1)Wang Peng-Ye1)2)Li Ming1)2)Xi Xu-Guang3)Dou Shuo-Xing1)2)?
1)(Beijing National Laboratory for Condensed Matter Physics,Key Laboratory of Soft Matter Physics,Institute of Physics,Chinese Academy of Sciences,Beijing 100190,China)
2)(School of Physical Sciences,University of Chinese Academy of Sciences,Beijing 100049,China)
3)(College of Life Sciences,Northwest A&F University,Yangling 712100,China)

8 December 2016;revised manuscript

14 March 2017)

The single-molecule fl uorescence resonance energy transfer(smFRET)technique plays an important role in the development of biophysics.Measuring the changes of the fl uorescence intensities of donor and acceptor and of the FRET efficiency can reveal the changes of distance between the labeling positions.The smFRET may be used to study conformational changes of DNA,proteins and other biomolecules.Traditional algorithm for smFRET data processing is highly dependent on manual operation,leading to high noise,low efficiency and low reliability of the outputs.In the present work,we propose an automatic and more accurate algorithm for smFRET data processing.It consists of three parts:algorithm for automatic pairing of donor and acceptor fl uorescence spots based on negative correlation between their intensities;algorithm for data screening by eliminating invalid fl uorescence spots sections;algorithm for global data fi tting based on Baum-Welch algorithm of hidden Markov model.

Based on the law of energy conservation,the light intensity of one pair of donor and acceptor shows a negative correlation.We can use this feature to fi nd the active smFRET pairs automatically.The algorithm will fi rst fi nd out three active smFRET pairs with correlation coefficient lower than the threshold we set.This three active smFRET pairs will provide enough coordinate data for the algorithm to calculate the pairing matrix in the rest of automatic pairing work.After obtaining all the smFRET pairs,the algorithm for data screening will check the correlation coefficient for each pair.The invalid pairs with correlation coefficient higher than the threshold value will be eliminated.The rest of smFRET pairs will be analyzed by the data fi tting algorithm.The Baum-Welch algorithm can be used for learning the global parameters.The global parameters we obtained will then be used to fi t each FRET-time curve with Viterbi algorithm.The global parameter learning part will help us fi nd the speci fi c FRET efficiency for each state and the curve fi tting part will provide more kinetic parameters.

The optimization algorithm signi fi cantly simpli fi es the procedures of manual operation in the traditional algorithm and eliminate several types of noises from the experimental data automatically.We apply the new optimization algorithm to the analyses of folding kinetics data for human telomere repeat sequence,the G-quadruplex DNA.It is demonstrated that the optimization algorithm is more efficient to produce data with higher S/N ratio than the traditional algorithm.The fi nal results reveal clearly the folding of G-quadruplex DNA in multiple states that are in fl uenced by the K+concentration.

smFRET,data processing algorithm,G-quadruplex DNA,folding kinetics

10.7498/aps.66.118701

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):11674383,11474346,11274374)、國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(批準(zhǔn)號(hào):2013CB837200)和國(guó)家重點(diǎn)研究發(fā)展計(jì)劃(批準(zhǔn)號(hào):2016YFA0301500)資助的課題.

?通信作者.E-mail:huili@iphy.ac.cn

?通信作者.E-mail:sxdou@iphy.ac.cn

?2017中國(guó)物理學(xué)會(huì)Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.11674383,11474346,11274374),the National Basic Research Program of China(Grant No.2013CB837200),and the National Key Research and Development Program of China(Grant No.2016YFA0301500).

?Corresponding author.E-mail:huili@iphy.ac.cn

?Corresponding author.E-mail:sxdou@iphy.ac.cn