李建勛，李 鑫，謝 康，徐德晨，李 杰

(哈爾濱工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 生物信息學(xué)研究中心，哈爾濱 150001)

2019年12月出現(xiàn)以干咳、發(fā)熱、渾身乏力等為典型癥狀的新型肺炎疫情，中國國家衛(wèi)生健康委員會(huì)暫時(shí)將其命名為新型冠狀病毒肺炎。2020年1月31日，世界衛(wèi)生組織(WHO)宣布將新型冠狀病毒肺炎列入“國際經(jīng)濟(jì)公共衛(wèi)生事件”，并于2月11日將之命名為COVID-19[1]。同日，國際病毒分類委員會(huì)聲明，將引發(fā)疾病的新型冠狀病毒命名為SARS-CoV-2[2]。

1月份，新型冠狀病毒肺炎擴(kuò)散到國內(nèi)多個(gè)省份和地區(qū)，后續(xù)世界各國也爆發(fā)了疫情[1]。截止到4月16日，中國累計(jì)確診COVID-19共83 799例，但疫情已得到基本的控制，現(xiàn)存確診1 891例，其中境外輸入1 534例，占絕大部分[注]數(shù)據(jù)來源于國家及各省市地區(qū)衛(wèi)健委，丁香醫(yī)生整理(https://ncov.dxy.cn/ncovh5/view/pneumonia?)。然而，與此同時(shí)，新型冠狀病毒已在世界多國相繼爆發(fā)，全球累計(jì)確診達(dá)2 047 104例，現(xiàn)存確診140 675例，世界整體形勢不容樂觀。另外，國內(nèi)外有多位相關(guān)領(lǐng)域的醫(yī)學(xué)專家表示，新型冠狀病毒短期內(nèi)不會(huì)消失，可能與人類長期共存。

冠狀病毒是一種包膜的單鏈RNA病毒，它的受體結(jié)合域能快速突變和重組，以適應(yīng)并感染大量物種[3]。自病毒疫情爆發(fā)以來，學(xué)術(shù)界對于新型冠狀病毒的基本傳染數(shù)R0(Basic reproduction number)有許多研究和討論，在多篇論文中給出的取值呈越來越大的趨勢。4月7日，美國疾病管制局期刊《新興傳染病》發(fā)表了關(guān)于新型冠狀病毒的最新研究結(jié)論[4]，研究人員通過兩種建模方法(到達(dá)模型和病例計(jì)數(shù)模型)，最終計(jì)算得出新型冠狀病毒肺炎R(shí)0中位數(shù)為5.7(95%CI 3.8-8.9)，也就是說，一名新型冠狀病毒肺炎患者可以傳染5.7人，是之前認(rèn)為的2～3倍。

在這種背景下，對新型冠狀病毒序列的變異情況進(jìn)行監(jiān)控、分析及可視化，就顯得尤為重要[5-6]。國內(nèi)外已有一些相關(guān)的軟件系統(tǒng)向用戶開放，但是這些已有系統(tǒng)沒有從時(shí)空角度進(jìn)行分析的功能或相關(guān)功能仍有進(jìn)一步完善的空間，例如北京志諾維思公司的“戰(zhàn)新冠”系統(tǒng)(https://fight-sars2.genowis.com/)，提供了對病毒序列的變異分析及可視化，但沒有從時(shí)空角度進(jìn)行分析的功能；國家生物信息中心開發(fā)的2019新型冠狀病毒信息庫線上系統(tǒng)(https://bigd.big.ac.cn/ncov/)，提供了包含時(shí)空角度的病毒序列變異分析功能，但是其時(shí)空變異分析功能是使用熱力圖來實(shí)現(xiàn)，完全是靜態(tài)圖的形式，也不支持用戶對基準(zhǔn)序列、時(shí)間區(qū)間和地域范圍進(jìn)行選擇，功能上仍然需要擴(kuò)展和進(jìn)一步完善。本文所設(shè)計(jì)和實(shí)現(xiàn)的新型冠狀病毒變異時(shí)空分析系統(tǒng)用于解決這一問題，能夠動(dòng)態(tài)顯示不同時(shí)間和地域病毒序列的變異境況(系統(tǒng)現(xiàn)已上線，網(wǎng)址http://bionet.org.cn/NCP-web/)。

1 系統(tǒng)的總體功能

新型冠狀病毒變異時(shí)空分析系統(tǒng)的主要功能有三個(gè)：病毒序列來源統(tǒng)計(jì)、病毒序列變異統(tǒng)計(jì)及病毒序列比對。

1.1 系統(tǒng)首頁

系統(tǒng)首頁主要包括系統(tǒng)導(dǎo)航欄、系統(tǒng)信息簡介、系統(tǒng)三個(gè)主要功能模塊的簡介。系統(tǒng)首頁的部分截圖如圖1所示。

圖1 系統(tǒng)首頁截圖展示Fig.1 Screenshot of the homepage

1.2 病毒序列來源統(tǒng)計(jì)

自新型冠狀病毒爆發(fā)以來，各國研究人員在SARS-CoV-2測序方面做了大量的工作[7-12]。全球共享流感數(shù)據(jù)倡議組織(GISAID)將各國研究人員提交的病毒序列數(shù)據(jù)進(jìn)行了整理和匯總，并在其官方網(wǎng)站(https://www.gisaid.org/)進(jìn)行公開。本文使用的源數(shù)據(jù)來自GISAID官網(wǎng)[13]。

新型冠狀病毒序列來源統(tǒng)計(jì)功能旨在對來自不同國家和地區(qū)的新型冠狀病毒序列的數(shù)量進(jìn)行統(tǒng)計(jì)、以地圖和餅狀圖的形式進(jìn)行可視化。

對新型冠狀病毒序列來源的統(tǒng)計(jì)和可視化有兩種形式——地圖和餅狀圖(見圖2)。

圖2 地圖形式的序列來源可視化Fig.2 Visualization of sequence sources in map form

圖2展示了中國地圖和世界地圖上的新型冠狀病毒序列來源統(tǒng)計(jì)，當(dāng)用戶將鼠標(biāo)懸浮在目標(biāo)國家或者地區(qū)上時(shí)，會(huì)以彈框的形式展示該國家或地區(qū)已提交的序列總數(shù)。

圖3展示了餅狀圖形式下世界各國提交序列數(shù)量統(tǒng)計(jì)和中國各省提交序列數(shù)量統(tǒng)計(jì)。此外，系統(tǒng)還提供將以上二者結(jié)合的餅狀圖可視化。同時(shí)，當(dāng)用戶將鼠標(biāo)懸浮在圖中的扇形上時(shí)，系統(tǒng)會(huì)以彈框的形式展示該國家或省份提交的序列數(shù)占總數(shù)的百分比。

根據(jù)系統(tǒng)給出的結(jié)果(截止到3月19日的數(shù)據(jù))，可以發(fā)現(xiàn)世界范圍內(nèi)共提交SARS-CoV-2序列934條，中國提交的病毒序列最多，共267條，占28.59%；國內(nèi)廣東省提交的序列最多，共90條，占33.71%，湖北省次之，共42條，占15.73%。希望世界各國及時(shí)共享更多的SARS-CoV-2序列，促進(jìn)COVID-19的研究。

1.3 病毒序列變異統(tǒng)計(jì)

系統(tǒng)將來自GISAID的序列數(shù)據(jù)進(jìn)行篩選和預(yù)處理之后，將所有序列轉(zhuǎn)換成等長的形式存入數(shù)據(jù)庫(具體的篩選與預(yù)處理過程及方法見3.1節(jié)病毒序列數(shù)據(jù)的處理)。在此基礎(chǔ)上，系統(tǒng)在第一次啟動(dòng)時(shí)，會(huì)比對數(shù)據(jù)庫中所有序列，將每一個(gè)位點(diǎn)上頻次最高的核苷酸作為該位點(diǎn)的標(biāo)準(zhǔn)核苷酸，從而得到一條統(tǒng)計(jì)意義上的標(biāo)準(zhǔn)序列(Standard sequence)，后續(xù)可以此標(biāo)準(zhǔn)序列為基準(zhǔn)進(jìn)行變異分析(也支持由用戶指定某一序列作為標(biāo)準(zhǔn)序列的變異分析)。

該模塊下又分為序列變異統(tǒng)計(jì)分析和序列變異時(shí)空分析兩個(gè)子功能模塊。

1.3.1 序列變異統(tǒng)計(jì)分析

圖4展示了系統(tǒng)為用戶提供的選擇入口，用戶可在該入口下指定標(biāo)準(zhǔn)序列、選擇待分析的目標(biāo)序列集合、指定要分析的起始位點(diǎn)和終止位點(diǎn)。通過圖4的設(shè)置面板完成分析基本參數(shù)選定并點(diǎn)擊分析按鈕，系統(tǒng)會(huì)給出病毒變異統(tǒng)計(jì)分析的兩個(gè)結(jié)果，如圖5和圖6所示。

圖4 選擇基準(zhǔn)序列與待分析序列Fig.4 Choice for base sequence and sequences to be analyzed

圖5 在選擇位點(diǎn)之間的變異病毒株數(shù)分布Fig.5 Distribution of variant virus strains between selected sites

圖5展示了用戶選定的目標(biāo)序列集合與基準(zhǔn)序列相比，具有不同變異位點(diǎn)數(shù)的病毒株數(shù)分布。除具有不同變異頻數(shù)的病毒株數(shù)分布統(tǒng)計(jì)(見圖5)之外，系統(tǒng)還會(huì)對病毒序列在不同時(shí)間和地區(qū)的平均變異率做統(tǒng)計(jì)(見圖6)。該功能模塊會(huì)將所有序列按照國家進(jìn)行分組，每一組用不同的顏色加以區(qū)分，組內(nèi)序列按照采樣時(shí)間進(jìn)行排序，給出所在地區(qū)在每一個(gè)采樣日的平均變異率(每一個(gè)采樣日內(nèi)，對所在地區(qū)當(dāng)天采樣的所有病毒序列的變異率求取平均值)。該柱狀圖可以直觀地展示某一地區(qū)內(nèi)病毒序列的變異率隨著時(shí)間遞增的變化情況。

根據(jù)系統(tǒng)給出的分析結(jié)果(截止到3月13日的病毒序列數(shù)據(jù))，由圖5可以發(fā)現(xiàn)，大部分序列的變異位點(diǎn)數(shù)少于10個(gè)，單條序列的最多變異數(shù)為31，有1條；由圖6可以發(fā)現(xiàn)，中國在2020年1月29日采樣的病毒序列的平均變異率最高，為0.09%，其余大部分序列變異率波動(dòng)不大，絕大部分序列的變異率在0.04%以下。

1.3.2 序列變異時(shí)空分析

隨著疫情的擴(kuò)散，病毒在不同地區(qū)的變異率差異、在不同時(shí)間點(diǎn)的變異率變化這兩個(gè)信息尤為重要，可以幫助科研人員對病毒進(jìn)行更加精準(zhǔn)的分析和研究，從而為疾病的研究和防治提供指導(dǎo)。該功能模塊可以實(shí)現(xiàn)這一目標(biāo)。

與病毒序列變異統(tǒng)計(jì)模塊的參數(shù)面板相比，圖7所示的參數(shù)面板具有更豐富的選項(xiàng)，包括標(biāo)準(zhǔn)序列、要分析的時(shí)間區(qū)間(以采樣時(shí)間為準(zhǔn))、要分析的地域(以采樣地域?yàn)闇?zhǔn)，可選擇全世界的不同國家，也可以選擇全國的不同省份)。在提交參數(shù)選項(xiàng)之后，系統(tǒng)會(huì)給出兩個(gè)分析結(jié)果圖，如圖8和圖9所示。

圖7 變異時(shí)空分析的基本參數(shù)選項(xiàng)Fig.7 Basic parameter options for spatiotemporal analysis of variation

圖8 不同時(shí)間下各地區(qū)的變異率Fig.8 Variation rates in different regions at different times

圖8展示了在2020年1月21日到2020年2月3日之間十五個(gè)不同國家的病毒序列平均變異率(對所在國家或地區(qū)的截止到觀測日期所采樣的所有病毒序列的變異率求取平均值)的變化情況。在系統(tǒng)中實(shí)際操作時(shí)，該條形圖是時(shí)間為變量進(jìn)行動(dòng)態(tài)播放的，在播放到最后一個(gè)時(shí)間節(jié)點(diǎn)時(shí)就會(huì)停止輪播。根據(jù)圖8我們可以發(fā)現(xiàn)，在1月21日到2月3日之間，中國的平均變異率為0.013%，大致在平均水平左右；平均變異率最高的國家是越南，為0.030%；平均變異率最低的國家是柬埔寨，為0.003%。

圖9展示的是序列變異時(shí)空分析模塊給出的第二個(gè)分析結(jié)果。與圖8相同，該折線圖在實(shí)際操作時(shí)也是隨時(shí)間變化動(dòng)態(tài)播放。該圖展示的是1月21日到2月3日之間，十五個(gè)不同國家整體的病毒序列變異情況。其中黃色折線為最大變異率、綠色折線為最小變異率、藍(lán)色折線為平均變異率(對全球所有國家或者國內(nèi)所有省份的截止到觀測日期所采樣的所有病毒序列的變異率求取平均值)，鼠標(biāo)在圖上懸浮時(shí)會(huì)給出詳細(xì)信息。根據(jù)圖9可以發(fā)現(xiàn)，1月21日到2月3日之間，病毒序列的最大變異率為0.09%，結(jié)合1.3.1節(jié)的圖6，可知該序列在1月29日采樣于中國廣東；最小變異率為0%，該序列在1月23日被采樣；平均變異率在0.01%左右波動(dòng)。

1.4 病毒序列比對

除病毒序列的來源統(tǒng)計(jì)和病毒序列變異統(tǒng)計(jì)兩個(gè)模塊外，系統(tǒng)還設(shè)有病毒序列比對模塊。該模塊下又細(xì)分為病毒序列基本信息展示和病毒序列差異比對兩個(gè)子模塊。

1.4.1 病毒序列基本信息展示

GISAID官網(wǎng)除了提供病毒序列之外，還給出了一些病毒序列相關(guān)的信息，其中包括病毒序列號(hào)、病毒名稱、采樣地點(diǎn)、采樣實(shí)驗(yàn)室、提交實(shí)驗(yàn)室、作者信息、采樣時(shí)間，系統(tǒng)以表格的形式對這些信息做了展示(見圖10)。表格的每一行對應(yīng)一個(gè)病毒序列，所有病毒序列按照國家分組并按照采樣時(shí)間升序排序。同時(shí)，在實(shí)際操作系統(tǒng)時(shí)，表中每一行前面都有一個(gè)單選框，用戶在感興趣的單選框內(nèi)打鉤，即可選中該行所對應(yīng)的序列，進(jìn)行病毒序列差異比對。

圖9 不同時(shí)間下所有地區(qū)的累計(jì)變異率Fig.9 Cumulative variation of all regions at different times

圖10 序列的基本信息表Fig.10 Sequence basic information table

1.4.2 病毒序列差異比對

圖11展示的是系統(tǒng)對EPI_ISL_403932、EPI_ISL_408484、EPI_ISL_410984、 EPI_ISL_415462 和EPI_ISL_415652共五條序列在1到30145這一位點(diǎn)區(qū)間上進(jìn)行差異比對的結(jié)果，最終以表格的形式將所有存在差異的位點(diǎn)一一列出。根據(jù)表格可以發(fā)現(xiàn)，上述五條序列在298、1 247、3 094、8 839、9 495、11 143、13 111、14 468、15 384、19 251、23 467、24 098、24 926、26 793、28 141、28 208、29 159這些位點(diǎn)上存在差異，即這些位點(diǎn)都是可能發(fā)生了變異的位點(diǎn)。

圖11 病毒序列差異比對Fig.11 Virus sequence difference alignment

2 系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)

2.1 數(shù)據(jù)庫

系統(tǒng)使用了關(guān)系型數(shù)據(jù)庫MySQL，以每一條序列作為一個(gè)實(shí)體，將所有數(shù)據(jù)都整合到一張表中。即表中綜合了每個(gè)序列的基本信息(序列名、病毒名、采樣地點(diǎn)、采樣時(shí)間和提交時(shí)間等)和核苷酸排布信息(以字符串形式存儲(chǔ)，在數(shù)據(jù)庫中字段的類型為不限長度的text類型)。表的主要字段如圖12所示。

圖12 表格主要字段Fig.12 Main fields of the data table

2.2 服務(wù)端

系統(tǒng)服務(wù)器端使用Spring、SpringMVC配合MyBatis框架作為主要的技術(shù)進(jìn)行實(shí)現(xiàn)[14]，并將項(xiàng)目部署在阿里云ECS服務(wù)器上提供服務(wù)。

2.3 客戶端

系統(tǒng)的客戶端主要采用jquery、bootstrap、echarts框架[15]進(jìn)行開發(fā)。jquery框架作為一個(gè)快速簡潔的輕量型框架，提供了完備的功能接口，支持豐富的插件便于擴(kuò)展，同時(shí)能夠較好的支持服務(wù)端與客戶端得數(shù)據(jù)交互。bootstrap框架具有響應(yīng)式布局設(shè)計(jì)，能讓網(wǎng)站兼容不同分辨率設(shè)備，給用戶較好的體驗(yàn)。而echarts框架能夠?qū)⑽覀兊臄?shù)據(jù)以圖表可視化的方式進(jìn)行呈現(xiàn)，有利于用戶對數(shù)據(jù)進(jìn)行閱讀和分析。

3 數(shù)據(jù)處理與分析

3.1 病毒序列數(shù)據(jù)的處理

為了尋找發(fā)生變異的位點(diǎn)，我們首先將新型冠狀病毒序列進(jìn)行了多重序列比對。本研究中使用的多序列對比工具是MUSCLE[16]，軟件版本為3.8.31。比起其它常用的序列對比工具，如ClustalW[17]、T-coffee[18]，MUSCLE在保證了較高的準(zhǔn)確度的同時(shí)，具有更快的多序列比對速度。在參數(shù)選擇上，我們設(shè)置了參數(shù)-diags來適用于高度同源的新型冠狀病毒序列，同時(shí)最大迭代次數(shù)(-maxiters)選擇16來保證多序列對比的準(zhǔn)確性，而其它參數(shù)選擇了默認(rèn)設(shè)置。

多重序列比對之后，多個(gè)新型冠狀病毒序列間同源序列的相同核苷酸將處在相同的位點(diǎn)上，而如果原始序列在該位點(diǎn)中不存在核苷酸，將用字符“-”來表示。經(jīng)過多重序列比對的序列都將擁有相同的序列長度。

對于對比后的序列，根據(jù)編碼區(qū)異常字符(N)的數(shù)量以及異常缺失位點(diǎn)數(shù)量，我們設(shè)置30個(gè)為病毒序列的篩選條件，異常位點(diǎn)大于篩選條件值的序列將被遺棄。

將篩選后的比對結(jié)果轉(zhuǎn)化為核苷酸矩陣，每條序列為行，每個(gè)位點(diǎn)為列，尋找每一列是否有超過一種核苷酸，如果結(jié)果為是，則該位點(diǎn)為有變異。

經(jīng)如上方法對來自GISAID的數(shù)據(jù)進(jìn)行預(yù)處理之后，即可將所有病毒序列轉(zhuǎn)換為等長的形式，在此基礎(chǔ)上即可按照2.1節(jié)所述內(nèi)容進(jìn)行數(shù)據(jù)的整合和數(shù)據(jù)庫表的設(shè)計(jì)與填充。

3.2 病毒變異的時(shí)空分析

隨著時(shí)間線的延長，各國提交的序列也逐漸增多，相應(yīng)的各地區(qū)病毒序列的平均變異率也在逐漸增加，如圖9所示。這也符合包膜的陽性單鏈RNA病毒在傳播過程中容易變異的特點(diǎn)。我們也基于圖8所對應(yīng)的功能模塊對截止到3月13日的不同國家或地區(qū)的平均變異率做了比較分析，找出了平均變異率較高和較低的部分國家以及國內(nèi)部分省份，分別如表1和表2所示?？梢钥闯?，到目前為止，病毒序列的平均變異率處于一個(gè)較小的波動(dòng)范圍內(nèi)，病毒仍然處于較為穩(wěn)定的狀態(tài)。

表1 病毒變異率較高和較低的部分國家Table 1 Countries with higher and lower virus variation rates

表2 病毒變異率較高和較低的中國部分省份Table 2 Provinces in China with higher and lower virus variation rates

結(jié)合在選擇位點(diǎn)之間的變異病毒株數(shù)分析結(jié)果(圖5所對應(yīng)的功能模塊)，大部分有變異發(fā)生的病毒序列的變異位點(diǎn)數(shù)量小于12個(gè)，這也印證了目前病毒變異率較小這一結(jié)果。

進(jìn)一步比較不同地區(qū)不同時(shí)間的病毒變異率，發(fā)現(xiàn)隨著地區(qū)與時(shí)間的不同，病毒平均變異率也有差異，但總體仍處于在較小范圍內(nèi)波動(dòng)的狀態(tài)。

對于病毒序列變異分析，用戶可以選擇感興趣的起始位點(diǎn)和終止位點(diǎn)對序列進(jìn)行分析，比如針對病毒的S蛋白編碼區(qū)，在本次多重序列對比中對應(yīng)起始與終止位置分別為21 627和25 448，在系統(tǒng)中選擇后，可以看到變異數(shù)量的信息，其中變異位點(diǎn)數(shù)量(SNP)為0和1的病毒數(shù)量最多，對應(yīng)病毒株數(shù)分別為302和194，其余SNP數(shù)分別為2、3和7，對應(yīng)病毒株數(shù)分別為28、6和1，大部分病毒序列均未發(fā)生明顯改變。

4 總結(jié)

本系統(tǒng)對來源于全球共享流感數(shù)據(jù)倡議組織(GISAID)的新型冠狀病毒數(shù)據(jù)進(jìn)行整合與分析，主要實(shí)現(xiàn)了以下功能：

1)對數(shù)據(jù)的來源進(jìn)行可視化分析。

2)對序列的基本信息進(jìn)行展示，并支持選定序列之間核苷酸比較。

3)以動(dòng)態(tài)圖的形式呈現(xiàn)不同時(shí)間下各地區(qū)的病毒變異率，以及所有地區(qū)的累計(jì)變異率，使用戶能夠直觀感受病毒的時(shí)空變異情況。

4)對選擇位點(diǎn)之間的變異病毒株數(shù)進(jìn)行統(tǒng)計(jì)，并以柱狀圖的形式對病毒的變異率進(jìn)行時(shí)間空間上的展示，方便用戶對病毒的變異情況進(jìn)行分析。

該系統(tǒng)以一種更為直觀的方式呈現(xiàn)病毒的具體信息，同時(shí)以圖表的方式，對病毒的變異情況進(jìn)行了展示，方便用戶對新型冠狀病毒的序列進(jìn)行分析，對新型冠狀病毒防治的研究具有重要的參考意義。