動(dòng)物病毒宏基因組數(shù)據(jù)分析平臺(tái)的建立及應(yīng)用

史智賓，王靖飛

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

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

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

本研究建立動(dòng)物病毒宏基因組數(shù)據(jù)分析平臺(tái)，能夠高效地進(jìn)行動(dòng)物宏病毒組數(shù)據(jù)分析，能夠?yàn)閯?dòng)物病毒流行監(jiān)測(cè)及疾病防控提供有效的技術(shù)手段。

1 材料與方法

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

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

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

1.1.2 讀長(zhǎng)序列分析讀長(zhǎng)序列分析方法(Assembly-free methods)，是不依賴于序列組裝直接利用測(cè)序后讀長(zhǎng)序列進(jìn)行注釋?zhuān)饕糜跇悠分胁《疚⑸锏姆N類(lèi)組成及多樣性研究。其主要利用注釋軟件BLASTn 對(duì)得到的讀長(zhǎng)序列集進(jìn)行注釋?zhuān)y(tǒng)計(jì)分析讀長(zhǎng)序列分配至各病毒分類(lèi)單元情況。

1.1.3 重疊序列分析重疊序列分析方法(Assembly-based methods)，是依賴于序列組裝，將小而短的讀長(zhǎng)序列拼接成大而長(zhǎng)的重疊序列，利用這些重疊序列進(jìn)行基因注釋?zhuān)饕糜跇悠分幸阎《炯拔粗《净蛐蛄醒芯俊Ｔ摲椒ɡ媒M裝軟件MEGAHIT 將短的讀長(zhǎng)序列拼接成長(zhǎng)的重疊序列，使用注釋軟件BLASTn 對(duì)得到的重疊序列集進(jìn)行注釋?zhuān)瑱z測(cè)重疊序列與各病毒科病毒序列的同源性。

1.2 臨床樣品宏基因組分析

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

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

2 結(jié) 果

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

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

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

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

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

表2 序列組裝拼接效果評(píng)價(jià)組裝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 討 論

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

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

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

病毒科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-

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