基于深度學習的冠狀病毒刺突蛋白拉曼特征峰的理論研究

倪 爽， 溫家星， 周民杰， 黃景林， 樂 瑋， 陳 果，何智兵， 李 波， 趙松楠， 趙宗清， 杜 凱

中國工程物理研究院激光聚變研究中心， 四川 綿陽 621900

引 言

新冠疫情對全球的經(jīng)濟造成了巨大破壞。 為了有效控制新冠疫情， 快速檢測新冠病毒(嚴重急性呼吸綜合癥冠狀病毒2, severe acute respiratory syndrome coronavirus 2, SARS-CoV-2)是一個急需解決的問題。 新冠病毒的刺突蛋白是病毒攻擊人體的鑰匙， 且在病毒表面大量分布， 因此成為拉曼光譜技術(shù)檢測新冠病毒的檢測點[1-3]。 要實現(xiàn)刺突蛋白的拉曼技術(shù)檢測， 關(guān)鍵之一在于構(gòu)建刺突蛋白的拉曼特征峰。 鑒于純刺突蛋白的拉曼譜圖較難獲得且成本較高， 需要從理論上快速構(gòu)建刺突蛋白拉曼特征峰。 此外， 人類已知可以感染的七種冠狀病毒刺突蛋白的結(jié)構(gòu)[4-9]相近[圖1(a)]， 是否可以通過拉曼技術(shù)區(qū)分他們對于新冠病毒的準確檢測是一個十分重要的問題。 基于理論構(gòu)建的刺突蛋白拉曼特征峰， 可以分析七種冠狀病毒刺突蛋白拉曼特征峰的不同， 為實驗提供指導。

(1)

式(1)中：εi為局域模頻率，βij為局域模之間的相互作用， 激子模型哈密頓量可以寫成矩陣形式

(2)

式(2)中：n為局域模的數(shù)量， 對角化H矩陣可以獲得激子的頻率及其對應特征矢量(局域模展開為簡正模的線性展開系數(shù))。 對于蛋白酰胺Ⅰ峰而言， 其局域模近似為組成蛋白骨架的每個酰胺單元的酰胺Ⅰ振動模。

本文關(guān)注的重點是定性分析七種冠狀病毒刺突蛋白拉曼特征峰的差異而非絕對值且在激子模型中H矩陣對角元的值遠大于非對角元的值。 因此， 本文采用簡化的激子模型， 僅考慮結(jié)構(gòu)變化對H矩陣對角元的校正而將非對角元舍去。

為了校正對角元， 需要對七種刺突蛋白結(jié)構(gòu)中的每個酰胺單元計算拉曼譜圖， 每個刺突蛋白有數(shù)千個酰胺單元， 逐個計算拉曼譜圖費時且不利于統(tǒng)一誤差。 江俊[15-16]等近期提出了基于機器學習構(gòu)建分子結(jié)構(gòu)和拉曼性質(zhì)(頻率、 強度)之間映射關(guān)系的策略。 其先通過分子動力學構(gòu)建分子的動力學結(jié)構(gòu)， 然后計算結(jié)構(gòu)的拉曼頻率以及強度， 最后使用機器學習擬合結(jié)構(gòu)和性質(zhì)的映射關(guān)系。 本文基于此策略， 構(gòu)建酰胺單元結(jié)構(gòu)和拉曼峰頻率以及強度之間的模型。

文獻中多數(shù)研究集中在酰胺Ⅰ峰， 對酰胺Ⅲ拉曼峰研究較少。 本文基于深度學習技術(shù)， 從理論上構(gòu)建蛋白酰胺單元結(jié)構(gòu)和酰胺Ⅰ、 Ⅲ拉曼特征峰的映射關(guān)系， 然后統(tǒng)計七種冠狀病毒刺突蛋白的結(jié)構(gòu)差異， 帶入模型中獲得拉曼特征峰。 最后通過洛倫茲線型展開譜線獲得譜圖， 并比較譜圖的差異。

1 模型的構(gòu)建和計算方法

使用N-甲基乙酰胺(NMA)分子來模擬蛋白的酰胺單元[圖1(b， c)]， 構(gòu)建酰胺Ⅰ、 Ⅲ峰的深度學習模型。 首先， 通過VASP軟件， 采用分子動力學的方法構(gòu)建10 000個NMA分子隨時間演化的結(jié)構(gòu)， 未考慮H2O的影響。 然后通過Gaussian09軟件計算這些結(jié)構(gòu)對應的拉曼譜峰， 最后通過深度學習技術(shù)構(gòu)建NMA分子結(jié)構(gòu)特征和酰胺Ⅰ、 Ⅲ峰頻率以及強度的映射關(guān)系。 VASP計算采用GGA-PBE泛函， 周期性的邊界條件為14.6×14.6×14.6 ?3， 平面波的截斷能設(shè)置為400 eV， 能量收斂標準為1×10-5eV， 分子動力學模擬選擇正則系綜(NVT)， 時間步長1 fs， 溫度為300 K， 總共模擬10 000步。 Gaussian09計算拉曼譜峰基于B3LYP雜化泛函水平， 基組使用6-311++G(d, p)。

深度學習模型由1個輸入層， 3個隱藏層以及1個輸出層組成， 對于每一個隱藏層， 使用線性修正單元(Rectified Linear Unit)激活函數(shù)。 對于NMA分子， 選用10個結(jié)構(gòu)特征來描述(4個鍵長， 4個鍵角， 2個二面角， 如圖2)。 為了增強模型的魯棒性， 對輸入特征以及輸出結(jié)果進行了歸一化處理。

圖2 用來預測NMA分子拉曼位移以及強度的描述符

2 結(jié)果與討論

2.1 酰胺Ⅰ、 Ⅲ峰模型

基于分子動力學構(gòu)建了10 000個NMA分子結(jié)構(gòu)， 計算了其均方根誤差， 為1.06 ?(圖3)， 這說明通過動力學演化得到的NMA分子結(jié)構(gòu)變化較大且相關(guān)性較小。

圖3 NMA分子結(jié)構(gòu)的均方根誤差

圖4 NMA描述符間的皮爾遜相關(guān)系數(shù)

圖5 酰胺Ⅰ峰的拉曼位移(a)、 強度 (b)以及酰胺Ⅲ峰的拉曼位移(c)、 強度(d)通過深度學習模型和DFT計算的結(jié)果比較

2.2 七種冠狀病毒刺突蛋白骨架的結(jié)構(gòu)特征

基于實驗上的冠狀病毒刺突蛋白結(jié)構(gòu)[4-9](PDB code: 6u7h, 5i08, 5szs, 6ohw, 5x59, 5x58, 6vsb)， 統(tǒng)計其骨架特征鍵長、 鍵角、 二面角的分布。

圖6 七種冠狀病毒刺突蛋白骨架的鍵長特征統(tǒng)計分布圖

從鍵角[圖7(a—d)]的分布可以看出： 鍵角的分布較為廣泛且均勻， 說明鍵角變化的力常數(shù)較小， 七種冠狀病毒的鍵角變化相差不大， 鍵角對酰胺特征峰的影響較小。

圖7 七種冠狀病毒刺突蛋白骨架的鍵角和二面角特征統(tǒng)計分布圖

從二面角的分布[圖7(e,f)]中可以看出， 七種冠狀病毒刺突蛋白的二面角均在180°(3.14弧度)附近， 這是由于酰胺平面共軛導致的。

2.3 七種冠狀病毒刺突蛋白的拉曼特征峰比較

根據(jù)前面獲得的酰胺Ⅰ、 Ⅲ峰的模型(酰胺單元的10個特征和酰胺Ⅰ、 Ⅲ振動峰的頻率以及強度的映射關(guān)系)， 從七種冠狀病毒刺突蛋白的實驗結(jié)構(gòu)中計算出每個酰胺單元的10個特征， 帶入到模型中獲得每個酰胺單元的酰胺Ⅰ、 Ⅲ拉曼峰(振動頻率和強度)， 最后通過洛倫茲線型將每個酰胺單元的Ⅰ、 Ⅲ拉曼峰展開獲得七種冠狀病毒刺突蛋白的酰胺Ⅰ、 Ⅲ譜帶(圖8)。 從圖中可以看出， 七種冠狀病毒刺突蛋白的酰胺Ⅰ、 Ⅲ譜帶根據(jù)最高峰頻率可以各分成三個組。 對于酰胺Ⅰ峰， SARS-CoV-2， SARS-CoV， MERS-CoV刺突蛋白頻率相近， 其頻率值在1 636～1 637 cm-1區(qū)間； HCoV-HKU1， HCoV-NL63刺突蛋白頻率相近， 其頻率值在1 657～1 658 cm-1區(qū)間； HCoV-229E， HCoV-OC43刺突蛋白頻率相近， 其頻率值在1 673～1 674 cm-1區(qū)間。 對于酰胺Ⅲ峰， SARS-CoV-2， SARS-CoV， MERS-CoV刺突蛋白頻率相近， 其頻率值在1 263～1 265 cm-1； HCoV-229E， HCoV-HKU1， HCoV-NL63刺突蛋白頻率相近， 其頻率值在1 272～1 275 cm-1； HCoV-OC43刺突蛋白單獨一個頻率， 其頻率值為1 285 cm-1。

圖8 七種冠狀病毒刺突蛋白的酰胺Ⅰ、 Ⅲ拉曼特征峰比較

根據(jù)上面的分析， 可以根據(jù)酰胺Ⅰ、 Ⅲ峰的頻率劃分七種冠狀病毒， 如圖9所示， 七種冠狀病毒分為四組： SARS-CoV-2， SARS-CoV， MERS-CoV在同一個組； HCoV-HKU1， HCoV-NL63為一組； HCoV-229E一組； HCoV-OC43一組。 不同組之間其刺突蛋白特征峰的數(shù)值差異較大， 可以區(qū)分開來。 同組中的刺突蛋白特征峰的數(shù)值差異較小， 較難區(qū)分。

圖9 通過酰胺Ⅰ、 Ⅲ特征峰頻率區(qū)分七種冠狀病毒

3 結(jié) 論

基于深度學習的技術(shù)構(gòu)建了刺突蛋白酰胺Ⅰ、 Ⅲ特征拉曼峰模型， 結(jié)合實驗上的冠狀病毒刺突蛋白結(jié)構(gòu)， 獲得了七種冠狀病毒刺突蛋白的酰胺Ⅰ、 Ⅲ拉曼特征峰， 通過比較七種冠狀病毒刺突蛋白的拉曼特征峰， 可以把七種冠狀病毒分為四組： SARS-CoV-2， SARS-CoV， MERS-CoV一組； HCoV-HKU1， HCoV-NL63一組； HCoV-229E一組； HCoV-OC43一組。 不同組的冠狀病毒特征峰差異較大， 可以區(qū)分開來； 同一組的冠狀病毒特征峰差異較小， 較難區(qū)分。