史智賓，王靖飛

(中國農(nóng)業(yè)科學院哈爾濱獸醫(yī)研究所獸醫(yī)生物技術(shù)國家重點實驗室，黑龍江 哈爾濱 150069)

病毒是自然界最豐富的生物體之一，其具有廣泛的宿主及復雜的遺傳多樣性。病毒在微生物生態(tài)系統(tǒng)中處于重要的地位，然而在某些環(huán)境中病毒微生物群落的全球多樣性研究還處于初級階段[1]。據(jù)估計目前對病毒世界的探知還不及1 %[2]。病毒因具有獨特的結(jié)構(gòu)及生物學特征，且絕大多數(shù)病毒沒有專門的宿主細胞系，使得傳統(tǒng)的病毒分離鑒定方法只能鑒定出很少一部分病毒。此外，病毒本身也在不斷進化以適應(yīng)宿主環(huán)境，增加其核酸多樣性，產(chǎn)生新的病毒物種或新的病毒亞型，擴大其宿主范圍，從而嚴重威脅人類的生命健康[3]。傳統(tǒng)的病毒分離鑒定技術(shù)已無法滿足病毒研究的需求，開發(fā)研究新的病毒發(fā)現(xiàn)技術(shù)是研究病毒首要解決的問題。

病毒宏基因組學(Viral metagenomics)技術(shù)的出現(xiàn)彌補了病毒研究方面的缺陷。該技術(shù)是從環(huán)境及其它生物樣品中濃縮病毒粒子遺傳物質(zhì)，利用生物信息學分析病毒遺傳信息。隨著病毒宏基因組學技術(shù)的發(fā)展，科學家們已經(jīng)利用該技術(shù)獲得了成千上萬個病毒基因組及大的基因片段[4]。2001年，病毒宏基因組學首次應(yīng)用于海洋病毒組的研究，在近20年的研究發(fā)展中，其應(yīng)用范圍已延伸至土壤、湖水、下水道等無機環(huán)境[5-6]以及動物組織、血液、呼吸道、消化道等有機環(huán)境[7-9]中。病毒宏基因組學的發(fā)展拓展了人類對病毒世界的認知，為新病毒的分離鑒定及未知病毒的發(fā)現(xiàn)提供了新的思路[10-12]。近幾年病毒宏基因組學技術(shù)在獸醫(yī)研究領(lǐng)域發(fā)展迅速，為動物病毒研究提供了便利條件。然而，目前絕大部分病毒宏基因組數(shù)據(jù)分析經(jīng)由測序公司完成，商業(yè)化分析雖可以保證數(shù)據(jù)分析的標準化，但其存在成本高昂、耗時長等問題，無法滿足及時監(jiān)控新發(fā)疫情時病毒的流行傳播情況。因此，低成本、高效率的病毒宏基因組數(shù)據(jù)分析是目前數(shù)據(jù)分析工作的重中之重。

本研究建立動物病毒宏基因組數(shù)據(jù)分析平臺，能夠高效地進行動物宏病毒組數(shù)據(jù)分析，能夠為動物病毒流行監(jiān)測及疾病防控提供有效的技術(shù)手段。

1 材料與方法

1.1 病毒宏基因組數(shù)據(jù)分析平臺的建立 利用實驗室本地linux 服務(wù)器建立病毒宏基因組數(shù)據(jù)分析平臺，目前硬件配置：一臺安裝Red Hat 系統(tǒng)(版本為4.8.5)的Linux 服務(wù)器(Linux 版本為3.10.0)，該服務(wù)器存儲容量為1 T，內(nèi)存48 G，速度266 mHZ。表1所示為數(shù)據(jù)分析過程中使用軟件，所有軟件均安裝于服務(wù)器平臺。

表1 主要數(shù)據(jù)分析軟件應(yīng)用Application數(shù)據(jù)質(zhì)量控制Data Quality Control數(shù)據(jù)預(yù)處理Data pre-processing讀長注釋Read annotation序列拼接及組裝效果評價Assembly and assessment重疊序列注釋Contig annotation后處理分析Post-processing軟件名稱Software name FastQC Cutadapt,Bowtie2 BLAST MEGAHIT,QUAST BLAST GenMark,MEGA,BEAST

1.1.1 數(shù)據(jù)的質(zhì)量控制及預(yù)處理高通量測序后下機得到的數(shù)據(jù)是長度為150 bp～300 bp 的單端讀長數(shù)據(jù)集或雙端讀長數(shù)據(jù)集，稱之為原始數(shù)據(jù)(Raw data)。這樣的結(jié)果往往會含有質(zhì)量低的序列或測序過程中插入的接頭序列，若不及時處理，會對后續(xù)數(shù)據(jù)分析造成嚴重的干擾。使用質(zhì)量控制軟件FastQC 對原始數(shù)據(jù)中接頭序列及低質(zhì)量序列進行檢測，然后使用數(shù)據(jù)處理軟件Cutadapt 清洗，同時利用Bowtie2 軟件對數(shù)據(jù)中宿主背景清除，以得到可以用于后續(xù)分析的清潔數(shù)據(jù)集(Clean data)。基于病毒多樣性及新病原發(fā)現(xiàn)的研究，建立兩種數(shù)據(jù)分析方法，即讀長序列分析和重疊序列分析。

1.1.2 讀長序列分析讀長序列分析方法(Assembly-free methods)，是不依賴于序列組裝直接利用測序后讀長序列進行注釋，主要用于樣品中病毒微生物的種類組成及多樣性研究。其主要利用注釋軟件BLASTn 對得到的讀長序列集進行注釋，統(tǒng)計分析讀長序列分配至各病毒分類單元情況。

1.1.3 重疊序列分析重疊序列分析方法(Assembly-based methods)，是依賴于序列組裝，將小而短的讀長序列拼接成大而長的重疊序列，利用這些重疊序列進行基因注釋，主要用于樣品中已知病毒及未知病毒基因序列研究。該方法利用組裝軟件MEGAHIT 將短的讀長序列拼接成長的重疊序列，使用注釋軟件BLASTn 對得到的重疊序列集進行注釋，檢測重疊序列與各病毒科病毒序列的同源性。

1.2 臨床樣品宏基因組分析

1.2.1 樣品信息本研究中樣品數(shù)據(jù)為本實驗室從安徽省部分豬場采集健康豬鼻拭子后處理的混合樣品，經(jīng)由上海派森諾生物科技有限公司Illumina 測序平臺測序所得。

1.2.2 樣品數(shù)據(jù)分析利用質(zhì)量控制軟件檢測樣品原始數(shù)據(jù)中存在的接頭序列及堿基質(zhì)量值低的讀長，利用數(shù)據(jù)處理軟件對樣品數(shù)據(jù)中接頭序列及低質(zhì)量序列過濾，完成數(shù)據(jù)過濾后重新進行質(zhì)量控制，比較過濾前后數(shù)據(jù)質(zhì)量，經(jīng)質(zhì)量控制檢測合格后經(jīng)下游數(shù)據(jù)處理分析。讀長序列分析直接利用過濾后的讀長序列進行注釋并統(tǒng)計；重疊序列分析使用序列組裝軟件對過濾后得到的讀長序列進行序列拼接，將拼接后的重疊序列進行基因注釋并統(tǒng)計。

2 結(jié) 果

2.1 病毒宏基因組數(shù)據(jù)分析平臺的建立 經(jīng)過長期的數(shù)據(jù)分析軟件篩選和分析方案的優(yōu)化完善，建立起病毒宏基因組數(shù)據(jù)分析流程(圖1)。該數(shù)據(jù)分析平臺部署了數(shù)據(jù)質(zhì)量控制、數(shù)據(jù)預(yù)處理、序列拼接、序列注釋等應(yīng)用程序。根據(jù)動物新病毒發(fā)現(xiàn)中對敏感性和特異性的不同需求，建立了兩種基因序列注釋策略：讀長注釋分析和重疊序列分析。前者能提高注釋的敏感性，有利于發(fā)現(xiàn)低拷貝的病毒基因組序列；后者可提高注釋的特異性，保障注釋結(jié)果的準確性?？梢愿鶕?jù)不同研究需要選擇不同的分析方法。

2.2 臨床樣品的宏基因組分析

2.2.1 樣品數(shù)據(jù)過濾及質(zhì)量控制利用數(shù)據(jù)質(zhì)量控制軟件FastQC 對原始測序數(shù)據(jù)進行質(zhì)量控制檢測。堿基質(zhì)量分布圖中藍色線表示各位置堿基質(zhì)量平均數(shù)，一般平均數(shù)高于25，所有位置的10 %中位數(shù)大于20 表示測序質(zhì)量結(jié)果較好。接頭序列分布圖中主要對4 種通用測序接頭序列進行檢測，根據(jù)測序數(shù)據(jù)中接頭序列的存在情況進行處理。原始讀長序列的3' 末端堿基質(zhì)量值偏低(圖2a)并包含有Illumina 測序通用接頭序列(圖2c)。利用數(shù)據(jù)處理軟件Cutadapt[13]進行數(shù)據(jù)接頭序列及低質(zhì)量序列的去除，過濾后重新進行質(zhì)量控制檢測。處理后所有讀長序列各位置堿基質(zhì)量均高于35，且所有位置10 %中位數(shù)大于26 (圖2b)，并且數(shù)據(jù)集中所含有的接頭序列已被去除干凈(圖2d)。上述結(jié)果表明經(jīng)處理后得到的數(shù)據(jù)可以進行下游分析。

2.2.2 序列組裝效果評價利用本地化安裝的組裝效果評價軟件對得到的重疊序列集進行評價，A 為本地化建立的病毒宏基因組數(shù)據(jù)分析方法組裝后得到重疊序列集的組裝效果評價結(jié)果，B 為經(jīng)公司測序分析后得到的重疊序列集的組裝效果評價結(jié)果(表2)。N50 作為組裝效果評價中的重要指標，其是將拼接得到的重疊序列由長到短進行排序并累加，當累加和達到重疊序列總長度的50 %時，最后參與加和的那一條重疊序列即為拼接后重疊序列數(shù)據(jù)集的N50 的長度。結(jié)果A 與結(jié)果B 的N50 長度分別為1 052 bp 和1 002 bp，表明二者之間在序列拼接方面無明顯差異。

2.2.3 基因注釋測序后共計產(chǎn)生32233976 條讀長，每條讀長序列長度為150 bp，經(jīng)過接頭序列及低質(zhì)量序列去除后得到31490033 條讀長。讀長序列分析注釋后發(fā)現(xiàn)共有742196 條讀長與病毒序列有較高的相似性，占總體讀長序列數(shù)據(jù)的2.36 %。重疊序列分析得到415 117 條重疊序列，其中307 條被注釋到病毒基因組包括指環(huán)病毒科(Anelloviridae)、沙粒病毒科(Arenaviridae)、星狀病毒科(Astroviridae)、環(huán)狀病毒科(Circoviridae)、冠狀病毒科(Coronaviridae)、細小病毒科(Parvoviridae)、小 RNA 病毒科(Picornaviridae)、痘病毒科(Poxviridae)、Smacoviridae、Genomoviridae、逆轉(zhuǎn)錄病毒科(Retroviridae)和副粘病毒科(Paramyxoviridae)等12 個病毒科的64 種病毒(表3)。從注釋結(jié)果看，該宏病毒組數(shù)據(jù)注釋到冠狀病毒科的多個種屬(同源性為70 %～100 %)，其中與豬德爾塔冠狀病毒(Porcine deltacoronavirus)的同源性最高(99.775 %)。讀長序列分析耗時1.5 d，重疊序列分析耗時3.5 d，表明本地化建立的動物病毒宏基因組數(shù)據(jù)分析平臺能夠高效的完成宏病毒組數(shù)據(jù)分析工作。

表2 序列組裝拼接效果評價組裝Assembly contigs (≥0)contigs (≥1,000 bp)contigs (≥5,000 bp)contigs (≥10,000 bp)contigs (≥25,000 bp)contigs (≥50,000 bp)Total length (≥0)Total length (≥1,000 bp)Total length (≥5,000 bp)Total length (≥10,000 bp)Total length (≥25,000 bp)Total length (≥50,000 bp)contigs Largest contig Total length GC (%)N50 N75 L50 L75 N's per 100 kbp結(jié)果A Result A 415117 53202 1148 219 53 10 275075706 99924084 10376714 4550115 2121335 657354 186544 88869 190152253 49.15 1052 693 48471 105387 0.00結(jié)果B Result B 503629 51875 1109 224 54 12 307875756 97195435 10317758 474527 2254595 789813 196074 88714 193816451 49.10 1002 673 51589 111890 0.00

3 討 論

本研究利用本地化建立的動物病毒宏基因組學數(shù)據(jù)分析平臺對豬源鼻拭子宏病毒組數(shù)據(jù)進行了分析，結(jié)果表明該樣品中包含大量指環(huán)病毒、環(huán)狀病毒、冠狀病毒和細小病毒等。有研究發(fā)現(xiàn)其它動物如熊貓[14]、蝙蝠[15]等體內(nèi)也同樣檢測到有上述病毒的存在，表明這些病毒可能在動物體內(nèi)普遍存在。同時，在該樣品中還發(fā)現(xiàn)有微量豬星狀病毒、捷申病毒和白血病病毒等，這也揭示了豬體內(nèi)豐富的病毒多樣性。因此，本地化建立的數(shù)據(jù)分析平臺能夠滿足動物病毒宏基因組數(shù)據(jù)的分析。

本研究建立的動物病毒宏基因組數(shù)據(jù)分析平臺融合了數(shù)據(jù)質(zhì)控、數(shù)據(jù)過濾、序列拼接和基因注釋等多種分析手段。與之前文獻及測序公司已有的分析平臺相比，該分析平臺所選擇分析軟件兼容性強，分析效率更高，分析流程相對簡單，在本地化平臺能夠高效的完成數(shù)據(jù)分析工作。同時，本研究對兩種序列分析方法進行比較，發(fā)現(xiàn)二者在數(shù)據(jù)分析中各存在優(yōu)劣勢。在微生物群落綜合性與復雜性分析方面，讀長序列分析擁有著更大的優(yōu)勢，其能夠?qū)Σ《疚⑸锶郝涔δ芘c結(jié)構(gòu)進行聚類分析，若能有足夠的測序深度和參考數(shù)據(jù)庫的覆蓋范圍，該方法可以對任意復雜的病毒群落進行分析；而重疊序列分析只有在基因組有足夠的豐度時才能夠完成多個基因組的構(gòu)建，且在復雜的病毒微生物群落中只有少部分基因組可以通過組裝得到分析。此外，在計算成本方面，重疊序列分析需要高昂的計算成本來進行基因組的組裝及注釋。但相比較而言，重疊序列分析在新病毒的發(fā)現(xiàn)及整合微生物基因組學方面更具優(yōu)勢。重疊序列分析能夠?qū)⒖蓟驇鞜o親緣關(guān)系的基因組序列進行解析，而讀長序列分析無法分辨未知基因組序列；重疊序列分析可以利用得到的基因組片段信息支持純培養(yǎng)分離基因組分析，而讀長序列分析無法完成[16]；重疊序列分析在新病毒的發(fā)現(xiàn)[17]及病毒的遺傳進化分析方面應(yīng)用較廣。此外，將原始讀長序列轉(zhuǎn)化為有意義的微生物特征分析工具也在持續(xù)不斷的更新改進，基于讀長序列的病毒株水平分析目前也已被應(yīng)用[18-20]。

值得注意的是，不論哪種序列分析方法，都依賴于病毒微生物群落的組成及復雜性、測序的深度、測序數(shù)據(jù)集的大小和計算分析資源的可用性。因此，在實際的樣品數(shù)據(jù)分析時需要考慮各種因素對數(shù)據(jù)分析及結(jié)果的影響，同時利用兩種序列分析方法可以做到相互驗證及補充，更能保證結(jié)果的可靠性及準確性。

病毒科Viral family Anelloviridae Viral Genus Alphatorquevirus Betatorquevirus Iotatorquevirus序列注釋條數(shù)Number of contigs--113 98 74 29 17 25 14 Arenaviridae Astroviridae Kappatorquevirus Mammarenavirus Mamastrovirus Circoviridae Circovirus Unclassified Number of reads 295 226 10273 11632 4700 879 331 594 188 104015 262862 248 270 593 13 393 235 24 14 217 354 102 1-111-51--1 2 28 Coronaviridae Alphacoronavirus Betacoronavirus 55 4--2 3 23 23 23 17 23 24 852 4 23 23 23 23 Parvoviridae Deltacoronavirus Torovirus Unclassified Bocaparvovirus Copiparvovirus Dependoparvovirus Protoparvovirus 2458 2483 2450 2415 1127 2450 2401 277 204 2857 2444 2490 2490 2 478 16 32 38 89 89 519 239 2564 4 189 10455 1574 6208-34662 6 Virus name Simian TTV TTMV TTSV-1 TTSV-1a TTSV-1b TTSV-2 Porcine TTV TTSV-k2a TTSV-k2b GuanaritoMammarenavirus Luna Mammarenavirus AstV MaAstV Porcine AstV PCV 1 PCV 2 PCV 3 BatCV BoCV FSfaCV Po-Circo-like virus PoSCV PRCV TGEV BetaCoV BoCoV Calf-giraffe CoV Canine respiratory CoV Equine CoV Giraffe CoV Human CoV Murine CoV Murine hepatitis virus PHEV CoV Sable antelope CoV Sambar deer CoV Waterbuck CoV White-tailed deer CoV BtRs-BetaCoV Sparrow CoV Porcine torovirus Porcine CoV Porcine deltaCoV Porcine bocavirus Rs-BtBoV PPV4 Goose parvovirus PPV2 PPV3 PPV5 PPV6 Canine parvovirus 2a Canine parvovirus 2b FPV UTV Chimpanzee parv4 Human parvovirus 4 Teschovirus A BeAn 58058 virus Vaccinia virus Camel associated porprismacovirus MiGyV MoMuLV Porcine respirovirus 17 12 999 Tetraparvovirus 181 22 Picornaviridae Poxviridae Smacoviridae Genomoviridae Retroviridae Paramyxoviridae Teschovirus Chordopoxvirinae Orthopoxvirus Porprismacovirus Gemycircularvirus Gammaretrovirus Respirovirus 151 15 13 3 33 41 18 165491116----21-2-

本地化建立的動物病毒宏基因組數(shù)據(jù)分析平臺，能夠高效的進行動物宏病毒組數(shù)據(jù)分析，這對于新疫情暴發(fā)時及時確定病因，提出解決措施意義重大。然而，目前本地化建立的病毒宏基因組學數(shù)據(jù)分析平臺還處于初步應(yīng)用階段，未來數(shù)據(jù)分析過程中在保證工作與時間效率的同時，還需優(yōu)化結(jié)果的報告形式及可視化展示，實現(xiàn)高標準化的數(shù)據(jù)分析，更好的應(yīng)用于動物病毒宏基因組數(shù)據(jù)分析。該病毒宏基因組數(shù)據(jù)分析平臺的建立為潛在致病原及新病原的發(fā)現(xiàn)提供便利工具，為全國流行病學分析及疾病防控研究提供技術(shù)支持。