基于LDA模型和MDS算法的多基因组的发展关联

　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. 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.
　　
　　[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.
　　. Information Theory， IEEE Transactions on， 2003， 49（7）： 1858-60.