Meta分析及其在动物遗传育种中的应用

殷诗舒，何 俊，黄生强*，吴晓林,2,3*

(1.湖南农业大学 动物科学技术学院，湖南 长沙 410128；2.纽勤公司生物统计和生物信息部，内布拉斯加州 林肯，美国 68504；3.威斯康星大学 动物科学系，威斯康星 麦迪逊，美国 53706)

Meta分析中文通常翻译为元分析、后设分析、整合分析、综合分析或荟萃分析，是一类可以综合多项研究结果的统计分析方法。该方法通过对以往若干项研究结果进行系统定量综合分析，从而得到一个新的、更为可靠的的研究结论。Meta分析方法的使用可追溯到17世纪的天文学研究[1]。卡尔·皮尔森(Karl Pearson)发表于1904年《英国医学杂志》的一篇论文，用Meta分析的方法综合了几个伤寒接种实验的临床研究结果[2-3]，这也是临床研究中第一个Meta的研究报告。1940年，杜克大学心理学家Pratt等[4]出版了一本关于Meta分析的论文集，其中包括了从1882年到1939年间共145篇有关的研究论文。Meta分析在医学研究上的应用始于1955年，从20世纪70年代开始，教育研究领域提出了一些复杂的Meta分析统计模型，Glass首先使用了Meta分析这个术语[5]。

Meta分析的基本思想是假设K个试验研究了某一特定的研究对象，每个研究都有一定的估计准确性，同时也存在一定的估计误差。因此，借助于一定的统计学方法，可以综合所有K个研究结果，从而获得对于该研究对象可靠性更高的研究结论。Meta分析的具体方法很多，但实质上这些方法都是计算各项研究结果的一个加权平均值，不同之处在于这些权重的计算方法以及对于未知参数的不确定性的评估方式上的一些差异。因此，Meta分析通过综合评估以前的研究结果，可对研究对象做出更为高效和可靠的统计估计(推论)。此外，Meta分析还可以对比不同的研究结果，发现它们之间的变异(差异)或者一些内在规律，研究他们之间可能存在的关系[6]。因此，Meta分析又是对特定研究问题的系统性评估方法[7]。

1 Meta分析的方法与步骤

1.1 前期预备

1.1.1 试验设计 Meta分析首先要确定研究问题和研究目标，然后确定下步研究方案、文献检索的方法与策略以及统计模型和统计方法的选择等。这是Meta分析的第一步，或相当于Meta分析的试验设计。所有的研究都必须有一个特定的问题，但研究问题本身并不等同于研究目的。例如，研究某个候选基因的效应，目的可能是通过多项研究结果来重新估计其效应，或是在多群体间评估其效应的一致性(或差异性)，或是在更大的群体中验证其效应，或是利用Meta分析的结果来指导设计新的研究。因为研究目的不同，Meta分析在方法设计、统计模型和实际实施上也是不一样的。

Meta分析不能只用已经发表论文的数据，因为学术期刊一般都倾向于选择性地发表统计分析上显著的研究结果，从而导致Meta分析的结果产生偏差，这情形被称为“档案抽屉问题”(the file drawer problem)，又称为“发表偏差”(publication bias)。因此，Meta分析若只考虑已经发表的论文的研究结果，将不足以代表所考察问题研究的总体。因此减少Meta分析偏差的一种实用的方法是在分析中考虑一些未发表的研究结果[8-9]。作为未发表的文献数据库，可以搜索各大学的硕博士研究论文以及一些公共书目数据库。例如，NTIS(国家技术信息服务)包括了由600多个美国联邦机构赞助的所有研究项目(http://www.ntis.gov/)。

1.1.2 文献检索 文献检索最常用的是公共文献检索库，例如Web of Science(https://login.webofknowledge.com/)或PubMed(https://www.ncbi.nlm.nih.gov/pubmed/)。手动搜索这些数据库需要指定搜索语句和搜索方法。为了方便文献检索，大多数图书馆都已经将文献检索的控制词汇添加到索引系统中，这些控制词汇包括书籍的分类代码和主题以及文章的描述符等。实际的检索往往是复式的，常包括两个以上的项目。例如，在PubMed中检索发表在Genetics(美国遗传学会的遗传学杂志)上的乳腺炎(mastitis)QTL，其搜索方式如下：(QTL AND mastitis AND Genetics[journal])。其中“QTL”和“mastitis”都是PubMed中使用的Mesh标题，“Genetics”是所要检索的学术期刊的名称，“[journal]”是学术期刊的名称的搜索字段标签。

1.1.3 检索文献的评估和搜索成功与否的评判 文献检索中要收集的内容很多。以QTL的Meta分析为例，通常包括研究报告的出版信息(例如作者，发表年份，出版物来源)，物种和品种(包括地理亚群类型)，统计分析方法和统计模型，QTL位置和效应大小等。然后需要评估所搜集的每一个原始研究的质量，以决定是否可以将该研究结果纳入到Meta分析中[10]。Cook等[11]提供了一些框架性的规则，如果违反了这些规则，就会影响到对于原始研究结果的解释和其研究质量。Chalmers等[12]提出了一个评估原始研究的评分方法，包括三个方面：研究设计，研究实施和数据分析。例如，两位读者在对作者和资料来源均不知情的情况下，从这三个方面对各研究报告的质量进行评分，然后计算出一个百分比分数，作为对检索的文献评估的依据。一个比较正式的方法是需要组成一个专家小组，来最终决定一项原始研究的结果是否可以纳入到Meta分析中。

现有文献数量巨大，不可能获得与研究课题有关的每一篇文献。文献检索一般采用检索的查出率(Recall)和精确率(Precision)这两个概念来描述搜索过程的成功率和相关文献检出的成功率。文献查出率代表了文献检索过程的成功率，计算方法如下：

文献查出率(R)=[(实际检出的相关文献数量)/((应该检索文献的总数量)]×100%

(1)

文献查出率越高，表明文献检索过程的成功率也越高。但应该检索的文献总数量往往是未知的，因此总数量(分母)通常按照某些众所周知的经验或规则设定。文献精确率反应了所检索文献的正确(或错误)的比率。

计算公式如下：

文献精确率(P)=(实际检出且符合主题的文献数量)/(实际检出文献总数量)

(2)

文献检索的目标是有尽可能高的文献查出率和精确率。实际检索中，这个目标可以通过使用一定搜索策略的组合来实现。

1.2 对以往研究结果的Meta分析

1.2.1 定义Meta分析的研究结果指标 文献检索收集了原始研究的结果之后，需要确定一个共有的研究结果指标作为Meta分析的数据。例如，QTL或候选基因的Meta分析可以用QTL(基因)位置，效应大小，或者是任何一个假设检验或统计分析的统计量(比如P值，或是LOD分数)。Meta分析通常对原始研究的指标的选择几乎没有控制，而是依赖于从原始研究中所能够获得的结果。并且，对于同一个主题的多项原始研究而言，其研究结果的指标也可能不尽一样，因此需要做一些换算。例如，Wu和Hu[13]介绍了用原始研究中的数据对于QTL位置和QTL效应做换算的一般性方法。

1.2.2 选择统计分析方法和统计模型 目前可以应用于Meta分析的统计模型很多。本文由于篇幅所限，简要介绍固定效应Meta分析模型(fixed-effect meta-analytical model)、随机效应Meta分析模型(random-effect meta-analytical model)以及两点混合分布Meta分析模型(two-point mixture meta-analytical model)。其它有关的Meta分析统计模型，例如随机效应的贝叶斯Meta分析模型和非参数的DPP-Meta分析模型，参见有关文献[14-15]。

Meta分析最常用的模型是固定效应模型和随机效应模型。设定收集了K个独立的研究结果，每个研究结果都估计了某个效应量(γi)及其标准差(Si2)。设θi为第i个独立研究中该效应量的真实值。而该真实效应未知，有待估计。就统计模型而言，不同Meta分析方法的主要区别在于对所估计效应量的假设上的差异。

固定效应Meta分析模型假设该效应量的真实值在所有K个独立研究中是相同的(同质性假设)。因此，各独立研究所估计效应量的差异是源于试验的抽样误差。设定所估计的效应量服从正态分布：

γi～N(θi,σi2)

(3)

式中：θ1=θ2=…=θn=θ，σi2为 γi的方差。在Meta分析中，γi的方差是未知的，但可以用样本方差(Si2)作为其近似估计值。固定效应Meta分析模型中的未知效应量θ可以用最大似然法(maximum likelihood，ML)估计[13]。

实际情况中，某效应量的真实值可能是不同的(异质性假设)。如在候选基因(包括QTL)的Meta分析中，常常会考虑到两种异质性情形[13]：位点异质性和效应异质性。位点异质性指一个基因位点可能在一些群体中有效应，但在另一些群体中没有效应。效应异质性指一个基因对于所有群体都有效应，但其效应大小因群体而异。存在效应异质性的情形下，所估计效应量的变异决定于其真实效应的变化以及各独立研究中的抽样误差。

随机效应的Meta分析模型假设每个独立研究中该效应量的真实值都不一样，后者服从一个均值为0、方差为σ2的正态分布：

θi～N(μ,σ2)

(4)

因此，每个独立研究中估计的效应量服从如下正态分布：

γi～N(μ,σ2+Si2)

(5)

随机效应Meta分析模型中的未知参数(μ和σ2)可以通过约束最大似然法(REML)估计[13]。从随机效应模型所估计的各独立研究的效应量又称为最佳线性无偏预估(BLUP)。显然，如果σ2=0，等价于所研究效应量的同质性假设成立。则随机Meta分析模型(4)就等同于固定效应Meta分析模型(3)。但是如果存在效应量的异质性(σ2>0)，固定效应Meta模型会低估其真实效应量的方差，使计算的P值偏小，从而夸大统计显著性，使假阳性率升高。两点混合分布的Meta分析模型描述的是位点异质性：所研究的效应量在一些群体中不为零(有效应)，而在另一些群体中为零(没有效应)。假设每个独立研究的效应量服从正态分布，就可以用式(4)来描述，所不同的是, 对于一些群体，θi=θ，而对于另一些群体，θi=0。γi～α×(θi～N(θ,Si2))+(1-α)×(θ≡0)。该模型中的未知参数可以用最大期望算法(expectation-maximum algorithm，EM)来计算[13]。

1.3 Meta分析的软件

Meta分析的软件很多，比如：CMA、MetAnalysis、MetaWin、MIX、RevMan、以及WEasyMA。这些软件在功能、适用性和价格方面存在一些差别。但是对于相同的模型，这些软件所计算的结果基本上是相同的[16]。Meta分析的R软件包也非常多。目前在CRAN (The Comprehensive R Archive Network；https://cran.r-project.org/) Meta分析(或相关)的R软件包就有93个(https://CRAN.R-project.org/view=MetaAnalysis)。这些R软件包可以用来处理不同学科和不同统计模型的Meta分析问题。以下简单介绍三个比较实用的Meta分析R软件包：

(1)Metaphor: 可以处理固定效应，随机效应和混合效应的Meta分析模型，创建各种类型的元分析图[17]。(2)Bayesmeta: 用于贝叶斯随机效应Meta分析，提供对于模型参数的联合和边缘后验概率分布，预测分布，收缩效应，后验预测P值等功能[18]。(3)MultiMeta：可以处理多变量全基因组关联研究(GWAS)的Meta分析。多变量设置意味着每个单核苷酸多态性(SNP)的结果包括几个效应大小(也被称为“β系数”，每个性状一个)，以及相应的方差值，还包括beta之间的协方差。该软件包可以绘制QQ图，曼哈顿图以及各种自定义的汇总图[19]。

3 Meta分析的应用

3.1 Meta分析与QTL定位研究

数量性状基因座(QTL)指与特定表型性状相关的DNA区域。QTL定位研究是用连锁分析的方法，通过与决定数量性状的基因相连锁的分子标记，来检测和估计这些基因的位置和效应。在过去几十年里，动植物的QTL定位研究结果层出不穷。由于研究所使用的参考群体各不相同，许多问题随之产生。例如，对于影响某个经济性状的QTL而言，其染色体上的位置和效应在群体间或者是物种间的关系如何，是否具有可比性或一致性。Meta分析就可以帮助回答这些问题，并解决不同研究间的不确定性和分歧。然而用Meta分析来综合评估多个QTL研究的结果具有较大的挑战性，因为这些研究在标记密度，连锁图谱，样本大小，种群类型，实验设计和统计方法等方面都是不同的。如Goffinet和Gerber[20]提出了一种基于模型选择的QTL数量及其位置的Meta分析策略来研究QTL的一致性。Etzel和Guerra[21]用Meta分析的方法来解决研究之间的异质性问题，从而可以更加精确估计QTL的位置和遗传效应。Veyrieras等[22]提出了一个两步法的Meta分析方法，以整合多个独立的QTL作图试验的结果。Wu等[13]对QTL定位的Meta分析方法进行了比较系统性的综述和回顾。QTL定位的Meta分析一般考察QTL效应或显著性检验的P值。此外，Wu等[15]还将非参数的DPP模型应用于QTL位置的Meta分析，该模型可以根据QTL位置分族(类)，每一类可能代表同一个QTL，或者是包括了功能相关或是在遗传图谱上连锁的一组QTL。

Meta分析应用于动物QTL定位的研究起步相对较晚。虽然Meta分析方法在动物QTL定位上的研究报道数量不多，但是已经应用到主要的家畜物种，如牛、猪、羊[23]和鸡[24]。Ben和Goddard[25]使用Meta分析法估计影响数量性状的基因效应的分布。Khatkar等[26]通过文献综合和Meta分析的方法，发现了若干个效应一致性的染色体区域，其中比较引人注目的是6号染色体上49cM和87cM这两个位置上的QTL，分别解释了牛奶产量的4.2%和3.6%的遗传变异。Silva等[27]用Meta分析的方法重新评估了猪4号染色体上的QTL定位结果，表明Meta分析所检出的显著QTL的数量要少于原始研究的QTL数量，并且与原始研究相比，Meta分析将QTL置信区间降低了85%，更加适合于QTL的精细定位研究。刘璐等[28]研究结果表明，Meta分析可以有效整合多个原始研究中的结果，最终得到的QTL定位的统计功效高于每个研究结果，还可以更加精确地估计其效应。

3.2 在全基因组关联性分析中的应用

家畜和家禽的全基因组关联分析(GWAS)现在已广泛应用于鉴定复杂性状的遗传变异，不仅大大丰富了畜禽标记辅助选择中可利用的分子标记，还为这些性状分子机理的探索提供了重要依据[29]。

利用Meta分析的方法可以重新整合和评估以往的GWAS研究结果。GWAS的Meta分析通常考察分子标记(如SNP)的相关效应或显著性检验的P值。Bolormaa等[30]发现多性状的Meta分析模型不仅可以提高QTL的检出功效(power)，而且还可区分对于多性状有一因多效(pleiotropic)的QTL和在染色体位置上连锁的影响不同性状的QTL。Pausch等[31]用Meta分析的方法重新评估了多个奶牛品种DNA序列数据的GWAS结果，共鉴定了25个QTL，其中6号和20号染色体上显著的QTL还包括ABCG2和GHR这两个基因中的两个错义突变。奶牛的肠道甲烷(CH4)产量与乳脂肪酸(FA)分布之间存在一定关系，Lingen等[32]由此利用Meta分析牛奶FA曲线具有一定预测单位饲料CH4产量的潜力。Rabiee等[33]用Meta分析方法研究了添加脂肪的奶牛日粮对牛奶产量和成分的影响。Laud[34]使用Meta分析模型和序列数据，提高了用多品种GWAS方法定位QTL的准确度。Silvia等[35]用系统文献综合和Meta分析方法评估了人畜共患病系统，通过对坦桑尼亚牛中布鲁氏菌属(Mycobacteriumspp)和锥虫属(Trypanosomaspp)的流行率估计，结果显示研究之间具有高度异质性。

众多研究表明，Meta分析能够将各独立的研究结果统一分析整理、验证和寻找新的研究点。如为了鉴定与猪肉性状相关的显著基因组区域，Bernal等[36]使用了3个猪数据集(美国肉类动物研究中心，商业和密歇根州立大学猪资源群体)的数据，证实了先前报道的pHu，PRL和CKL(PRKAG3)和SF(CAPN1和CAST)的候选基因，并鉴定出CKL、IMF(ACSF3)和CIE a*(GYS1和FTL)的新候选基因。Meta分析还具有强大的识别基因组位点的功能，Guo等[37]对猪四肢长度的GWAS-meta分析中共鉴定出16个位点，其中包括4个新位点，还找出两个与猪肢骨长度相关的候选基因HMGA1和PLAG1。Minozzi等[38]采用GRAMMAR-CG方法使用两种不同定义的对照种群进行全基因组关联分析，证实了12号染色体上的关联SNP，以及仅在意大利群体中发现的未分配SNP与ParaTB；此外，当ELISA和组织阳性动物与组织阴性样品进行比较时，还发现了几个额外的基因组区域与ParaTB相关。在过去十余年里，基因组预测和基因组选择日益成为动物遗传改良的主流方法[39]。Bernal Rubio等[40]研究表明，Meta分析增加了关联分析的统计效率，还可以在统计模型中考虑种群结构和不同种群间方差组分的异质性。

3.3 在候选基因效应中的应用

候选基因(candidate gene，CG)为鉴定和表征复杂表型的主要基因提供了一种策略，如生产性状对疾病的易感性。候选基因大体上分为两类，一类是在功能上相关，或者是假定参与所研究性状的调节和表达的基因，另一类是在遗传连锁图谱上位于所发现QTL的区间内的基因。前者称为功能性候选基因(functional candidate gene)，后者称为位置性候选基因(positional candidate gene)。CG分析的优点是它既不需要开发像QTL定位研究的遗传分离群体(参考群体)，也不需要像GWAS那样做高密度分子标记的遗传分型。因此CG研究在试验方法上简便易行，试验成本也相对较低。CG分析仅仅考察所研究基因的多态性与目标性状变异间的关联性。或者对于某个疾病性状，发病和健康动物之间该基因的频率是否存在显著差异。过去的几十年里，有大量候选基因的研究报道。但对于某一个特定性状和基因而言，它们之间的关联性往往在不同的研究、不同的群体，甚至是同一个群体的不同重复试验中也可能存在差异。造成这些差异的原因可能是遗传性的，也可能是环境因素，甚至于可能与试验设计和实验方法、统计模型和统计方法，以及人为误差有关。Meta分析为重新整合和评估以往的研究结果提供了一个有效的方法。候选基因的Meta分析，通常考察候选基因的效应大小或者其显著性检验的P值。所采用的统计模型一般为固定效应模型和随机效应模型。例如Wu等[15]将贝叶斯参数模型和非参数模型应用于候选基因的Meta分析。

Meta分析在畜禽候选基因上应用比较广泛。在已发表的研究中，Alfoso[41]使用固定效应模型和随机效应模型的Meta分析方法，研究了ESR PvuII基因的多态性与母猪产仔数之间的关系，发现“AA”基因型组母猪的窝产仔数显著降低，并且该基因的效应在研究群体间具有明显的异质性。Grigoryev等[42]对人类、猪、大鼠和小鼠模型中肾脏IRI的基因表达谱进行Meta分析，分别确定了10个和11个新的候选基因，在验证了这种新方法的同时也可以把实验结果应用到人类的疾病中。LU等[43]利用Meta分析与基因组富集分析(GSEA)对不同物种肝细胞癌(HCC)基因表达的芯片数据进行分析，结果显示在Meta分析下，包括cdc25a在内的25个基因在细胞周期途径中有差异表达，结合GSEA综合分析，确定了HCC发病机理的关键途径和基因。Purfield等[44]为了了解奶牛和肉牛产犊性能，使用高密度基因芯片进行了全基因组关联研究，最后对每个产犊性能特征进行三个品种(荷斯坦，夏洛来和利穆赞)的Meta分析，结果显示产犊性能的一部分遗传变异是所有三个品种共同的。Tetens等[45]使用Illumina EquineSNP50芯片对782匹德国血统马进行了影响马肩隆高度的基因座全基因组关联扫描，结果证实了与其他马群体最近报道一致的ECA3，并且找到了一个代表该QTL的强候选基因座LCORL/NCAPG。Casas等[46]对肉牛循环中25-羟基维生素D(25OHD)浓度相关的基因组区域的基因进行研究，鉴定出了CYP2J2为与牛血清中25OHD浓度相关的候选基因。这些研究结果均表明Meta分析的确可以综合以往的研究结果，确定和鉴定出候选基因，验证基因的功能性。

4 小 结

Meta分析是一种功能强大的系统研究方法，能够为科学研究及时间决策提供有效的帮助。该方法可使用统计学方法来整合多项研究结果，提高统计分析的效率和统计推理的可靠性和准确性，解决不同研究之间的不一致和不确定性，同时也是对多项研究做系统评价的统计方法。Meta分析的结果实质上是对研究结果的加权平均，该方法的优点主要有：分析结果可以推广到更大的群体；随着更多数据的使用，可提高估算的精确性和准确性，又可增加统计功效提高检测效果；可以量化和分析不同研究结果的不一致性。如是因为抽样误差不一致，还是研究结果(部分)受研究间异质性的影响；可对总结性统计量(群体均值)做假设检验；包括利用间接变量(moderator)来解释研究之间的差异；可以用于研究发表偏差。

但是，Meta分析中存在的问题也是客观存在的。除了前面提到的“档案抽屉问题”(发表偏差)外，其他的问题包括与不报告非统计学显著效应的研究有关的问题(缺乏对于不显著结果的统计资料)和与统计方法有关的问题(对于统计模型的选择具有很大的主观性，以及各种统计模型的内在假设可能是不符合实际情况的)等。此外，对于Meta分析的批评还包括一些Meta分析的研究结论也可能被一定的社会经济和政治原因所影响[47]。

Meta分析统计模型的发展，包括参数模型、半参数模型和非参数模型[13]。且统计方法仍然在发展之中，新的统计方法和统计计算技术也将继续出现。如殷诗舒等[48]提出了一个并行计算的DPM-Meta分析模型，并且通过模拟数据应用于候选基因的研究。该研究从大数据并行计算的角度，推导出可以准确估计总体数据SNP效应的公式，证明了用加权平均的方法，并行计算各子集数据估计的SNP效应，如果该权重为估计SNP效应的方差的倒数(即精确度)，等同于一个固定效应的Meta分析模型。