基于大數(shù)據(jù)挖掘與可視化的HPV病毒基因組分析研究

陳翠霞，曹宗富，李天君，于磊，喻浴飛，蔡瑞琨，羅敏娜，李乾，沈玥，陸超，高華方，馬旭*

(1.國(guó)家衛(wèi)生健康委科學(xué)技術(shù)研究所，北京 100081；2.國(guó)家人類遺傳資源中心，北京 102206)

宮頸癌是中國(guó)女性最常見的惡性腫瘤之一，僅次于乳腺癌[1]。宮頸癌主要病因是高危人乳頭瘤病毒(Human papilloma virus，HPV)的持續(xù)性感染[2]。目前，宮頸癌治療手段以毀損性手術(shù)為主，缺乏特效無(wú)創(chuàng)的阻斷方法，臨床檢測(cè)鑒定方法僅覆蓋23種亞型，占現(xiàn)存亞型10%，無(wú)法避免假陽(yáng)性率和漏診率問題，因此有必要對(duì)所有HPV進(jìn)行快速準(zhǔn)確的分型檢測(cè)[3-6]。最新的疫苗只覆蓋9種亞型的HPV病毒，保護(hù)71%的患者[3-6]。隨著高通量測(cè)序技術(shù)的成熟，病原體的溯源、進(jìn)化關(guān)系、基因組比較、致病危險(xiǎn)性預(yù)測(cè)等研究工作已成為傳染病的防控中重點(diǎn)研究方向[7-10]。有研究工作通過相同亞型HPV病毒的某段基因核苷酸的變異或基因整合熱點(diǎn)位置，來研究病毒進(jìn)化與疾病的相互關(guān)系以及病毒基因整合致病機(jī)制[10-13]，也有文獻(xiàn)對(duì)13種高危亞型的HPV病毒感染后其致病危險(xiǎn)性發(fā)展的自然史進(jìn)行監(jiān)測(cè)研究[14]。但仍舊缺乏一種簡(jiǎn)單易操作的自動(dòng)化流程來對(duì)HPV基因組大數(shù)據(jù)進(jìn)行全面深入的挖掘分析。鑒于此，本文基于比較基因組學(xué)研究方法，針對(duì)HPV病毒全基因組數(shù)據(jù)，設(shè)計(jì)和集成一系列的算法工具包，構(gòu)建了一種基于大數(shù)據(jù)挖掘技術(shù)的HPV基因組信息可視化分析流程框架，不僅覆蓋迄今為止發(fā)現(xiàn)的298種亞型的HPV病毒，且具備對(duì)HPV病毒基因組數(shù)據(jù)的深度比對(duì)挖掘能力。

方 法

HPV基因組數(shù)據(jù)挖掘分析包含3個(gè)步驟(圖1)：(1)數(shù)據(jù)預(yù)處理：從搜集HPV基因組數(shù)據(jù)最全面的genBank(genBank[http：//www.ncbi.nih.gov/nucleotide/SRA/genBank]數(shù)據(jù)庫(kù)下載所有病毒類基因組數(shù)據(jù)作為原始數(shù)據(jù)。采用Perl 5.0語(yǔ)言腳本從原始數(shù)據(jù)中抽提298個(gè)HPV病毒亞型的全基因組數(shù)據(jù)，及其關(guān)鍵核心基因(E6、E7)的核酸和氨基酸序列，采用Perl腳本將關(guān)鍵基因按照染色體上的順序連接組裝，建立關(guān)鍵核心基因組，全部整理為fasta格式。(2)系統(tǒng)進(jìn)化分析：采用Dnasp[15]軟件評(píng)估全基因組序列的替換飽和度，對(duì)于通過評(píng)估篩選的HPV全基因組，采用MAFFT[16]軟件進(jìn)行多序列比對(duì)，并采用FastTree[17]軟件的最大似然法完成進(jìn)化樹的構(gòu)建，使用Dendroscope 3.0[18]軟件實(shí)現(xiàn)全基因組進(jìn)化樹的展示、編輯和導(dǎo)出，至此完成全基因組進(jìn)化樹的構(gòu)建。另外，關(guān)鍵核心基因組進(jìn)化樹的構(gòu)建方法與全基因組進(jìn)化樹構(gòu)建方法相同。(3)氨基酸分類比較分析：基于序列組分和分類比對(duì)進(jìn)行氨基酸表達(dá)偏好模式分析。通過Perl腳本統(tǒng)計(jì)在HPV關(guān)鍵核心基因組的氨基酸序列中，20種氨基酸的占比，然后根據(jù)氨基酸的親疏水性、極性、酸堿性分類，匯總為HPV病毒氨基酸分類的占比，用R語(yǔ)言的Scale包對(duì)結(jié)果進(jìn)行歸一化校準(zhǔn)后得到氨基酸分類矩陣，最后用R語(yǔ)言的Heatmap包繪制氨基酸分類矩陣的熱圖(Heatmap)，進(jìn)而研究基因組分子結(jié)構(gòu)和氨基酸表達(dá)偏好模式與毒株生物學(xué)特性表型的關(guān)系。

圖1 HPV基因組數(shù)據(jù)挖掘分析框架流程圖

結(jié) 果

核酸和氨基酸分子水平結(jié)構(gòu)、氨基酸分類導(dǎo)致物種系統(tǒng)發(fā)育進(jìn)化的不同，進(jìn)而構(gòu)成HPV致病危險(xiǎn)性和主要感染侵襲部位不同的基礎(chǔ)條件。所以本文從基于分子結(jié)構(gòu)差異的系統(tǒng)進(jìn)化分析和氨基酸表達(dá)偏好模式比較分析兩個(gè)方面，來研究他們與亞型分化和致病危險(xiǎn)性之間的關(guān)系。

1.數(shù)據(jù)預(yù)處理：抽提出上述298種HPV病毒的全基因組和核心關(guān)鍵基因(E6、E7)核酸和氨基酸序列，采用Perl腳本將關(guān)鍵核心基因按照染色體上的順序連接組裝，建立關(guān)鍵核心基因組，全基因組序列與關(guān)鍵核心基因組序列均為fasta格式。從這298種HPV中選出致病危險(xiǎn)性較為明確的毒株37個(gè)和5個(gè)侵襲感染類型或致病危險(xiǎn)性不太明確的毒株[1-6，10，12-14，19-34]，共42個(gè)毒株，將其標(biāo)記感染類型和危險(xiǎn)性分類(字母S表示主要侵襲感染皮膚Skin；M表示主要侵襲感染粘膜Mucosa；N表示不明確；HR表示致病性為高危類型High-Risk；LR表示致病性為低危類型Low-Risk)后，與剩余的其他HPV病毒基因組作為參比毒株進(jìn)入后續(xù)分析。統(tǒng)計(jì)得知參比的HPV全基因組長(zhǎng)度范圍為7 080～8 104 bp。

圖2 42個(gè)HPV參比毒株在目前現(xiàn)存的298株HPV全基因組進(jìn)化樹上的分布情況

熱圖的紅色塊表示對(duì)該類型氨基酸較為偏好，占比較高；灰色塊表示該類型氨基酸占比較少；熱圖第一列是 非極性/疏水性(Non_polar-hydrophobic)類型氨基酸，第二列是極性/中性(Polarity-neutral)類型氨基酸， 第三列是堿性(Alkaline)類型氨基酸，第四列是酸性(Acidic)類型氨基酸圖3 42個(gè)參比HPV毒株的關(guān)鍵核心基因組進(jìn)化樹與其氨基酸分類熱圖對(duì)應(yīng)關(guān)系圖

3.氨基酸分類比較分析：基于進(jìn)化樹上毒株的位置，對(duì)42個(gè)參比毒株的關(guān)鍵核心基因組進(jìn)行氨基酸序列組分和分類比對(duì)分析，進(jìn)而發(fā)現(xiàn)氨基酸表達(dá)偏好模式與毒株生物學(xué)特性表型之間的關(guān)系(圖3)。通過Perl腳本分別統(tǒng)計(jì)HPV關(guān)鍵核心基因組的氨基酸序列中20種氨基酸各自在整條序列中所占百分比，然后將氨基酸分為非極性/疏水性(Non_polar-hydrophobic)、極性/中性(Polarity-neutral)、堿性(Alkaline)，酸性(Acidic)四大類，匯總HPV病毒20種氨基酸的百分比為四類氨基酸的占比，然后用R語(yǔ)言的scale包對(duì)結(jié)果進(jìn)行歸一化標(biāo)準(zhǔn)化后，得到氨基酸分類矩陣，最后用R語(yǔ)言的heatmap包繪制氨基酸分類矩陣的熱圖(Heatmap)，從而可以展示氨基酸表達(dá)偏好模式的差異與毒株生物學(xué)特性表型的關(guān)系。由圖3可見，皮膚高危型(SHR，)毒株比較偏好表達(dá)非極性/疏水性和酸性氨基酸產(chǎn)物；黏膜高危型(MHR，)毒株比較偏好表達(dá)極性/中性和堿性氨基酸產(chǎn)物，其中毒性較高的HPV16，HPV31，HPV18，HPV35，HPV73，HPV56[14]毒株也同時(shí)偏好表達(dá)酸性氨基酸；黏膜低危型(MLR，)毒株比較偏好表達(dá)非極性/疏水性和極性/中性氨基酸產(chǎn)物，有趣的是HPV40_MLR和HPV43_MLR的氨基酸表達(dá)種類模式與MHR的模式很類似，也偏好極性/中性和堿性氨基酸產(chǎn)物，HPV61_MLR和HPV81_NLR的表達(dá)模式與SHR類似，也偏好表達(dá)非極性/疏水性和酸性氨基酸產(chǎn)物，這也許正是HPV各亞型之間轉(zhuǎn)化以及低危和高危類型之間轉(zhuǎn)化的分子水平的物質(zhì)基礎(chǔ)；皮膚低危型()比較偏好表達(dá)非極性/疏水性和堿性氨基酸產(chǎn)物。

討 論

HPV病毒的親緣性或危險(xiǎn)性與進(jìn)化樹上位置密切相關(guān)，而進(jìn)化樹是基于基因組分子組分和結(jié)構(gòu)之間的遺傳距離繪制的。由圖2、圖3可見，HPV病毒根據(jù)親緣關(guān)系遠(yuǎn)近匯聚到不同的簇，聚集位置相近，說明可能是由同一祖先共同進(jìn)化而來，致病危險(xiǎn)性也相似。如果毒株呈現(xiàn)獨(dú)立分支，可以預(yù)測(cè)該毒株是新亞型。本方法曾成功應(yīng)用于鏈球菌新菌的分離鑒定工作[36-37]，而本文將該方法推廣應(yīng)用到HPV基因組分析，結(jié)果顯示全基因組進(jìn)化樹(參見圖2)的Clade分布趨勢(shì)與關(guān)鍵核心基因組進(jìn)化樹(參見圖3)結(jié)果相吻合，一方面表明關(guān)鍵核心基因E6/E7主導(dǎo)了HPV物種的進(jìn)化和分化，從分子進(jìn)化的角度驗(yàn)證E6/E7蛋白對(duì)宮頸癌的發(fā)生發(fā)展起決定性的作用。同時(shí)證明了本流程框架從細(xì)菌基因組分型研究[36-37]推廣應(yīng)用于病毒分型研究仍然是可行的，可以輔助HPV病毒新亞型鑒定、亞型間親緣性和危險(xiǎn)性的預(yù)測(cè)，為HPV病毒的防控提供依據(jù)。例如HPV82_NHR，已知屬于高危亞型，但是侵襲部位不很明確，進(jìn)化樹上與粘膜高危型聚集到一個(gè)clade(參見圖2、圖3左)，且其氨基酸偏好表達(dá)模式與MHR相同(參見圖3右)，因此預(yù)測(cè)其侵襲感染部位為粘膜類型，即為粘膜高危型(HPV82_MHR)；HPV8_SNR，已知其屬于侵襲皮膚類型，進(jìn)化樹上與皮膚高危型聚集到一個(gè)clade(參見圖2、圖3左)，且其氨基酸偏好表達(dá)模式與SHR相同(參見圖3右)，因此可以預(yù)測(cè)其為皮膚高危型(HPV8_SHR)；同樣，HPV54_NLR、HPV72_NLR均與粘膜低危型聚集到一個(gè)clade(參見圖2、圖3左)，同時(shí)兩個(gè)亞型毒株的氨基酸偏好表達(dá)模式與MLR相似，可以預(yù)測(cè)他們?yōu)檎衬さ臀Ｐ?，即分別為HPV54_MLR、HPV72_MLR。因此可以說，圖2、圖3從分子進(jìn)化和氨基酸分類偏好模式角度驗(yàn)證了文獻(xiàn)[19-22]結(jié)果的正確性，致病危險(xiǎn)性類似的毒株，其基因組結(jié)構(gòu)特點(diǎn)也類似，找到了HPV病毒的生物學(xué)表型特性的差異與核酸、氨基酸分子水平差異有緊密相關(guān)性的證據(jù)。

基于不同的氨基酸類型在基因組序列中的表達(dá)比率，來研究HPV基因組中氨基酸的表達(dá)偏好模式，進(jìn)而發(fā)現(xiàn)氨基酸表達(dá)偏好模式與病毒侵襲皮膚的特性之間的關(guān)系。由圖3可見，致病危險(xiǎn)性相似的毒株，氨基酸表達(dá)偏好模式也類似，從而保證了其物種基本的生物學(xué)特性。MHR和MLR共同偏好表達(dá)極性/中性(Polarity-neutral)類型氨基酸，可能說明極性/中性氨基酸的偏好表達(dá)與病毒侵襲粘膜的特性有密切關(guān)系，而MHR因?yàn)楦悠帽磉_(dá)堿性和酸性氨基酸而成為高危亞型，MLR因?yàn)楦悠帽磉_(dá)非極性/疏水性氨基酸產(chǎn)物而成為低危亞型；SLR和SHR共同偏好非極性/疏水性，可能表明非極性/疏水性氨基酸的高表達(dá)與病毒侵襲皮膚的特性有密切關(guān)系，SHR因?yàn)楦帽磉_(dá)酸性氨基酸而成為高危亞型，SLR因?yàn)楦帽磉_(dá)堿性氨基酸而成為低危亞型?；蚪M核酸分子結(jié)構(gòu)和基因結(jié)構(gòu)的差異，氨基酸偏好表達(dá)模式的不同，均導(dǎo)致了HPV病毒侵襲特性的不同和致病危險(xiǎn)性的多樣化。

綜上所述，本研究具有以下特點(diǎn)：(1)目前研究最多的HPV亞型數(shù)量是131個(gè)，本文囊括目前發(fā)現(xiàn)的298種HPV亞型毒株代表株，并整理獲得HPV全基因組序列和病毒所有基因(E1/E2/E4/E5/E6/E7/L1/L2)的序列、起止位點(diǎn)、長(zhǎng)度、分離病毒的亞型和名稱等詳細(xì)信息。(2)通過全基因組比較和系統(tǒng)進(jìn)化關(guān)系分析，可以根據(jù)新發(fā)病毒在已知致病危險(xiǎn)性毒株的進(jìn)化樹上的分布情況，實(shí)現(xiàn)對(duì)未知或新發(fā)病毒的亞型、感染侵襲特性及其致病危險(xiǎn)性的快速預(yù)測(cè)，為病毒的預(yù)防和控制提供證據(jù)支持和技術(shù)補(bǔ)充。(3)本文也嘗試探究氨基酸表達(dá)偏好模式與致病危險(xiǎn)性之間的關(guān)系，流程框架能夠迅速發(fā)現(xiàn)其氨基酸表達(dá)偏好模式，利用該表達(dá)模式推斷HPV毒株的侵襲特性和致病性危險(xiǎn)性，嘗試將基因組結(jié)構(gòu)差異和氨基酸偏好表達(dá)模式與生物表型特性做關(guān)聯(lián)性研究，為指導(dǎo)臨床用藥和控制感染提供依據(jù)。