人源甲型H1N1流感HA基因進化特征分析

金夢哲，綦朝暉，李素麗

（石家莊鐵道大學 信息科學與技術(shù)學院，河北 石家莊 050043）

過去的三次流感大爆發(fā)（1918年由H1N1引發(fā)的西班牙流感，1957年由H2N2病毒引發(fā)的亞洲流感和1968年由H3N2引發(fā)的香港流感）都并非起源于人類，源頭實際為禽流感病毒的變異［1］。最近一次甲型H1N1流感疫情在2009年4月大規(guī)模爆發(fā)于墨西哥，研究人員從墨西哥和美國的感染人群中收集到的流感病毒未曾有過記錄。在短短不到一個月的時間里，四十多個國家證實了甲型H1N1流感的大規(guī)模爆發(fā)［2］。2009年6月11日，WHO宣布將甲型H1N1流感大流行警告級別提升為6級，全球進入流感大流行階段，H1N1流感已席卷全球160多個國家和地區(qū)。直到2010年，WTO宣布甲型H1N1流感的大流行已結(jié)束，但據(jù)統(tǒng)計目前流感高峰期的主流病毒仍為甲型H1N1流感病毒［3］。2013年底H1N1流感再度流行，美國聯(lián)邦疾病防治中心12月27日宣布，流感目前已經(jīng)擴散到全美十個州，引起四名兒童死亡，上萬成年人入院治療。所幸的是，這一病毒尚未發(fā)生大的抗原性變異。

隨著生物信息學在生命科學的研究中的興起，使得以計算機為工具對生物信息進行分析和計算的方法成為生物特征信息獲取的重要途徑。利用生物信息學方法對生物序列進行分析已成為人們的研究熱點。在目前對高維度空間中的生物序列的研究中，數(shù)據(jù)降維方法被成功地用來獲取低維數(shù)據(jù)以便進一步分析。利用生物信息學方法對病毒序列進行表達、分析病毒的進化特征對于預防流感的大爆發(fā)有著重要作用。

1 甲型H1N1病毒的HA基因樣本序列

流感大爆發(fā)與流感病毒編碼蛋白中的血凝素（HA）基因的變異密切相關(guān)，對于流感病毒的變異和傳播而言，最關(guān)鍵的是HA抗原所引起的人體免疫應答［4］，因此本文選擇以HA基因樣本序列來研究H1N1流感病毒的變異規(guī)律。

1.1 病毒樣本

本文所應用的數(shù)據(jù)集為從美國國家生物技術(shù)信息中心（NCBI）的流感病毒資源庫（Influenza Virus Resource，http：//www.ncbi.nlm.nih.gov/genomes/FLU/Database/nph－select.cgi？go＝database）獲取的H1N1病毒HA基因序列源數(shù)據(jù)。該流感病毒數(shù)據(jù)庫開放性地提供了流感宿主，感染國家、地區(qū)，蛋白質(zhì)類型，病毒亞型，序列長度和收集時間為分類的流感病毒的核苷酸序列和氨基酸序列數(shù)據(jù)?？紤]到2009年4月開始的甲型H1N1流感大爆發(fā)，具體的索引數(shù)據(jù)集設(shè)定為2009年3月到8月的長度不限的人源甲型H1N1流感病毒HA基因的蛋白質(zhì)全序列，如圖1所示。

圖1 索引收集時間為2009年3月1日到8月31日的H1N1病毒HA基因序列

滿足條件的甲型H1N1流感病毒HA基因序列共有3602條，按數(shù)據(jù)保存于FASTA文件proteinHA.fas中。

1.2 MUSCLE對位

圖2是proteinHA.fas中三條流感病毒原始序列的數(shù)據(jù)形式。每條H1N1流感病毒HA基因序列由兩行數(shù)據(jù)組成：第一行為該序列的描述，依次為該病毒HA序列的GenBank索引號，收集地點和時間等；第二行為蛋白質(zhì)序列的氨基酸排列。經(jīng)過查詢后，滿足條件的3602條序列有長有短，序列的長短不同在序列的聚類分析過程中造成很大的干擾。為了使聚類分析不受序列長度的影響，采用經(jīng)典的蛋白質(zhì)水平上多序列比對 MUSCLE（Multiple Protein Sequence Alignment）［5］方法。在算法層面上，MUSCLE是一個漸進比對并結(jié)合迭代優(yōu)化的綜合算法，該對位不但要求對位后序列長度均相等，而且不同序列在相同位置盡可能多地存在相同片段［6］。

圖2 H1N1流感病毒序列的部分原始數(shù)據(jù)

圖3 H1N1流感病毒序列MUSCLE對位后的部分數(shù)據(jù)

MUSCLE對位排列后的流感病毒序列如圖3所示，對某些序列氨基酸不足的位置用“－”占位，這樣造成了序列的錯位使得在相同位置更多的同種氨基酸片段得以產(chǎn)生。MUSCLE對位后，該流感病毒序列樣本中3602條甲型H1N1流感病毒HA基因序列長度均為582。

2 基于PCA的聚類分析

經(jīng)過MUSCLE對位排列后的流感病毒序列，本身具有對位病毒序列的相同位置越多相似程度越高的特征。一次流感大爆發(fā)簡單地說是由流感病毒結(jié)構(gòu)發(fā)生的變異引起的，對于大爆發(fā)期間獲取到的數(shù)據(jù)量相當大的流感病毒樣本來說，其中既有大爆發(fā)前通過變異產(chǎn)生的新型病毒菌株，又有舊的常年隨季節(jié)變化的病毒菌株。在研究流感大爆發(fā)的傳播特征時舊病毒是沒有產(chǎn)生作用的，因此需要挑選新型病毒菌株，這就產(chǎn)生了對H1N1病毒數(shù)據(jù)樣本中的序列的進行分類的思考。

2.1 病毒序列的數(shù)字化表達

MUSCLE對位后的序列本身有等長的性質(zhì)，但是要進行對位后的序列數(shù)字化聚類分析，首先應該對病毒的HA基因樣本序列數(shù)字化表達，使其轉(zhuǎn)化為數(shù)值序列。

構(gòu)成蛋白質(zhì)的氨基酸一共有二十種，簡寫分別為：P，L，Q，H，R，S，F(xiàn)，Y，W，C，T，I，M，K，N，A，V，D，E，G。為了能夠使得序列中所有的氨基酸排列信息都保存在數(shù)值序列中，必須要求病毒序列的數(shù)字化表達后形成的數(shù)值序列與原氨基酸字母序列能夠形成一一映射關(guān)系。本文引入Xiao等人氨基酸編碼的二進制符號語言的方法［7］，將氨基酸符號序列通過二進制編碼數(shù)字化表達為數(shù)值序列，具體氨基酸字母對應的編碼如表1所示。

表1 二十種氨基酸的二進制編碼

文獻［7］中數(shù)字化表達病毒序列的方法，是通過相似原則、互補原則和分子識別理論制定的，可以完整地捕捉到氨基酸序列的物理和化學特性。此外，由于1.2節(jié)MUSCLE對位后的序列中會產(chǎn)生“－”占位符，本文在文獻［7］數(shù)字表達序列的方法基礎(chǔ)之上利用00000表示“－”。在原始序列中會有一些位置由于統(tǒng)計失誤等原因產(chǎn)生其他字符，這些字符實際上沒有意義，統(tǒng)一編碼為00000。原始序列中還會使用“X”代表任意氨基酸，考慮到任意氨基酸位與占位符、無意義位的嚴格區(qū)別，我們將“X”編碼為11111。

這樣，之前得到的MUSCLE對位后的HA基因序列樣本就表達成由0、1組成的數(shù)值序列樣本。由于采用了5位二進制編碼，數(shù)值序列的長為對位后字母序列長度的5倍，即數(shù)值序列長度為2910。

2.2 PCA降維方法

主成分分析（Principal Component Analysis，PCA）從統(tǒng)計學的角度來說是一種多元統(tǒng)計方法，它可以將多個變量通過線性變換選出較少的重要變量。具體地說，PCA通過線性變換求出原數(shù)據(jù)矩陣到主元得分空間的映射矩陣，選取最重要的部分，將其余的維數(shù)省去，達到降維的目的［8］。PCA方法間接地對數(shù)據(jù)進行了壓縮處理，同時很大程度上保留了原數(shù)據(jù)的信息。

經(jīng)過2.1節(jié)病毒序列的數(shù)字化表達后，3602條甲型H1N1的HA蛋白序列數(shù)據(jù)轉(zhuǎn)化為3602×2910的高維數(shù)值矩陣，樣本數(shù)量為3602，每個樣本有2910個特征點。將矩陣簡化記為：

MATLAB提供了所需的主成分分析，只需在軟件中調(diào)用PCA命令princomp，即可實現(xiàn)對數(shù)字化表達后的高維數(shù)據(jù)的降維，具體的調(diào)用格式如下：

其中m為選取的score矩陣的列數(shù)，也就是選取的維數(shù)。隨著m的增大，選取的score矩陣的列數(shù)增加，ρ也逐漸趨近于1。為了保證選取的主成分能夠忽略掉不重要的成分，同時不會丟失太多的重要成分，這里要求ρ≥0.85［9］。

通過計算，m＝1時，ρ＝0.8322，m＝2時，ρ＝0.8769≥0.85，于是選取score的前兩列作為原矩陣X的主成分進行提取，得到一個n×2的降維矩陣。這樣，甲型H1N1流感病毒HA蛋白全序列數(shù)據(jù)集組成的3602×2910的原數(shù)值矩陣通過主成分分析降維方法得到了3602×2的數(shù)值矩陣。通過PCA降維方法，高維數(shù)值矩陣在保存了其數(shù)值特性的前提下降維成二維數(shù)值矩陣，即一條序列可以簡單地用兩個主成分得分來表示。

2.3 病毒樣本分類

將數(shù)值序列經(jīng)PCA得到的第一主成分作為x軸坐標，第二主成分作為y軸坐標，即可將3602條甲型H1N1流感病毒HA基因序列樣本表達在平面上，如圖4（a）所示。

圖4 病毒樣本的二維圖像表達

由于PCA中score矩陣是原數(shù)值矩陣通過線性變換得到的，因此圖4（a）上的點和原病毒序列之間存在著映射的關(guān)系。圖4（a）中，雖然3602條病毒序列樣本均屬于甲型流感病毒H1N1亞型，但是根據(jù)主元分析后的第一主成分得分，這些序列大致可以劃分成了兩個特征顯著不同的病毒簇：第一主成分在－14附近的病毒簇和第一主成分在0到1的病毒簇。

3 H1N1流感HA基因進化特征分析（2009.3－2009.8）

3.1 病毒的新舊區(qū)分

在2009年的甲型H1N1流感大爆發(fā)期間采集到的病毒菌株中，既有流感大爆發(fā)之前許多年里一直存在的沒有高致病性的舊病毒，也有在大爆發(fā)前或者大爆發(fā)期間產(chǎn)生的新變異的高傳染性的流感病毒。對病毒的新舊區(qū)分能夠?qū)ρ芯坎《镜倪M化特征找到正確的靶標，使得研究結(jié)果更加準確。

利用PCA方法對2009年3月到8月采集到的3602個病毒數(shù)據(jù)樣本進行了降維，并對降維得到的低維數(shù)據(jù)進行了兩類病毒簇的分類。經(jīng)過分類后的病毒序列的氨基酸構(gòu)成相似，因而具有相近的生物性質(zhì)。但是，到目前為止還無法確定圖4（a）中哪一區(qū)域是新型病毒簇，哪一區(qū)域是舊型病毒簇。

世界衛(wèi)生組織推薦以A/California/07/2009 2009/04/09作為甲型H1N1流感疫苗的參考菌株［10］，它代表著2009年流感大爆發(fā)中的新型病毒群體。因此，A/California/07/2009 2009/04/09病毒菌株一定屬于新型流感病毒簇。圖4（b）中，右下角的“＊”號位置標示出了該病毒樣本在二維平面圖中的位置，坐標為（0.8556，－1.2513）。在圖像上即可判斷出矩形區(qū)域表示新型流感病毒簇，圓形區(qū)域?qū)獮榕f流感病毒簇。除此之外，可以看出圓形區(qū)域集中在一個相對小的范圍內(nèi)，說明舊病毒的變異程度很低，然而矩形區(qū)域的范圍相對來說較大，這就代表了新型流感病毒的變異程度高，發(fā)生了頻繁的抗原飄移。

3.2 病毒傳播分析

通過PCA降維方法，原本高維的病毒序列樣本轉(zhuǎn)換成了二維數(shù)據(jù)，該樣本數(shù)據(jù)在平面上可以清晰地表現(xiàn)出來。我們將病毒樣本按菌株的收集時間順序排序，把排序序號作為z軸，即可得出流感病毒數(shù)據(jù)集的三維表示，如圖5：

圖5 病毒樣本的三維圖像表達

圖5顯示了第一、二主成分在時間序列上的變化，箭頭指向了在3月到8月收集到的H1N1流感病毒菌株群。通過觀察，所有3月份收集到的菌株均屬于舊型病毒簇，正與3月流感大爆發(fā)尚未發(fā)生這一事實相符。而4月開始，流感大爆發(fā)于墨西哥，新型H1N1病毒開始廣泛傳染，也與4月收集到的流感菌株一部分來自舊型病毒簇一部分來自新型病毒簇的事實相符。5、6、7、8四個月流感大爆發(fā)，收集到的流感菌株數(shù)量很大，且絕大部分屬于新型病毒簇。圖4（b）中，新型變異病毒由于抗原漂移造成的小幅度變異也分為幾個簇，但由于差異很小，只體現(xiàn)在了第二主成分的差異。在圖5中也可以看到，從4月開始數(shù)據(jù)樣本的第二主成分較多地在－1附近，而隨著時間推移到8月，數(shù)據(jù)樣本的第二主成分更多地分布在0.5到1，說明這些新型病毒隨著時間也在發(fā)生著細微的變化，在圖4（b）中，逐漸由I區(qū)向II區(qū)變化，但是這些微小變異并未造成抗原轉(zhuǎn)變而形成新的亞型。

4 結(jié)論

以2009年3月到8月的3602條甲型H1N1病毒HA基因全序列為數(shù)據(jù)樣本，將經(jīng)MUSCLE比對的符號序列數(shù)字化表達成一組等長的高維數(shù)值序列。而后以高維數(shù)值序列為數(shù)據(jù)矩陣，利用主成分分析（PCA）映射到得分空間，選取其中兩個主要成分進行降維，得到了病毒樣本的二維平面表達圖，利用H1N1流感疫苗制作的參考菌株在圖中確定并區(qū)分了新舊病毒類型，其中新型病毒對于研究流感大爆發(fā)中病毒的傳播有關(guān)鍵性作用。添加收集時間為第三維度，得到三維圖，通過第一、二主成分隨時間的變化分析了甲型H1N1流感病毒隨時間的變化，得到了甲型H1N1流感病毒HA基因在2009年3月到8月期間的進化特征。

［1］GARTEN R J，DAVIS C T，RUSSELL C A，et al.Antigenic and genetic characteristics of swine－origin 2009A （H1N1）influenza viruses circulating in humans［J］.Science，2009，325（5937）：197－201.

［2］Malik Peiris J S，POON L L M，GUAN Y.Emergence of a novel swine－origin influenza A virus（S－OIV）H1N1virus in humans［J］.Journal of Clinical Virology，2009，45（3）：169－173.

［3］GUNSON R N，CARMAN W F.During the summer 2009outbreak of"swine flu"in Scotland what respiratory pathogens were diagnosed as H1N1/2009？［J］.BMC infectious diseases，2011，11（1）：192.

［4］李國強.流感病毒H1亞型血凝素單抗庫的構(gòu)建及其抗原性變異分析［D］.福建：廈門大學，2009.

［5］EDGAR R C.MUSCLE：multiple sequence alignment with high accuracy and high throughput［J］.Nucleic acids research，2004，32（5）：1792－1797.

［6］楊凡.生物序列分析中若干問題的研究［D］.電子科技大學，2011.

［7］XIAO X，SHAO S，DING Y，et al.Using cellular automata to generate image representation for biological sequences［J］.Amino Acids，2005，28（1）：29－35.

［8］JACKSON J E.A user＇s guide to principal components［M］.John Wiley ＆ Sons，2005.

［9］QI Zhao Hui，WEI Ruo Yan.A combination dimensionality reduction approach to codon position patterns of eubacteria based on their complete genomes［J］.Journal of theoretical biology，2011，272（1）：26－34.

［10］World Health Organization.Recommended viruses for influenza vaccines for use in the 2010－2011northern hemisphere influenza season［J］.Wkly Epidemiol Rec，2010，85（10）：81－92.