基于LDA模型和MDS算法的多基因組可視化

雋立然　王亞東

摘 要：面向多基因組的研究，以建模多個體關系和比較個體差異為主要研究內容。多基因組可視化可以幫助研究者依據多個體關系，有目的地分析、比較多基因組之間的差異。多個基因組遺傳變異層面的比較，因為變異數量巨大、并且絕大部分變異并無信息性，故而很難在有限的顯示空間內可視化。本文根據多基因組可視化的需求，分析了多基因組可視化的數據降維策略，提出了基于LDA模型及KL散度的多基因組相似度求解方法，建立了基于MDS算法的多基因組可視化降維方法，并使用千人基因組第三階段的基因組變異數據，驗證上述方法的可靠性。

關鍵詞：可視化；多基因組；LDA模型；MDS算法

中圖分類號：TP18 文獻標識號：A 文章編號：2095-2163（2015）01-

Abstract： Multiple genome research focuses on modelling multiple individual relationships and comparing individual differences. Visualization of multiple genomes can help researchers to analyze and compare differences among multiple individual genomes intuitively. However， multiple genomes are difficult to be visualized in limited displaying space owing to the huge number of variants. And most of genome variants are less informative. This paper analyzed the data dimension reducing strategy for multiple genome visualization， proposes multiple genome similarity solving algorithm based on LDA model and KL-divergence， and designs multiple genome visualization dimention reducing method. The phase3 datasets of 1000 genomes project are used to verify the effectiveness and reliability of the above methods.

Keywords： Visualization； Multiple Genomes； LDA model； MDS Algorithm

0 引言

近年來隨著個體基因組測序的普及，測序的個體基因組數量大大增加，基于多個個體基因組數據的研究也日漸增多。遺傳學、人類學、社會學等許多學科都高度重視對人群的研究，早期的HapMap計劃[1]就有意識地搜集世界各地不同人群的基因組數據，作為其延續，2014年年中，千人基因組計劃公開發布了第三階段的數據，共包括屬于5個超級人群、26個人群的2 504個個體的基因組變異數據。在遺傳疾病的研究中，對多個疾病樣本與多個正常樣本的基因組進行比照分析、對多個疾病亞型的樣本基因組進行測試分析，均為常見的科學手段。因此，同時對多個個體的基因組進行比較、分析即已成為生命科學和醫學研究中的重要需求。多基因組可視化能夠顯著提升多個個體基因組的比較和分析效率，也是重要的研究課題。

多基因組可視化并非多個個人基因組可視化的簡單集成，特別是當需要可視化的個人基因組數量較多時，簡單集成的方法無法直觀地表達多個基因組之間的異同。多基因組可視化關注的是多個基因組之間的關系，也并非是基因組的一般性特征，這又不同于一般性的基因組可視化。多個基因組遺傳變異層面的比較，因為變異數量巨大、并且絕大部分變異并無信息性，故而很難在有限的顯示空間內可視化，也即使研究者很難從大量變異數據中篩選出重要的變異。通過幫助研究者們查看多個基因組在遺傳變異層面的比較結果，并且尋找多個基因組中對研究有用的變異，則是多基因組可視化的主要目的。

本文根據多基因組可視化的需求，探討了多基因組可視化面臨的主要問題，分析了多基因組可視化的數據降維策略；提出了基于LDA模型及KL散度的多基因組相似度求解方法，其中，LDA模型由于可以給出相似基因組之間的共同潛在特征相關的變異列表和概率分布，將更加有利于對研究者所關心的變異進行識別；并且建立了基于MDS算法的多基因組可視化降維方法；最后，本文使用千人基因組第三階段的基因組變異數據，分析和測試了上述方法的有效性。

1 LDA模型的基本理論

Latent Dirichlet Allocation（LDA）模型[3]是無監督學習的概率主題模型，該模型假設每個文檔有不同概率的多個主題，而文檔中的詞匯則通過這些主題以一定概率具體選擇而生成。因此，通過學習語料庫中的文檔，LDA模型可以獲取語料庫中潛在的主題，并得到每個文檔的主題混合分布，以及每個主題中的主題-詞匯概率。由于是無監督學習算法，使得LDA并不需要輸入標注后的語料庫，并且對于每一個主題，都可以找出一個詞的集合對其進行描述。LDA模型基于bag-of-words假設，即不考慮詞在文檔中的順序，并且也不考慮文檔之間的順序。但LDA并不假設詞匯或文檔是獨立同分布的。LDA模型可以用于文本主題識別，文本分類和文本相似度計算等問題。

LDA模型的主要任務是尋找使語料庫文檔具有較高概率的模型，并使語料庫以外的其他類似文檔也具有較高概率，以識別非語料庫文檔的主題。因此，LDA模型的基本策略是用一組隨機混合的隱含主題分布表示文檔，并使用詞匯的概率分布來描述每個主題。通過觀測到的語料庫內文檔中詞匯的分布估計文檔-主題向量和主題-詞匯矩陣等參數，從而獲得主題的詞匯描述、文本的主題相似度等信息。

LDA模型最早通過變分貝葉斯期望最大化算法（Variational Bayesian Expectation Maximization，VBEM）估計參數[3]，也可以使用較快的吉布斯采樣（Gibbs Sampling）方法估計參數[4]。在LDA模型基礎上，D.Blei等和D. Ramage等隨后又提出了有監督的LDA模型sLDA[5]和L-LDA[6]。

2基于LDA模型的多基因組相似度計算方法

基于此，一般化地考慮LDA模型，LDA模型試圖從可觀測到的離散數據單元和離散數據單元的無序集合的關聯關系中，為這些集合學習到一個有意義的隱含屬性，該屬性是集合包含其內容數據的標志，也是集合與集合之間進行語義性比較的基礎，并且該屬性還可以使用一些離散數據單元所描述或定義。但是LDA模型并不限定該屬性和數據之間或屬性和集合之間有因果關系。

綜上理解可知，LDA模型可以應用于多基因組的相似度計算和比較研究中。人類基因組之間有高達99.9%的相似性，個體基因組一般被表示為相對于一個標準的參考基因組的一組變異信息。因此基因組可以表示為許多變異的集合，對于基因組而言，這些變異只有分子位置上的順序關系，語義關聯上的順序關系可以被忽略。而根據不同的問題背景，該基因組可能具有不同的屬性，如不同人群、超級人群，或者疾病-正常、疾病的不同亞型等等。本文以屬于不同人群的多個基因組為例，應用LDA模型計算多基因組的相似程度，本例中，個體的基因組事實上可能是多個人群的混合，如混血。圖1顯示了多基因組相似度計算問題的數據與LDA模型術語間的映射關系。

本文根據先驗知識的變異預篩選策略能夠使多基因組相似度計算問題的規模降低到LDA模型的求解算法能夠求解的范圍內，并更好地識別有意義的人群-變異關系和基因組相似特征變異。

與一般的LDA模型解類似，人群（多基因組的子類別）的相似度以及多個個體基因組之間的相似度可以使用Kullback-Leibler Divergence，即KL散度（KL距離）[7]來刻畫，由于KL散度的不對稱性，也可以使用對稱KL散度，即KL散度的算術平均數、幾何平均數、調和平均數，或者JS散度（Jensen-Shannon divergence）及其平方根[8， 9]，本文將使用JS散度的平方根作為兩個個體基因組之間的相似性度量。

根據LDA模型的基本理論，使用變分貝葉斯期望最大化（VBEM）算法，可以迭代求解人群分布向量。

VBEM算法引入變分參數γ和?，簡化了原來由于θ、z和w的條件關系而難以求解的概率圖模型。指定了簡化的可優化下界的函數后，即需尋找使下界函數和真實聯合后驗分布的KL散度極小化的變分參數γ和?，具體公式為：

3 基于MDS算法的多基因組可視化降維算法

MDS方法的基本流程為：

（1）給定M個樣本的K維數據，計算每對樣本之間的相似度/距離，并存入M×M的矩陣Δ。

（2）把數據投射到低維（r維，r << K）空間，為樣本在r維空間中隨機初始位置，使用一個M×r 的矩陣X存放投影后每個樣本在r維空間的坐標。

（3）根據樣本在低維空間的坐標，計算每對樣本之間的距離，一般為歐氏距離，并存入M×M的矩陣D。

（4）測量Δ與D的差別，差別使用應力值衡量，計算公式如：

（5）如果應力值大于閾值，即低維空間的樣本距離關系還沒有足夠近似高維空間中的樣本相似關系，則移動矩陣X中的樣本坐標，使高維空間相似度高的樣本之間的距離減小，以減小總體應力值。

（6）重復（3）-（5），如果應力值小于閾值，或多次循環閾值差別不大（收斂到局部最優），算法停止。

本文使用LDA模型求解了多基因組之間的相似性，相似度量是KL散度。KL散度是不對稱的，但是在MDS算法中，作為輸入的高維空間樣本距離應是對稱的。因此在實踐中，常常使用KL散度的算數平均數、幾何平均數、調和平均數，或者JS散度及其平方根作為相似性的度量。

JS散度（Jensen-Shannon Divergence）是對稱和平滑版本的KL散度，數學定義如下：

MDS算法的另一個關鍵問題是應力函數的優化方法。本文采用SMACOF算法[10]最小化應力函數，應力函數定義為：

（5）

也就是說，最小化應力函數實際是盡量令。具體地，dij為r維空間上樣本i和j的歐氏距離。

4 實驗結果與結論分析

研究采用千人基因組第三階段數據作為本文方法的測試數據。

在遺傳過程中，子代個體將繼承兩個親代個體的部分變異，并產生少量（<100）新的變異。因此，對于子代個體來說，性狀和表型主要由親代個體遺傳信息的重新組合決定，這是子代與親代相似性的遺傳基礎，同時也將使親代的遺傳特征以變異為表現形式而保留在后代的基因組中。由于地理因素和社會因素，人類在漫長的進化和發展過程中，總是在一定的人群（population）范圍內通婚，這就使得每個人群中廣泛存在某些從祖先獲得的較穩定遺傳特征，而不見或少見于其他人群。典型的特征如膚色、瞳孔、發色等等。在基因組層面，這些遺傳特征可以用一個變異或多個變異的組合進行描述，而且這些變異在不同人群中，則呈現為高內聚、低耦合的特點。

圖2就是所有個體的基因組相似度計算和可視化結果。從圖2（a）中可以看到，本文的多基因組可視化方法，盡管采用了無監督的算法，但卻完美重現了全部5個超級人群的劃分：歐洲人群（European，EUR），東亞人群（East Asian，EAS），混血美洲人群（Ad Mixed American，AMR），南亞人群（South Asian，SAS）和非洲人群（African，AFR）。

圖2（a）還準確地定位了混血美洲人群的位置，即歐洲人群和非洲人群之間，但也延伸向東亞人群和非洲人群之間。這與中南美洲是歐洲殖民者、非裔奴隸和當地原住民長期混血的事實十分吻合。特別地，其中延伸向非洲人群和東亞人群中間部分的趨勢，也與中南美洲原住民是冰川期從歐亞大陸沿白令海峽遷往美洲大陸的理論自洽。而南亞人群與歐洲人群同屬白色人種，兩者的距離比非洲人群和東亞人群更加接近。

由于各超級人群的區別十分明顯，為了進一步討論超級人群內部各人群的相似與區別，圖2（b）用不同形狀代表各超級人群，用顏色區分超級人群內部的各個人群。超級人群下的各人群之間劃分也十分顯著。圖2（b）還顯示了在歐洲人群中，伊比利亞人群（IBS，圖b-1）最接近混血美洲人群中的波多黎各人群和哥倫比亞人群（PUR，CLM，圖b-2），這與中南美洲的主要殖民者是西班牙人和葡萄牙人也十分吻合。而南美洲太平洋沿岸國家秘魯的人群（PEL，圖b-3）則保留了較多的原住民血統。盡管混血美洲人群血統較為復雜，但本文的可視化方法仍然較為清晰地展示了混血美洲人群內部各個人群的區別，以及這些人群與其他超級人群的聯系。

4 結束語

通過對圖2的分析，作為無監督方法，本文的多基因組相似度計算和可視化降維算法完好地還原了千人基因組各超級人群的分布和內部各人群的區分。可視化所展示的各人群遠近關系符合生物學和歷史、地理上的各種證據，由此進一步證實了本文多基因組可視化方法的能力和準確性。

參考文獻：

[1] INTERNATIONAL H C. The International HapMap Project [J]. Nature， 2003， 426（6968）： 789-796.

[2] GENOMES P C， ABECASIS G R， AUTON A， et al. An integrated map of genetic variation from 1，092 human genomes [J]. Nature， 2012， 491（7422）： 56-65.

[3] BLEI D M， NG A Y， JORDAN M I. Latent Dirichlet allocation [J]. J Mach Learn Res， 2003， 3（4-5）： 993-1022.

[4] GRIFFITHS T L， STEYVERS M. Finding scientific topics [J]. Proc Natl Acad Sci U S A， 2004， 101 Suppl 1（5228-35.

[5] MCAULIFFE J D， BLEI D M. Supervised topic models； proceedings of the Advances in neural information processing systems， F， 2008 [C]// NIPS in Whistler， British Columbia， Canada，2007.

[6] RAMAGE D， HALL D， NALLAPATI R， et al. Labeled LDA： A supervised topic model for credit attribution in multi-labeled corpora； proceedings of the Proceedings of the 2009 Conference on Empirical Methods in Natural Language Processing： Volume 1-Volume 1， F， 2009 [C]// Association for Computational Linguistics. EMNLP in Singapore， 2009.

[7] KULLBACK S， LEIBLER R A. On information and sufficiency [J]. The Annals of Mathematical Statistics， 1951：79-86.

[8] JOHNSON D H， SINANOVIC S. Symmetrizing the kullback-leibler distance [J]. IEEE Transactions on Information Theory， 2001， 1（1）： 1-10.

[9] ENDRES D M， SCHINDELIN J E. A new metric for probability distributions [J]. Information Theory， IEEE Transactions on， 2003， 49（7）： 1858-60.

[10] De LEEUW J， MAIR P. Multidimensional scaling using majorization： SMACOF in R [J]. Department of Statistics， UCLA， 2011，